LCOV - code coverage report
Current view: top level - vgen - od_mkgxyz3.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 127 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             :       MODULE m_od_mkgxyz3
       8             :       CONTAINS
       9           0 :       SUBROUTINE od_mkgxyz3(
      10           0 :      >     ndm,jsdm,ng3,jspins,vl,r,
      11           0 :      >     dvz,dvp,dvr,dvzz,dvpp,dvrr,dvpr,dvrz,dvzp,
      12           0 :      <     agrt,agru,agrd, g2rt,g2ru,g2rd, gggrt,gggru,gggrd,
      13           0 :      <     gzgr)
      14             : 
      15             : c     by use of cylindrical z,phi,r components of charge density gradients,
      16             : c     make the quantities
      17             : c     agrt,agru,agrd,g2rt,g2ru,g2rd,gggrt,gggru,gggrd,gzgr
      18             : c     used to calculate gradient contribution to xc potential and
      19             : c     energy.
      20             : c     Have no time so writing after this crazy subroutines of t.asada
      21             : c     Y.Mokrousov, 15-16 october 2002
      22             : 
      23             :       IMPLICIT NONE
      24             : 
      25             :       integer, intent (in) ::  ndm,ng3,jsdm,jspins
      26             :       real,    intent (in) ::  vl(ndm,jsdm),r
      27             :       real,    intent (in) ::  dvz(ndm,jsdm),dvp(ndm,jsdm)
      28             :       real,    intent (in) ::  dvr(ndm,jsdm),dvzz(ndm,jsdm)
      29             :       real,    intent (in) ::  dvpp(ndm,jsdm),dvrr(ndm,jsdm)
      30             :       real,    intent (in) ::  dvpr(ndm,jsdm),dvrz(ndm,jsdm)
      31             :       real,    intent (in) ::  dvzp(ndm,jsdm)
      32             : 
      33             :       real,    intent (out) :: agrt(ndm),agru(ndm),agrd(ndm),g2rt(ndm)
      34             :       real,    intent (out) :: g2ru(ndm),g2rd(ndm),gggrt(ndm)
      35             :       real,    intent (out) :: gggru(ndm),gggrd(ndm)
      36             :       real,    intent (out) :: gzgr(ndm)
      37             : 
      38             : c     local
      39             : 
      40             :       real    :: vlt,dvzt,dvpt,dvrt,dvzzt,dvppt,dvrrt,dvprt,dvrzt,dvzpt
      41             :       real    :: vlu,dvzu,dvpu,dvru,dvzzu,dvppu,dvrru,dvpru,dvrzu,dvzpu
      42             :       real    :: vld,dvzd,dvpd,dvrd,dvzzd,dvppd,dvrrd,dvprd,dvrzd,dvzpd
      43             :       real    :: dagrzt,dagrzd,dagrzu,dagrpt,dagrpd,dagrpu,dagrrt,dagrrd
      44             :       real    :: dagrru,dzdz,dzdp,dzdr
      45             :       real    :: sml
      46             :       integer :: i
      47             : 
      48             : c     so what do all these hieroglyphs mean?
      49             : c     for the description of input parameters go to vvacxcg_1.F 
      50             : c     as far as in a cylindrical case r and phi variables are somehow 
      51             : c     connected, all the staff will depend on the r-coordinate of the point
      52             : c     at which it is calculated, which is noted by r
      53             : c     I use three types of densities: rho(up),rho(down),rho(total)
      54             : c     so that all the corresponding variables are noted by 'u','d', and 't'
      55             : c     in the end, special case is variable zeta which is an expression:
      56             : c     zeta = (rho(up)-rho(down))/rho(total)
      57             : c     different types of differential expressions are calculated here:
      58             : c     agrt,agru,agrd:
      59             : c                  --->  abs(grad(rho))
      60             : c     g2rt,g2ru,g2rd:
      61             : c                  --->  laplasian(rho)
      62             : c     gggrt,gggru,gggrd :
      63             : c                  --->  grad(rho)*grad(abs(grad(rho)))
      64             : c     gzgr:
      65             : c                  --->  grad(zeta)*grad(rho(total))  
      66             : 
      67           0 :       sml = 1.e-14
      68             : 
      69           0 :       do i = 1,ndm
      70           0 :           agrt(i) = 0.0
      71           0 :           agru(i) = 0.0
      72           0 :           agrd(i) = 0.0
      73           0 :           gggrt(i) = 0.0
      74           0 :           gggru(i) = 0.0
      75           0 :           gggrd(i) = 0.0
      76           0 :           gzgr(i) = 0.0
      77           0 :           g2rt(i) = 0.0
      78           0 :           g2ru(i) = 0.0
      79           0 :           g2rd(i) = 0.0
      80             :       end do
      81             : 
      82           0 :          if (jspins.eq.1) then
      83             :             
      84             :             
      85           0 :             do 10 i = 1,ng3
      86           0 :                vlu=max(vl(i,1)/2.,sml)
      87           0 :                dvzu=dvz(i,1)/2.
      88           0 :                dvpu=dvp(i,1)/2.
      89           0 :                dvru=dvr(i,1)/2.
      90           0 :                dvzzu=dvzz(i,1)/2.
      91           0 :                dvppu=dvpp(i,1)/2.
      92           0 :                dvrru=dvrr(i,1)/2.
      93           0 :                dvpru=dvpr(i,1)/2.
      94           0 :                dvrzu=dvrz(i,1)/2.
      95           0 :                dvzpu=dvzp(i,1)/2.
      96             :                
      97           0 :                vld=vlu
      98           0 :                dvzd=dvzu
      99           0 :                dvpd=dvpu
     100           0 :                dvrd=dvru
     101           0 :                dvzzd=dvzzu
     102           0 :                dvppd=dvppu
     103           0 :                dvrrd=dvrru
     104           0 :                dvprd=dvpru
     105           0 :                dvrzd=dvrzu
     106           0 :                dvzpd=dvzpu
     107             :                
     108             :                
     109           0 :                vlt = vlu + vld
     110             :                
     111           0 :                dvzt = dvzu + dvzd
     112           0 :                dvpt = dvpu + dvpd
     113           0 :                dvrt = dvru + dvrd
     114           0 :                dvzzt = dvzzu + dvzzd
     115           0 :                dvppt = dvppu + dvppd
     116           0 :                dvrrt = dvrru + dvrrd
     117           0 :                dvprt = dvpru + dvprd
     118           0 :                dvrzt = dvrzu + dvrzd
     119           0 :                dvzpt = dvzpu + dvzpd
     120             :                
     121             : c     agr: abs(grad(ro)), t,u,d for total, up and down.
     122           0 :            agrt(i) = max(sqrt(dvzt**2+(dvpt**2)/(r**2)+dvrt**2),sml)
     123           0 :            agru(i) = max(sqrt(dvzu**2+(dvpu**2)/(r**2)+dvru**2),sml)
     124           0 :            agrd(i) = max(sqrt(dvzd**2+(dvpd**2)/(r**2)+dvrd**2),sml)
     125             : c     d(abs(grad(rho)))/dz
     126             :                dagrzt = (dvzt*dvzzt+(dvpt*dvzpt)/(r**2)+dvrt*dvrzt)/
     127           0 :      /              agrt(i)
     128             :                dagrzu = (dvzu*dvzzu+(dvpu*dvzpu)/(r**2)+dvru*dvrzu)/
     129           0 :      /              agru(i)
     130             :                dagrzd = (dvzd*dvzzd+(dvpd*dvzpd)/(r**2)+dvrd*dvrzd)/
     131           0 :      /              agrd(i)
     132             : c     d(abs(grad(rho)))/d(phi)
     133             :                dagrpt = (dvzt*dvzpt+(dvpt*dvppt)/(r**2)+dvrt*dvprt)/
     134           0 :      /              agrt(i)
     135             :                dagrpu = (dvzu*dvzpu+(dvpu*dvppu)/(r**2)+dvru*dvpru)/
     136           0 :      /              agru(i)
     137             :                dagrpd = (dvzd*dvzpd+(dvpd*dvppd)/(r**2)+dvrd*dvprd)/
     138           0 :      /              agrd(i)
     139             : c     d(abs(grad(rho)))/dr
     140             :                dagrrt = (dvzt*dvrzt+dvrt*dvrrt+
     141           0 :      +              (dvpt*dvprt)/(r**2)-(dvpt**2)/(r**3) )/agrt(i)
     142             :                dagrru = (dvzu*dvrzu+dvru*dvrru+
     143           0 :      +              (dvpu*dvpru)/(r**2)-(dvpu**2)/(r**3) )/agru(i)
     144             :                dagrrd = (dvzd*dvrzd+dvrd*dvrrd+
     145           0 :      +              (dvpd*dvprd)/(r**2)-(dvpd**2)/(r**3) )/agrd(i)
     146             : c     grad(rho)*grad(abs(grad(rho)))
     147             :                gggrt(i) = dvzt*dagrzt + (dvpt*dagrpt)/(r**2) + dvrt*
     148           0 :      *              dagrrt
     149             :                gggru(i) = dvzu*dagrzu + (dvpu*dagrpu)/(r**2) + dvru*
     150           0 :      *              dagrru
     151             :                gggrd(i) = dvzd*dagrzd + (dvpd*dagrpd)/(r**2) + dvrd*
     152           0 :      *              dagrrd
     153             : c     dzdz=d(zeta)/dz,dzdp=d(zeta)/dp,dzdr=d(zeta)/dr
     154           0 :                dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/(vlt**2)
     155           0 :                dzdp = (dvpu-dvpd)/vlt - (vlu-vld)*dvpt/(vlt**2)
     156           0 :                dzdr = (dvru-dvrd)/vlt - (vlu-vld)*dvrt/(vlt**2)
     157             : c     gzgr=grad(zeta)*grad(ro)
     158           0 :                gzgr(i) = dzdz*dvzt + (dzdp*dvpt)/(r**2) + dzdr*dvrt
     159             : c     g2r: laplasian(rho)
     160           0 :                g2rt(i) = dvzzt + (dvppt)/(r**2) + dvrrt + dvrt/r
     161           0 :                g2ru(i) = dvzzu + (dvppu)/(r**2) + dvrru + dvru/r
     162           0 :                g2rd(i) = dvzzd + (dvppd)/(r**2) + dvrrd + dvrd/r
     163           0 :  10         end do
     164             :             
     165             :          else
     166             :             
     167           0 :             do 20 i = 1,ng3
     168           0 :                vlu=max(vl(i,1),sml)
     169           0 :                dvzu=dvz(i,1)
     170           0 :                dvpu=dvp(i,1)
     171           0 :                dvru=dvr(i,1)
     172           0 :                dvzzu=dvzz(i,1)
     173           0 :                dvppu=dvpp(i,1)
     174           0 :                dvrru=dvrr(i,1)
     175           0 :                dvpru=dvpr(i,1)
     176           0 :                dvrzu=dvrz(i,1)
     177           0 :                dvzpu=dvzp(i,1)
     178             :                
     179           0 :                vld=max(vl(i,jspins),sml)
     180           0 :                dvzd=dvz(i,jspins)
     181           0 :                dvpd=dvp(i,jspins)
     182           0 :                dvrd=dvr(i,jspins)
     183           0 :                dvzzd=dvzz(i,jspins)
     184           0 :                dvppd=dvpp(i,jspins)
     185           0 :                dvrrd=dvrr(i,jspins)
     186           0 :                dvprd=dvpr(i,jspins)
     187           0 :                dvrzd=dvrz(i,jspins)
     188           0 :                dvzpd=dvzp(i,jspins)
     189             :                
     190           0 :                vlt = vlu + vld
     191             :                
     192           0 :                dvzt = dvzu + dvzd
     193           0 :                dvpt = dvpu + dvpd
     194           0 :                dvrt = dvru + dvrd
     195           0 :                dvzzt = dvzzu + dvzzd
     196           0 :                dvppt = dvppu + dvppd
     197           0 :                dvrrt = dvrru + dvrrd
     198           0 :                dvprt = dvpru + dvprd
     199           0 :                dvrzt = dvrzu + dvrzd
     200           0 :                dvzpt = dvzpu + dvzpd
     201             : c     agr: abs(grad(ro)), t,u,d for total, up and down.
     202           0 :                agrt(i) = max(sqrt(dvzt**2+(dvpt**2)/(r**2)+dvrt**2),sml)
     203           0 :                agru(i) = max(sqrt(dvzu**2+(dvpu**2)/(r**2)+dvru**2),sml)
     204           0 :                agrd(i) = max(sqrt(dvzd**2+(dvpd**2)/(r**2)+dvrd**2),sml)
     205             : c     d(abs(grad(rho)))/dz
     206             :                dagrzt = (dvzt*dvzzt+(dvpt*dvzpt)/(r**2)+dvrt*dvrzt)/
     207           0 :      /              agrt(i)
     208             :                dagrzu = (dvzu*dvzzu+(dvpu*dvzpu)/(r**2)+dvru*dvrzu)/
     209           0 :      /              agru(i)
     210             :                dagrzd = (dvzd*dvzzd+(dvpd*dvzpd)/(r**2)+dvrd*dvrzd)/
     211           0 :      /              agrd(i)
     212             : c     d(abs(grad(rho)))/d(phi)
     213             :                dagrpt = (dvzt*dvzpt+(dvpt*dvppt)/(r**2)+dvrt*dvprt)/
     214           0 :      /              agrt(i)
     215             :                dagrpu = (dvzu*dvzpu+(dvpu*dvppu)/(r**2)+dvru*dvpru)/
     216           0 :      /              agru(i)
     217             :                dagrpd = (dvzd*dvzpd+(dvpd*dvppd)/(r**2)+dvrd*dvprd)/
     218           0 :      /              agrd(i)
     219             : c     d(abs(grad(rho)))/dr
     220             :                dagrrt = (dvzt*dvrzt+dvrt*dvrrt+
     221           0 :      +              (dvpt*dvprt)/(r**2)-(dvpt**2)/(r**3) )/agrt(i)
     222             :                dagrru = (dvzu*dvrzu+dvru*dvrru+
     223           0 :      +              (dvpu*dvpru)/(r**2)-(dvpu**2)/(r**3) )/agru(i)
     224             :                dagrrd = (dvzd*dvrzd+dvrd*dvrrd+
     225           0 :      +              (dvpd*dvprd)/(r**2)-(dvpd**2)/(r**3) )/agrd(i)
     226             : c     grad(rho)*grad(abs(grad(rho)))
     227             :                gggrt(i) = dvzt*dagrzt + (dvpt*dagrpt)/(r**2) + dvrt*
     228           0 :      *              dagrrt
     229             :                gggru(i) = dvzu*dagrzu + (dvpu*dagrpu)/(r**2) + dvru*
     230           0 :      *              dagrru
     231             :                gggrd(i) = dvzd*dagrzd + (dvpd*dagrpd)/(r**2) + dvrd*
     232           0 :      *              dagrrd
     233             : c     dzdz=d(zeta)/dz,dzdp=d(zeta)/dp,dzdr=d(zeta)/dr
     234           0 :                dzdz = (dvzu-dvzd)/vlt - (vlu-vld)*dvzt/(vlt**2)
     235           0 :                dzdp = (dvpu-dvpd)/vlt - (vlu-vld)*dvpt/(vlt**2)
     236           0 :                dzdr = (dvru-dvrd)/vlt - (vlu-vld)*dvrt/(vlt**2)
     237             : c     gzgr=grad(zeta)*grad(ro)
     238           0 :                gzgr(i) = dzdz*dvzt + (dzdp*dvpt)/(r**2) + dzdr*dvrt
     239             : c     g2r: laplasian(rho)
     240           0 :                g2rt(i) = dvzzt + (dvppt)/(r**2) + dvrrt + dvrt/r
     241           0 :                g2ru(i) = dvzzu + (dvppu)/(r**2) + dvrru + dvru/r
     242           0 :                g2rd(i) = dvzzd + (dvppd)/(r**2) + dvrrd + dvrd/r
     243             :                
     244           0 :  20         end do
     245             :             
     246             :          end if
     247             :       
     248           0 :       END SUBROUTINE od_mkgxyz3
     249             :       END MODULE m_od_mkgxyz3

Generated by: LCOV version 1.13