rhosphnlo.f90 5.97 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! 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.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15
MODULE m_rhosphnlo
  !***********************************************************************
  ! Add the local orbital contributions to the charge density. The 
  ! corresponding summation of the pure apw contribuions is done in
  ! cdnval.
  ! Philipp Kurz 99/04
  !***********************************************************************
CONTAINS
  SUBROUTINE rhosphnlo(itype,atoms,sphhar, uloulopn,dulon,uulon,&
16
       ello,vr, aclo,bclo,cclo,acnmt,bcnmt,ccnmt,f,g, rho,qmtllo,rhoLRes)
17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

    USE m_constants, ONLY : c_light,sfp_const
    USE m_radsra
    USE m_radsrdn
    USE m_types
    IMPLICIT NONE
    TYPE(t_sphhar),INTENT(IN)   :: sphhar
    TYPE(t_atoms),INTENT(IN)    :: atoms
    !     ..
    !     .. Scalar Arguments ..
    INTEGER,    INTENT (IN) :: itype 
    !     ..
    !     .. Array Arguments ..
    REAL,    INTENT (IN) :: aclo(atoms%nlod),bclo(atoms%nlod),cclo(atoms%nlod,atoms%nlod)
    REAL,    INTENT (IN) :: acnmt(0:atoms%lmaxd,atoms%nlod,sphhar%nlhd)
    REAL,    INTENT (IN) :: bcnmt(0:atoms%lmaxd,atoms%nlod,sphhar%nlhd)
    REAL,    INTENT (IN) :: ccnmt(atoms%nlod,atoms%nlod,sphhar%nlhd)
    REAL,    INTENT (IN) :: dulon(atoms%nlod),uulon(atoms%nlod),vr(atoms%jmtd)
    REAL,    INTENT (IN) :: uloulopn(atoms%nlod,atoms%nlod),ello(atoms%nlod)
    REAL,    INTENT (IN) :: f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
    REAL,    INTENT (INOUT) :: qmtllo(0:atoms%lmaxd)
    REAL,    INTENT (INOUT) :: rho(atoms%jmtd,0:sphhar%nlhd)
39
    REAL,    INTENT (INOUT) :: rhoLRes(atoms%jmtd,0:sphhar%nlhd,0:(atoms%lmaxd*(atoms%lmaxd+1))/2+atoms%lmaxd)
40 41 42
    !     ..
    !     .. Local Scalars ..
    REAL dsdum,usdum ,c_1,c_2
43 44
    INTEGER j,l,lh,lo,lop,lp,nodedum,llp
    REAL dus,ddn,c,temp
45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93
    !     ..
    !     .. Local Arrays ..
    REAL,    ALLOCATABLE :: flo(:,:,:),glo(:,:)
    REAL filo(atoms%jmtd,2)
    !     ..
    c = c_light(1.0)
    c_1 = 1.0 / atoms%neq(itype)
    c_2 = 1.0 /(atoms%neq(itype)*sfp_const)
    !
    DO lo = 1,atoms%nlo(itype)
       l = atoms%llo(lo,itype)
       qmtllo(l) = qmtllo(l) + (aclo(lo)*uulon(lo) +bclo(lo)*dulon(lo)) * c_1
       DO lop = 1,atoms%nlo(itype)
          IF (atoms%llo(lop,itype).EQ.l) THEN
             qmtllo(l) = qmtllo(l) + (cclo(lop,lo) *uloulopn(lop,lo)) * c_1
          END IF
       END DO
    END DO
    ALLOCATE ( flo(atoms%jmtd,2,atoms%nlod),glo(atoms%jmtd,2) )

    !---> calculate the local ortital radial functions

    DO lo = 1,atoms%nlo(itype)
       l = atoms%llo(lo,itype)
       CALL radsra(ello(lo),l,vr,atoms%rmsh(1,itype),atoms%dx(itype),atoms%jri(itype),c,&
            usdum,dus,nodedum,flo(:,1,lo),flo(:,2,lo))
       !+apw+lo
       IF (atoms%l_dulo(lo,itype).or.atoms%ulo_der(lo,itype).ge.1) THEN
          !--->    calculate orthogonal energy derivative at e
          j = atoms%ulo_der(lo,itype)
          IF(atoms%l_dulo(lo,itype)) j = 1
          CALL radsrdn(ello(lo),l,vr,atoms%rmsh(1,itype),atoms%dx(itype),atoms%jri(itype),c,&
               usdum,dsdum,ddn,nodedum,glo,filo,flo(:,:,lo),dus,j) ! filo is a dummy array&
          DO j=1,atoms%jri(itype)
             flo(j,1,lo) = glo(j,1)
             flo(j,2,lo) = glo(j,2)
          ENDDO
          ddn = sqrt(ddn)
          IF(atoms%l_dulo(lo,itype)) ddn=1.0
          flo(:,:,lo) = flo(:,:,lo)/ddn ! Normalize ulo (flo) if APW+lo is not used
       ENDIF
       !-apw+lo
    END DO

    !---> add the contribution of the local orbitals and flapw - lo cross-
    !---> terms to the spherical chargedensity inside the muffin tins.

    DO lo = 1,atoms%nlo(itype)
       l = atoms%llo(lo,itype)
