LCOV - code coverage report
Current view: top level - cdn - pwden.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 224 281 79.7 %
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_pwden
       8             : CONTAINS
       9        1184 :   SUBROUTINE pwden(stars,kpts,banddos,oneD, input,mpi,noco,cell,atoms,sym, &
      10        1184 :        ikpt,jspin,lapw,ne,we,eig,den,results,f_b8,zMat,dos)
      11             :     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      12             :     !     In this subroutine the star function expansion coefficients of
      13             :     !     the plane wave charge density is determined.
      14             :     !
      15             :     !     This subroutine is called for each k-point and each spin.
      16             :     !
      17             :     !
      18             :     !     Two methods are implemented to calculate the charge density
      19             :     !     1) which uses the FFT. The effort in calculating the charge
      20             :     !        density is proportional to M * N * log(N) , M being number of
      21             :     !        states and N being number of plane waves. This is the method
      22             :     !        which we use for production runs
      23             :     !     2) the traditional method for calculating the charge density
      24             :     !        using the double summation. In this case the effort scales as
      25             :     !        M * N * N. The method is only used for test purposes or for
      26             :     !        special cases.
      27             :     !
      28             :     !
      29             :     !     INPUT:    eigen vectors
      30             :     !               reciprocal lattice information
      31             :     !               Brillouine zone sampling
      32             :     !               FFT information
      33             :     !
      34             :     !     OUTPUT:   den%pw(s)
      35             :     !               1) using FFT
      36             :     !
      37             :     !                2) traditional method
      38             :     !
      39             :     !                             -1             ef
      40             :     !               den%pw(g) = vol * sum{ sum{ sum{ sum{ w(k) * f(nu) *
      41             :     !                                  sp   k    nu   g'
      42             :     !                                     *
      43             :     !                                    c(g'-g,nu,k) * c(g',nu,k) } } } }
      44             :     !                or :
      45             :     !                             -1             ef
      46             :     !               den%pw(g) = vol * sum{ sum{ sum{ sum{ w(k) * f(nu) *
      47             :     !                                  sp   k    nu   g'
      48             :     !                                     *
      49             :     !                                    c(g',nu,k) * c(g'+g,nu,k) } } } }
      50             :     !
      51             :     !                den%pw(g) are actuall 
      52             :     ! 
      53             :     !                the weights w(k) are normalized: sum{w(k)} = 1
      54             :     !                                                  k                -6
      55             :     !                         a) 1                           for kT < 10
      56             :     !                f(nu) = {                           -1             -6
      57             :     !                         b){ 1+exp(e(k,nu) -ef)/kt) }   for kt >=10
      58             :     !
      59             :     !
      60             :     !                                      Stefan Bl"ugel, JRCAT, Feb. 1997
      61             :     !                                      Gustav Bihlmayer, UniWien       
      62             :     !
      63             :     !     In non-collinear calculations the density becomes a hermitian 2x2
      64             :     !     matrix. This subroutine generates this density matrix in the 
      65             :     !     interstitial region. The diagonal elements of this matrix 
      66             :     !     (n_11 & n_22) are stored in den%pw, while the real and imaginary part
      67             :     !     of the off-diagonal element are store in den%pw(:,3). 
      68             :     !
      69             :     !     Philipp Kurz 99/07
      70             :     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
      71             :     !
      72             :     !
      73             : !DEC$ NOOPTIMIZE
      74             : #include"cpp_double.h"
      75             :     USE m_forceb8
      76             :     USE m_pwint
      77             :     USE m_juDFT
      78             :     USE m_rfft
      79             :     USE m_cfft
      80             :     USE m_types
      81             :     USE m_fft_interface
      82             :     IMPLICIT NONE
      83             :     TYPE(t_lapw),INTENT(IN)       :: lapw
      84             :     TYPE(t_mpi),INTENT(IN)        :: mpi
      85             :     TYPE(t_oneD),INTENT(IN)       :: oneD
      86             :     TYPE(t_banddos),INTENT(IN)    :: banddos
      87             :     TYPE(t_input),INTENT(IN)      :: input
      88             :     TYPE(t_noco),INTENT(IN)       :: noco
      89             :     TYPE(t_sym),INTENT(IN)        :: sym
      90             :     TYPE(t_stars),INTENT(IN)      :: stars
      91             :     TYPE(t_cell),INTENT(IN)       :: cell
      92             :     TYPE(t_kpts),INTENT(IN)       :: kpts
      93             :     TYPE(t_atoms),INTENT(IN)      :: atoms
      94             :     TYPE(t_mat),INTENT(IN)        :: zMat
      95             :     TYPE(t_potden),INTENT(INOUT)  :: den
      96             :     TYPE(t_results),INTENT(INOUT) :: results
      97             :     TYPE(t_dos), INTENT(INOUT)    :: dos
      98             : 
      99             :     REAL,INTENT(IN)   :: we(:) !(nobd) 
     100             :     REAL,INTENT(IN)   :: eig(:)!(dimension%neigd)
     101             :     !----->  BASIS FUNCTION INFORMATION
     102             :     INTEGER,INTENT(IN):: ne
     103             :     !----->  CHARGE DENSITY INFORMATION
     104             :     INTEGER,INTENT(IN)    :: ikpt,jspin
     105             :     COMPLEX, INTENT (INOUT) ::  f_b8(3,atoms%ntype)
     106             : 
     107             :     !-----> LOCAL VARIABLES
     108             : 
     109             :     !----->  FFT  INFORMATION
     110             :     INTEGER :: ifftq2d,ifftq3d
     111             : 
     112             :     INTEGER  isn,nu,iv,ir,ik,il,im,in,istr,nw1,nw2,nw3,i,j
     113             :     INTEGER  ifftq1,ifftq2,ifftq3
     114             :     INTEGER  idens,ndens,ispin,jkpt,jsp_start,jsp_end
     115             :     REAL     q0,q0_11,q0_22,scale,xk(3)
     116             :     REAL     s
     117             :     COMPLEX  x
     118             :     INTEGER,PARAMETER::  ist(-1:1)=(/1,0,0/)
     119             :     REAL,PARAMETER:: zero   = 0.00,  tol_3=1.0e-3 
     120             :     !
     121        2368 :     INTEGER  iv1d(SIZE(lapw%gvec,2),input%jspins)
     122        2368 :     REAL wtf(ne),wsave(stars%kq3_fft+15)
     123        2368 :     REAL,    ALLOCATABLE :: psir(:),psii(:),rhon(:)
     124        1184 :     REAL,    ALLOCATABLE :: psi1r(:),psi1i(:),psi2r(:),psi2i(:)
     125        1184 :     REAL,    ALLOCATABLE :: rhomat(:,:)
     126        1184 :     REAL,    ALLOCATABLE :: kpsir(:),kpsii(:)
     127        1184 :     REAL,    ALLOCATABLE :: ekin(:)
     128        1184 :     COMPLEX, ALLOCATABLE :: cwk(:),ecwk(:)
     129             :     !
     130             :     LOGICAL l_real
     131             :     REAL     CPP_BLAS_sdot
     132             :     EXTERNAL CPP_BLAS_sdot
     133             :     COMPLEX  CPP_BLAS_cdotc
     134             :     EXTERNAL CPP_BLAS_cdotc
     135             : 
     136             :     LOGICAL forw
     137             :     INTEGER length_zfft(3)
     138        1184 :     COMPLEX, ALLOCATABLE :: zfft(:)
     139             : 
     140             :     
     141             :     !------->          ABBREVIATIONS
     142             :     !
     143             :     !     rhon  : charge density in real space
     144             :     !     ne    : number of occupied states
     145             :     !     nv    : number of g-components in eigenstate
     146             :     !     cv=z  : wavefunction in g-space (reciprocal space)
     147             :     !     psir   : wavefunction in r-space (real-space)
     148             :     !     cwk   : complex work array: charge density in g-space (as stars)
     149             :     !     den%pw : charge density stored as stars
     150             :     !     trdchg: logical key, determines the mode of charge density
     151             :     !             calculation: false (default) : fft
     152             :     !                          true            : double sum over stars
     153             :     !     we    : weights for the BZ-integration for a particular k-point
     154             :     !     omtil : volume (slab) unit cell, between -.5*D_tilde and +.5*D_tilde
     155             :     !     k1   : reciprocal lattice vectors G=G(k1,k2,k3) for wavefunction
     156             :     !     k2   :                             =k1*a_1 + k2*a_2 + k3*a_3
     157             :     !     k3   : where a_i= Bravais lattice vectors in reciprocal space
     158             :     !             kwi, depend on k-point.                            
     159             :     !     kq1d  : dimension of the charge density FFT box in the pos. domain
     160             :     !     kq2d  : defined in dimens.f program (subroutine apws).1,2,3 indicate
     161             :     !     kq3d  ; a_1, a_2, a_3 directions.
     162             :     !     kq(i) : i=1,2,3 actual length of the fft-box for which FFT is done.
     163             :     !     nstr  : number of members (arms) of reciprocal lattice (g) vector
     164             :     !             of each star.
     165             :     !     ng3_fft: number of stars in the  charge density  FFT-box
     166             :     !     ng3   : number of 3 dim. stars in the charge density sphere defined
     167             :     !             by gmax
     168             :     !     kmxq_fft: number of g-vectors forming the ng3_fft stars in the
     169             :     !               charge density sphere 
     170             :     !     kimax : number of g-vectors forming the ng3 stars in the gmax-sphere
     171             :     !     iv1d  : maps vector (k1,k2,k3) of wave function into one
     172             :     !             dimensional vector of cdn-fft box in positive domain.
     173             :     !     ifftq3d: elements (g-vectors) in the charge density  FFT-box
     174             :     !     igfft : pointer from the g-sphere (stored as stars) to fft-grid 
     175             :     !             and     from fft-grid to g-sphere (stored as stars)
     176             :     !     pgfft : contains the phases of the g-vectors of sph.     
     177             :     !     isn   : isn = +1, FFT transform for g-space to r-space
     178             :     !             isn = -1, vice versa
     179             : 
     180        1184 :     CALL timestart("pwden")
     181             : 
     182        1184 :     ALLOCATE(cwk(stars%ng3),ecwk(stars%ng3))
     183             : 
     184        1184 :     IF (noco%l_noco) THEN
     185             :        ALLOCATE ( psi1r(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     186             :             psi1i(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     187             :             psi2r(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     188             :             psi2i(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     189         664 :             rhomat(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1,4) )
     190             :     ELSE
     191         520 :        IF (zmat%l_real) THEN
     192             :           ALLOCATE ( psir(-stars%kq1_fft*stars%kq2_fft:2*stars%kq1_fft*stars%kq2_fft*(stars%kq3_fft+1)-1),&
     193             :                psii(1),&
     194         350 :                rhon(-stars%kq1_fft*stars%kq2_fft:stars%kq1_fft*stars%kq2_fft*(stars%kq3_fft+1)-1) )
     195         350 :           IF (input%l_f) ALLOCATE ( kpsii(1),&
     196             :                kpsir(-stars%kq1_fft*stars%kq2_fft:2*stars%kq1_fft*stars%kq2_fft*(stars%kq3_fft+1)-1),&
     197           0 :                ekin(-stars%kq1_fft*stars%kq2_fft:2*stars%kq1_fft*stars%kq2_fft*(stars%kq3_fft+1)-1))
     198             :        ELSE
     199             :           ALLOCATE ( psir(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     200             :                psii(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     201             :                zfft(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     202         170 :                rhon(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1) )
     203         170 :           IF (input%l_f) ALLOCATE ( kpsir(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     204             :                kpsii(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1),&
     205          12 :                ekin(0:stars%kq1_fft*stars%kq2_fft*stars%kq3_fft-1) )
     206             :        ENDIF
     207             :     ENDIF
     208             :     !
     209             :     !=======>  CALCULATE CHARGE DENSITY USING FFT
     210             :     ! 
     211             :     !
     212             :     !------> setup FFT
     213             :     !
     214        1184 :     ifftq1  = stars%kq1_fft
     215        1184 :     ifftq2  = stars%kq1_fft*stars%kq2_fft
     216        1184 :     ifftq3  = stars%kq1_fft*stars%kq2_fft*stars%kq3_fft
     217        1184 :     ifftq3d = stars%kq1_fft*stars%kq2_fft*stars%kq3_fft
     218        1184 :     ifftq2d = stars%kq1_fft*stars%kq2_fft
     219             :     !
     220        1184 :     nw1=NINT(stars%kq1_fft/4.+0.3)
     221        1184 :     nw2=NINT(stars%kq2_fft/4.+0.3)
     222        1184 :     nw3=NINT(stars%kq3_fft/4.+0.3)
     223             :     !
     224             :     !------> g=0 star: calculate the charge for this k-point and spin
     225             :     !                  analytically to test the quality of FFT
     226             :     !
     227        1184 :     q0 = zero
     228        1184 :     q0_11 = zero
     229        1184 :     q0_22 = zero
     230        1184 :     IF (noco%l_noco) THEN
     231         664 :        q0_11 = zero
     232         664 :        q0_22 = zero
     233         664 :        IF (.NOT.zmat%l_real ) THEN
     234        6168 :           DO nu = 1 , ne
     235        5504 :              q0_11 = q0_11 + we(nu) * CPP_BLAS_cdotc(lapw%nv(1),zMat%data_c(1,nu),1,zMat%data_c(1,nu),1)
     236        6168 :              q0_22 = q0_22 + we(nu) * CPP_BLAS_cdotc(lapw%nv(2),zMat%data_c(lapw%nv(1)+atoms%nlotot+1,nu),1, zMat%data_c(lapw%nv(1)+atoms%nlotot+1,nu),1)
     237             :           ENDDO
     238             :        ENDIF
     239         664 :        q0_11 = q0_11/cell%omtil
     240         664 :        q0_22 = q0_22/cell%omtil
     241             :     ELSE
     242         520 :        IF (zmat%l_real) THEN
     243        9458 :           DO nu = 1 , ne
     244        9458 :              q0=q0+we(nu)*CPP_BLAS_sdot(lapw%nv(jspin),zMat%data_r(1,nu),1,zMat%data_r(1,nu),1)
     245             :           ENDDO
     246             :        ELSE
     247        3434 :           DO nu = 1 , ne
     248         170 :              q0=q0+we(nu) *REAL(CPP_BLAS_cdotc(lapw%nv(jspin),zMat%data_c(1,nu),1,zMat%data_c(1,nu),1))
     249             :           ENDDO
     250             :        ENDIF
     251         520 :        q0 = q0/cell%omtil
     252             :     ENDIF
     253             :     !
     254             :     !--------> initialize charge density with zero
     255             :     !
     256        1184 :     IF (noco%l_noco) THEN
     257        3320 :        rhomat = 0.0
     258         664 :        IF (ikpt.LE.mpi%isize) THEN
     259         456 :           dos%qis=0.0
     260             :        ENDIF
     261             :     ELSE
     262     4587352 :        rhon=0.0
     263         520 :        IF (input%l_f) ekin=0.0
     264             :     ENDIF
     265             :     !
     266             :     !------> calculate:  wtf(nu,k) =  w(k)*f(nu,k)/vol
     267             :     !
     268        1184 :     wtf(:ne) = we(:ne)/cell%omtil
     269             :     !
     270             :     !------> prepare mapping from wave function box to cdn FFT box
     271             :     !
     272        1184 :     IF (noco%l_ss) THEN
     273             :        jsp_start = 1
     274             :        jsp_end   = 2
     275             :     ELSE
     276        1144 :        jsp_start = jspin
     277        1144 :        jsp_end   = jspin
     278             :     ENDIF
     279        2408 :     DO ispin = jsp_start,jsp_end
     280      240207 :        DO iv = 1 , lapw%nv(ispin)
     281             :           !                                              -k1d <= L <= k1d
     282             :           !                                              -k2d <= M <= k2d
     283             :           !                                              -k3d <= N <= k3d
     284      237799 :           il = lapw%gvec(1,iv,ispin)
     285      237799 :           im = lapw%gvec(2,iv,ispin)
     286      237799 :           in = lapw%gvec(3,iv,ispin)
     287             :           !
     288             :           !------>  L,M,N LATTICE POINTS OF G-VECTOR IN POSITIVE DOMAIN
     289             :           !         (since charge density box = two times charge density box
     290             :           !          wrap arround error should not occur )
     291             :           !                                           0<= L <=2*k1-1 = kq1_fft-1
     292             :           !                                           0<= M <=2*k2-1 = kq2_fft-1
     293             :           !                                           0<= N <=2*k3-1 = kq3_fft-1
     294             :           !
     295      237799 :           il = il  +  stars%kq1_fft * ist( isign(1,il) )
     296      237799 :           im = im  +  stars%kq2_fft * ist( isign(1,im) )
     297      237799 :           in = in  +  stars%kq3_fft * ist( isign(1,in) )
     298             :           !
     299      239023 :           iv1d(iv,ispin) =  in*ifftq2 + im*ifftq1 + il
     300             :        ENDDO
     301             :     ENDDO
     302             : 
     303             :     !
     304             :     !------------> LOOP OVER OCCUPIED STATES
     305             :     !
     306       36936 :     DO  nu = 1 , ne
     307             :        !
     308             :        !---> FFT transform c_nu,k(g) --> psi_nu,k(r), for each k-point
     309             :        !                                              and each nu-state
     310       17876 :        IF (noco%l_noco) THEN
     311    12973440 :           psi1r=0.0
     312    12973440 :           psi1i=0.0
     313    12973440 :           psi2r=0.0
     314    12973440 :           psi2i=0.0
     315             :           !------> map WF into FFTbox
     316        5504 :           IF (noco%l_ss) THEN
     317        9440 :              DO iv = 1 , lapw%nv(1)
     318        9280 :                 psi1r( iv1d(iv,1) )   = REAL( zMat%data_c(iv,nu) )
     319        9440 :                 psi1i( iv1d(iv,1) )   = AIMAG( zMat%data_c(iv,nu) )
     320             :              ENDDO
     321        9440 :              DO iv = 1 , lapw%nv(2)
     322        9280 :                 psi2r( iv1d(iv,2) ) =  REAL(zMat%data_c(lapw%nv(1)+atoms%nlotot+iv,nu))
     323        9440 :                 psi2i( iv1d(iv,2) ) = AIMAG(zMat%data_c(lapw%nv(1)+atoms%nlotot+iv,nu))
     324             :              ENDDO
     325             :           ELSE
     326      556192 :              DO iv = 1 , lapw%nv(jspin)
     327      550848 :                 psi1r( iv1d(iv,jspin) ) = REAL( zMat%data_c(iv,nu) )
     328      550848 :                 psi1i( iv1d(iv,jspin) ) = AIMAG( zMat%data_c(iv,nu) )
     329      550848 :                 psi2r(iv1d(iv,jspin))=REAL( zMat%data_c(lapw%nv(1)+atoms%nlotot+iv,nu))
     330      556192 :                 psi2i(iv1d(iv,jspin))=AIMAG(zMat%data_c(lapw%nv(1)+atoms%nlotot+iv,nu))
     331             :              ENDDO
     332             :           ENDIF
     333             : 
     334             :        ELSE
     335   251967824 :           psir=0.0
     336    15111240 :           psii=0.0
     337             :           !------> map WF into FFTbox
     338       12372 :           IF (zmat%l_real) THEN
     339     4119864 :              DO iv = 1 , lapw%nv(jspin)
     340        9108 :                 psir( iv1d(iv,jspin) ) = zMat%data_r(iv,nu)
     341             :              ENDDO
     342             :           ELSE
     343      655784 :              DO iv = 1 , lapw%nv(jspin)
     344      652520 :                 psir( iv1d(iv,jspin) ) =  REAL(zMat%data_c(iv,nu))
     345      655784 :                 psii( iv1d(iv,jspin) ) = AIMAG(zMat%data_c(iv,nu))
     346             :              ENDDO
     347             :           ENDIF
     348             :        ENDIF
     349             :        !
     350             :        !------> do (real) inverse FFT; notice that the array psir is filled from
     351             :        !        0 to ifftq3-1, but starts at -ifftq2 to give work space for rfft
     352             :        !
     353       17876 :        IF (noco%l_noco) THEN
     354        5504 :           isn = 1
     355             : 
     356        5504 :           CALL cfft(psi1r,psi1i,ifftq3,stars%kq1_fft,ifftq1,isn)
     357        5504 :           CALL cfft(psi1r,psi1i,ifftq3,stars%kq2_fft,ifftq2,isn)
     358        5504 :           CALL cfft(psi1r,psi1i,ifftq3,stars%kq3_fft,ifftq3,isn)
     359             : 
     360        5504 :           CALL cfft(psi2r,psi2i,ifftq3,stars%kq1_fft,ifftq1,isn)
     361        5504 :           CALL cfft(psi2r,psi2i,ifftq3,stars%kq2_fft,ifftq2,isn)
     362        5504 :           CALL cfft(psi2r,psi2i,ifftq3,stars%kq3_fft,ifftq3,isn)
     363             :        ELSE
     364       12372 :           isn = 1
     365       12372 :           IF (zmat%l_real) THEN
     366             :              CALL rfft(isn,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft+1,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,&
     367        9108 :                   nw1,nw2,nw3,wsave,psir(ifftq3d), psir(-ifftq2))
     368             : 
     369             :              ! GM forces part
     370        9108 :              IF (input%l_f) THEN
     371           0 :                 DO in=-1,stars%kq3_fft,2
     372           0 :                    DO im=0,ifftq2-1
     373           0 :                       ir = ifftq2 * in + im
     374           0 :                       ekin(ir) = ekin(ir) - wtf(nu) * eig(nu) * (psir(ir)**2 + psir(ir+ifftq2)**2)
     375             :                    ENDDO
     376             :                 ENDDO
     377             : 
     378           0 :                 DO j = 1,3
     379           0 :                    kpsir(ifftq3d:)=0.0
     380           0 :                    kpsir(-ifftq2d:ifftq3d)=0.0
     381           0 :                    DO iv = 1 , lapw%nv(jspin)
     382           0 :                       xk=lapw%gvec(:,iv,jspin)+lapw%bkpt
     383             :                       s = 0.0
     384           0 :                       DO i = 1,3
     385           0 :                          s = s + xk(i)*cell%bmat(i,j)
     386             :                       ENDDO
     387           0 :                       kpsir( iv1d(iv,jspin) ) = s * zMat%data_r(iv,nu)
     388             :                    ENDDO
     389             :                    CALL rfft(isn,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft+1,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,&
     390           0 :                         nw1,nw2,nw3,wsave,kpsir(ifftq3d), kpsir(-ifftq2))
     391           0 :                    DO in=-1,stars%kq3_fft,2
     392           0 :                       DO im=0,ifftq2-1
     393           0 :                          ir = ifftq2 * in + im
     394           0 :                          ekin(ir) = ekin(ir) + wtf(nu) * 0.5 * (kpsir(ir)**2 + kpsir(ir+ifftq2)**2)
     395             :                       ENDDO
     396             :                    ENDDO
     397             :                 ENDDO
     398             :              ENDIF
     399             :           ELSE
     400             :              !--------------------------------
     401             :              ! FFT transform
     402        3264 :              zfft = cmplx(psir,psii)
     403             :              if (isn == -1) then
     404             :                 forw = .true.
     405             :              else
     406        3264 :                 forw = .false.
     407             :              end if
     408        3264 :              length_zfft(1) = stars%kq1_fft
     409        3264 :              length_zfft(2) = stars%kq2_fft
     410        3264 :              length_zfft(3) = stars%kq3_fft
     411        3264 :              call fft_interface(3,length_zfft,zfft,forw)
     412        3264 :              psir = real(zfft)
     413        3264 :              psii = aimag(zfft)
     414             :              !--------------------------------
     415             :              ! GM forces part
     416        3264 :              IF (input%l_f) THEN
     417     1959720 :                 DO ir = 0,ifftq3d-1
     418      979944 :                    ekin(ir) = ekin(ir) - wtf(nu)*eig(nu)* (psir(ir)**2+psii(ir)**2)
     419             :                 ENDDO
     420             : 
     421        1176 :                 DO j = 1,3
     422     2939832 :                    kpsir=0.0
     423     2939832 :                    kpsii=0.0
     424      120876 :                    DO iv = 1 , lapw%nv(jspin)
     425      120372 :                       xk=lapw%gvec(:,iv,jspin)+lapw%bkpt
     426             :                       s = 0.0
     427      842604 :                       DO i = 1,3
     428      481488 :                          s = s + xk(i)*cell%bmat(i,j)
     429             :                       ENDDO
     430      120372 :                       kpsir( iv1d(iv,jspin) ) = s *  REAL(zMat%data_c(iv,nu))
     431      120876 :                       kpsii( iv1d(iv,jspin) ) = s * AIMAG(zMat%data_c(iv,nu))
     432             :                    ENDDO
     433             : 
     434             :                    !--------------------------------
     435             :                    ! FFT transform
     436         504 :                    zfft = cmplx(kpsir,kpsii)
     437             :                    if (isn == -1) then
     438             :                       forw = .true.
     439             :                    else
     440         504 :                       forw = .false.
     441             :                    end if
     442         504 :                    length_zfft(1) = stars%kq1_fft
     443         504 :                    length_zfft(2) = stars%kq2_fft
     444         504 :                    length_zfft(3) = stars%kq3_fft
     445         504 :                    call fft_interface(3,length_zfft,zfft,forw)
     446         504 :                    kpsir = real(zfft)
     447         504 :                    kpsii = aimag(zfft)
     448             :                    !--------------------------------
     449             : 
     450     5879328 :                    DO ir = 0,ifftq3d-1
     451     2939832 :                       ekin(ir) = ekin(ir) + wtf(nu) * 0.5 * (kpsir(ir)**2+kpsii(ir)**2)
     452             :                    ENDDO
     453             :                 ENDDO
     454             :              ENDIF
     455             :           ENDIF
     456             :        ENDIF
     457             :        !----> calculate rho(r) = sum w(k)*f(nu)*conjg(psi_nu,k(r))*psi_nu,k(r)
     458             :        !                         k,nu
     459             :        !      again, we fill rhon() from -ifftq2 to ifftq3-1 for the rfft
     460             :        !
     461       19060 :        IF (noco%l_noco) THEN
     462             :           !--->             in the non-collinear case rho becomes a hermitian 2x2
     463             :           !--->             matrix (rhomat).
     464    25941376 :           DO ir = 0,ifftq3d-1
     465    12967936 :              rhomat(ir,1) = rhomat(ir,1) + wtf(nu)*( psi1r(ir)**2 + psi1i(ir)**2 )
     466    12967936 :              rhomat(ir,2) = rhomat(ir,2) + wtf(nu)*( psi2r(ir)**2 + psi2i(ir)**2 )
     467    12967936 :              rhomat(ir,3) = rhomat(ir,3) + wtf(nu)* (psi2r(ir)*psi1r(ir)+psi2i(ir)*psi1i(ir))
     468    12973440 :              rhomat(ir,4) = rhomat(ir,4) + wtf(nu)* (psi2r(ir)*psi1i(ir)-psi2i(ir)*psi1r(ir))
     469             :           ENDDO
     470             :           !--->             in a non-collinear calculation the interstitial charge
     471             :           !--->             cannot be calculated by a simple substraction if the
     472             :           !--->             muffin-tin (and vacuum) charge is know, because the
     473             :           !--->             total charge does not need to be one in each spin-
     474             :           !--->             channel. Thus it has to be calculated explicitly, if
     475             :           !--->             it is needed.
     476        5504 :           IF ((banddos%dos.OR.banddos%vacdos.OR.input%cdinf)) THEN
     477           0 :              DO ir = 0,ifftq3d-1
     478           0 :                 psi1r(ir) = (psi1r(ir)**2 + psi1i(ir)**2)
     479           0 :                 psi2r(ir) = (psi2r(ir)**2 + psi2i(ir)**2)
     480             :              ENDDO
     481           0 :              isn = -1
     482           0 :              psi1i=0.0
     483           0 :              CALL cfft(psi1r,psi1i,ifftq3,stars%kq1_fft,ifftq1,isn)
     484           0 :              CALL cfft(psi1r,psi1i,ifftq3,stars%kq2_fft,ifftq2,isn)
     485           0 :              CALL cfft(psi1r,psi1i,ifftq3,stars%kq3_fft,ifftq3,isn)
     486           0 :              psi2i=0.0
     487           0 :              CALL cfft(psi2r,psi2i,ifftq3,stars%kq1_fft,ifftq1,isn)
     488           0 :              CALL cfft(psi2r,psi2i,ifftq3,stars%kq2_fft,ifftq2,isn)
     489           0 :              CALL cfft(psi2r,psi2i,ifftq3,stars%kq3_fft,ifftq3,isn)
     490           0 :              cwk=0.0
     491           0 :              DO ik = 0 , stars%kmxq_fft - 1
     492             :                 cwk(stars%igfft(ik,1))=cwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     493           0 :                      CMPLX(psi1r(stars%igq_fft(ik)),psi1i(stars%igq_fft(ik)))
     494             :              ENDDO
     495           0 :              DO istr = 1,stars%ng3_fft
     496           0 :                 CALL pwint(stars,atoms,sym, oneD,cell,istr,x)
     497           0 :                 dos%qis(nu,ikpt,1) = dos%qis(nu,ikpt,1) + REAL(cwk(istr)*x)/cell%omtil/REAL(ifftq3)
     498             :              ENDDO
     499             : 
     500           0 :              cwk=0.0
     501           0 :              DO ik = 0 , stars%kmxq_fft - 1
     502             :                 cwk(stars%igfft(ik,1))=cwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     503           0 :                                                               CMPLX(psi2r(stars%igq_fft(ik)),psi2i(stars%igq_fft(ik)))
     504             :              ENDDO
     505           0 :              DO istr = 1,stars%ng3_fft
     506           0 :                 CALL pwint(stars,atoms,sym, oneD,cell, istr, x)
     507           0 :                 dos%qis(nu,ikpt,input%jspins) = dos%qis(nu,ikpt,input%jspins) + REAL(cwk(istr)*x)/cell%omtil/REAL(ifftq3)
     508             :              ENDDO
     509             :           ENDIF
     510             :        ELSE
     511       12372 :           IF (zmat%l_real) THEN
     512      250705 :              DO in=-1,stars%kq3_fft,2
     513    59748252 :                 DO im=0,ifftq2-1
     514    59739144 :                    ir = ifftq2 * in + im
     515    59980741 :                    rhon(ir) = rhon(ir) + wtf(nu) * ( psir(ir)**2 + psir(ir+ifftq2)**2 )
     516             :                 ENDDO
     517             :              ENDDO
     518             :           ELSE
     519    30182784 :              DO ir = 0,ifftq3d-1
     520        3264 :                 rhon(ir)=rhon(ir)+wtf(nu)*(psir(ir)**2+psii(ir)**2)
     521             :              ENDDO
     522             :           ENDIF
     523             :        ENDIF
     524             :        !              DO ir = -ifftq2 , ifftq3-1
     525             :        !     +                      + wtf(nu)*(psi(ir+ifftq3d) * psi(ir+ifftq3d)
     526             :        !     +                               + psi(ir  ) * psi(ir  )
     527             :        !     +                                 )
     528             :        !              ENDDO
     529             : 
     530             :     ENDDO
     531             :     !
     532             :     !<<<<<<<<<<<<<< END OUTER LOOP OVER STATES NU  >>>>>>>>>>>>>>>>>>
     533             :     !
     534             :     !
     535             :     !----> perform back  FFT transform: rho(r) --> chgn(star)
     536             :     !        ( do direct FFT)                    = cwk(star)
     537             : 
     538             :     !--->  In a collinear calculation pwden is calles once per spin.
     539             :     !--->  However in a non-collinear calculation pwden is only called once
     540             :     !--->  and all four components of the density matrix (n_11 n_22 n_12
     541             :     !--->  n_21) have to be calculated at once.
     542        1184 :     ndens = 1
     543        1184 :     IF (noco%l_noco) ndens = 4
     544        4360 :     DO idens = 1,ndens
     545        3176 :        IF (noco%l_noco) THEN
     546     6408800 :           psi1r=0.0
     547        2656 :           isn = -1
     548        2656 :           CALL cfft(rhomat(0,idens),psi1r,ifftq3,stars%kq1_fft,ifftq1,isn)
     549        2656 :           CALL cfft(rhomat(0,idens),psi1r,ifftq3,stars%kq2_fft,ifftq2,isn)
     550        2656 :           CALL cfft(rhomat(0,idens),psi1r,ifftq3,stars%kq3_fft,ifftq3,isn)
     551             :        ELSE
     552             :           !--->  psir is used here as work array, charge is real ,but fft complex
     553         520 :           IF (zmat%l_real) THEN
     554         350 :              psir(ifftq3d:)=0.0
     555         350 :              IF (input%l_f) kpsir(ifftq3d:)=0.0
     556             :           ELSE
     557      895610 :              psir=0.0
     558      895610 :              psii=0.0
     559         170 :              IF (input%l_f) kpsir=0.0
     560         170 :              IF (input%l_f) kpsii=0.0
     561             :           ENDIF
     562             : 
     563         520 :           isn = -1
     564         520 :           IF (zmat%l_real) THEN
     565             :              CALL rfft(isn,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft+1,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,&
     566         350 :                   stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,wsave,psir(ifftq3d), rhon(-ifftq2))
     567         350 :              IF (input%l_f) CALL rfft(isn,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft+1,stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,&
     568           0 :                   stars%kq1_fft,stars%kq2_fft,stars%kq3_fft,wsave,kpsir(ifftq3d), ekin(-ifftq2))
     569             :           ELSE
     570             :              !--------------------------------
     571             :              ! FFT transform
     572         170 :              zfft = cmplx(rhon,psir)
     573             :              if (isn == -1) then
     574         170 :                 forw = .true.
     575             :              else
     576             :                 forw = .false.
     577             :              end if
     578         170 :              length_zfft(1) = stars%kq1_fft
     579         170 :              length_zfft(2) = stars%kq2_fft
     580         170 :              length_zfft(3) = stars%kq3_fft
     581         170 :              call fft_interface(3,length_zfft,zfft,forw)
     582         170 :              rhon = real(zfft)
     583         170 :              psir = aimag(zfft)
     584             :              !--------------------------------
     585             :              !+apw
     586         170 :              IF (input%l_f) THEN 
     587             :                 !--------------------------------
     588             :                 ! FFT transform
     589          12 :                 zfft = cmplx(ekin,psii)
     590             :                 if (isn == -1) then
     591          12 :                    forw = .true.
     592             :                 else
     593             :                    forw = .false.
     594             :                 end if
     595             :                 length_zfft(1) = stars%kq1_fft
     596             :                 length_zfft(2) = stars%kq2_fft
     597             :                 length_zfft(3) = stars%kq3_fft
     598          12 :                 call fft_interface(3,length_zfft,zfft,forw)
     599          12 :                 ekin = real(zfft)
     600          12 :                 psii = aimag(zfft)
     601             :                 !--------------------------------
     602             :              ENDIF
     603             :           ENDIF
     604             :        ENDIF
     605             :        !  ---> collect stars from the fft-grid
     606             :        !
     607     1686724 :        cwk=0.0
     608     1686724 :        ecwk=0.0
     609        3176 :        IF (noco%l_noco) THEN
     610     2202496 :           DO ik = 0 , stars%kmxq_fft - 1
     611             :              cwk(stars%igfft(ik,1))=cwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     612     2202496 :                                                            CMPLX(rhomat(stars%igq_fft(ik),idens),psi1r(stars%igq_fft(ik)))
     613             :           ENDDO
     614             :        ELSE
     615         520 :           IF (zmat%l_real) THEN
     616     1035208 :              DO ik = 0 , stars%kmxq_fft - 1
     617             :                 cwk(stars%igfft(ik,1))=cwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     618     1035208 :                                                               CMPLX(rhon(stars%igq_fft(ik)),zero)
     619             :              ENDDO
     620             :           ELSE
     621      313156 :              DO ik = 0 , stars%kmxq_fft - 1
     622             :                 cwk(stars%igfft(ik,1))=cwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     623         170 :                                                               CMPLX(rhon(stars%igq_fft(ik)),psir(stars%igq_fft(ik)))
     624             :              ENDDO
     625             :           ENDIF
     626             :           !+apw
     627         520 :           IF (input%l_f) THEN 
     628          12 :              IF (zmat%l_real) THEN
     629           0 :                 DO ik = 0 , stars%kmxq_fft - 1
     630             :                    ecwk(stars%igfft(ik,1))=ecwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     631           0 :                                                                    CMPLX(ekin(stars%igq_fft(ik)),zero)
     632             :                 ENDDO
     633             :              ELSE
     634       23016 :                 DO ik = 0 , stars%kmxq_fft - 1
     635             :                    ecwk(stars%igfft(ik,1))=ecwk(stars%igfft(ik,1))+CONJG(stars%pgfft(ik))*&
     636          12 :                                                                    CMPLX(ekin(stars%igq_fft(ik)),psii(stars%igq_fft(ik)))
     637             :                 ENDDO
     638             :              ENDIF
     639             :           ENDIF
     640             :           !-apw
     641             :        ENDIF
     642             :        !
     643        3176 :        scale=1.0/ifftq3
     644      496912 :        DO istr = 1 , stars%ng3_fft
     645      496912 :           cwk(istr) = scale * cwk(istr) / REAL( stars%nstr(istr) )
     646             :        ENDDO
     647        3176 :        IF (input%l_useapw) THEN
     648             : 
     649           0 :           IF (input%l_f) THEN
     650           0 :              DO istr = 1 , stars%ng3_fft
     651           0 :                 ecwk(istr) = scale * ecwk(istr) / REAL( stars%nstr(istr) )
     652             :              ENDDO
     653           0 :              CALL force_b8(atoms,ecwk,stars, sym,cell, jspin, results%force,f_b8)
     654             :           ENDIF
     655             :        ENDIF
     656             :        !
     657             :        !---> check charge neutralilty
     658             :        !
     659        3176 :        IF ((idens.EQ.1).OR.(idens.EQ.2)) THEN
     660        1848 :           IF (noco%l_noco) THEN
     661        1328 :              IF (idens.EQ.1) THEN
     662         664 :                 q0 = q0_11
     663             :              ELSE
     664         664 :                 q0 = q0_22
     665             :              ENDIF
     666             :           ENDIF
     667        1848 :           IF ( ABS( q0 ) .GT. 1.0e-9) THEN
     668        1720 :              IF ( ABS( q0 - REAL(cwk(1)) )/q0 .GT. tol_3 ) THEN
     669           0 :                 WRITE(99,*) "XX:",ne,lapw%nv
     670           0 :                 IF (zmat%l_real) THEN
     671           0 :                    DO istr=1,SIZE(zMat%data_r,2)
     672           0 :                       WRITE(99,*) "X:",istr,zMat%data_r(:,istr)
     673             :                    ENDDO
     674             :                 ELSE
     675           0 :                    DO istr=1,SIZE(zMat%data_c,2)
     676           0 :                       WRITE(99,*) "X:",istr,zMat%data_c(:,istr)
     677             :                    ENDDO
     678             :                 ENDIF
     679           0 :                 WRITE ( 6,'(''bad quality of charge density'',2f13.8)')q0, REAL( cwk(1) )
     680           0 :                 CALL juDFT_warn('pwden: bad quality of charge')
     681             :              ENDIF
     682             :           ENDIF
     683             :        ENDIF
     684             :        !
     685             :        !---> add charge density to existing one
     686             :        !
     687        4360 :        IF(idens.LE.2) THEN
     688             :           !--->       add to spin-up or -down density (collinear & non-collinear)
     689        1848 :           ispin = jspin
     690        1848 :           IF (noco%l_noco) ispin = idens
     691      349136 :           DO istr = 1 , stars%ng3_fft
     692      349136 :              den%pw(istr,ispin) = den%pw(istr,ispin) + cwk(istr)
     693             :           ENDDO
     694        1328 :        ELSE IF (idens.EQ.3) THEN
     695             :           !--->       add to off-diag. part of density matrix (only non-collinear)
     696       73888 :           DO istr = 1 , stars%ng3_fft
     697         664 :              den%pw(istr,3) = den%pw(istr,3) + cwk(istr)
     698             :           ENDDO
     699             :        ELSE
     700             :           !--->       add to off-diag. part of density matrix (only non-collinear)
     701       73888 :           DO istr = 1 , stars%ng3_fft
     702         664 :              den%pw(istr,3) = den%pw(istr,3) + CMPLX(0.0,1.0)*cwk(istr)
     703             :           ENDDO
     704             :        ENDIF
     705             : 
     706             :     ENDDO
     707             : 
     708        1184 :     DEALLOCATE(cwk,ecwk)
     709             : 
     710        1184 :     IF (noco%l_noco) THEN
     711         664 :        DEALLOCATE ( psi1r,psi1i,psi2r,psi2i,rhomat )
     712             :     ELSE
     713         520 :        DEALLOCATE ( psir,psii,rhon )
     714         520 :        IF (input%l_f) DEALLOCATE ( kpsir,kpsii,ekin)
     715             :     ENDIF
     716             : 
     717        1184 :     CALL timestop("pwden")
     718             : 
     719        3552 :   END SUBROUTINE pwden
     720             : END MODULE m_pwden

Generated by: LCOV version 1.13