LCOV - code coverage report
Current view: top level - vgen - vvacxcg.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 150 240 62.5 %
Date: 2019-09-08 04:53:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : MODULE m_vvacxcg
       2             :   use m_juDFT
       3             :   !-----------------------------------------------------------------------
       4             :   !     calculates 2-d star function coefficients of exchange-correlation*
       5             :   !     potential in the vacuum regions and adds them to the corresponding
       6             :   !     coeffs of the coulomb potential            c.l.fu, r.podloucky   *
       7             :   !     for the gradient contribution.   t.a. 1996
       8             :   !-----------------------------------------------------------------------
       9             : CONTAINS
      10           5 :   SUBROUTINE vvacxcg(ifftd2,stars,vacuum,noco,oneD,cell,xcpot,input,obsolete,den, vxc,exc)
      11             : 
      12             :     !-----------------------------------------------------------------------
      13             :     !     instead of vvacxcor.f: the different exchange-correlation
      14             :     !     potentials defined through the key icorr are called through
      15             :     !     the driver subroutine vxcallg.f, subroutines vectorized
      16             :     !     in case of total = .true. calculates the ex-corr. energy
      17             :     !     density through the driver subroutine excallg.f
      18             :     !     ** r.pentcheva 08.05.96
      19             :     !-----------------------------------------------------------------------
      20             :     use m_constants
      21             :     USE m_grdrsvac
      22             :     USE m_grdchlh
      23             :     USE m_mkgz
      24             :     USE m_mkgxyz3
      25             :     USE m_od_mkgxyz3
      26             :     USE m_od_mkgz
      27             :     USE m_fft2d
      28             :     USE m_types
      29             :     USE m_types_xcpot_libxc
      30             :     IMPLICIT NONE
      31             :     CLASS(t_xcpot),INTENT(IN)    :: xcpot
      32             :     TYPE(t_oneD),INTENT(IN)      :: oneD
      33             :     TYPE(t_obsolete),INTENT(IN)  :: obsolete
      34             :     TYPE(t_input),INTENT(IN)     :: input
      35             :     TYPE(t_vacuum),INTENT(IN)    :: vacuum
      36             :     TYPE(t_noco),INTENT(IN)      :: noco
      37             :     TYPE(t_stars),INTENT(IN)     :: stars
      38             :     TYPE(t_cell),INTENT(IN)      :: cell
      39             :     TYPE(t_potden),INTENT(IN)    :: den
      40             :     TYPE(t_potden),INTENT(INOUT) :: vxc
      41             :     TYPE(t_potden),INTENT(INOUT) :: exc
      42             :     !     ..
      43             :     !     .. Scalar Arguments ..
      44             :     INTEGER, INTENT (IN) :: ifftd2
      45             :    
      46             :     !     ..
      47             :     !     .. Local Scalars ..
      48             :     INTEGER :: js,nt,i,iq,irec2,nmz0,nmzdiff,ivac,ip
      49             :     REAL    :: rhti,zro,fgz,rhmnv,d_15,bmat1(3,3),rd
      50             :     !     ..
      51             :     !     .. Local Arrays ..
      52          10 :     REAL, ALLOCATABLE :: af2(:,:),bf2(:)
      53          15 :     REAL, ALLOCATABLE :: rhdx(:,:),rhdy(:,:),rhdz(:,:)
      54          15 :     REAL, ALLOCATABLE :: rhdxx(:,:),rhdyy(:,:),rhtdz(:,:),rhtdzz(:,:)
      55          15 :     REAL, ALLOCATABLE :: rhdzz(:,:),rhdyz(:,:),rhdzx(:,:),rhdxy(:,:)
      56          10 :     REAL, ALLOCATABLE :: v_x(:,:),v_xc(:,:),e_xc(:),vxz(:,:),vxcz(:,:)
      57           5 :     REAL, ALLOCATABLE :: rxydzr(:),rxydzi(:)
      58          10 :     REAL, ALLOCATABLE :: rxydzzr(:),rxydzzi(:),rhtxyr(:),rhtxyi(:)
      59          15 :     REAL, ALLOCATABLE :: rhtxc(:,:),rhtz(:,:),dummy(:)
      60          20 :     COMPLEX, ALLOCATABLE :: fgxy(:),rxydz(:,:,:),rxydzz(:,:,:),cqpw(:)
      61             : 
      62           5 :     TYPE(t_gradients)::grad
      63             :     !     ..
      64             :     !     for the noco-case only 
      65             :     REAL :: chdens
      66          15 :     REAL, ALLOCATABLE :: magmom(:,:), dxmagmom(:),ddxmagmom(:,:)
      67          10 :     REAL, ALLOCATABLE :: dymagmom(:),ddymagmom(:,:), dzmagmom(:,:),ddzmagmom(:,:)
      68           5 :     REAL, ALLOCATABLE :: mx(:),my(:) 
      69           5 :     REAL, ALLOCATABLE :: af2noco(:,:,:)
      70             : 
      71             :     !     .. unused input (needed for other noco GGA-implementations) ..
      72             : 
      73             : 
      74             :     SELECT TYPE(xcpot)
      75             :     TYPE IS (t_xcpot_libxc)
      76           0 :        CALL judft_error("libxc GGA functionals not implemented in film setups")
      77             :     END SELECT
      78             :    
      79             : 
      80           5 :     d_15     = 1.e-15
      81           5 :     zro      = 0.0
      82           5 :     nt       = ifftd2
      83           5 :     if(oneD%odi%d1)then
      84           0 :        bmat1(:,:) = 0.
      85           0 :        bmat1(1,1) = cell%bmat(3,3)
      86           0 :        bmat1(2,2) = 1.
      87             :     else
      88          20 :        bmat1(:,:) = cell%bmat(:,:)
      89             :     endif
      90             : 
      91           5 :     WRITE (6,'(/'' ifftd2,vacuum%nmz='',2i7)') ifftd2,vacuum%nmz
      92           5 :     WRITE(6,'('' 9990nmzxy='',2i5)') vacuum%nmzxy
      93             : 
      94           5 :     ALLOCATE ( rxydz(vacuum%nmzxy,stars%ng2-1,input%jspins),rxydzz(vacuum%nmzxyd,stars%ng2-1,input%jspins) )
      95           5 :     ALLOCATE ( rhtdz(vacuum%nmzd,input%jspins),rhtdzz(vacuum%nmzd,input%jspins) )
      96             : 
      97          10 :     DO ivac=1,vacuum%nvac 
      98             : 
      99             :        ! the charge density in vacuum is expanded in 2-dim stars on a mesh
     100             :        ! in z-direction. the g||.ne.zero-components expand from 1 to nmzxy
     101             :        ! the g||.eq.zero-components expand from 1 to nmz
     102             :        ! first we calculate vxc in the warping region 
     103             : 
     104           5 :        IF (noco%l_noco) THEN
     105             : 
     106           0 :           ALLOCATE ( magmom(0:ifftd2-1,vacuum%nmzxy) ) 
     107           0 :           ALLOCATE ( dzmagmom(0:ifftd2-1,vacuum%nmzxy) ) 
     108           0 :           ALLOCATE ( ddzmagmom(0:ifftd2-1,vacuum%nmzxy) ) 
     109           0 :           ALLOCATE ( mx(0:ifftd2-1),my(0:ifftd2-1) )
     110           0 :           ALLOCATE ( dxmagmom(0:ifftd2-1),dymagmom(0:ifftd2-1) )
     111           0 :           ALLOCATE ( ddxmagmom(0:ifftd2-1,2),ddymagmom(0:ifftd2-1,2) ) 
     112           0 :           ALLOCATE ( af2noco(0:ifftd2-1,vacuum%nmzxy,input%jspins),bf2(0:ifftd2-1) ) 
     113             : 
     114             :           ! Transform charge and magnetization to real-space.
     115             :           ! In the collinear case that is done later within
     116             :           ! another loop over the vacuum-layers in order to 
     117             :           ! save memory.
     118             : 
     119           0 :           DO ip=1,vacuum%nmzxy
     120             : 
     121           0 :              DO js=1,input%jspins
     122             :                 CALL fft2d(stars, af2noco(0,ip,js),bf2, den%vacz(ip,ivac,js),0.,&
     123           0 :                      den%vacxy(ip,1,ivac,js), vacuum%nmzxy,+1)
     124             :              END DO
     125             :              CALL fft2d(stars, mx,my, den%vacz(ip,ivac,3),den%vacz(ip,ivac,4), &
     126           0 :                   den%vacxy(ip,1,ivac,3), vacuum%nmzxy,+1)
     127             : 
     128           0 :              DO i=0,9*stars%mx1*stars%mx2-1
     129           0 :                 magmom(i,ip)= mx(i)**2 + my(i)**2 + ((af2noco(i,ip,1)-af2noco(i,ip,2))/2.)**2
     130           0 :                 magmom(i,ip)= SQRT(magmom(i,ip))
     131           0 :                 chdens= af2noco(i,ip,1)/2.+af2noco(i,ip,2)/2.
     132           0 :                 af2noco(i,ip,1)= chdens + magmom(i,ip)
     133           0 :                 af2noco(i,ip,2)= chdens - magmom(i,ip)
     134             :              END DO
     135             : 
     136             :           END DO ! ip=1,vacuum%nmzxy 
     137           0 :           DEALLOCATE ( bf2 )
     138             :        END IF ! noco%l_noco 
     139             : 
     140             :        !      ENDDO    ! ivac
     141             :        !      DO ivac = 1,nvac
     142             : 
     143           5 :        IF (xcpot%needs_grad()) THEN
     144          13 :           DO js=1,input%jspins
     145             :              !
     146             :              ! calculate first (rhtdz) & second (rhtdzz) derivative of den%vacz(1:nmz)
     147             :              !
     148           8 :              ALLOCATE ( dummy(vacuum%nmz) )
     149             :              CALL grdchlh(0,1,vacuum%nmz,vacuum%delz,dummy,den%vacz(1,ivac,js),obsolete%ndvgrd,&
     150           8 :                   rhtdz(1,js),rhtdzz(1,js))
     151           8 :              DEALLOCATE ( dummy )
     152           8 :              ALLOCATE ( rhtxyr(vacuum%nmzxy), rhtxyi(vacuum%nmzxy),dummy(vacuum%nmzxy) )
     153           8 :              ALLOCATE ( rxydzr(vacuum%nmzxy), rxydzi(vacuum%nmzxy) )
     154           8 :              ALLOCATE ( rxydzzr(vacuum%nmzxy),rxydzzi(vacuum%nmzxy) )
     155             : 
     156        2080 :              DO iq = 1, stars%ng2-1
     157             :                 !
     158             :                 ! calculate first (rxydz) & second (rxydzz) derivative of den%vacxy:
     159             :                 !
     160      416472 :                 DO ip=1,vacuum%nmzxy
     161      209272 :                    rhtxyr(ip)=den%vacxy(ip,iq,ivac,js)
     162             :                 ENDDO
     163        2072 :                 CALL grdchlh(0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyr,obsolete%ndvgrd, rxydzr,rxydzzr) 
     164             : 
     165      209272 :                 DO ip=1,vacuum%nmzxy
     166      209272 :                    rhtxyi(ip)=aimag(den%vacxy(ip,iq,ivac,js))
     167             :                 ENDDO
     168        2072 :                 CALL grdchlh(0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyi,obsolete%ndvgrd, rxydzi,rxydzzi)
     169             : 
     170      209280 :                 DO ip=1,vacuum%nmzxy
     171      207200 :                    rxydz(ip,iq,js)=cmplx(rxydzr(ip),rxydzi(ip))
     172      209272 :                    rxydzz(ip,iq,js)=cmplx(rxydzzr(ip),rxydzzi(ip))
     173             :                 ENDDO
     174             : 
     175             :              ENDDO ! loop over 2D stars (iq)
     176             : 
     177           8 :              DEALLOCATE ( rhtxyr,rhtxyi,rxydzr,rxydzi,rxydzzr,rxydzzi )
     178          13 :              DEALLOCATE ( dummy )
     179             : 
     180             :           ENDDO ! input%jspins
     181             : 
     182           5 :           IF (noco%l_noco) THEN 
     183             :              !  calculate  dzmagmom = d magmom / d z  and ddzmagmom= d dmagmom / d z 
     184             : 
     185           0 :              ALLOCATE ( rhtxyr(vacuum%nmzxy),dummy(vacuum%nmzxy)   )
     186           0 :              ALLOCATE ( rxydzr(vacuum%nmzxy),rxydzzr(vacuum%nmzxy) )
     187           0 :              DO i=0,9*stars%mx1*stars%mx2-1 
     188           0 :                 DO ip=1,vacuum%nmzxy
     189           0 :                    rhtxyr(ip)=magmom(i,ip)
     190             :                 ENDDO
     191           0 :                 CALL grdchlh(0,1,vacuum%nmzxy,vacuum%delz,dummy,rhtxyr,obsolete%ndvgrd, rxydzr,rxydzzr)
     192           0 :                 DO ip=1,vacuum%nmzxy
     193           0 :                    dzmagmom(i,ip)= rxydzr(ip)
     194           0 :                    ddzmagmom(i,ip)= rxydzzr(ip)
     195             :                 ENDDO
     196             :              END DO
     197           0 :              DEALLOCATE ( rhtxyr,rxydzr,rxydzzr,dummy )
     198             :           END IF ! noco%l_noco 
     199             : 
     200             :        ENDIF   ! xcpot%igrd.GT.0
     201             : 
     202             :        !       WRITE(6,'('' 9990nmzxy='',2i5)') nmzxy
     203           5 :        ALLOCATE ( rhdx(0:ifftd2-1,input%jspins),rhdy(0:ifftd2-1,input%jspins) )
     204           5 :        ALLOCATE ( rhdz(0:ifftd2-1,input%jspins),rhdxx(0:ifftd2-1,input%jspins) )
     205           5 :        ALLOCATE ( rhdyy(0:ifftd2-1,input%jspins),rhdzz(0:ifftd2-1,input%jspins) )
     206           5 :        ALLOCATE ( rhdyz(0:ifftd2-1,input%jspins),rhdzx(0:ifftd2-1,input%jspins) )
     207           5 :        ALLOCATE ( rhdxy(0:ifftd2-1,input%jspins))
     208             : 
     209           5 :        CALL xcpot%alloc_gradients(ifftd2,input%jspins,grad)
     210           5 :        ALLOCATE ( v_xc(0:ifftd2-1,input%jspins),e_xc(0:ifftd2-1) )
     211           5 :        ALLOCATE ( cqpw(stars%ng2-1),af2(0:ifftd2-1,input%jspins) )
     212           5 :        ALLOCATE ( fgxy(stars%ng2-1),bf2(0:ifftd2-1) )
     213             : 
     214           5 :        ALLOCATE( v_x(0:ifftd2-1,input%jspins) )
     215             : 
     216             :  
     217           5 :        rd = 0.0
     218          13 :        af2=0.0
     219         505 :        DO ip = 1,vacuum%nmzxy
     220             :           ! loop over warping region
     221             : 
     222         500 :           IF (.not. noco%l_noco) THEN
     223             :              ! Transform charge and magnetization to real-space.
     224             : 
     225        1300 :              DO js=1,input%jspins
     226             :                 CALL fft2d(stars, af2(0,js),bf2, den%vacz(ip,ivac,js),0.,&
     227        1300 :                      den%vacxy(ip,1,ivac,js), vacuum%nmzxyd,+1)
     228             :              END DO
     229             : 
     230             :           ELSE
     231             : 
     232           0 :              DO i=0,9*stars%mx1*stars%mx2-1
     233           0 :                 af2(i,1)= af2noco(i,ip,1) 
     234           0 :                 af2(i,2)= af2noco(i,ip,2)
     235             :              END DO
     236             : 
     237             :           END IF
     238             : 
     239         500 :           IF (xcpot%needs_grad()) THEN 
     240             :              ! calculate derivatives with respect to x,y in g-space 
     241             :              ! and transform them to real-space.  
     242             : 
     243        1300 :              DO js = 1,input%jspins
     244             : 
     245      208000 :                 DO iq=1,stars%ng2-1
     246      208000 :                    cqpw(iq)=ImagUnit*den%vacxy(ip,iq,ivac,js)
     247             :                 ENDDO
     248             : 
     249         800 :                 rhti = 0.0                    ! d(rho)/atoms%dx is obtained by a FFT of i*gx*den%vacxy
     250             :                 ! (den%vacz is set to zero and gx is included in 
     251             :                 !    dn/atoms =  FFT(0,i*gx*den%vacxy)
     252             : 
     253             :               
     254             : 
     255         800 :                 CALL fft2d(stars, rhdx(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfx)
     256             :                 !TODO    &                 pgft2x) 
     257             : 
     258         800 :                 rhti = 0.0
     259             :                 CALL fft2d(    &               ! dn/dy =  FFT(0,i*gy*den%vacxy)&
     260         800 :                       stars, rhdy(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy)
     261             : 
     262         800 :                 rhti = 0.0
     263             :                 CALL fft2d(     &              ! dn/dz = FFT(rhtdz,rxydz)&
     264         800 :                         stars, rhdz(0,js),bf2, rhtdz(ip,js),rhti,rxydz(ip,1,js), vacuum%nmzxyd,+1) 
     265             : 
     266      208000 :                 DO iq=1,stars%ng2-1
     267      208000 :                    cqpw(iq)=-den%vacxy(ip,iq,ivac,js)
     268             :                 ENDDO
     269             : 
     270         800 :                 rhti = 0.0
     271             :                 CALL fft2d(      &          ! d2n/dx2 = FFT(0,-gx^2*den%vacxy)&
     272         800 :                        stars, rhdxx(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfx*stars%ft2_gfx)
     273             : 
     274         800 :                 rhti = 0.0
     275             :                 CALL fft2d(       &          ! d2n/dy2 = FFT(0,-gy^2*den%vacxy)&
     276         800 :                       stars, rhdyy(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy*stars%ft2_gfy)
     277             : 
     278         800 :                 rhti = 0.0
     279             :                 CALL fft2d(        &         ! d2n/dz2 = FFT(rhtdzz,rxydzz)&
     280         800 :                        stars, rhdzz(0,js),bf2, rhtdzz(ip,js),rhti,rxydzz(ip,1,js), vacuum%nmzxyd,+1)
     281             : 
     282             : 
     283      208000 :                 DO iq=1,stars%ng2-1
     284      208000 :                    cqpw(iq)=ImagUnit*rxydz(ip,iq,js)
     285             :                 ENDDO
     286             : 
     287         800 :                 rhti = 0.0
     288             :                 CALL fft2d(         &         ! d2n/dyz = FFT(0,i*gy*rxydz)&
     289         800 :                        stars, rhdyz(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy)
     290             : 
     291         800 :                 rhti = 0.0
     292             :                 CALL fft2d(          &        ! d2n/dzx = FFT(0,i*gx*rxydz)&
     293         800 :                        stars, rhdzx(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfx)
     294             : 
     295      208000 :                 DO iq=1,stars%ng2-1
     296      208000 :                    cqpw(iq)=-den%vacxy(ip,iq,ivac,js)
     297             :                 ENDDO
     298             : 
     299         800 :                 rhti = 0.0
     300             :                 CALL fft2d(           &    ! d2n/dxy = FFT(0,-gx*gy*den%vacxy)&
     301        1300 :                       stars, rhdxy(0,js),bf2, zro,rhti,cqpw, 1,+1,stars%ft2_gfy*stars%ft2_gfx)
     302             : 
     303             :              END DO ! js=1,input%jspins
     304             : 
     305             : 
     306         500 :              IF (noco%l_noco) THEN
     307             :                 ! ! In non-collinear calculations the derivatives of |m| are calculated
     308             :                 ! ! in real-space. The derivatives of the charge density, that are 
     309             :                 ! ! already calculated in g-space, will be used. 
     310             : 
     311           0 :                 CALL grdrsvac(magmom(0,ip),bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, dxmagmom,dymagmom) 
     312           0 :                 DO i=0,9*stars%mx1*stars%mx2-1 
     313           0 :                    chdens= rhdx(i,1)/2.+rhdx(i,2)/2.
     314           0 :                    rhdx(i,1)= chdens + dxmagmom(i) 
     315           0 :                    rhdx(i,2)= chdens - dxmagmom(i) 
     316           0 :                    chdens= rhdy(i,1)/2.+rhdy(i,2)/2.
     317           0 :                    rhdy(i,1)= chdens + dymagmom(i) 
     318           0 :                    rhdy(i,2)= chdens - dymagmom(i) 
     319           0 :                    chdens= rhdz(i,1)/2.+rhdz(i,2)/2.
     320           0 :                    rhdz(i,1)= chdens + dzmagmom(i,ip) 
     321           0 :                    rhdz(i,2)= chdens - dzmagmom(i,ip) 
     322             :                 END DO
     323             : 
     324             :                 CALL grdrsvac(dxmagmom,bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, &
     325           0 :                      ddxmagmom(0,1),ddymagmom(0,1))
     326             :                 CALL grdrsvac(&
     327           0 :                      dymagmom,bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd,ddxmagmom(0,2),ddymagmom(0,2))
     328           0 :                 DO i=0,9*stars%mx1*stars%mx2-1
     329           0 :                    chdens= rhdxx(i,1)/2.+rhdxx(i,2)/2. 
     330           0 :                    rhdxx(i,1)= chdens + ddxmagmom(i,1) 
     331           0 :                    rhdxx(i,2)= chdens - ddxmagmom(i,1) 
     332           0 :                    chdens= rhdyy(i,1)/2.+rhdyy(i,2)/2. 
     333           0 :                    rhdyy(i,1)= chdens + ddymagmom(i,2) 
     334           0 :                    rhdyy(i,2)= chdens - ddymagmom(i,2) 
     335           0 :                    chdens= rhdxy(i,1)/2.+rhdxy(i,2)/2. 
     336           0 :                    rhdxy(i,1)= chdens + ( ddxmagmom(i,2) + ddymagmom(i,1) )/2.
     337           0 :                    rhdxy(i,2)= chdens - ( ddxmagmom(i,2) + ddymagmom(i,1) )/2.
     338             :                 END DO
     339             :                 CALL grdrsvac(dzmagmom(0,ip),bmat1,3*stars%mx1,3*stars%mx2,obsolete%ndvgrd, &
     340           0 :                      ddxmagmom(0,1),ddymagmom(0,1))
     341           0 :                 DO i=0,9*stars%mx1*stars%mx2-1
     342           0 :                    chdens= rhdzx(i,1)/2.+rhdzx(i,2)/2. 
     343           0 :                    rhdzx(i,1)= chdens + ddxmagmom(i,1) 
     344           0 :                    rhdzx(i,2)= chdens - ddxmagmom(i,1) 
     345           0 :                    chdens= rhdyz(i,1)/2.+rhdyz(i,2)/2. 
     346           0 :                    rhdyz(i,1)= chdens + ddymagmom(i,1) 
     347           0 :                    rhdyz(i,2)= chdens - ddymagmom(i,1) 
     348           0 :                    chdens= rhdzz(i,1)/2.+rhdzz(i,2)/2. 
     349           0 :                    rhdzz(i,1)= chdens + ddzmagmom(i,ip) 
     350           0 :                    rhdzz(i,2)= chdens - ddzmagmom(i,ip) 
     351             :                 END DO
     352             : 
     353             :              END IF ! noco%l_noco  
     354             : 
     355             :           END IF ! vxc_is_gga 
     356             :           !
     357             :           ! set minimal value of af2 to 1.0e-13
     358             :           !
     359             : 
     360         500 :           af2=max(af2,10e-13)
     361             : 
     362             : 
     363             :           !   calculate the quantities such as abs(grad(rho)),.. used in
     364             :           !c  evaluating the gradient contributions to potential and energy.
     365             : 
     366         500 :           rd = cell%z1 + vacuum%delz*(ip-1)
     367             : 
     368         500 :           if(oneD%odi%d1)then
     369             : !!$             CALL od_mkgxyz3(&
     370             : !!$                  &           ifftd2,input%jspins,ifftd2,input%jspins,&
     371             : !!$                  &           af2,rd,rhdx,rhdy,rhdz,rhdxx,rhdyy,&
     372             : !!$                  &           rhdzz,rhdyz,rhdzx,rhdxy,&
     373             : !!$                  &           agr,agru,agrd,g2r,g2ru,g2rd,&
     374             : !!$                  &           gggr,gggru,gggrd,gzgr)
     375           0 :              CALL judft_error("OneD not implemented")
     376             :           ELSE
     377         500 :              CALL mkgxyz3(af2,rhdx,rhdy, rhdz,rhdxx,rhdyy,rhdzz,rhdyz,rhdzx,rhdxy, grad)
     378             :           endif
     379             : 
     380             :           !     rhmnv: rho_minimum_vacuum.
     381             : 
     382         500 :           rhmnv=10.e+10
     383        1300 :           DO js=1,input%jspins
     384      280100 :              DO i=0,stars%kimax2
     385      278800 :                 af2(i,js)=max(af2(i,js),d_15)
     386      279600 :                 rhmnv=min(rhmnv,af2(i,js))
     387             :              ENDDO
     388             :           ENDDO
     389             : 
     390         500 :           IF (rhmnv.LT.obsolete%chng) THEN
     391           0 :              WRITE(6,'(/'' rhmn.lt.obsolete%chng. rhmn,obsolete%chng='',2d9.2)') rhmnv,obsolete%chng
     392             :              !             CALL juDFT_error("vvacxcg: rhmn.lt.chng",calledby="vvacxcg")
     393             :           ENDIF
     394             : 
     395             :           !         calculate the exchange-correlation potential in  real space
     396             :           !
     397         500 :           CALL xcpot%get_vxc(input%jspins,af2,v_xc,v_x,grad)
     398             :        
     399             : 
     400        1300 :           DO js = 1,input%jspins
     401             :              !
     402             :              !           ----> 2-d back fft to g space
     403             :              !
     404      765800 :              bf2=0.0
     405         800 :              CALL fft2d(stars, v_xc(0,js),bf2, fgz,rhti,fgxy, 1,-1)
     406             : 
     407             :              !            ----> and add vxc to coulomb potential
     408             :              !                  the g||.eq.zero component is added to vxc%vacz
     409             :              !
     410         800 :              vxc%vacz(ip,ivac,js) = fgz + vxc%vacz(ip,ivac,js)
     411             :              !
     412             :              !            the g||.ne.zero components are added to vxc%vacxy
     413             :              !
     414      208500 :              DO irec2 = 1,stars%ng2-1
     415      208000 :                 vxc%vacxy(ip,irec2,ivac,js)=vxc%vacxy(ip,irec2,ivac,js)+fgxy(irec2)
     416             :              ENDDO
     417             : 
     418             :           END DO
     419             : 
     420             : 
     421             :           !         calculate the exchange-correlation energy density in  real space
     422             :           !
     423         505 :           IF (ALLOCATED(exc%vacz)) THEN
     424             : 
     425         500 :              CALL xcpot%get_exc(input%jspins,af2,e_xc,grad, mt_call=.False.)
     426             :    
     427             :              !           ----> 2-d back fft to g space
     428             :              !
     429      512600 :              bf2=0.0
     430         500 :              CALL fft2d(stars, e_xc,bf2, exc%vacz(ip,ivac,1),rhti,exc%vacxy(ip,1,ivac,1), vacuum%nmzxyd,-1)
     431             : 
     432             :           ENDIF
     433             : 
     434             :        END DO ! ip=1,vacuum%nmzxy 
     435           5 :        DEALLOCATE ( rhdx,rhdy,rhdz,rhdxx,rhdyy,rhdzz )
     436           5 :        DEALLOCATE ( cqpw,fgxy,     rhdyz,rhdzx,rhdxy )
     437             : 
     438           5 :        IF (noco%l_noco) THEN
     439           0 :           DEALLOCATE ( dzmagmom,ddzmagmom,dxmagmom,af2noco )
     440           0 :           DEALLOCATE ( dymagmom,ddxmagmom,ddymagmom )
     441             :        END IF
     442             : 
     443             :        ! now treat the non-warping region 
     444             : 
     445           5 :        nmzdiff = vacuum%nmz - vacuum%nmzxy
     446             :        !       WRITE(6,'(/'' 9992excz''/(8f15.7))') (excz(ip,1),ip=1,nmz)
     447           5 :        WRITE(6,'(/'' 9992nmzdiff='',i5)') nmzdiff
     448             : 
     449             :        ! The non-warping region runs from nmzxy+1 to nmz.
     450             :        ! The values from nmz0 to nmzxy are taken into account in order
     451             :        ! to get the real-space derivative smooth around nmzxy+1. 
     452           5 :        nmz0= vacuum%nmzxy+1+(obsolete%ndvgrd/2)-obsolete%ndvgrd
     453           5 :        IF (nmz0 <= 0) THEN ! usually vacuum%nmzxy>obsolete%ndvgrd 
     454           0 :           nmz0= 1
     455             :        END IF
     456             : 
     457           5 :        ALLOCATE ( rhtz(vacuum%nmzd,input%jspins) )
     458             : 
     459        1535 :        DO ip=nmz0,vacuum%nmz 
     460         770 :           IF (.not. noco%l_noco) THEN
     461        1989 :              DO js=1,input%jspins 
     462         765 :                 rhtz(ip,js)= den%vacz(ip,ivac,js)
     463             :              END DO
     464             :           ELSE
     465           0 :              af2(0,1) = den%vacz(ip,ivac,1)
     466           0 :              af2(0,2) = den%vacz(ip,ivac,2)
     467           0 :              mx(0) = den%vacz(ip,ivac,3)
     468           0 :              my(0) = den%vacz(ip,ivac,4)
     469           0 :              chdens= (af2(0,1)+af2(0,2))/2.
     470           0 :              magmom(0,1)= mx(0)**2 + my(0)**2 + ((af2(0,1)-af2(0,2))/2.)**2
     471           0 :              magmom(0,1)= SQRT(magmom(0,1))
     472           0 :              rhtz(ip,1)= chdens + magmom(0,1)
     473           0 :              rhtz(ip,2)= chdens - magmom(0,1) 
     474             :           END IF
     475             :        END DO
     476             : 
     477           5 :        IF (noco%l_noco) THEN 
     478           0 :           DEALLOCATE ( magmom,mx,my ) 
     479           0 :           ALLOCATE ( dummy(vacuum%nmz) )
     480           0 :           DO js=1,input%jspins
     481             :              CALL grdchlh(0,1,vacuum%nmz-nmz0+1,vacuum%delz,dummy,rhtz(nmz0,js),obsolete%ndvgrd,&
     482           0 :                   rhtdz(nmz0,js),rhtdzz(nmz0,js))
     483             :           END DO
     484           0 :           DEALLOCATE ( dummy )
     485             : 
     486             :        END IF
     487             : 
     488             :        !       calculate the quantities such as abs(grad(rho)),.. used in
     489             :        !c      evaluating the gradient contributions to potential and
     490             :        !c      energy.
     491             : 
     492           5 :        CALL xcpot%alloc_gradients(nmzdiff,input%jspins,grad)
     493           5 :        IF (xcpot%needs_grad())  THEN
     494           5 :           if(oneD%odi%d1)then
     495             : !             CALL od_mkgz(&
     496             : !                              cell%z1,vacuum%nmzxy,vacuum%delz,&
     497             : !                              nmzdiff,input%jspins,&
     498             : !                              rhtz(vacuum%nmzxy+1,1),rhtz(vacuum%nmzxy+1,input%jspins),&
     499             : !                              rhtdz(vacuum%nmzxy+1,1), rhtdz(vacuum%nmzxy+1,input%jspins),&
     500             : !                              rhtdzz(vacuum%nmzxy+1,1),rhtdzz(vacuum%nmzxy+1,input%jspins),&
     501             : !                              agr,agru,agrd,g2r,g2ru,g2rd,gggr,gggru,gggrd,&
     502             : !                              gzgr)
     503           0 :              CALL judft_error("OneD not implemented")
     504             :           ELSE
     505             :              CALL mkgz(nmzdiff,input%jspins, rhtz(vacuum%nmzxy+1,1),rhtz(vacuum%nmzxy+1,input%jspins),&
     506             :                   rhtdz(vacuum%nmzxy+1,1),rhtdz(vacuum%nmzxy+1,input%jspins),rhtdzz(vacuum%nmzxy+1,1),&
     507           5 :                   rhtdzz(vacuum%nmzxy+1,input%jspins),grad)
     508             :           endif
     509             :        ENDIF
     510             : 
     511             :        !       calculate vxc for z now beyond warping region
     512           5 :        DEALLOCATE ( af2)
     513           5 :        ALLOCATE ( rhtxc(vacuum%nmzd,input%jspins) )
     514           5 :        ALLOCATE ( vxcz(vacuum%nmzd,input%jspins) )
     515           5 :        ALLOCATE ( vxz(vacuum%nmzd,input%jspins) )
     516             : 
     517          21 :        DO js=1,input%jspins
     518             :           !DO i=0,ifftd2-1
     519             :           !  af2(i,js)=0.0
     520             :           !ENDDO
     521        1213 :           DO ip=vacuum%nmzxy+1,vacuum%nmz
     522        1200 :              rhtxc(ip-vacuum%nmzxy,js) = max(rhtz(ip,js),d_15) !+gb
     523        1208 :              rhmnv=min(rhmnv,rhtxc(ip-vacuum%nmzxy,js))
     524             :           ENDDO
     525             :        ENDDO
     526             : 
     527             :    
     528           5 :        IF (rhmnv.lt.obsolete%chng) THEN
     529           0 :           WRITE (6,'('' rhmn.lt.obsolete%chng. rhmn,obsolete%chng='',2d9.2)') rhmnv,obsolete%chng
     530             :           !           CALL juDFT_error("vvacxcg: rhmn.lt.chng",calledby="vvacxcg")
     531             :        ENDIF
     532             : 
     533           5 :        CALL xcpot%get_vxc(input%jspins,rhtxc(:nmzdiff,:),vxcz(:nmzdiff,:),vxz(:nmzdiff,:),grad)
     534             :    
     535          13 :        DO js = 1,input%jspins
     536        1213 :           DO ip = vacuum%nmzxy + 1,vacuum%nmz
     537        1208 :              vxc%vacz(ip,ivac,js) = vxc%vacz(ip,ivac,js) + vxcz(ip-vacuum%nmzxy,js)
     538             :           ENDDO
     539             :        ENDDO
     540             : 
     541             :        !
     542           5 :        WRITE (6,fmt=8020) ivac, (vxc%vacz(vacuum%nmz,ivac,js),js=1,input%jspins)
     543             : 8020   FORMAT(/,5x,'vacuum zero for vacuum',i3,' = ',2f14.10)
     544             :        !
     545             :        !     calculate the ex-corr. energy density now beyond warping region
     546             :        !
     547           5 :        IF (ALLOCATED(exc%vacz)) THEN
     548           5 :           CALL xcpot%get_exc(input%jspins,rhtxc(:nmzdiff,:),exc%vacz(vacuum%nmzxy+1:,ivac,1),grad, mt_call=.False.)
     549             :        ENDIF
     550             : 
     551             : 
     552           5 :        DEALLOCATE ( bf2)
     553          10 :        DEALLOCATE ( v_x,v_xc,e_xc,vxz,vxcz,rhtz,rhtxc )
     554             : 
     555             :     ENDDO    ! loop over vacua (ivac)
     556             :     
     557             : 
     558             : 
     559           5 :   END SUBROUTINE vvacxcg
     560           5 : END MODULE m_vvacxcg

Generated by: LCOV version 1.13