94
       llp = (l* (l+1))/2 + l
95
       DO j = 1,atoms%jri(itype)
96 97 98 99 100
          temp = c_2 *&
                 (aclo(lo) * ( f(j,1,l)*flo(j,1,lo) +f(j,2,l)*flo(j,2,lo) ) +&
                 bclo(lo) * ( g(j,1,l)*flo(j,1,lo) +g(j,2,l)*flo(j,2,lo) ) )
          rho(j,0) = rho(j,0) + temp
          rhoLRes(j,0,llp) = rhoLRes(j,0,llp) + temp
101 102 103 104
       END DO
       DO lop = 1,atoms%nlo(itype)
          IF (atoms%llo(lop,itype).EQ.l) THEN
             DO j = 1,atoms%jri(itype)
105
                temp = c_2 * cclo(lop,lo) *&
106
                     ( flo(j,1,lop)*flo(j,1,lo) +flo(j,2,lop)*flo(j,2,lo) )
107 108
                rho(j,0) = rho(j,0) + temp
                rhoLRes(j,0,llp) = rhoLRes(j,0,llp) + temp
109 110 111 112 113 114 115 116 117 118 119
             END DO
          END IF
       END DO
    END DO

    !---> add the contribution of the local orbitals and flapw - lo cross-
    !---> terms to the non-spherical chargedensity inside the muffin tins.

    DO lh = 1,sphhar%nlh(atoms%ntypsy(atoms%nat))
       DO lp = 0,atoms%lmax(itype)
          DO lo = 1,atoms%nlo(itype)
120 121
             l = atoms%llo(lo,itype)
             llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
122
             DO j = 1,atoms%jri(itype)
123
                temp = c_1 * (&
124 125
                     acnmt(lp,lo,lh) * (f(j,1,lp)*flo(j,1,lo) +f(j,2,lp)*flo(j,2,lo) ) +&
                     bcnmt(lp,lo,lh) * (g(j,1,lp)*flo(j,1,lo) +g(j,2,lp)*flo(j,2,lo) ) )
126 127
                rho(j,lh) = rho(j,lh) + temp
                rhoLRes(j,lh,llp) = rhoLRes(j,lh,llp) + temp
128 129 130 131
             END DO
          END DO
       END DO
       DO lo = 1,atoms%nlo(itype)
132
          l = atoms%llo(lo,itype)
133
          DO lop = 1,atoms%nlo(itype)
134 135
             lp = atoms%llo(lop,itype)
             llp = (MAX(l,lp)* (MAX(l,lp)+1))/2 + MIN(l,lp)
136
             DO j = 1,atoms%jri(itype)
137
                temp = c_1 * ccnmt(lop,lo,lh) *&
138
                     ( flo(j,1,lop)*flo(j,1,lo) +flo(j,2,lop)*flo(j,2,lo) )
139 140
                rho(j,lh) = rho(j,lh) + temp
                rhoLRes(j,lh,llp) = rhoLRes(j,lh,llp) + temp
141 142 143 144 145 146 147 148
             END DO
          END DO
       END DO
    END DO
    DEALLOCATE (flo,glo)

  END SUBROUTINE rhosphnlo
END MODULE m_rhosphnlo