LCOV - code coverage report
Current view: top level - vgen - mkgylm.f90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 100.0 % 189 189
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_mkgylm
       7              :    !-----------------------------------------------------------------------------
       8              :    ! Using the components and derivatives of a charge density rho in spherical
       9              :    ! coordinates on the real space grid, make the following quantaties:
      10              :    ! 
      11              :    ! gr(js):      [partial_r (rho),partial_theta (rho),partial_phi (rho)]
      12              :    ! sigma:       |grad(rho)|^2 for jspins==1, otherwise three components, namely
      13              :    !              |grad(rho_up)|^2,grad(rho_up)*grad(rho_down),|grad(rho_down)|^2 
      14              :    ! laplace(js): laplace(rho_js)
      15              :    ! 
      16              :    ! and these older components:
      17              :    ! agrt/u/d:    |grad(rho)| for total density/spin-up/spin-down
      18              :    ! g2rt/u/d:    laplace(rho)
      19              :    ! gggrt/u/d:   (grad(rho))*(grad(|grad(rho)|)) [scalar product]
      20              :    ! gzgr:        (grad(zeta))*(grad(rho)) for zeta=(rho_up-rho_down)/rho 
      21              :    ! 
      22              :    ! which are used to calculate gradient contributions to the xc potential and
      23              :    ! energy.
      24              :    ! 
      25              :    ! rh is rho, rhd[i][j] are the partial derivatives along one/two directions.
      26              :    ! rv and thet are the radius and polar angles and t/f in derivatives stand for
      27              :    ! theta/phi.
      28              :    ! 
      29              :    ! Modified so only allocated old quantities require calculations. A.N. 2019
      30              :    ! 
      31              :    ! Quantities fo libxc are calculated as well. D.W./M.R. 2018
      32              :    ! 
      33              :    ! Original script by T.A. 1996
      34              :    !-----------------------------------------------------------------------------
      35              : CONTAINS
      36       338811 :    SUBROUTINE mkgylm(jspins,rv,thet,nsp, rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff, rhdtf,rhdrt,rhdrf,grad,kt)
      37              :       USE m_types
      38              :       IMPLICIT NONE
      39              :       REAL, INTENT (IN)                :: rv
      40              :       REAL, INTENT (IN)                :: thet(:)
      41              :       REAL, INTENT (IN)                :: rh(:,:)
      42              :       REAL, INTENT (IN)                :: rhdr(:,:),rhdt(:,:),rhdf(:,:)
      43              :       REAL, INTENT (IN)                :: rhdrr(:,:),rhdtt(:,:),rhdff(:,:)
      44              :       REAL, INTENT (IN)                :: rhdtf(:,:),rhdrf(:,:),rhdrt(:,:)
      45              :       TYPE(t_gradients), INTENT(INOUT) :: grad
      46              :       INTEGER, INTENT(IN)              :: jspins,nsp
      47              :       INTEGER, INTENT(IN)              :: kt !index of first point to use in gradients
      48              : 
      49              :       REAL dagrf,dagrfd,dagrfu,dagrr,dagrrd,dagrru,dagrt, &
      50              :            dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr, &
      51              :            grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu, &
      52              :            rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu, &
      53              :            rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu, &
      54              :            ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1
      55              :       REAL, PARAMETER  :: chsml = 1.e-10
      56              :       INTEGER i,js,fac
      57              : 
      58       338811 :       rv1 = rv
      59       338811 :       rv2 = rv1**2
      60       338811 :       rv3 = rv1**3
      61              : 
      62              :       ! Gradients for libxc calculations.
      63       338811 :       IF (ALLOCATED(grad%gr)) THEN
      64        42252 :          DO js=1,jspins
      65      5521002 :             DO i=1,nsp
      66     21944601 :                grad%gr(:,kt+i,js)=[rhdr(i,js),rhdt(i,js)/rv1,rhdf(i,js)/(rv1*sin(thet(i)))]
      67              :             END DO
      68              :          END DO
      69              :          ! Use contracted gradients only for libxc.
      70        12651 :          IF (ALLOCATED(grad%sigma)) THEN
      71         4217 :             IF (jspins==1) THEN
      72       290277 :                DO i=1,nsp
      73      1158627 :                   grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
      74              :                END DO
      75              :             ELSE
      76       579690 :                DO i=1,nsp
      77      2305200 :                   grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
      78      2305200 :                   grad%sigma(2,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,2))
      79      2308590 :                   grad%sigma(3,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,2),grad%gr(:,kt+i,2))
      80              :                END DO
      81              :             END IF
      82              :          END IF
      83        12651 :          IF (ALLOCATED(grad%laplace)) THEN
      84              :             !Lapacian of density
      85              :             !fac=MERGE(2,1,jspins==1)
      86              :             fac=1
      87        11824 :             DO js=1,jspins
      88      1453874 :                DO i=1,nsp
      89              :                   grad%laplace(kt+i,js) = (rhdrr(i,js) + 2.e0*rhdr(i,js)/rv1 +&
      90      1449657 :                                            (rhdtt(i,js)+rhdt(i,js)/TAN(thet(i))+rhdff(i,js)/SIN(thet(i))**2)/rv2)/fac
      91              :                ENDDO
      92              :             ENDDO
      93              :          ENDIF
      94              :          ! Cartesian components of the gradient for sourcefree calculations. TODO: Need phi here as well.
      95              :          !IF (ALLOCATED(grad%grxyz)) THEN
      96              :          !   grad%grxyz=SIN(th)*COS(ph)*grad%gr(1,kt+k,1) + COS(th)*COS(ph)*grad%gr(2,kt+k,1)/r - SIN(ph)*grad%gr(3,kt+k,1)/(r*SIN(th))
      97              :          !END IF
      98        12651 :          RETURN ! Do not calculate arrays for in-build GGA.
      99              :       END IF
     100              : 
     101              :       IF (ALLOCATED(grad%gr)) THEN
     102              :          
     103              :       END IF
     104              :       !     Old code for in-build xcpots
     105       326160 :       IF(allocated(grad%agrt)) THEN
     106       326160 :          IF (jspins==1) THEN
     107              : 
     108     34433706 :             points_1 : DO i = 1,nsp
     109              : 
     110     34262388 :                grad%agrt(kt+i)  = 0.0
     111     34262388 :                grad%agru(kt+i)  = 0.0
     112     34262388 :                grad%agrd(kt+i)  = 0.0
     113     34262388 :                grad%g2rt(kt+i)  = 0.0
     114     34262388 :                grad%g2ru(kt+i)  = 0.0
     115     34262388 :                grad%g2rd(kt+i)  = 0.0
     116     34262388 :                grad%gggrt(kt+i) = 0.0
     117     34262388 :                grad%gggru(kt+i) = 0.0
     118     34262388 :                grad%gggrd(kt+i) = 0.0
     119     34262388 :                grad%grgru(kt+i) = 0.0
     120     34262388 :                grad%grgrd(kt+i) = 0.0
     121     34262388 :                grad%gzgr(kt+i)  = 0.0
     122              : 
     123     34262388 :                ro = rh(i,1)
     124              : 
     125     34262388 :                IF (ro<chsml) CYCLE points_1
     126     34262388 :                sint1 = sin(thet(i))
     127     34262388 :                sint2 = sint1**2
     128     34262388 :                tant1 = tan(thet(i))
     129     34262388 :                rvsin1 = rv1*sint1
     130              : 
     131     34262388 :                rou   = ro/2
     132     34262388 :                rdru  = rhdr(i,1)/2
     133     34262388 :                rdtu  = rhdt(i,1)/2
     134     34262388 :                rdfu  = rhdf(i,1)/2
     135     34262388 :                rdrru = rhdrr(i,1)/2
     136     34262388 :                rdttu = rhdtt(i,1)/2
     137     34262388 :                rdffu = rhdff(i,1)/2
     138     34262388 :                rdtfu = rhdtf(i,1)/2
     139     34262388 :                rdrtu = rhdrt(i,1)/2
     140     34262388 :                rdrfu = rhdrf(i,1)/2
     141              : 
     142     34262388 :                rod   = rou
     143     34262388 :                rdrd  = rdru
     144     34262388 :                rdtd  = rdtu
     145     34262388 :                rdfd  = rdfu
     146     34262388 :                rdrrd = rdrru
     147     34262388 :                rdttd = rdttu
     148     34262388 :                rdffd = rdffu
     149     34262388 :                rdtfd = rdtfu
     150     34262388 :                rdrtd = rdrtu
     151     34262388 :                rdrfd = rdrfu
     152              : 
     153     34262388 :                rdr  = rdru  + rdrd
     154     34262388 :                rdt  = rdtu  + rdtd
     155     34262388 :                rdf  = rdfu  + rdfd
     156     34262388 :                drdr = rdrru + rdrrd
     157     34262388 :                rdtt = rdttu + rdttd
     158     34262388 :                rdff = rdffu + rdffd
     159     34262388 :                rdtf = rdtfu + rdtfd
     160     34262388 :                rdrt = rdrtu + rdrtd
     161     34262388 :                rdrf = rdrfu + rdrfd
     162              : 
     163     34262388 :                ro2 = ro**2
     164              : 
     165     34262388 :                grr = rdr
     166     34262388 :                grt = rdt/rv1
     167     34262388 :                grf = rdf/rvsin1
     168              : 
     169     34262388 :                grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
     170              : 
     171     34262388 :                IF (grad%agrt(kt+i)<chsml) CYCLE points_1
     172              : 
     173              :                dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ &
     174     34262388 :                         rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
     175              : 
     176              :                dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/&
     177     34262388 :                        (grad%agrt(kt+i)*rv3)
     178              : 
     179              :                dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/&
     180     34262388 :                        (grad%agrt(kt+i)*rv3*sint1)
     181              : 
     182              :                grad%g2rt(kt+i)=drdr+2.0*rdr/rv1+&
     183     34262388 :                                 (rdtt+rdt/tant1+rdff/sint2)/rv2
     184              : 
     185     34262388 :                dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
     186              : 
     187              :                !--   >      dzdtr,dzdfs vanish by definition.
     188              : 
     189     34262388 :                dzdtr = 0.0
     190     34262388 :                dzdfs = 0.0
     191              : 
     192     34262388 :                grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
     193              : 
     194     34262388 :                grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
     195              : 
     196     34262388 :                grru = rdru
     197     34262388 :                grtu = rdtu/rv1
     198     34262388 :                grfu = rdfu/rvsin1
     199              : 
     200     34262388 :                grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
     201              : 
     202              :                dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
     203     34262388 :                          rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
     204              : 
     205              :                dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+&
     206     34262388 :                          rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
     207              : 
     208              :                dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/&
     209     34262388 :                         (grad%agru(kt+i)*rv3*sint1)
     210              : 
     211              :                grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 +&
     212     34262388 :                                  (rdttu+rdtu/tant1+rdffu/sint2)/rv2
     213              : 
     214     34262388 :                grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
     215              : 
     216     34262388 :                grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
     217              : 
     218     34262388 :                grrd = rdrd
     219     34262388 :                grtd = rdtd/rv1
     220     34262388 :                grfd = rdfd/rvsin1
     221              : 
     222     34262388 :                grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
     223              : 
     224              :                dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+&
     225     34262388 :                          rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
     226              : 
     227              :                dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+&
     228     34262388 :                          rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
     229              : 
     230              :                dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/&
     231     34262388 :                         (grad%agrd(kt+i)*rv3*sint1)
     232              : 
     233              :                grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 +&
     234     34262388 :                                  (rdttd+rdtd/tant1+rdffd/sint2)/rv2
     235              : 
     236     34262388 :                grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
     237              : 
     238     34433706 :                grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
     239              : 
     240              :             ENDDO points_1
     241              : 
     242              :          ELSE
     243              : 
     244     29781844 :             points_2 : DO i = 1,nsp
     245              : 
     246     29627002 :                grad%agrt(kt+i)  = 0.0
     247     29627002 :                grad%agru(kt+i)  = 0.0
     248     29627002 :                grad%agrd(kt+i)  = 0.0
     249     29627002 :                grad%g2rt(kt+i)  = 0.0
     250     29627002 :                grad%g2ru(kt+i)  = 0.0
     251     29627002 :                grad%g2rd(kt+i)  = 0.0
     252     29627002 :                grad%gggrt(kt+i) = 0.0
     253     29627002 :                grad%gggru(kt+i) = 0.0
     254     29627002 :                grad%gggrd(kt+i) = 0.0
     255     29627002 :                grad%grgru(kt+i) = 0.0
     256     29627002 :                grad%grgrd(kt+i) = 0.0
     257     29627002 :                grad%gzgr(kt+i)  = 0.0
     258              : 
     259     29627002 :                ro = rh(i,1) + rh(i,jspins)
     260              : 
     261     29627002 :                IF (ro<chsml) CYCLE points_2
     262              : 
     263     29627002 :                sint1 = sin(thet(i))
     264     29627002 :                sint2 = sint1**2
     265     29627002 :                tant1 = tan(thet(i))
     266     29627002 :                rvsin1 = rv1*sint1
     267              : 
     268     29627002 :                rou   = rh(i,1)
     269     29627002 :                rdru  = rhdr(i,1)
     270     29627002 :                rdtu  = rhdt(i,1)
     271     29627002 :                rdfu  = rhdf(i,1)
     272     29627002 :                rdrru = rhdrr(i,1)
     273     29627002 :                rdttu = rhdtt(i,1)
     274     29627002 :                rdffu = rhdff(i,1)
     275     29627002 :                rdtfu = rhdtf(i,1)
     276     29627002 :                rdrtu = rhdrt(i,1)
     277     29627002 :                rdrfu = rhdrf(i,1)
     278              : 
     279     29627002 :                rod   = rh(i,jspins)
     280     29627002 :                rdrd  = rhdr(i,jspins)
     281     29627002 :                rdtd  = rhdt(i,jspins)
     282     29627002 :                rdfd  = rhdf(i,jspins)
     283     29627002 :                rdrrd = rhdrr(i,jspins)
     284     29627002 :                rdttd = rhdtt(i,jspins)
     285     29627002 :                rdffd = rhdff(i,jspins)
     286     29627002 :                rdtfd = rhdtf(i,jspins)
     287     29627002 :                rdrtd = rhdrt(i,jspins)
     288     29627002 :                rdrfd = rhdrf(i,jspins)
     289              : 
     290     29627002 :                rdr   = rdru  + rdrd
     291     29627002 :                rdt   = rdtu  + rdtd
     292     29627002 :                rdf   = rdfu  + rdfd
     293     29627002 :                drdr  = rdrru + rdrrd
     294     29627002 :                rdtt  = rdttu + rdttd
     295     29627002 :                rdff  = rdffu + rdffd
     296     29627002 :                rdtf  = rdtfu + rdtfd
     297     29627002 :                rdrt  = rdrtu + rdrtd
     298     29627002 :                rdrf  = rdrfu + rdrfd
     299              : 
     300     29627002 :                ro2 = ro**2
     301              : 
     302     29627002 :                grr = rdr
     303     29627002 :                grt = rdt/rv1
     304     29627002 :                grf = rdf/rvsin1
     305              : 
     306     29627002 :                grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
     307              : 
     308     29627002 :                IF (grad%agrt(kt+i)<chsml) CYCLE points_2
     309              : 
     310     29627002 :                dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
     311              : 
     312     29627002 :                dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ (grad%agrt(kt+i)*rv3)
     313              : 
     314     29627002 :                dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/(grad%agrt(kt+i)*rv3*sint1)
     315              : 
     316     29627002 :                grad%g2rt(kt+i)= drdr+2.0*rdr/rv1+(rdtt+rdt/tant1+rdff/sint2)/rv2
     317              : 
     318     29627002 :                dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
     319              : 
     320              :                !     dzdtr,dzdfs vanish by definition.
     321     29627002 :                dzdtr = 0.0
     322     29627002 :                dzdfs = 0.0
     323              : 
     324     29627002 :                grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
     325              : 
     326     29627002 :                grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
     327              : 
     328     29627002 :                grru = rdru
     329     29627002 :                grtu = rdtu/rv1
     330     29627002 :                grfu = rdfu/rvsin1
     331              : 
     332     29627002 :                grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
     333              : 
     334              :                dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
     335     29627002 :                          rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
     336              : 
     337     29627002 :                dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
     338              : 
     339     29627002 :                dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ (grad%agru(kt+i)*rv3*sint1)
     340              : 
     341     29627002 :                grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 + (rdttu+rdtu/tant1+rdffu/sint2)/rv2
     342              : 
     343     29627002 :                grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
     344              : 
     345     29627002 :                grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
     346              : 
     347     29627002 :                grrd = rdrd
     348     29627002 :                grtd = rdtd/rv1
     349     29627002 :                grfd = rdfd/rvsin1
     350              : 
     351     29627002 :                grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
     352              : 
     353     29627002 :                dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
     354              : 
     355     29627002 :                dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
     356              : 
     357     29627002 :                dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ (grad%agrd(kt+i)*rv3*sint1)
     358              : 
     359     29627002 :                grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 + (rdttd+rdtd/tant1+rdffd/sint2)/rv2
     360              : 
     361     29627002 :                grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
     362              : 
     363     29781844 :                grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
     364              : 
     365              :             ENDDO points_2
     366              : 
     367              :          ENDIF
     368              :       ENDIF
     369              :    END SUBROUTINE mkgylm
     370              : END MODULE m_mkgylm
        

Generated by: LCOV version 2.0-1