LCOV - code coverage report
Current view: top level - vgen - mkgxyz3.f90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 100.0 % 140 140
Test Date: 2025-06-14 04:34:23 Functions: 100.0 % 1 1

            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_mkgxyz3
       7              :    USE m_judft
       8              :    !-----------------------------------------------------------------------------
       9              :    ! Using the cartesian components and derivatives of a charge density rho on
      10              :    ! the real space grid, make the following quantaties:
      11              :    !
      12              :    ! gr(js):      grad(rho_js)
      13              :    ! sigma:       |grad(rho)|^2 for jspins==1, otherwise three components, namely
      14              :    !              |grad(rho_up)|^2,grad(rho_up)*grad(rho_down),|grad(rho_down)|^2
      15              :    ! laplace(js): laplace(rho_js)
      16              :    !
      17              :    ! and these older components:
      18              :    ! agrt/u/d:    |grad(rho)| for total density/spin-up/spin-down
      19              :    ! g2rt/u/d:    laplace(rho)
      20              :    ! gggrt/u/d:   (grad(rho))*(grad(|grad(rho)|)) [scalar product]
      21              :    ! gzgr:        (grad(zeta))*(grad(rho)) for zeta=(rho_up-rho_down)/rho
      22              :    !
      23              :    ! which are used to calculate gradient contributions to the xc potential and
      24              :    ! energy.
      25              :    !
      26              :    ! vl is rho, dv[i][j] are the partial derivatives along one/two directions.
      27              :    !
      28              :    ! Modified so only allocated old quantities require calculations. A.N. 2019
      29              :    !
      30              :    ! Quantities fo libxc are calculated as well. D.W./M.R. 2018
      31              :    !
      32              :    ! Original script by T.A. 1996
      33              :    !-----------------------------------------------------------------------------
      34              : CONTAINS
      35         2984 :    SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvxz,dvxy,idx,grad)
      36              :       USE m_types
      37              :       IMPLICIT NONE
      38              :       REAL, INTENT (IN)                :: vl(:,:)
      39              :       REAL, INTENT (IN)                :: dvx(:,:),dvy(:,:),dvz(:,:)
      40              :       REAL, INTENT (IN)                :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
      41              :       REAL, INTENT (IN)                :: dvyz(:,:),dvxz(:,:),dvxy(:,:)
      42              :       INTEGER ,intent(in)              :: idx
      43              :       TYPE(t_gradients), INTENT(INOUT) :: grad
      44              : 
      45              :       REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvxzt,dvxyt, &
      46              :            vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvxzu,dvxyu, &
      47              :            vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvxzd,dvxyd, &
      48              :            dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
      49              :            dagrzu,dzdx,dzdy,dzdz
      50              :       REAL, PARAMETER  :: sml = 1.e-14
      51              :       INTEGER i,js,jspins,nsp
      52              : 
      53         2984 :       nsp=SIZE(dvx,1)
      54         2984 :       jspins=SIZE(dvx,2)
      55              : 
      56              :       ! Gradients for libxc and sourcefree calculations.
      57         2984 :       IF (ALLOCATED(grad%gr)) THEN
      58           39 :          DO js=1,jspins
      59       135591 :             DO i=1,nsp
      60       542235 :                grad%gr(:,i+idx,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
      61              :             END DO
      62              :          END DO
      63              :          ! Use contracted gradients only for libxc.
      64           12 :          IF(ALLOCATED(grad%sigma)) THEN
      65           12 :             IF (jspins==1) THEN
      66        41475 :                DO i=1,nsp
      67        41475 :                   grad%sigma(1,i+idx) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
      68              :                END DO
      69              :             ELSE
      70        35289 :                DO i=1,nsp
      71        35280 :                   grad%sigma(1,i+idx) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
      72        35280 :                   grad%sigma(2,i+idx) = dvx(i,1)*dvx(i,2) + dvy(i,1)*dvy(i,2) + dvz(i,1)*dvz(i,2)
      73        35289 :                   grad%sigma(3,i+idx) = dvx(i,2)*dvx(i,2) + dvy(i,2)*dvy(i,2) + dvz(i,2)*dvz(i,2)
      74              :                END DO
      75              :             END IF
      76              :          END IF
      77           12 :          IF(ALLOCATED(grad%laplace)) THEN
      78           39 :             DO js=1,jspins
      79       135591 :                DO i=1,nsp
      80       135579 :                   grad%laplace(i+idx,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
      81              :                END DO
      82              :             END DO
      83              :          END IF
      84           12 :          RETURN ! Do not calculate arrays for in-build GGA.
      85              :       END IF
      86              : 
      87         8916 :       IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
      88              : 
      89         2972 :            IF(ALLOCATED(grad%agrt)) THEN
      90      9964713 :          DO i = 1,nsp
      91      9961741 :             grad%agrt(idx+i)  = 0.0
      92      9961741 :             grad%agru(idx+i)  = 0.0
      93      9961741 :             grad%agrd(idx+i)  = 0.0
      94      9961741 :             grad%gggrt(idx+i) = 0.0
      95      9961741 :             grad%gggru(idx+i) = 0.0
      96      9961741 :             grad%gggrd(idx+i) = 0.0
      97      9961741 :             grad%gzgr(idx+i)  = 0.0
      98      9961741 :             grad%g2rt(idx+i)  = 0.0
      99      9961741 :             grad%g2ru(idx+i)  = 0.0
     100      9964713 :             grad%g2rd(idx+i)  = 0.0
     101              :          END DO
     102              : 
     103         2972 :          IF (jspins.eq.1) THEN
     104              : 
     105      8641275 :             DO i = 1,nsp
     106              : 
     107      8638845 :                vlu   = max(vl(i,1)/2,sml)
     108      8638845 :                dvxu  = dvx(i,1)/2
     109      8638845 :                dvyu  = dvy(i,1)/2
     110      8638845 :                dvzu  = dvz(i,1)/2
     111      8638845 :                dvxxu = dvxx(i,1)/2
     112      8638845 :                dvyyu = dvyy(i,1)/2
     113      8638845 :                dvzzu = dvzz(i,1)/2
     114      8638845 :                dvyzu = dvyz(i,1)/2
     115      8638845 :                dvxzu = dvxz(i,1)/2
     116      8638845 :                dvxyu = dvxy(i,1)/2
     117              : 
     118      8638845 :                vld   = vlu
     119      8638845 :                dvxd  = dvxu
     120      8638845 :                dvyd  = dvyu
     121      8638845 :                dvzd  = dvzu
     122      8638845 :                dvxxd = dvxxu
     123      8638845 :                dvyyd = dvyyu
     124      8638845 :                dvzzd = dvzzu
     125      8638845 :                dvyzd = dvyzu
     126      8638845 :                dvxzd = dvxzu
     127      8638845 :                dvxyd = dvxyu
     128              : 
     129      8638845 :                vlt = vlu + vld
     130      8638845 :                dvxt  = dvxu  + dvxd
     131      8638845 :                dvyt  = dvyu  + dvyd
     132      8638845 :                dvzt  = dvzu  + dvzd
     133      8638845 :                dvxxt = dvxxu + dvxxd
     134      8638845 :                dvyyt = dvyyu + dvyyd
     135      8638845 :                dvzzt = dvzzu + dvzzd
     136      8638845 :                dvyzt = dvyzu + dvyzd
     137      8638845 :                dvxzt = dvxzu + dvxzd
     138      8638845 :                dvxyt = dvxyu + dvxyd
     139              : 
     140      8638845 :                grad%agrt(idx+i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
     141      8638845 :                grad%agru(idx+i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
     142      8638845 :                grad%agrd(idx+i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
     143              : 
     144      8638845 :                dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvxzt) / grad%agrt(idx+i)
     145      8638845 :                dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvxzu) / grad%agru(idx+i)
     146      8638845 :                dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvxzd) / grad%agrd(idx+i)
     147              : 
     148      8638845 :                dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(idx+i)
     149      8638845 :                dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(idx+i)
     150      8638845 :                dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(idx+i)
     151              : 
     152      8638845 :                dagrzt = (dvxt*dvxzt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(idx+i)
     153      8638845 :                dagrzu = (dvxu*dvxzu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(idx+i)
     154      8638845 :                dagrzd = (dvxd*dvxzd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(idx+i)
     155              : 
     156      8638845 :                grad%gggrt(idx+i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
     157      8638845 :                grad%gggru(idx+i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
     158      8638845 :                grad%gggrd(idx+i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
     159              : 
     160      8638845 :                dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
     161      8638845 :                dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
     162      8638845 :                dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
     163              : 
     164      8638845 :                grad%gzgr(idx+i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
     165              : 
     166      8638845 :                grad%g2rt(idx+i) = dvxxt + dvyyt + dvzzt
     167      8638845 :                grad%g2ru(idx+i) = dvxxu + dvyyu + dvzzu
     168      8641275 :                grad%g2rd(idx+i) = dvxxd + dvyyd + dvzzd
     169              :             ENDDO
     170              :          ELSE
     171      1323438 :             DO i = 1,nsp
     172              : 
     173      1322896 :                vlu   = max(vl(i,1),sml)
     174      1322896 :                dvxu  = dvx(i,1)
     175      1322896 :                dvyu  = dvy(i,1)
     176      1322896 :                dvzu  = dvz(i,1)
     177      1322896 :                dvxxu = dvxx(i,1)
     178      1322896 :                dvyyu = dvyy(i,1)
     179      1322896 :                dvzzu = dvzz(i,1)
     180      1322896 :                dvyzu = dvyz(i,1)
     181      1322896 :                dvxzu = dvxz(i,1)
     182      1322896 :                dvxyu = dvxy(i,1)
     183              : 
     184      1322896 :                vld   = max(vl(i,jspins),sml)
     185      1322896 :                dvxd  = dvx(i,jspins)
     186      1322896 :                dvyd  = dvy(i,jspins)
     187      1322896 :                dvzd  = dvz(i,jspins)
     188      1322896 :                dvxxd = dvxx(i,jspins)
     189      1322896 :                dvyyd = dvyy(i,jspins)
     190      1322896 :                dvzzd = dvzz(i,jspins)
     191      1322896 :                dvyzd = dvyz(i,jspins)
     192      1322896 :                dvxzd = dvxz(i,jspins)
     193      1322896 :                dvxyd = dvxy(i,jspins)
     194              : 
     195      1322896 :                vlt = vlu + vld
     196              : 
     197      1322896 :                dvxt  = dvxu  + dvxd
     198      1322896 :                dvyt  = dvyu  + dvyd
     199      1322896 :                dvzt  = dvzu  + dvzd
     200      1322896 :                dvxxt = dvxxu + dvxxd
     201      1322896 :                dvyyt = dvyyu + dvyyd
     202      1322896 :                dvzzt = dvzzu + dvzzd
     203      1322896 :                dvyzt = dvyzu + dvyzd
     204      1322896 :                dvxzt = dvxzu + dvxzd
     205      1322896 :                dvxyt = dvxyu + dvxyd
     206              : 
     207      1322896 :                grad%agrt(idx+i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
     208      1322896 :                grad%agru(idx+i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
     209      1322896 :                grad%agrd(idx+i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
     210              : 
     211      1322896 :                dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvxzt) / grad%agrt(idx+i)
     212      1322896 :                dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvxzu) / grad%agru(idx+i)
     213      1322896 :                dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvxzd) / grad%agrd(idx+i)
     214              : 
     215      1322896 :                dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(idx+i)
     216      1322896 :                dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(idx+i)
     217      1322896 :                dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(idx+i)
     218              : 
     219      1322896 :                dagrzt = (dvxt*dvxzt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(idx+i)
     220      1322896 :                dagrzu = (dvxu*dvxzu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(idx+i)
     221      1322896 :                dagrzd = (dvxd*dvxzd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(idx+i)
     222              : 
     223      1322896 :                grad%gggrt(idx+i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
     224      1322896 :                grad%gggru(idx+i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
     225      1322896 :                grad%gggrd(idx+i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
     226              : 
     227      1322896 :                dzdx = (dvxu-dvxd)/vlt -  (vlu-vld)*dvxt/vlt**2
     228      1322896 :                dzdy = (dvyu-dvyd)/vlt -  (vlu-vld)*dvyt/vlt**2
     229      1322896 :                dzdz = (dvzu-dvzd)/vlt -  (vlu-vld)*dvzt/vlt**2
     230              : 
     231      1322896 :                grad%gzgr(idx+i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
     232              : 
     233      1322896 :                grad%g2rt(idx+i) = dvxxt + dvyyt + dvzzt
     234      1322896 :                grad%g2ru(idx+i) = dvxxu + dvyyu + dvzzu
     235      1323438 :                grad%g2rd(idx+i) = dvxxd + dvyyd + dvzzd
     236              : 
     237              :             END DO
     238              :          END IF
     239              :       END IF
     240              : 
     241              :    END SUBROUTINE mkgxyz3
     242              : END MODULE m_mkgxyz3
        

Generated by: LCOV version 2.0-1