mt_tofrom_grid.F90 7.02 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------
MODULE m_mt_tofrom_grid
  USE m_types
  PRIVATE
  REAL,PARAMETER    :: d_15 = 1.e-15
  INTEGER,PARAMETER :: ndvgrd=6 ! this should be consistent across GGA derivative routines
  REAL, ALLOCATABLE :: ylh(:,:,:),ylht(:,:,:),ylhtt(:,:,:)
  REAL, ALLOCATABLE :: ylhf(:,:,:),ylhff(:,:,:),ylhtf(:,:,:)
  REAL, ALLOCATABLE :: wt(:),rx(:,:),thet(:)
  PUBLIC :: init_mt_grid,mt_to_grid,mt_from_grid,finish_mt_grid
CONTAINS
16
  SUBROUTINE init_mt_grid(nsp,jspins,atoms,sphhar,xcpot,sym)
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
    USE m_gaussp
    USE m_lhglptg
    USE m_lhglpts
    IMPLICIT NONE
    INTEGER,INTENT(IN)          :: nsp,jspins
    TYPE(t_atoms),INTENT(IN)    :: atoms
    TYPE(t_sphhar),INTENT(IN)   :: sphhar
    CLASS(t_xcpot),INTENT(IN)   :: xcpot
    TYPE(t_sym),INTENT(IN)      :: sym
    
    ! generate nspd points on a sherical shell with radius 1.0
    ! angular mesh equidistant in phi,
    ! theta are zeros of the legendre polynomials
    ALLOCATE(wt(nsp),rx(3,nsp),thet(nsp))
    CALL gaussp(atoms%lmaxd, rx,wt)
    ! generate the lattice harmonics on the angular mesh
    ALLOCATE ( ylh(nsp,0:sphhar%nlhd,sphhar%ntypsd))
34 35 36 37 38 39 40 41 42
    IF (xcpot%needs_grad()) THEN
      ALLOCATE(ylht,MOLD=ylh )
      ALLOCATE(ylhtt,MOLD=ylh )
      ALLOCATE(ylhf,MOLD=ylh )
      ALLOCATE(ylhff,MOLD=ylh )
      ALLOCATE(ylhtf,MOLD=ylh )
      
      CALL lhglptg(sphhar,atoms,rx,nsp,xcpot,sym,&
           ylh,thet,ylht,ylhtt,ylhf,ylhff,ylhtf)
43 44 45 46 47
    ELSE
       CALL lhglpts( sphhar,atoms, rx,nsp, sym, ylh)
    END IF
  END SUBROUTINE init_mt_grid
  
48 49
  SUBROUTINE mt_to_grid(xcpot,jspins,atoms,sphhar,den_mt,nsp,n,grad,ch)
!  SUBROUTINE pw_to_grid(xcpot,jspins,l_noco,stars,cell,den_pw,grad,rho)
50 51 52
    USE m_grdchlh
    USE m_mkgylm
    IMPLICIT NONE
53
    CLASS(t_xcpot),INTENT(IN)   :: xcpot
54 55
    TYPE(t_atoms),INTENT(IN)    :: atoms
    TYPE(t_sphhar),INTENT(IN)   :: sphhar
56
    REAL,INTENT(IN)             :: den_mt(:,0:,:)
57
    INTEGER,INTENT(IN)          :: n,jspins,nsp
58 59
    REAL,INTENT(OUT),OPTIONAL      :: ch(:,:)
    TYPE(t_gradients),INTENT(INOUT):: grad
60 61
    
    REAL, ALLOCATABLE :: chlh(:,:,:),chlhdr(:,:,:),chlhdrr(:,:,:)
62
    REAL, ALLOCATABLE :: chdr(:,:),chdt(:,:),chdf(:,:),ch_tmp(:,:)
63 64 65 66 67 68 69
    REAL, ALLOCATABLE :: chdrr(:,:),chdtt(:,:),chdff(:,:),chdtf(:,:)
    REAL, ALLOCATABLE :: chdrt(:,:),chdrf(:,:)
    INTEGER:: nd,lh,js,jr,kt,k

    nd = atoms%ntypsy(SUM(atoms%neq(:n-1))+1)

    ALLOCATE ( chlh(atoms%jmtd,0:sphhar%nlhd,jspins))
70
    ALLOCATE ( ch_tmp(nsp,jspins) )
71
    IF (xcpot%needs_grad())  THEN
72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
       ALLOCATE(chdr(nsp,jspins),chdt(nsp,jspins),chdf(nsp,jspins),chdrr(nsp,jspins),&
            chdtt(nsp,jspins),chdff(nsp,jspins),chdtf(nsp,jspins),chdrt(nsp,jspins),&
            chdrf(nsp,jspins) )
       ALLOCATE (chlhdr(atoms%jmtd,0:sphhar%nlhd,jspins))
       ALLOCATE (chlhdrr(atoms%jmtd,0:sphhar%nlhd,jspins))
    ENDIF
    
    DO lh = 0,sphhar%nlh(nd)

       !         calculates gradients of radial charge densities of l=> 0.
       !         rho*ylh/r**2 is charge density. chlh=rho/r**2.
       !         charge density=sum(chlh*ylh).
       !         chlhdr=d(chlh)/dr, chlhdrr=dd(chlh)/drr.
       
       DO js = 1,jspins      
          DO jr = 1,atoms%jri(n)
88
             chlh(jr,lh,js) = den_mt(jr,lh,js)/(atoms%rmsh(jr,n)*atoms%rmsh(jr,n))
