LCOV - code coverage report
Current view: top level - vgen - mkgxyz3.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 136 138 98.6 %
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             : MODULE m_mkgxyz3
       7             :   USE m_judft
       8             :   !.....------------------------------------------------------------------
       9             :   !c     by use of cartesian x,y,z components of charge density gradients,
      10             :   !c     make the quantities
      11             :   !cc      agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
      12             :   !cc    used to calculate gradient contribution to xc potential and
      13             :   !cc    energy.
      14             :   !c.....------------------------------------------------------------------
      15             : CONTAINS
      16         681 :   SUBROUTINE mkgxyz3(vl,dvx,dvy,dvz,dvxx,dvyy,dvzz,dvyz,dvzx,dvxy,grad)
      17             :     USE m_types
      18             :     IMPLICIT NONE
      19             :     REAL,    INTENT (IN) :: vl(:,:)
      20             :     REAL,    INTENT (IN) :: dvx(:,:),dvy(:,:),dvz(:,:)
      21             :     REAL, INTENT (IN) :: dvxx(:,:),dvyy(:,:),dvzz(:,:)
      22             :     REAL, INTENT (IN) :: dvyz(:,:),dvzx(:,:),dvxy(:,:)
      23             : 
      24             :     TYPE(t_gradients),INTENT(INOUT)::grad
      25             : 
      26             :     REAL vlt,dvxt,dvyt,dvzt,dvxxt,dvyyt,dvzzt,dvyzt,dvzxt,dvxyt,&
      27             :          vlu,dvxu,dvyu,dvzu,dvxxu,dvyyu,dvzzu,dvyzu,dvzxu,dvxyu,&
      28             :          vld,dvxd,dvyd,dvzd,dvxxd,dvyyd,dvzzd,dvyzd,dvzxd,dvxyd,&
      29             :          dagrxt,dagrxd,dagrxu,dagryt,dagryd,dagryu,dagrzt,dagrzd,&
      30             :          dagrzu,dzdx,dzdy,dzdz
      31             :     REAL, PARAMETER  :: sml = 1.e-14
      32             :     INTEGER i,js,jspins,nsp
      33             : 
      34         681 :     nsp=SIZE(dvx,1)
      35         681 :     jspins=SIZE(dvx,2)
      36             : 
      37         681 :     IF (ALLOCATED(grad%gr)) THEN
      38             :        !      Gradients for libxc
      39         240 :        DO js=1,jspins
      40     1024140 :           DO i=1,nsp
      41      512100 :              grad%gr(:,i,js)=(/dvx(i,js),dvy(i,js),dvz(i,js)/)
      42             :           ENDDO
      43             :        END DO
      44          40 :        IF(ALLOCATED(grad%sigma)) THEN
      45             :           !Use only contracted gradients for libxc
      46          40 :           IF (jspins==1) THEN
      47           0 :              DO i=1,nsp
      48           0 :                 grad%sigma(1,i) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
      49             :              ENDDO
      50             :           ELSE
      51      409640 :              DO i=1,nsp
      52      204800 :                 grad%sigma(1,i) = dvx(i,1)*dvx(i,1) + dvy(i,1)*dvy(i,1) + dvz(i,1)*dvz(i,1)
      53      204800 :                 grad%sigma(2,i) = dvx(i,1)*dvx(i,2) + dvy(i,1)*dvy(i,2) + dvz(i,1)*dvz(i,2)
      54      204840 :                 grad%sigma(3,i) = dvx(i,2)*dvx(i,2) + dvy(i,2)*dvy(i,2) + dvz(i,2)*dvz(i,2)
      55             :              ENDDO
      56             :           ENDIF
      57             :        END IF
      58          40 :        IF(ALLOCATED(grad%laplace)) THEN
      59         240 :           DO js=1,jspins
      60     1024140 :              DO i=1,nsp
      61      512100 :                 grad%laplace(i,js)= dvxx(i,js)+dvyy(i,js)+dvzz(i,js)
      62             :              ENDDO
      63             :           ENDDO
      64             :        ENDIF
      65             :        RETURN
      66             :     ENDIF
      67             : 
      68         641 :     IF (ANY(SHAPE(vl).NE.SHAPE(dvx))) CALL judft_error("Gradients for internal GGA called with inconsistent sizes",hint="This is a bug")
      69             :     
      70     1653082 :     DO i = 1,size(grad%agrt)
      71     1652441 :        grad%agrt(i)  = 0.0
      72     1652441 :        grad%agru(i)  = 0.0
      73     1652441 :        grad%agrd(i)  = 0.0
      74     1652441 :        grad%gggrt(i) = 0.0
      75     1652441 :        grad%gggru(i) = 0.0
      76     1652441 :        grad%gggrd(i) = 0.0
      77     1652441 :        grad%gzgr(i)  = 0.0
      78     1652441 :        grad%g2rt(i)  = 0.0
      79     1652441 :        grad%g2ru(i)  = 0.0
      80     1653082 :        grad%g2rd(i)  = 0.0
      81             :     ENDDO
      82             : 
      83         641 :     IF (jspins.eq.1) THEN
      84             : 
      85     1371107 :        DO i = 1,nsp
      86             : 
      87      685440 :           vlu   = max(vl(i,1)/2,sml)
      88      685440 :           dvxu  = dvx(i,1)/2
      89      685440 :           dvyu  = dvy(i,1)/2
      90      685440 :           dvzu  = dvz(i,1)/2
      91      685440 :           dvxxu = dvxx(i,1)/2
      92      685440 :           dvyyu = dvyy(i,1)/2
      93      685440 :           dvzzu = dvzz(i,1)/2
      94      685440 :           dvyzu = dvyz(i,1)/2
      95      685440 :           dvzxu = dvzx(i,1)/2
      96      685440 :           dvxyu = dvxy(i,1)/2
      97             : 
      98      685440 :           vld   = vlu
      99      685440 :           dvxd  = dvxu
     100      685440 :           dvyd  = dvyu
     101      685440 :           dvzd  = dvzu
     102      685440 :           dvxxd = dvxxu
     103      685440 :           dvyyd = dvyyu
     104      685440 :           dvzzd = dvzzu
     105      685440 :           dvyzd = dvyzu
     106      685440 :           dvzxd = dvzxu
     107      685440 :           dvxyd = dvxyu
     108             : 
     109             : 
     110      685440 :           vlt = vlu + vld
     111             : 
     112      685440 :           dvxt  = dvxu  + dvxd
     113      685440 :           dvyt  = dvyu  + dvyd
     114      685440 :           dvzt  = dvzu  + dvzd
     115      685440 :           dvxxt = dvxxu + dvxxd
     116      685440 :           dvyyt = dvyyu + dvyyd
     117      685440 :           dvzzt = dvzzu + dvzzd
     118      685440 :           dvyzt = dvyzu + dvyzd
     119      685440 :           dvzxt = dvzxu + dvzxd
     120      685440 :           dvxyt = dvxyu + dvxyd
     121             : 
     122             :           !         agr: abs(grad(ro)), t,u,d for total, up and down.
     123             : 
     124      685440 :           grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
     125      685440 :           grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
     126      685440 :           grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
     127             : 
     128      685440 :           dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
     129      685440 :           dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
     130      685440 :           dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
     131             : 
     132      685440 :           dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
     133      685440 :           dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
     134      685440 :           dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
     135             : 
     136      685440 :           dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
     137      685440 :           dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
     138      685440 :           dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
     139             : 
     140      685440 :           grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
     141      685440 :           grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
     142      685440 :           grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
     143             : 
     144             :           !         dzdx=d(zeta)/dx,..
     145             : 
     146      685440 :           dzdx = (dvxu-dvxd)/vlt - (vlu-vld)*dvxt/vlt**2
     147      685440 :           dzdy = (dvyu-dvyd)/vlt - (vlu-vld)*dvyt/vlt**2
     148      685440 :           dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/vlt**2
     149             : 
     150             :           !         gzgr=grad(zeta)*grad(ro).
     151             : 
     152      685440 :           grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
     153             : 
     154             :           !         g2r: grad(grad(ro))
     155             : 
     156      685440 :           grad%g2rt(i) = dvxxt + dvyyt + dvzzt
     157      685440 :           grad%g2ru(i) = dvxxu + dvyyu + dvzzu
     158      685667 :           grad%g2rd(i) = dvxxd + dvyyd + dvzzd
     159             : 
     160             : 
     161             :        ENDDO
     162             : 
     163             :     ELSE
     164             : 
     165     1934416 :        DO i = 1,nsp
     166             : 
     167      967001 :           vlu   = max(vl(i,1),sml)
     168      967001 :           dvxu  = dvx(i,1)
     169      967001 :           dvyu  = dvy(i,1)
     170      967001 :           dvzu  = dvz(i,1)
     171      967001 :           dvxxu = dvxx(i,1)
     172      967001 :           dvyyu = dvyy(i,1)
     173      967001 :           dvzzu = dvzz(i,1)
     174      967001 :           dvyzu = dvyz(i,1)
     175      967001 :           dvzxu = dvzx(i,1)
     176      967001 :           dvxyu = dvxy(i,1)
     177             : 
     178      967001 :           vld   = max(vl(i,jspins),sml)
     179      967001 :           dvxd  = dvx(i,jspins)
     180      967001 :           dvyd  = dvy(i,jspins)
     181      967001 :           dvzd  = dvz(i,jspins)
     182      967001 :           dvxxd = dvxx(i,jspins)
     183      967001 :           dvyyd = dvyy(i,jspins)
     184      967001 :           dvzzd = dvzz(i,jspins)
     185      967001 :           dvyzd = dvyz(i,jspins)
     186      967001 :           dvzxd = dvzx(i,jspins)
     187      967001 :           dvxyd = dvxy(i,jspins)
     188             : 
     189      967001 :           vlt = vlu + vld
     190             : 
     191      967001 :           dvxt  = dvxu  + dvxd
     192      967001 :           dvyt  = dvyu  + dvyd
     193      967001 :           dvzt  = dvzu  + dvzd
     194      967001 :           dvxxt = dvxxu + dvxxd
     195      967001 :           dvyyt = dvyyu + dvyyd
     196      967001 :           dvzzt = dvzzu + dvzzd
     197      967001 :           dvyzt = dvyzu + dvyzd
     198      967001 :           dvzxt = dvzxu + dvzxd
     199      967001 :           dvxyt = dvxyu + dvxyd
     200             : 
     201             :           !c         agr: abs(grad(ro)), t,u,d for total, up and down.
     202             : 
     203      967001 :           grad%agrt(i) = max(sqrt(dvxt**2 + dvyt**2 + dvzt**2),sml)
     204      967001 :           grad%agru(i) = max(sqrt(dvxu**2 + dvyu**2 + dvzu**2),sml)
     205      967001 :           grad%agrd(i) = max(sqrt(dvxd**2 + dvyd**2 + dvzd**2),sml)
     206             : 
     207      967001 :           dagrxt = (dvxt*dvxxt + dvyt*dvxyt + dvzt*dvzxt) / grad%agrt(i)
     208      967001 :           dagrxu = (dvxu*dvxxu + dvyu*dvxyu + dvzu*dvzxu) / grad%agru(i)
     209      967001 :           dagrxd = (dvxd*dvxxd + dvyd*dvxyd + dvzd*dvzxd) / grad%agrd(i)
     210             : 
     211      967001 :           dagryt = (dvxt*dvxyt + dvyt*dvyyt + dvzt*dvyzt) / grad%agrt(i)
     212      967001 :           dagryu = (dvxu*dvxyu + dvyu*dvyyu + dvzu*dvyzu) / grad%agru(i)
     213      967001 :           dagryd = (dvxd*dvxyd + dvyd*dvyyd + dvzd*dvyzd) / grad%agrd(i)
     214             : 
     215      967001 :           dagrzt = (dvxt*dvzxt + dvyt*dvyzt + dvzt*dvzzt) / grad%agrt(i)
     216      967001 :           dagrzu = (dvxu*dvzxu + dvyu*dvyzu + dvzu*dvzzu) / grad%agru(i)
     217      967001 :           dagrzd = (dvxd*dvzxd + dvyd*dvyzd + dvzd*dvzzd) / grad%agrd(i)
     218             : 
     219      967001 :           grad%gggrt(i) = dvxt*dagrxt + dvyt*dagryt + dvzt*dagrzt
     220      967001 :           grad%gggru(i) = dvxu*dagrxu + dvyu*dagryu + dvzu*dagrzu
     221      967001 :           grad%gggrd(i) = dvxd*dagrxd + dvyd*dagryd + dvzd*dagrzd
     222             : 
     223             :           !c         dzdx=d(zeta)/dx,..
     224             : 
     225      967001 :           dzdx = (dvxu-dvxd)/vlt -  (vlu-vld)*dvxt/vlt**2
     226      967001 :           dzdy = (dvyu-dvyd)/vlt -  (vlu-vld)*dvyt/vlt**2
     227      967001 :           dzdz = (dvzu-dvzd)/vlt -  (vlu-vld)*dvzt/vlt**2
     228             : 
     229             :           !c         gzgr=grad(zeta)*grad(ro).
     230             : 
     231      967001 :           grad%gzgr(i) = dzdx*dvxt + dzdy*dvyt + dzdz*dvzt
     232             : 
     233             :           !c         g2r: grad(grad(ro))
     234             : 
     235      967001 :           grad%g2rt(i) = dvxxt + dvyyt + dvzzt
     236      967001 :           grad%g2ru(i) = dvxxu + dvyyu + dvzzu
     237      967415 :           grad%g2rd(i) = dvxxd + dvyyd + dvzzd
     238             : 
     239             :        ENDDO
     240             : 
     241             :     ENDIF
     242             : 
     243             :   END SUBROUTINE mkgxyz3
     244             : END MODULE m_mkgxyz3

Generated by: LCOV version 1.13