LCOV - code coverage report
Current view: top level - vgen - mkgylm.f90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 187 189 98.9 %
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_mkgylm
       7             : CONTAINS
       8      206921 :    SUBROUTINE mkgylm(jspins,rv,thet,nsp, rh,rhdr,rhdt,rhdf,rhdrr,rhdtt,rhdff, rhdtf,rhdrt,rhdrf,grad,kt)
       9             :       !c.....------------------------------------------------------------------
      10             :       !c     by use of charge density and its polar coord. gradient components
      11             :       !c     c    calculate agr and others used to evaluate gradient
      12             :       !c     c    contributions to potential and energy. t.a. 1996.
      13             :       !c.....------------------------------------------------------------------
      14             :       !c     ro=sum(ro*ylh), rdr=sum(drr*ylh), drdr=sum(ddrr*ylh),
      15             :       !c     c    rdt=sum(ro*ylht1), rdtt=sum(ro*ylht2), ...
      16             :       !c     c    rdf=sum(ro*ylhf1), rdff=sum(ro*ylhf2), ...
      17             :       !c     c    rdtf=sum(ro*ylhtf), rdff=sum(ro*ylhf2), ...
      18             :       !c     c    rdrt=sum(drr*ylht1),rdrf=sum(drr*ylhf1),!
      19             :       !c     agr: abs(grad(ro)), g2r: laplacian(ro),
      20             :       !c     c    gggr: grad(ro)*grad(agr),
      21             :       !c     c    grgru,d: grad(ro)*grad(rou),for rod., gzgr: grad(zeta)*grad(ro).
      22             :       !
      23             :       !c     dagrr,-t,-f: d(agr)/dr, d(agr)/dth/r, d(agr)/dfi/r/sint.
      24             :       !c.....------------------------------------------------------------------
      25             :       USE m_types
      26             :       IMPLICIT NONE
      27             :       !C     ..
      28             :       !C     .. Arguments ..
      29             :       REAL,    INTENT (IN) :: rv
      30             :       REAL,    INTENT (IN) :: thet(:)
      31             :       REAL,    INTENT (IN) :: rh(:,:),rhdr(:,:)
      32             :       REAL,    INTENT (IN) :: rhdf(:,:),rhdrr(:,:)
      33             :       REAL,    INTENT (IN) :: rhdtt(:,:),rhdff(:,:)
      34             :       REAL,    INTENT (IN) :: rhdtf(:,:),rhdrt(:,:)
      35             :       REAL,    INTENT (IN) :: rhdrf(:,:),rhdt(:,:)
      36             :       TYPE(t_gradients),INTENT(INOUT) :: grad
      37             :       INTEGER,INTENT(IN)              :: jspins,nsp
      38             :       INTEGER,INTENT(IN)              :: kt !index of first point to use in gradients
      39             : 
      40             :       !C     ..
      41             :       !C     .. Locals ..
      42             :       INTEGER i,js,fac
      43             :       REAL    chsml, dagrf,dagrfd,dagrfu,dagrr,dagrrd,dagrru,dagrt,&
      44             :          dagrtd,dagrtu,drdr,dzdfs,dzdr,dzdtr,grf,grfd,grfu,grr,&
      45             :          grrd,grru,grt,grtd,grtu,rdf,rdfd,rdff,rdffd,rdffu,rdfu,&
      46             :          rdr,rdrd,rdrf,rdrfd,rdrfu,rdrrd,rdrru,rdrt,rdrtd,rdrtu,&
      47             :          rdru,rdt,rdtd,rdtf,rdtfd,rdtfu,rdtt,rdttd,rdttu,rdtu,&
      48             :          ro,ro2,rod,rou,rv1,rv2,rv3,rvsin1,sint1,sint2,tant1
      49             : 
      50      206921 :       chsml = 1.e-10
      51             : 
      52      206921 :       rv1 = rv
      53      206921 :       rv2 = rv1**2
      54      206921 :       rv3 = rv1**3
      55             : 
      56      206921 :       IF (ALLOCATED(grad%gr)) THEN
      57             :          !      Gradients for libxc
      58      158200 :          DO js=1,jspins
      59    19368200 :             DO i=1,nsp
      60    19323000 :                grad%gr(:,kt+i,js)=[rhdr(i,js),rhdt(i,js),rhdf(i,js)]
      61             :             ENDDO
      62             :          ENDDO
      63             :          !contracted gradients for libxc
      64       45200 :          IF (ALLOCATED(grad%sigma)) THEN
      65             :             !     Contracted gradients for libxc
      66       22600 :             IF (jspins==1) THEN
      67           0 :                DO i=1,nsp
      68           0 :                   grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
      69             :                ENDDO
      70             :             ELSE
      71     3864600 :                DO i=1,nsp
      72     3842000 :                   grad%sigma(1,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,1))
      73     3842000 :                   grad%sigma(2,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,1),grad%gr(:,kt+i,2))
      74    15390600 :                   grad%sigma(3,kt+i)= dot_PRODUCT(grad%gr(:,kt+i,2),grad%gr(:,kt+i,2))
      75             :                ENDDO
      76             :             ENDIF
      77             :          END IF
      78       45200 :          IF (ALLOCATED(grad%laplace)) THEN
      79             :             !Lapacian of density
      80             :             !fac=MERGE(2,1,jspins==1)
      81             :             fac=1
      82      113000 :             DO js=1,jspins
      83     7751800 :                DO i=1,nsp
      84             :                   grad%laplace(kt+i,js) = (rhdrr(i,js) + 2.e0*rhdr(i,js)/rv1 +&
      85     7729200 :                                            (rhdtt(i,js)+rhdt(i,js)/TAN(thet(i))+rhdff(i,js)/SIN(thet(i))**2)/rv2)/fac
      86             :                ENDDO
      87             :             ENDDO
      88             :          ENDIF
      89             :          RETURN !Do not calculate arrays for in-build GGA
      90             :       END IF
      91             :       !     Old code for in-build xcpots
      92      161721 :       IF(allocated(grad%agrt)) THEN
      93      161721 :          IF (jspins==1) THEN
      94             : 
      95     7766442 :             points_1 : DO i = 1,nsp
      96             : 
      97     7729260 :                grad%agrt(kt+i)  = 0.0
      98     7729260 :                grad%agru(kt+i)  = 0.0
      99     7729260 :                grad%agrd(kt+i)  = 0.0
     100     7729260 :                grad%g2rt(kt+i)  = 0.0
     101     7729260 :                grad%g2ru(kt+i)  = 0.0
     102     7729260 :                grad%g2rd(kt+i)  = 0.0
     103     7729260 :                grad%gggrt(kt+i) = 0.0
     104     7729260 :                grad%gggru(kt+i) = 0.0
     105     7729260 :                grad%gggrd(kt+i) = 0.0
     106     7729260 :                grad%grgru(kt+i) = 0.0
     107     7729260 :                grad%grgrd(kt+i) = 0.0
     108     7729260 :                grad%gzgr(kt+i)  = 0.0
     109             : 
     110     7729260 :                ro = rh(i,1)
     111             : 
     112     7729260 :                IF (ro<chsml) CYCLE points_1
     113     7729260 :                sint1 = sin(thet(i))
     114     7729260 :                sint2 = sint1**2
     115     7729260 :                tant1 = tan(thet(i))
     116     7729260 :                rvsin1 = rv1*sint1
     117             : 
     118     7729260 :                rou   = ro/2
     119     7729260 :                rdru  = rhdr(i,1)/2
     120     7729260 :                rdtu  = rhdt(i,1)/2
     121     7729260 :                rdfu  = rhdf(i,1)/2
     122     7729260 :                rdrru = rhdrr(i,1)/2
     123     7729260 :                rdttu = rhdtt(i,1)/2
     124     7729260 :                rdffu = rhdff(i,1)/2
     125     7729260 :                rdtfu = rhdtf(i,1)/2
     126     7729260 :                rdrtu = rhdrt(i,1)/2
     127     7729260 :                rdrfu = rhdrf(i,1)/2
     128             : 
     129     7729260 :                rod   = rou
     130     7729260 :                rdrd  = rdru
     131     7729260 :                rdtd  = rdtu
     132     7729260 :                rdfd  = rdfu
     133     7729260 :                rdrrd = rdrru
     134     7729260 :                rdttd = rdttu
     135     7729260 :                rdffd = rdffu
     136     7729260 :                rdtfd = rdtfu
     137     7729260 :                rdrtd = rdrtu
     138     7729260 :                rdrfd = rdrfu
     139             : 
     140     7729260 :                rdr  = rdru  + rdrd
     141     7729260 :                rdt  = rdtu  + rdtd
     142     7729260 :                rdf  = rdfu  + rdfd
     143     7729260 :                drdr = rdrru + rdrrd
     144     7729260 :                rdtt = rdttu + rdttd
     145     7729260 :                rdff = rdffu + rdffd
     146     7729260 :                rdtf = rdtfu + rdtfd
     147     7729260 :                rdrt = rdrtu + rdrtd
     148     7729260 :                rdrf = rdrfu + rdrfd
     149             : 
     150     7729260 :                ro2 = ro**2
     151             : 
     152     7729260 :                grr = rdr
     153     7729260 :                grt = rdt/rv1
     154     7729260 :                grf = rdf/rvsin1
     155             : 
     156     7729260 :                grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
     157             : 
     158     7729260 :                IF (grad%agrt(kt+i)<chsml) CYCLE points_1
     159             : 
     160             :                dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ &
     161     7729260 :                         rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
     162             : 
     163             :                dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/&
     164     7729260 :                        (grad%agrt(kt+i)*rv3)
     165             : 
     166             :                dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/&
     167     7729260 :                        (grad%agrt(kt+i)*rv3*sint1)
     168             : 
     169             :                grad%g2rt(kt+i)=drdr+2.0*rdr/rv1+&
     170     7729260 :                                 (rdtt+rdt/tant1+rdff/sint2)/rv2
     171             : 
     172     7729260 :                dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
     173             : 
     174             :                !--   >      dzdtr,dzdfs vanish by definition.
     175             : 
     176     7729260 :                dzdtr = 0.0
     177     7729260 :                dzdfs = 0.0
     178             : 
     179     7729260 :                grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
     180             : 
     181     7729260 :                grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
     182             : 
     183     7729260 :                grru = rdru
     184     7729260 :                grtu = rdtu/rv1
     185     7729260 :                grfu = rdfu/rvsin1
     186             : 
     187     7729260 :                grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
     188             : 
     189             :                dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
     190     7729260 :                          rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
     191             : 
     192             :                dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+&
     193     7729260 :                          rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
     194             : 
     195             :                dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/&
     196     7729260 :                         (grad%agru(kt+i)*rv3*sint1)
     197             : 
     198             :                grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 +&
     199     7729260 :                                  (rdttu+rdtu/tant1+rdffu/sint2)/rv2
     200             : 
     201     7729260 :                grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
     202             : 
     203     7729260 :                grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
     204             : 
     205     7729260 :                grrd = rdrd
     206     7729260 :                grtd = rdtd/rv1
     207     7729260 :                grfd = rdfd/rvsin1
     208             : 
     209     7729260 :                grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
     210             : 
     211             :                dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+&
     212     7729260 :                          rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
     213             : 
     214             :                dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+&
     215     7729260 :                          rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
     216             : 
     217             :                dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/&
     218     7729260 :                         (grad%agrd(kt+i)*rv3*sint1)
     219             : 
     220             :                grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 +&
     221     7729260 :                                  (rdttd+rdtd/tant1+rdffd/sint2)/rv2
     222             : 
     223     7729260 :                grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
     224             : 
     225     7766442 :                grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
     226             : 
     227             :             ENDDO points_1
     228             : 
     229             :          ELSE
     230             : 
     231    21449427 :             points_2 : DO i = 1,nsp
     232             : 
     233    21324888 :                grad%agrt(kt+i)  = 0.0
     234    21324888 :                grad%agru(kt+i)  = 0.0
     235    21324888 :                grad%agrd(kt+i)  = 0.0
     236    21324888 :                grad%g2rt(kt+i)  = 0.0
     237    21324888 :                grad%g2ru(kt+i)  = 0.0
     238    21324888 :                grad%g2rd(kt+i)  = 0.0
     239    21324888 :                grad%gggrt(kt+i) = 0.0
     240    21324888 :                grad%gggru(kt+i) = 0.0
     241    21324888 :                grad%gggrd(kt+i) = 0.0
     242    21324888 :                grad%grgru(kt+i) = 0.0
     243    21324888 :                grad%grgrd(kt+i) = 0.0
     244    21324888 :                grad%gzgr(kt+i)  = 0.0
     245             : 
     246    21324888 :                ro = rh(i,1) + rh(i,jspins)
     247             : 
     248    21324888 :                IF (ro<chsml) CYCLE points_2
     249             : 
     250    21324888 :                sint1 = sin(thet(i))
     251    21324888 :                sint2 = sint1**2
     252    21324888 :                tant1 = tan(thet(i))
     253    21324888 :                rvsin1 = rv1*sint1
     254             : 
     255    21324888 :                rou   = rh(i,1)
     256    21324888 :                rdru  = rhdr(i,1)
     257    21324888 :                rdtu  = rhdt(i,1)
     258    21324888 :                rdfu  = rhdf(i,1)
     259    21324888 :                rdrru = rhdrr(i,1)
     260    21324888 :                rdttu = rhdtt(i,1)
     261    21324888 :                rdffu = rhdff(i,1)
     262    21324888 :                rdtfu = rhdtf(i,1)
     263    21324888 :                rdrtu = rhdrt(i,1)
     264    21324888 :                rdrfu = rhdrf(i,1)
     265             : 
     266    21324888 :                rod   = rh(i,jspins)
     267    21324888 :                rdrd  = rhdr(i,jspins)
     268    21324888 :                rdtd  = rhdt(i,jspins)
     269    21324888 :                rdfd  = rhdf(i,jspins)
     270    21324888 :                rdrrd = rhdrr(i,jspins)
     271    21324888 :                rdttd = rhdtt(i,jspins)
     272    21324888 :                rdffd = rhdff(i,jspins)
     273    21324888 :                rdtfd = rhdtf(i,jspins)
     274    21324888 :                rdrtd = rhdrt(i,jspins)
     275    21324888 :                rdrfd = rhdrf(i,jspins)
     276             : 
     277    21324888 :                rdr   = rdru  + rdrd
     278    21324888 :                rdt   = rdtu  + rdtd
     279    21324888 :                rdf   = rdfu  + rdfd
     280    21324888 :                drdr  = rdrru + rdrrd
     281    21324888 :                rdtt  = rdttu + rdttd
     282    21324888 :                rdff  = rdffu + rdffd
     283    21324888 :                rdtf  = rdtfu + rdtfd
     284    21324888 :                rdrt  = rdrtu + rdrtd
     285    21324888 :                rdrf  = rdrfu + rdrfd
     286             : 
     287    21324888 :                ro2 = ro**2
     288             : 
     289    21324888 :                grr = rdr
     290    21324888 :                grt = rdt/rv1
     291    21324888 :                grf = rdf/rvsin1
     292             : 
     293    21324888 :                grad%agrt(kt+i) = sqrt(grr**2+grt**2+grf**2)
     294             : 
     295    21324888 :                IF (grad%agrt(kt+i)<chsml) CYCLE points_2
     296             : 
     297    21324888 :                dagrr = (rdr*drdr*rv3+rdt* (rdrt*rv1-rdt)+ rdf* (rdrf*rv1-rdf)/sint2)/grad%agrt(kt+i)/rv3
     298             : 
     299    21324888 :                dagrt =(rdr*rdrt*rv2+rdt*rdtt+rdf* (-rdf/tant1+rdtf)/sint2)/ (grad%agrt(kt+i)*rv3)
     300             : 
     301    21324888 :                dagrf = (rdr*rdrf*rv2+rdt*rdtf+rdf*rdff/sint2)/(grad%agrt(kt+i)*rv3*sint1)
     302             : 
     303    21324888 :                grad%g2rt(kt+i)= drdr+2.0*rdr/rv1+(rdtt+rdt/tant1+rdff/sint2)/rv2
     304             : 
     305    21324888 :                dzdr = ((rdru-rdrd)*ro- (rou-rod)*rdr)/ro2
     306             : 
     307             :                !     dzdtr,dzdfs vanish by definition.
     308    21324888 :                dzdtr = 0.0
     309    21324888 :                dzdfs = 0.0
     310             : 
     311    21324888 :                grad%gggrt(kt+i) = grr*dagrr + grt*dagrt + grf*dagrf
     312             : 
     313    21324888 :                grad%gzgr(kt+i) = dzdr*grr + dzdtr*grt + dzdfs*grf
     314             : 
     315    21324888 :                grru = rdru
     316    21324888 :                grtu = rdtu/rv1
     317    21324888 :                grfu = rdfu/rvsin1
     318             : 
     319    21324888 :                grad%agru(kt+i) = sqrt(grru**2+grtu**2+grfu**2)
     320             : 
     321             :                dagrru = (rdru*rdrru*rv3+rdtu* (rdrtu*rv1-rdtu)+&
     322    21324888 :                          rdfu* (rdrfu*rv1-rdfu)/sint2)/grad%agru(kt+i)/rv3
     323             : 
     324    21324888 :                dagrtu = (rdru*rdrtu*rv2+rdtu*rdttu+ rdfu* (-rdfu/tant1+rdtfu)/sint2)/ (grad%agru(kt+i)*rv3)
     325             : 
     326    21324888 :                dagrfu = (rdru*rdrfu*rv2+rdtu*rdtfu+rdfu*rdffu/sint2)/ (grad%agru(kt+i)*rv3*sint1)
     327             : 
     328    21324888 :                grad%g2ru(kt+i) = rdrru + 2.e0*rdru/rv1 + (rdttu+rdtu/tant1+rdffu/sint2)/rv2
     329             : 
     330    21324888 :                grad%gggru(kt+i) = grru*dagrru + grtu*dagrtu + grfu*dagrfu
     331             : 
     332    21324888 :                grad%grgru(kt+i) = grr*grru + grt*grtu + grf*grfu
     333             : 
     334    21324888 :                grrd = rdrd
     335    21324888 :                grtd = rdtd/rv1
     336    21324888 :                grfd = rdfd/rvsin1
     337             : 
     338    21324888 :                grad%agrd(kt+i) = sqrt(grrd**2+grtd**2+grfd**2)
     339             : 
     340    21324888 :                dagrrd = (rdrd*rdrrd*rv3+rdtd* (rdrtd*rv1-rdtd)+ rdfd* (rdrfd*rv1-rdfd)/sint2)/grad%agrd(kt+i)/rv3
     341             : 
     342    21324888 :                dagrtd = (rdrd*rdrtd*rv2+rdtd*rdttd+ rdfd* (-rdfd/tant1+rdtfd)/sint2)/ (grad%agrd(kt+i)*rv3)
     343             : 
     344    21324888 :                dagrfd = (rdrd*rdrfd*rv2+rdtd*rdtfd+rdfd*rdffd/sint2)/ (grad%agrd(kt+i)*rv3*sint1)
     345             : 
     346    21324888 :                grad%g2rd(kt+i) = rdrrd + 2*rdrd/rv1 + (rdttd+rdtd/tant1+rdffd/sint2)/rv2
     347             : 
     348    21324888 :                grad%gggrd(kt+i) = grrd*dagrrd + grtd*dagrtd + grfd*dagrfd
     349             : 
     350    21449427 :                grad%grgrd(kt+i) = grr*grrd + grt*grtd + grf*grfd
     351             : 
     352             :             ENDDO points_2
     353             : 
     354             :          ENDIF
     355             :       ENDIF
     356             :    END SUBROUTINE mkgylm
     357             : END MODULE m_mkgylm

Generated by: LCOV version 1.13