89
          ENDDO
90
          IF (xcpot%needs_grad()) CALL grdchlh(1,1,atoms%jri(n),atoms%dx(n),atoms%rmsh(1,n),&
91 92 93 94
               chlh(1,lh,js),ndvgrd, chlhdr(1,lh,js),chlhdrr(1,lh,js))
          
       ENDDO ! js
    ENDDO   ! lh
95
    
96 97
    kt=0
    DO jr = 1,atoms%jri(n)
98
       ch_tmp(:,:)    = 0.0     ! charge density (on extended grid for all jr)
99 100 101 102 103
       !         following are at points on jr-th sphere.
       !  generate the densities on an angular mesh
       DO js = 1,jspins
          DO lh = 0,sphhar%nlh(nd)
             DO k = 1,nsp
104
                ch_tmp(k,js) = ch_tmp(k,js) + ylh(k,lh,nd)*chlh(jr,lh,js)
105 106 107
             ENDDO
          ENDDO
       ENDDO
108
       IF (xcpot%needs_grad()) THEN
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
          chdr(:,:)  = 0.0     ! d(ch)/dr
          chdt(:,:)  = 0.0     ! d(ch)/dtheta
          chdf(:,:)  = 0.0     ! d(ch)/dfai
          chdrr(:,:) = 0.0     ! dd(ch)/drr
          chdtt(:,:) = 0.0     ! dd(ch)/dtt
          chdff(:,:) = 0.0     ! dd(ch)/dff
          chdtf(:,:) = 0.0     ! dd(ch)/dtf
          chdrt(:,:) = 0.0     ! d(d(ch)/dr)dt
          chdrf(:,:) = 0.0     ! d(d(ch)/dr)df
          !  generate the derivatives on an angular mesh
          DO js = 1,jspins
             DO lh = 0,sphhar%nlh(nd)
                ! 
                DO k = 1,nsp
                   chdr(k,js) =chdr(k,js)+ ylh(k,lh,nd)*chlhdr(jr,lh,js)
                   chdrr(k,js)=chdrr(k,js)+ylh(k,lh,nd)*chlhdrr(jr,lh,js)
                ENDDO
                
                DO k = 1,nsp
                   chdrt(k,js)=chdrt(k,js)+ylht(k,lh,nd)*chlhdr(jr,lh,js)
                   chdrf(k,js)=chdrf(k,js)+ylhf(k,lh,nd)*chlhdr(jr,lh,js)
                   chdt(k,js) =chdt(k,js) +ylht(k,lh,nd)*chlh(jr,lh,js)
                   chdf(k,js) =chdf(k,js) +ylhf(k,lh,nd)*chlh(jr,lh,js)
                   chdtt(k,js)=chdtt(k,js)+ylhtt(k,lh,nd)*chlh(jr,lh,js)
                   chdff(k,js)=chdff(k,js)+ylhff(k,lh,nd)*chlh(jr,lh,js)
                   chdtf(k,js)=chdtf(k,js)+ylhtf(k,lh,nd)*chlh(jr,lh,js)
                ENDDO
             ENDDO ! lh
          ENDDO   ! js
          
       
          CALL mkgylm(jspins,atoms%rmsh(jr,n),thet,nsp,&
141
               ch_tmp,chdr,chdt,chdf,chdrr,chdtt,chdff,chdtf,chdrt,chdrf,grad,kt)
142
       ENDIF
143 144
       !Set charge to minimum value
       IF (PRESENT(ch)) ch(kt+1:kt+nsp,:)=MAX(ch_tmp(:nsp,:),d_15)
145 146 147 148 149 150 151 152 153 154 155 156
       kt=kt+nsp
    END DO
    
  END SUBROUTINE mt_to_grid


  SUBROUTINE mt_from_grid(atoms,sphhar,nsp,n,jspins,v_in,vr)
    IMPLICIT NONE
    TYPE(t_atoms),INTENT(IN) :: atoms
    TYPE(t_sphhar),INTENT(IN):: sphhar
    INTEGER,INTENT(IN)       :: nsp,jspins,n
    REAL,INTENT(IN)          :: v_in(:,:)
157
    REAL,INTENT(INOUT)       :: vr(:,0:,:)
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
    
    REAL    :: vpot(nsp),vlh
    INTEGER :: js,kt,lh,jr,nd
    nd=atoms%ntypsy(SUM(atoms%neq(:n-1))+1)
    DO js = 1,jspins
       !
       kt=0
       DO jr=1,atoms%jri(n)
          vpot=v_in(kt+1:kt+nsp,js)*wt(:)!  multiplicate v_in with the weights of the k-points
          
          DO lh = 0,sphhar%nlh(nd)
             !
             ! --->        determine the corresponding potential number
             !c            through gauss integration
             !
             vlh=dot_PRODUCT(vpot(:),ylh(:nsp,lh,nd))
174
             vr(jr,lh,js) = vr(jr,lh,js) + vlh
175
          ENDDO ! lh
176
          kt=kt+nsp
177 178 179 180 181 182
       ENDDO   ! jr
    ENDDO

  END SUBROUTINE mt_from_grid
    
  SUBROUTINE finish_mt_grid()
Daniel Wortmann's avatar
Daniel Wortmann committed
183 184
    DEALLOCATE(ylh,wt,rx,thet)
    IF (ALLOCATED(ylht)) DEALLOCATE(ylht,ylhtt,ylhf,ylhff,ylhtf)
185 186 187
  END SUBROUTINE finish_mt_grid

END MODULE m_mt_tofrom_grid