LCOV - code coverage report
Current view: top level - vgen - pw_tofrom_grid.F90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 94.3 % 176 166
Test Date: 2025-06-14 04:34:23 Functions: 100.0 % 4 4

            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              :    USE m_types_fftGrid
       9              :    PRIVATE
      10              :    REAL,PARAMETER:: d_15=1.e-15
      11              : 
      12              :    TYPE(t_fftgrid) :: fftgrid
      13              :    INTEGER         :: griddim
      14              :    REAL            :: gmax
      15              :    REAL            :: gmaxz
      16              : 
      17              :    PUBLIC :: init_pw_grid, pw_to_grid, pw_from_grid, finish_pw_grid
      18              : CONTAINS
      19          339 :   SUBROUTINE init_pw_grid(stars,sym,cell,xcpot)
      20              :     USE m_types
      21              :     use m_types_xcpot
      22              :     IMPLICIT NONE
      23              :     TYPE(t_stars),INTENT(IN)      :: stars
      24              :     TYPE(t_sym),INTENT(IN)        :: sym
      25              :     TYPE(t_cell),INTENT(IN)       :: cell
      26              :     CLASS(t_xcpot),INTENT(IN),OPTIONAL :: xcpot
      27              : 
      28              :       !---> set up pointer for backtransformation of from g-vector in
      29              :       !     positive domain of xc density fftbox into stars.
      30              :       !     also the x,y,z components of the g-vectors are set up to calculate
      31              :       !     derivatives.
      32              :       !     in principle this can also be done in main program once.
      33              :       !     it is done here to save memory.
      34              : 
      35          339 :     if (present(xcpot)) THEN
      36          337 :       gmax=xcpot%gmaxxc
      37          337 :       gmaxz=stars%gmaxz
      38          337 :       if (xcpot%needs_grad().and.gmax>0.0) then
      39          276 :          if (gmaxz>0.0) then
      40            0 :             call fftgrid%init(cell,sym,gmax,gmaxz)
      41              :          else
      42          276 :             call fftgrid%init(cell,sym,gmax)
      43              :          end if
      44              :       else
      45          244 :          call fftgrid%init((/3*stars%mx1,3*stars%mx2,3*stars%mx3/))
      46           61 :          gmax=stars%gmax
      47              :       endif
      48              :     else
      49            2 :       gmax=stars%gmax
      50            2 :       gmaxz=stars%gmaxz
      51            2 :       if (gmaxz>0.0) then
      52            0 :          call fftgrid%init(cell,sym,gmax,gmaxz)
      53              :       else
      54            2 :          call fftgrid%init(cell,sym,gmax)
      55              :       end if
      56              :    endif
      57          339 :     griddim=size(fftgrid%grid)
      58              : 
      59          339 :   END SUBROUTINE init_pw_grid
      60              : 
      61          702 :   SUBROUTINE pw_to_grid(dograds,jspins,l_noco,stars,cell,den_pw,grad,xcpot,rho,rhoim)
      62              :     !.....------------------------------------------------------------------
      63              :     !------->          abbreviations
      64              :     !
      65              :     !     ph_wrk: work array containing phase * g_x,gy......
      66              :     !     den%pw: charge density stored as stars
      67              :     !     rho   : charge density stored in real space
      68              :     !     v_xc   : exchange-correlation potential in real space
      69              :     !     exc   : exchange-correlation energy density in real space
      70              :     !     kxc1d  : dimension of the charge density fft box in the pos. domain
      71              :     !     kxc2d  : defined in dimens.f program (subroutine apws).1,2,3 indic
      72              :     !     kxc3d  ; a_1, a_2, a_3 directions.
      73              :     !     kq(i) : i=1,2,3 actual length of the fft-box for which fft is done
      74              :     !     nstr  : number of members (arms) of reciprocal lattice (g) vector
      75              :     !             of each star.
      76              :     !     nxc3_fft: number of stars in the  charge density  fft-box
      77              :     !     ng3   : number of 3 dim. stars in the charge density sphere define
      78              :     !             by gmax
      79              :     !     kmxxc_fft: number of g-vectors forming the nxc3_fft stars in the
      80              :     !               charge density or xc-density sphere
      81              :     !     kimax : number of g-vectors forming the ng3 stars in the gmax-sphe
      82              :     !     igfft : pointer from the g-sphere (stored as stars) to fft-grid
      83              :     !             and     from fft-grid to g-sphere (stored as stars)
      84              :     !     pgfft : contains the phases of the g-vectors of sph.
      85              :     !     isn   : isn = +1, fft transform for g-space to r-space
      86              :     !             isn = -1, vice versa
      87              :     !     l_rdm : if true, calculate noco gradients in local frame of density matrix
      88              :     !             improves convergence in case of GGA + noco (best with large G_max_xc)
      89              :     !
      90              :     !-------------------------------------------------------------------
      91              :     USE m_grdrsis
      92              :     USE m_polangle
      93              :     USE m_mkgxyz3
      94              :     USE m_types
      95              :     USE m_constants
      96              :     IMPLICIT NONE
      97              : 
      98              :     LOGICAL,INTENT(IN)                    :: dograds
      99              :     INTEGER,INTENT(IN)                    :: jspins
     100              :     LOGICAL,INTENT(IN)                    :: l_noco
     101              :     TYPE(t_stars),INTENT(IN)              :: stars
     102              :     TYPE(t_cell),INTENT(IN)               :: cell
     103              :     COMPLEX,INTENT(IN)                    :: den_pw(:,:)
     104              :     TYPE(t_gradients),INTENT(OUT)         :: grad
     105              :     CLASS(t_xcpot), INTENT(IN),OPTIONAL   :: xcpot
     106              :     REAL,ALLOCATABLE,INTENT(OUT),OPTIONAL :: rho(:,:),rhoim(:,:)
     107              : 
     108              : 
     109              :     INTEGER      :: js,i,idm,ig,ndm,jdm,j
     110              :     REAL         :: rhotot,mmx,mmy,mmz,theta,phi,fd(3),sd(3), rhoCutGrad
     111              :     COMPLEX      :: rho21
     112              :     !     .. Local Arrays ..
     113          351 :     COMPLEX, ALLOCATABLE :: cqpw(:,:),ph_wrk(:)
     114          351 :     REAL,    ALLOCATABLE :: bf3(:)
     115          351 :     REAL,    ALLOCATABLE :: rhd1(:,:,:),rhd2(:,:,:)
     116          351 :     REAL,    ALLOCATABLE :: mx(:),my(:)
     117          351 :     REAL,    ALLOCATABLE :: magmom(:),dmagmom(:,:),ddmagmom(:,:,:)
     118          351 :     REAL,    ALLOCATABLE :: rhodiag(:,:),der(:,:,:),dder(:,:,:,:)
     119          351 :     REAL,    ALLOCATABLE :: sinsqu(:),cossqu(:),sincos(:),rhdd(:,:,:,:)
     120          351 :     COMPLEX, ALLOCATABLE :: exi(:)
     121              : 
     122              :     LOGICAL, PARAMETER :: l_rdm=.true.
     123              : 
     124              :     ! Allocate arrays
     125          351 :     ALLOCATE( bf3(0:griddim-1))
     126          351 :     IF (dograds) THEN
     127          560 :        IF (PRESENT(rho)) ALLOCATE(rho(0:griddim-1,jspins))
     128          284 :        IF (PRESENT(rhoim)) ALLOCATE(rhoim(0:griddim-1,jspins))
     129          852 :        ALLOCATE( ph_wrk(0:griddim-1),rhd1(0:griddim-1,jspins,3))
     130          852 :        ALLOCATE( rhd2(0:griddim-1,jspins,6) )
     131              :      ELSE
     132          134 :         IF (PRESENT(rho)) ALLOCATE(rho(0:griddim-1,jspins))
     133           67 :         IF (PRESENT(rhoim)) ALLOCATE(rhoim(0:griddim-1,jspins))
     134              :      ENDIF
     135          351 :     IF (l_noco)  THEN
     136           96 :        IF (dograds) THEN
     137           67 :           ALLOCATE( mx(0:griddim-1),my(0:griddim-1),magmom(0:griddim-1))
     138              :           IF (l_rdm) THEN
     139          134 :            ALLOCATE( rhodiag(0:griddim-1,jspins),der(0:griddim-1,3,4),dder(0:griddim-1,3,3,4),rhdd(0:griddim-1,2,3,3) )
     140           67 :            ALLOCATE( sinsqu(0:griddim-1),cossqu(0:griddim-1),sincos(0:griddim-1),exi(0:griddim-1) )
     141              :           ELSE
     142              :             ALLOCATE( dmagmom(0:griddim-1,3),ddmagmom(0:griddim-1,3,3) )
     143              :           ENDIF
     144              :        ELSE
     145           29 :           ALLOCATE( mx(0:griddim-1),my(0:griddim-1),magmom(0:griddim-1))
     146              :        ENDIF
     147              :     END IF
     148              : !    if (gmaxz>0.0) write(*,*) "Present in pw_to_grid"
     149          351 :     IF (PRESENT(rho)) THEN
     150              :     !Put den_pw on grid and store into rho(:,1:2)
     151          885 :         DO js=1,jspins
     152          542 :             if (gmaxz>0.0) then
     153            0 :                call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz)
     154              :             else
     155          542 :                call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax)
     156              :             end if
     157          542 :             call fftgrid%perform_fft(forward=.false.)
     158              :             ! TODO: grid is technically still complex right? The REAL cast happens here:
     159              :             !rho(0:,js)=fftgrid%grid
     160     10085863 :             rho(0:,js)   =  REAL(fftgrid%grid)
     161          885 :             IF (PRESENT(rhoim)) rhoim(0:,js) = AIMAG(fftgrid%grid)
     162              :          END DO
     163              : 
     164          343 :        IF (l_noco) THEN
     165              :           !  Get mx,my on real space grid and recalculate rho and magmom
     166           96 :           if (gmaxz>0.0) then
     167            0 :              call fftgrid%putFieldOnGrid(stars,den_pw(:,3),cell,gmax,gmaxz)
     168              :           else
     169           96 :              call fftgrid%putFieldOnGrid(stars,den_pw(:,3),cell,gmax)
     170              :           end if
     171           96 :           call fftgrid%perform_fft(forward=.false.)
     172      2394138 :           mx=real(fftgrid%grid)
     173      2394138 :           my=aimag(fftgrid%grid) ! TODO: There is NO magic minus here. This is so scuffed....
     174              : 
     175      2394042 :           DO i=0,MIN(SIZE(rho,1),size(mx))-1
     176      2393946 :              rhotot= 0.5*( rho(i,1) + rho(i,2) )
     177      2393946 :              magmom(i)= SQRT(  (0.5*(rho(i,1)-rho(i,2)))**2 + mx(i)**2 + my(i)**2 )
     178      3040794 :              IF (l_rdm.AND.dograds) rhodiag(i,1:2) = rho(i,1:2)
     179      2393946 :              rho(i,1)= rhotot+magmom(i)
     180      2393946 :              rho(i,2)= rhotot-magmom(i)
     181      2394042 :              IF (l_rdm.AND.dograds) THEN              ! prepare rotation matrix
     182       323424 :                 mmx = 2.0 * mx(i)
     183       323424 :                 mmy = 2.0 * my(i)
     184       323424 :                 mmz = rhodiag(i,1) - rhodiag(i,2)
     185       323424 :                 CALL pol_angle(mmx,mmy,mmz,theta,phi)
     186       323424 :                 cossqu(i) = cos(0.5*theta)**2
     187       323424 :                 sinsqu(i) = sin(0.5*theta)**2
     188       323424 :                 sincos(i) = 2.0*sin(0.5*theta)*cos(0.5*theta)
     189       323424 :                 exi(i)    = exp(-ImagUnit*phi)
     190              :              ENDIF
     191              :           END DO
     192              :        ENDIF
     193              :     ENDIF
     194          351 :     IF (dograds) THEN
     195          284 :       IF (PRESENT(xcpot)) THEN
     196          284 :          CALL xcpot%alloc_gradients(griddim,jspins,grad)
     197              :       END IF
     198              : 
     199              : 
     200              :     ! in non-collinear calculations the derivatives of |m| are calculated in real space.
     201              : 
     202              : 
     203          284 :        IF (.not.l_noco) THEN
     204              : 
     205              :       ! In collinear calculations all derivatives are calculated in g-spce,
     206              :        ndm = 0
     207          868 :        DO idm = 1,3
     208          651 :          fd=0.0;fd(idm)=1
     209         1572 :          DO js=1,jspins
     210          921 :             if (gmaxz>0.0) then
     211            0 :                call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz,firstderiv=fd)
     212              :             else
     213          921 :                call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,firstderiv=fd)
     214              :             end if
     215          921 :             call fftgrid%perform_fft(forward=.false.)
     216      9179895 :             rhd1(0:,js,idm)=fftgrid%grid
     217              :          END DO
     218          868 :          IF (allocated(grad%laplace).or.allocated(grad%agrt)) THEN
     219              :            !Higher derivatives needed
     220         1953 :            DO jdm = 1,idm
     221         1302 :              sd=0;sd(jdm)=1
     222         1302 :              ndm = ndm + 1
     223         3795 :              DO js=1,jspins
     224         1842 :                 if (gmaxz>0.0) then
     225            0 :                    call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,gmaxz,firstderiv=fd,secondderiv=sd)
     226              :                 else
     227         1842 :                    call fftgrid%putFieldOnGrid(stars,den_pw(:,js),cell,gmax,firstderiv=fd,secondderiv=sd)
     228              :                 end if
     229         1842 :                 call fftgrid%perform_fft(forward=.false.)
     230     18359790 :                 rhd2(0:,js,ndm)=fftgrid%grid
     231              :              END DO
     232              :            END DO ! jdm
     233              :          ENDIF
     234              :        END DO   ! idm
     235              : 
     236              :        ELSE !noco case
     237              : 
     238              :           IF (l_rdm) THEN
     239              : 
     240           67 :              CALL grdrsis( rhodiag(:,1),cell,fftgrid%dimensions,der(:,:,1) )
     241           67 :              CALL grdrsis( rhodiag(:,2),cell,fftgrid%dimensions,der(:,:,2) )
     242           67 :              CALL grdrsis( mx,cell,fftgrid%dimensions,der(:,:,3) )
     243           67 :              CALL grdrsis( my,cell,fftgrid%dimensions,der(:,:,4) )
     244              : 
     245       323491 :              DO i=0,griddim-1 ! project on magnetization axis
     246      1293763 :                 DO idm=1,3
     247       970272 :                    rho21 = der(i,idm,3) + ImagUnit * der(i,idm,4)
     248              :                    rhd1(i,1,idm) = cossqu(i) * der(i,idm,1) + &
     249              :                                    sincos(i) * real( exi(i)*rho21 ) + &
     250       970272 :                                    sinsqu(i) * der(i,idm,2)
     251              :                    rhd1(i,2,idm) = sinsqu(i) * der(i,idm,1) - &
     252              :                                    sincos(i) * real( exi(i)*rho21 ) + &
     253      1293696 :                                    cossqu(i) * der(i,idm,2)
     254              :                 ENDDO
     255              :              ENDDO
     256              :              !Now second derivatives
     257          268 :              DO idm = 1,3
     258          201 :                 CALL grdrsis( der(:,idm,1),cell,fftgrid%dimensions, dder(:,:,idm,1) )
     259          201 :                 CALL grdrsis( der(:,idm,2),cell,fftgrid%dimensions, dder(:,:,idm,2) )
     260          201 :                 CALL grdrsis( der(:,idm,3),cell,fftgrid%dimensions, dder(:,:,idm,3) )
     261          201 :                 CALL grdrsis( der(:,idm,4),cell,fftgrid%dimensions, dder(:,:,idm,4) )
     262       970540 :                 DO i=0,griddim-1 ! project on magnetization axis
     263      3881289 :                   DO jdm=1,3
     264      2910816 :                     rho21 = dder(i,jdm,idm,3) + ImagUnit * dder(i,jdm,idm,4)
     265              :                     rhdd(i,1,jdm,idm) = cossqu(i) * dder(i,jdm,idm,1) + &
     266              :                                     sincos(i) * real( exi(i)*rho21 ) + &
     267      2910816 :                                     sinsqu(i) * dder(i,jdm,idm,2)
     268              :                     rhdd(i,2,jdm,idm) = sinsqu(i) * dder(i,jdm,idm,1) - &
     269              :                                     sincos(i) * real( exi(i)*rho21 ) + &
     270      3881088 :                                     cossqu(i) * dder(i,jdm,idm,2)
     271              :                   ENDDO
     272              :                 ENDDO
     273              :               ENDDO
     274          201 :               DO j=1,2
     275       647049 :                 DO i=0,griddim-1
     276       646848 :                   rhd2(i,j,1) = rhdd(i,j,1,1)
     277       646848 :                   rhd2(i,j,2) = 0.5*(rhdd(i,j,1,2)+rhdd(i,j,2,1)) ! xy - averaging should be unneccessary
     278       646848 :                   rhd2(i,j,3) = rhdd(i,j,2,2)
     279       646848 :                   rhd2(i,j,4) = 0.5*(rhdd(i,j,1,3)+rhdd(i,j,3,1)) ! zx
     280       646848 :                   rhd2(i,j,5) = 0.5*(rhdd(i,j,2,3)+rhdd(i,j,3,2)) ! yz
     281       646982 :                   rhd2(i,j,6) = rhdd(i,j,3,3)
     282              :                 ENDDO
     283              :               ENDDO
     284           67 :               DEALLOCATE (rhodiag,der,rhdd,dder,sinsqu,cossqu,sincos,exi)
     285              : 
     286              :           ELSE
     287              : 
     288              :              CALL grdrsis(magmom,cell,fftgrid%dimensions,dmagmom )
     289              : 
     290              :              DO i=0,griddim-1
     291              :                 DO idm=1,3
     292              :                    rhotot= rhd1(i,1,idm)/2.+rhd1(i,2,idm)/2.
     293              :                    rhd1(i,1,idm)= rhotot+dmagmom(i,idm)
     294              :                    rhd1(i,2,idm)= rhotot-dmagmom(i,idm)
     295              :                 END DO
     296              :              END DO
     297              :              !Second derivatives
     298              :              DO idm = 1,3
     299              :                 CALL grdrsis(dmagmom(:,idm),cell,fftgrid%dimensions,ddmagmom(:,:,idm) )
     300              :              END DO
     301              :              ndm= 0
     302              :              DO idm = 1,3
     303              :                 DO jdm = 1,idm
     304              :                    ndm = ndm + 1
     305              :                    DO i=0,griddim-1
     306              :                       rhotot= rhd2(i,1,ndm)/2.+rhd2(i,2,ndm)/2.
     307              :                       rhd2(i,1,ndm)= rhotot + ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2.
     308              :                       rhd2(i,2,ndm)= rhotot - ( ddmagmom(i,jdm,idm) + ddmagmom(i,idm,jdm) )/2.
     309              :                    END DO
     310              :                 ENDDO !jdm
     311              :              ENDDO   !idm
     312              :              DEALLOCATE(dmagmom,ddmagmom)
     313              :           END IF
     314              : 
     315              :        END IF
     316              : 
     317              : 
     318              : 
     319              : 
     320              :        !
     321              :        !     calculate the quantities such as abs(grad(rho)),.. used in
     322              :        !     evaluating the gradient contributions to potential and energy.
     323              :        !
     324          284 :        IF (PRESENT(rho)) THEN
     325              :           CALL mkgxyz3 (rho,rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
     326          276 :                rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),0,grad)
     327              :        ELSE
     328              :           !Dummy rho (only possible if grad is used for libxc mode)
     329              :           !CALL mkgxyz3 (RESHAPE((/0.0/),(/1,1/)),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
     330              :           !     rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),grad)
     331              : 
     332            8 :           IF (dograds.and.(.not.PRESENT(xcpot))) THEN
     333            0 :              ALLOCATE(grad%gr(3,griddim,1))
     334              :           END IF
     335              : 
     336              :           CALL mkgxyz3 (0*rhd1(0:,:,1),rhd1(0:,:,1),rhd1(0:,:,2),rhd1(0:,:,3),&
     337        98236 :                rhd2(0:,:,1),rhd2(0:,:,3),rhd2(0:,:,6), rhd2(0:,:,5),rhd2(0:,:,4),rhd2(0:,:,2),0,grad)
     338              :        END IF
     339              : 
     340              :     ENDIF
     341          351 :     rhoCutGrad = 3.0e-6
     342          351 :     IF (PRESENT(rho)) THEN
     343     10086206 :        WHERE(ABS(rho) < d_15) rho = d_15
     344          343 :        IF(PRESENT(xcpot)) THEN
     345              :           SELECT TYPE(xcpot)
     346              :           TYPE IS (t_xcpot_inbuild)
     347          332 :              IF (dograds) THEN
     348              :                 ! In the following, gradients at places with very small densities are set to zero. This is done for stability reasons.
     349      2770113 :                 DO i = 1, griddim
     350      2769841 :                    IF(rho(i-1,1).LT.rhoCutGrad) THEN
     351         7060 :                       grad%agru(i) = 0.0
     352         7060 :                       grad%g2ru(i) = 0.0
     353         7060 :                       grad%gggru(i) = 0.0
     354              : 
     355         7060 :                       grad%agrt(i) = grad%agrd(i)
     356         7060 :                       grad%g2rt(i) = grad%g2rd(i)
     357         7060 :                       grad%gggrt(i) = grad%gggrd(i)
     358              :                    END IF
     359              : 
     360      2769841 :                    IF(rho(i-1,jspins).LT.rhoCutGrad) THEN
     361         7256 :                       grad%agrd(i) = 0.0
     362         7256 :                       grad%g2rd(i) = 0.0
     363         7256 :                       grad%gggrd(i) = 0.0
     364              : 
     365         7256 :                       grad%agrt(i) = grad%agru(i)
     366         7256 :                       grad%g2rt(i) = grad%g2ru(i)
     367         7256 :                       grad%gggrt(i) = grad%gggru(i)
     368              :                    END IF
     369              : 
     370      2770113 :                    IF((rho(i-1,1).LT.rhoCutGrad).AND.(rho(i-1,jspins).LT.rhoCutGrad)) THEN
     371         7060 :                       grad%gzgr(i) = 0.0
     372              :                    END IF
     373              :                 END DO
     374              :              END IF
     375              :           END SELECT
     376              :        END IF
     377              :     ENDIF
     378          351 :     IF (PRESENT(rhoim)) THEN
     379            0 :        WHERE(ABS(rhoim) < d_15) rhoim = d_15
     380              :     ENDIF
     381              : 
     382          351 :   END SUBROUTINE pw_to_grid
     383              : 
     384              : 
     385         1362 :   SUBROUTINE pw_from_grid(stars,v_in,v_out_pw,v_out_pw_w)
     386              :     USE m_convol
     387              :     USE m_types
     388              :     IMPLICIT NONE
     389              :     TYPE(t_stars),INTENT(IN)      :: stars
     390              :     REAL,INTENT(INOUT)            :: v_in(0:,:)
     391              :     COMPLEX,INTENT(INOUT)         :: v_out_pw(:,:)
     392              :     COMPLEX,INTENT(INOUT),OPTIONAL:: v_out_pw_w(:,:)
     393              : 
     394              : 
     395              :     INTEGER              :: js,k,i
     396         1362 :     REAL,ALLOCATABLE     :: bf3(:),vcon(:)
     397         1362 :     COMPLEX, ALLOCATABLE :: fg3(:)
     398         1362 :     if (present(v_out_pw_w)) ALLOCATE( bf3(size(stars%ufft)),vcon(size(stars%ufft)))
     399         1362 :     ALLOCATE ( fg3(stars%ng3) )
     400         3333 :     DO js = 1,SIZE(v_in,2)
     401     36536920 :        fftgrid%grid=v_in(0:,js)
     402         1971 :        call fftgrid%perform_fft(forward=.true.)
     403         1971 :        if (gmaxz>0.0) then
     404            0 :           call fftgrid%takeFieldFromGrid(stars,fg3,gmax,gmaxz)
     405              :        else
     406         1971 :           call fftgrid%takeFieldFromGrid(stars,fg3,gmax)
     407              :        END IF
     408      4346987 :        v_out_pw(:,js) = v_out_pw(:,js) + fg3(:)
     409              : 
     410              :        !----> add to warped coulomb potential
     411         3333 :        IF (present(v_out_pw_w)) THEN
     412         1415 :          if (size(fftgrid%grid)==size(stars%ufft)) THEN
     413     16338444 :             fftgrid%grid=v_in(0:,js)*stars%ufft
     414          291 :             call fftgrid%perform_fft(forward=.true.)
     415          291 :             if (gmaxz>0.0) then
     416            0 :                call fftgrid%takeFieldFromGrid(stars,fg3,gmax,gmaxz)
     417              :             else
     418          291 :                call fftgrid%takeFieldFromGrid(stars,fg3,gmax)
     419              :             end if
     420      1545476 :             fg3 = fg3*stars%nstr
     421              :          else
     422         1124 :           call convol(stars,fg3)
     423              :          ENDIF
     424      3148886 :           v_out_pw_w(:,js) = v_out_pw_w(:,js) + fg3
     425              :        ENDIF
     426              :     END DO
     427         1362 :   END SUBROUTINE pw_from_grid
     428              : 
     429          339 :   SUBROUTINE finish_pw_grid()
     430              :     IMPLICIT NONE
     431              :     !IF (ALLOCATED(igxc_fft)) DEALLOCATE(igxc_fft,gxc_fft)
     432          339 :   END SUBROUTINE finish_pw_grid
     433              : 
     434          337 : END MODULE m_pw_tofrom_grid
        

Generated by: LCOV version 2.0-1