vmtxc.f90 4.92 KB
Newer Older
1 2 3 4 5
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------
6 7
MODULE m_vmtxc
CONTAINS
8
  SUBROUTINE vmtxc( DIMENSION,sphhar,atoms,den,xcpot,input,sym, vxc, exc,vx)
9
    !     ********************************************************************
10
    !     Calculate the LDA xc-potential in the MT-spheres
11 12 13 14 15 16
    !     *********************************************************

    USE m_lhglpts
    USE m_gaussp
    USE m_types
    IMPLICIT NONE
17
    CLASS(t_xcpot),INTENT(IN)      :: xcpot
18 19 20 21 22
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_input),INTENT(IN)       :: input
    TYPE(t_sym),INTENT(IN)         :: sym
    TYPE(t_sphhar),INTENT(IN)      :: sphhar
    TYPE(t_atoms),INTENT(IN)       :: atoms
23
    TYPE(t_potden),INTENT(IN)      :: den
24
    TYPE(t_potden),INTENT(INOUT)   :: vxc,exc
25 26
    TYPE(t_potden),INTENT(INOUT),optional     :: vx
    
27 28 29 30 31 32
    !     ..
    !     .. Local Scalars ..
    REAL elh,rr2,r2rho,vlh
    INTEGER jr,js,k,lh,n,nd,i,nsp,nat
    !     ..
    !     .. Local Arrays ..
33 34
    REAL v_x(DIMENSION%nspd,DIMENSION%jspd),v_xc(DIMENSION%nspd,DIMENSION%jspd),pos0(3)
    REAL rhoxc(DIMENSION%nspd,DIMENSION%jspd),e_ex(DIMENSION%nspd),wt(DIMENSION%nspd),rx(3,DIMENSION%nspd)
35 36 37 38 39 40 41 42 43 44 45 46
    REAL ylh(DIMENSION%nspd,0:sphhar%nlhd,sphhar%ntypsd)
    !     ..
    !     generates nspd points on a sherical shell with radius 1.0
    !     angular mesh equidistant in phi, 
    !     theta are zeros of the legendre polynomials
    !     
    CALL gaussp(atoms%lmaxd,rx,wt)

    nsp = DIMENSION%nspd
    !
    !     generates the lattice harmonics on the angular mesh 
    !     
Daniel Wortmann's avatar
Daniel Wortmann committed
47
    CALL lhglpts( sphhar,atoms, rx,nsp, sym, ylh)
48 49 50 51
    !
    !
    !     loop over topologically non-equivalent atoms
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
52
    !$OMP PARALLEL DO DEFAULT(NONE)&
53
    !$OMP SHARED(atoms,input,sphhar,xcpot,den,ylh,nsp,wt,vxc,vx,exc)&
Daniel Wortmann's avatar
Daniel Wortmann committed
54
    !$OMP PRIVATE(elh,rr2,r2rho,vlh,jr,js,k,lh,n,nd,i,nat)&
55
    !$OMP PRIVATE(v_x,v_xc,pos0,rhoxc,e_ex)
56 57 58 59 60 61 62 63 64 65 66
    DO n = 1,atoms%ntype
       nat=SUM(atoms%neq(:n-1))+1
       nd = atoms%ntypsy(nat)
       !
       !     loop over radial mesh
       !
       DO  jr = 1,atoms%jri(n)
          rr2 = 1.e0/ (atoms%rmsh(jr,n)*atoms%rmsh(jr,n))
          DO  js = 1,input%jspins
             rhoxc(:,js) = 0.
             DO  lh = 0,sphhar%nlh(nd)
67
                r2rho = rr2*den%mt(jr,lh,n,js)
68 69 70 71
                !
                !     generate the densities on an angular mesh
                !
                DO  k = 1,nsp
Daniel Wortmann's avatar
Daniel Wortmann committed
72
                   rhoxc(k,js) = rhoxc(k,js) + ylh(k,lh,nd)*r2rho
73 74 75 76 77
                ENDDO
             ENDDO
          ENDDO
          !     calculate the ex.-cor. potential
          !
78
          CALL xcpot%get_vxc(input%jspins,rhoxc(:nsp,:), v_xc,v_x) 
79 80 81 82 83 84
          !     ----> now determine the corresponding potential number 
          DO js = 1,input%jspins
             !
             ! ---> multiplikate vxc with the weights of the k-points    
             !
             DO  k = 1,nsp
85
                v_xc(k,js) = v_xc(k,js)*wt(k)
86

87
                v_x (k,js) = v_x (k,js)*wt(k)
88 89 90 91 92 93
             ENDDO
             DO  lh = 0,sphhar%nlh(nd)

                !
                ! ---> determine the corresponding potential number through gauss integration
                !
94
                vlh=DOT_PRODUCT( v_xc(:nsp,js),ylh(:nsp,lh,nd))
95 96 97 98

                !
                ! ---> add to the given potential

99
                vxc%mt(jr,lh,n,js) = vxc%mt(jr,lh,n,js) + vlh
100

101
                IF (PRESENT(vx) ) THEN
102 103
                   vlh = 0
                   DO  k = 1,nsp
104
                      vlh = vlh + v_x(k,js)*ylh(k,lh,nd)
105 106
                   END DO

107
                   vx%mt(jr,lh,n,js) = vx%mt(jr,lh,n,js) + vlh
108 109 110
                ENDIF
             ENDDO
          ENDDO
111
          IF (ALLOCATED(exc%mt)) THEN 
112 113 114
             !
             !     calculate the ex.-cor energy density
             !
115 116
             CALL xcpot%get_exc(input%jspins,rhoxc(:nsp,:), e_ex) 
             !     ----> now determine the corresponding energy density number 
117
             !
118
             ! ---> multiplikate exc with the weights of the k-points    
119
             !
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136
             DO  k = 1,nsp
                e_ex(k) = e_ex(k)*wt(k)
             ENDDO
             DO  lh = 0,sphhar%nlh(nd)
                
                !
                ! ---> determine the corresponding potential number through gauss integration
                !
                elh=DOT_PRODUCT(e_ex(:nsp),ylh(:nsp,lh,nd))
                
                !
                ! ---> add to the given potential
                
                exc%mt(jr,lh,n,1) = elh
             ENDDO
          END IF
          
137 138 139 140 141 142 143 144 145
       ENDDO

    ENDDO
    !$OMP END PARALLEL DO
    !

    RETURN
  END SUBROUTINE vmtxc
END MODULE m_vmtxc