LCOV - code coverage report
Current view: top level - vgen - mt_tofrom_grid.F90 (source / functions) Coverage Total Hit
Test: FLEUR test coverage Lines: 92.0 % 87 80
Test Date: 2025-06-14 04:34:23 Functions: 100.0 % 4 4

            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_mt_tofrom_grid
       7              :    USE m_types
       8              :    PRIVATE
       9              :    REAL, PARAMETER    :: d_15 = 1.e-15
      10              :    REAL, ALLOCATABLE :: ylh(:, :, :), ylht(:, :, :), ylhtt(:, :, :)
      11              :    REAL, ALLOCATABLE :: ylhf(:, :, :), ylhff(:, :, :), ylhtf(:, :, :)
      12              :    REAL, ALLOCATABLE :: wt(:), rx(:, :), thet(:), phi(:)
      13              :    PUBLIC :: init_mt_grid, mt_to_grid, mt_from_grid, finish_mt_grid
      14              : CONTAINS
      15          718 :    SUBROUTINE init_mt_grid(jspins, atoms, sphhar, dograds, sym, thout, phout)
      16              :       USE m_gaussp
      17              :       USE m_lhglptg
      18              :       USE m_lhglpts
      19              :       IMPLICIT NONE
      20              :       INTEGER, INTENT(IN)          :: jspins
      21              :       TYPE(t_atoms), INTENT(IN)    :: atoms
      22              :       TYPE(t_sphhar), INTENT(IN)   :: sphhar
      23              :       LOGICAL, INTENT(IN)          :: dograds
      24              :       TYPE(t_sym), INTENT(IN)      :: sym
      25              :       REAL, INTENT(OUT), OPTIONAL  :: thout(:)
      26              :       REAL, INTENT(OUT), OPTIONAL  :: phout(:)
      27              : 
      28          718 :       call timestart("init_mt_grid")
      29              :       ! generate nspd points on a sherical shell with radius 1.0
      30              :       ! angular mesh equidistant in phi,
      31              :       ! theta are zeros of the legendre polynomials
      32          718 :       ALLOCATE (wt(atoms%nsp()), rx(3, atoms%nsp()), thet(atoms%nsp()), phi(atoms%nsp()))
      33          718 :       CALL gaussp(atoms%lmaxd, rx, wt)
      34              :       ! generate the lattice harmonics on the angular mesh
      35         2154 :       ALLOCATE (ylh(atoms%nsp(), 0:sphhar%nlhd, sphhar%ntypsd))
      36          718 :       IF (dograds) THEN
      37         1656 :          ALLOCATE (ylht, MOLD=ylh)
      38         1104 :          ALLOCATE (ylhtt, MOLD=ylh)
      39         1104 :          ALLOCATE (ylhf, MOLD=ylh)
      40         1104 :          ALLOCATE (ylhff, MOLD=ylh)
      41         1104 :          ALLOCATE (ylhtf, MOLD=ylh)
      42              : 
      43              :          CALL lhglptg(sphhar, atoms, rx, atoms%nsp(), dograds, sym, &
      44          552 :                       ylh, thet, phi, ylht, ylhtt, ylhf, ylhff, ylhtf)
      45          552 :          IF (PRESENT(thout)) THEN
      46            0 :             thout=thet
      47            0 :             phout=phi
      48              :          END IF
      49              : 
      50              :       ELSE
      51          166 :          CALL lhglpts(sphhar, atoms, rx, atoms%nsp(), sym, ylh)
      52              :       END IF
      53          718 :       call timestop("init_mt_grid")
      54          718 :    END SUBROUTINE init_mt_grid
      55              : 
      56          678 :    SUBROUTINE mt_to_grid(dograds, jspins, atoms, sym,sphhar,rotch, den_mt, n, noco ,grad, ch)
      57              :       USE m_grdchlh
      58              :       USE m_mkgylm
      59              :       IMPLICIT NONE
      60              :       LOGICAL, INTENT(IN)          :: dograds
      61              :       TYPE(t_atoms), INTENT(IN)    :: atoms
      62              :       TYPE(t_sym), INTENT(IN)      :: sym
      63              :       LOGICAL, INTENT(IN)          :: rotch
      64              :       TYPE(t_sphhar), INTENT(IN)   :: sphhar
      65              :       REAL, INTENT(IN)             :: den_mt(:, 0:, :)
      66              :       INTEGER, INTENT(IN)          :: n, jspins
      67              :       REAL, INTENT(OUT), OPTIONAL  :: ch(:, :)
      68              :       TYPE(t_gradients), INTENT(INOUT):: grad
      69              :       TYPE(t_noco), INTENT(IN)     :: noco
      70              :       REAL                         :: dentot
      71              : 
      72              : 
      73              : 
      74              :       REAL    :: rho_11,rho_22,rho_21r,rho_21i,mx,my,mz,magmom
      75              :       REAL    :: rhotot,rho_up,rho_down
      76          678 :       REAL, ALLOCATABLE :: chlh(:, :, :), chlhdr(:, :, :), chlhdrr(:, :, :)
      77          678 :       REAL, ALLOCATABLE :: chdr(:, :), chdt(:, :), chdf(:, :), ch_tmp(:, :),ch_calc(:,:)
      78          678 :       REAL, ALLOCATABLE :: chdrr(:, :), chdtt(:, :), chdff(:, :), chdtf(:, :)
      79          678 :       REAL, ALLOCATABLE :: chdrt(:, :), chdrf(:, :)
      80          678 :       REAL, ALLOCATABLE :: drm(:,:), drrm(:,:), mm(:,:)
      81          678 :       REAL, ALLOCATABLE :: chlhtot(:),chlhdrtot(:),chlhdrrtot(:)
      82              :       INTEGER:: nd, lh, js, jr, kt, k, nsp,j,i,jspV
      83              : 
      84          678 :       call timestart("mt_to_grid")
      85              :       !This snippet is crucial to determine over which spins (Only diagonals in colinear case or also off diags in non colin case.)
      86         1945 :       IF (any(noco%l_unrestrictMT)) THEN
      87              :          jspV=4
      88              :       ELSE
      89          593 :          jspV=jspins
      90              :       END IF
      91              : 
      92          678 :       nd = sym%ntypsy(atoms%firstAtom(n))
      93          678 :       nsp = atoms%nsp()
      94              : 
      95              :       !General Allocations
      96         2034 :       ALLOCATE (chlh(atoms%jmtd, 0:sphhar%nlhd, jspV))
      97         2034 :       ALLOCATE (ch_tmp(nsp, jspV),ch_calc(nsp*atoms%jmtd, jspV))
      98              : 
      99              :       !Allocations in dograds case
     100          678 :       IF (dograds) THEN
     101            0 :          ALLOCATE (chdr(nsp, jspV), chdt(nsp, jspV), chdf(nsp, jspV), chdrr(nsp, jspV), &
     102            0 :                    chdtt(nsp, jspV), chdff(nsp, jspV), chdtf(nsp, jspV), chdrt(nsp, jspV), &
     103          501 :                    chdrf(nsp, jspV))
     104         1002 :          ALLOCATE (chlhdr(atoms%jmtd, 0:sphhar%nlhd, jspV))
     105         1002 :          ALLOCATE (chlhdrr(atoms%jmtd, 0:sphhar%nlhd, jspV))
     106              :       ENDIF
     107              : 
     108              :       !Allocations in mtNoco case
     109         1945 :       IF (any(noco%l_unrestrictMT)) THEN
     110              :          !General Noco Allocations
     111          170 :          ALLOCATE(mm(atoms%jmtd, 0:sphhar%nlhd))
     112              : 
     113              :          !Allocations in case one uses e.g. GGA with mtNoco
     114           85 :          IF (dograds) THEN
     115            0 :             ALLOCATE(drm(atoms%jmtd,0:sphhar%nlhd),drrm(atoms%jmtd, 0:sphhar%nlhd))
     116            0 :             ALLOCATE(chlhtot(0:sphhar%nlhd),chlhdrtot(0:sphhar%nlhd),chlhdrrtot(0:sphhar%nlhd))
     117              :          END IF
     118              :       END IF
     119              : 
     120              :       !Calc magnetization (This is necessary only in mtNoco case)
     121      5160625 :       IF(any(noco%l_unrestrictMT)) mm(:,:)=SQRT((0.5*(den_mt(:,:,1)-den_mt(:,:,2)))**2+4*den_mt(:,:,3)**2+4*den_mt(:,:,4)**2)
     122              : 
     123              : 
     124              :       !Loop to calculate chlh and necessary gradients (if needed)
     125        21402 :       DO lh = 0, sphhar%nlh(nd)
     126              :          !         calculates gradients of radial charge densities of l=> 0.
     127              :          !         rho*ylh/r**2 is charge density. chlh=rho/r**2.
     128              :          !         charge density=sum(chlh*ylh).
     129              :          !         chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr.
     130              : 
     131              :          !Scaling of the magnetic moments in the same way the charge density is scaled in chlh.
     132      5197007 :          IF(any(noco%l_unrestrictMT)) mm(:,lh)=0.5*mm(:,lh)/(atoms%rmsh(:, n)*atoms%rmsh(:, n))
     133              : 
     134        71616 :          DO js = 1, jspV
     135              :             chlh(1:atoms%jri(n), lh, js) = den_mt(1:atoms%jri(n), lh, js) / &
     136     37051960 :                                            (atoms%rmsh(1:atoms%jri(n), n)*atoms%rmsh(1:atoms%jri(n), n))
     137              : 
     138              :             !Necessary gradients
     139        70938 :             IF (dograds) THEN
     140              :                !Colinear case only needs radial derivatives of chlh
     141        15512 :                CALL grdchlh(atoms%dx(n), chlh(1:atoms%jri(n), lh, js),  chlhdr(1:, lh, js), chlhdrr(1:, lh,js),atoms%rmsh(:,n))
     142        41039 :                IF (any(noco%l_unrestrictMT)) THEN
     143              :                !Noco case also needs radial derivatives of mm
     144            0 :                   CALL grdchlh(atoms%dx(n), mm(:atoms%jri(n),lh),  drm(:,lh), drrm(:,lh), atoms%rmsh(:, n))
     145              :                END IF
     146              :             END IF
     147              :          END DO ! js
     148              :       END DO   ! lh
     149              : 
     150              :       !The following Loop maps chlh on the k-Grid using the lattice harmonics ylh
     151              :       !$OMP parallel do default( NONE ) &
     152              :       !$OMP SHARED(atoms,sphhar,noco,n,nsp,jspV,nd,ylh,chlh,dograds,chlhdr,chlhdrr,mm,drm,drrm)&
     153              :       !$OMP SHARED(ylht,ylhf,ylhtt,ylhff,ylhtf,jspins,thet,grad,ch,ch_calc) &
     154              :       !$OMP private(kt,ch_tmp,js,lh,k,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf) &
     155          678 :       !$OMP private(chlhtot,chlhdrtot,chlhdrrtot)
     156              :       DO jr = 1, atoms%jri(n)
     157              :          kt = (jr-1)*nsp
     158              :          ! charge density (on extended grid for all jr)
     159              :          ! following are at points on jr-th sphere.
     160              :          ch_tmp(:, :) = 0.0
     161              :          !  generate the densities on an angular mesh (ch_tmp is needed in mkgylm call later on)
     162              :          DO js = 1, jspV
     163              :             DO lh = 0, sphhar%nlh(nd)
     164              :                DO k = 1, nsp
     165              :                      ch_tmp(k, js) = ch_tmp(k, js) + ylh(k, lh, nd)*chlh(jr, lh, js)
     166              :                ENDDO
     167              :             ENDDO
     168              :          ENDDO
     169              : 
     170              :          !Initialize derivatives of ch on grid if needed.
     171              :          IF (dograds) THEN
     172              :             chdr(:, :) = 0.0     ! d(ch)/dr
     173              :             chdt(:, :) = 0.0     ! d(ch)/dtheta
     174              :             chdf(:, :) = 0.0     ! d(ch)/dfai
     175              :             chdrr(:, :) = 0.0     ! dd(ch)/drr
     176              :             chdtt(:, :) = 0.0     ! dd(ch)/dtt
     177              :             chdff(:, :) = 0.0     ! dd(ch)/dff
     178              :             chdtf(:, :) = 0.0     ! dd(ch)/dtf
     179              :             chdrt(:, :) = 0.0     ! d(d(ch)/dr)dt
     180              :             chdrf(:, :) = 0.0     ! d(d(ch)/dr)df
     181              :             !  generate the derivatives on an angular mesh
     182              :             DO js = 1, jspV
     183              :                DO lh = 0, sphhar%nlh(nd)
     184              : 
     185              :                   !The following snippet maps chlh and its radial derivatives on a colinear system
     186              :                   !using mm and its radial derivatives.
     187              :                   IF (any(noco%l_unrestrictMT)) THEN
     188              :                       IF (js.EQ.1) THEN
     189              :                          chlhtot(lh)=0.5*(chlh(jr, lh, 1)+chlh(jr, lh, 2))
     190              :                          chlhdrtot(lh)=0.5*(chlhdr(jr, lh, 1)+chlhdr(jr, lh, 2))
     191              :                          chlhdrrtot(lh)=0.5*(chlhdrr(jr, lh, 1)+chlhdrr(jr, lh, 2))
     192              :                          chlh(jr,lh,js)=chlhtot(lh)+mm(jr,lh)
     193              :                          chlhdr(jr, lh, js)=chlhdrtot(lh)+drm(jr,lh)
     194              :                          chlhdrr(jr, lh, js)=chlhdrrtot(lh)+drrm(jr,lh)
     195              :                       ELSE IF (js.EQ.2) THEN
     196              :                          chlh(jr,lh,js)=chlhtot(lh)-mm(jr,lh)
     197              :                          chlhdr(jr, lh, js)=chlhdrtot(lh)-drm(jr,lh)
     198              :                          chlhdrr(jr, lh, js)=chlhdrrtot(lh)-drrm(jr,lh)
     199              :                       ELSE
     200              :                          chlh(jr,lh,js)=0
     201              :                          chlhdr(jr, lh, js)=0
     202              :                          chlhdrr(jr, lh, js)=0
     203              :                      END IF
     204              :                   END IF
     205              : 
     206              :                   !The following loop brings chlhdr and chlhdrr on the k-grid.
     207              :                   DO k = 1, nsp
     208              :                      chdr(k, js) = chdr(k, js) + ylh(k, lh, nd)*chlhdr(jr, lh, js)
     209              :                      chdrr(k, js) = chdrr(k, js) + ylh(k, lh, nd)*chlhdrr(jr, lh, js)
     210              :                   ENDDO
     211              : 
     212              :                   !This loop calculates the other derviatives of ch (Angular terms) on the k-grid
     213              :                   !by using the lattice harmonics derivatives and chlh with its derivatives.
     214              :                   DO k = 1, nsp
     215              :                      chdrt(k, js) = chdrt(k, js) + ylht(k, lh, nd)*chlhdr(jr, lh, js)
     216              :                      chdrf(k, js) = chdrf(k, js) + ylhf(k, lh, nd)*chlhdr(jr, lh, js)
     217              :                      chdt(k, js) = chdt(k, js) + ylht(k, lh, nd)*chlh(jr, lh, js)
     218              :                      chdf(k, js) = chdf(k, js) + ylhf(k, lh, nd)*chlh(jr, lh, js)
     219              :                      chdtt(k, js) = chdtt(k, js) + ylhtt(k, lh, nd)*chlh(jr, lh, js)
     220              :                      chdff(k, js) = chdff(k, js) + ylhff(k, lh, nd)*chlh(jr, lh, js)
     221              :                      chdtf(k, js) = chdtf(k, js) + ylhtf(k, lh, nd)*chlh(jr, lh, js)
     222              :                   ENDDO
     223              :                ENDDO ! lh
     224              :             ENDDO   ! js
     225              :          !Rotation to local if needed (Indicated by rotch)
     226              :             !Makegradients
     227              :             !IF(jspins>2) CALL mkgylm(2, atoms%rmsh(jr, n), thet, nsp, &
     228              :             !            ch_tmp, chdr, chdt, chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt)
     229              :             !IF(jspins.LE.2)
     230              :             CALL mkgylm(jspins, atoms%rmsh(jr, n), thet, nsp, &
     231              :                         ch_tmp, chdr, chdt, chdf, chdrr, chdtt, chdff, chdtf, chdrt, chdrf, grad, kt)
     232              :          END IF
     233              :          !Set charge to minimum value
     234              :          IF (PRESENT(ch)) THEN
     235              :             WHERE (ABS(ch_tmp(:nsp,:)) < d_15) ch_tmp(:nsp,:) = d_15
     236              :             ch_calc(kt + 1:kt + nsp, :) = ch_tmp(:nsp, :)
     237              :          ENDIF
     238              :       END DO
     239              :       !$OMP END PARALLEL DO
     240              : 
     241              : 
     242          678 :       IF (PRESENT(ch)) THEN
     243              :       !Rotation to local if needed (Indicated by rotch)
     244         1905 :          IF (rotch.AND.any(noco%l_unrestrictMT)) THEN
     245      2247301 :             DO jr = 1,nsp*atoms%jri(n)
     246      2247284 :                rho_11  = ch_calc(jr,1)
     247      2247284 :                rho_22  = ch_calc(jr,2)
     248      2247284 :                rho_21r = ch_calc(jr,3)
     249      2247284 :                rho_21i = ch_calc(jr,4)
     250      2247284 :                mx      =  2*rho_21r
     251      2247284 :                my      = -2*rho_21i
     252      2247284 :                mz      = (rho_11-rho_22)
     253      2247284 :                magmom  = SQRT(mx**2 + my**2 + mz**2)
     254      2247284 :                rhotot  = rho_11 + rho_22
     255      2247284 :                rho_up  = (rhotot + magmom)/2
     256      2247284 :                rho_down= (rhotot - magmom)/2
     257      2247284 :                ch(jr,1) = rho_up
     258      2247301 :                ch(jr,2) = rho_down
     259              :             END DO
     260              :          ELSE
     261    148448480 :             ch(:nsp*atoms%jri(n),1:jspins)=ch_calc(:nsp*atoms%jri(n),1:jspins)
     262              : 
     263              :          END IF
     264              :       END IF
     265          678 :       call timestop("mt_to_grid")
     266          678 :    END SUBROUTINE mt_to_grid
     267              : 
     268         2466 :    SUBROUTINE mt_from_grid(atoms, sym, sphhar, n, jspins, v_in, vr)
     269              :       IMPLICIT NONE
     270              :       TYPE(t_atoms), INTENT(IN) :: atoms
     271              :       TYPE(t_sym), INTENT(IN)   :: sym
     272              :       TYPE(t_sphhar), INTENT(IN):: sphhar
     273              :       INTEGER, INTENT(IN)       :: jspins, n
     274              :       REAL, INTENT(IN)          :: v_in(:, :)
     275              :       REAL, INTENT(INOUT)       :: vr(:, 0:, :)
     276              : 
     277         2466 :       REAL    :: vpot(atoms%nsp()), vlh
     278              :       INTEGER :: js, kt, lh, jr, nd, nsp
     279              : 
     280         2466 :       call timestart("mt_from_grid")
     281              : 
     282         2466 :       nsp = atoms%nsp()
     283         2466 :       nd = sym%ntypsy(atoms%firstAtom(n))
     284              : 
     285         6103 :       DO js = 1, jspins
     286              :          !$OMP parallel do default( NONE ) &
     287              :          !$OMP SHARED(atoms,sphhar,n,v_in,nsp,wt,nd,ylh,vr,js)&
     288         6103 :          !$OMP private(vpot,kt,lh,vlh)
     289              :          DO jr = 1, atoms%jri(n)
     290              :             kt = (jr-1)*nsp
     291              :             vpot = v_in(kt + 1:kt + nsp, js)*wt(:)!  multiplicate v_in with the weights of the k-points
     292              : 
     293              :             DO lh = 0, sphhar%nlh(nd)
     294              :                !
     295              :                ! --->        determine the corresponding potential number
     296              :                !c            through gauss integration
     297              :                !
     298              :                vlh = dot_PRODUCT(vpot(:), ylh(:nsp, lh, nd))
     299              :                vr(jr, lh, js) = vr(jr, lh, js) + vlh
     300              :             ENDDO ! lh
     301              :          ENDDO   ! jr
     302              :          !$OMP END PARALLEL DO
     303              :       ENDDO
     304         2466 :       call timestop("mt_from_grid")
     305         2466 :    END SUBROUTINE mt_from_grid
     306              : 
     307          718 :    SUBROUTINE finish_mt_grid()
     308              :       implicit NONE 
     309          718 :       call timestart("finish_mt_grid")
     310          718 :       DEALLOCATE (ylh, wt, rx, thet, phi)
     311          718 :       IF (ALLOCATED(ylht)) DEALLOCATE (ylht, ylhtt, ylhf, ylhff, ylhtf)
     312          718 :       call timestop("finish_mt_grid")
     313          718 :    END SUBROUTINE finish_mt_grid
     314              : 
     315              : END MODULE m_mt_tofrom_grid
        

Generated by: LCOV version 2.0-1