abccoflo.f90 6.63 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
MODULE m_abccoflo
  USE m_juDFT
  !*********************************************************************
  ! Calculates the (upper case) A, B and C coefficients for the local
  ! orbitals.
  ! Philipp Kurz 99/04
  !*********************************************************************
CONTAINS
15
  SUBROUTINE abccoflo(atoms,sym, con1,rph,cph,ylm,ntyp,na,k,nv, l_lo1,alo1,blo1,&
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
       clo1, nkvec, enough,alo,blo,clo,kvec)
    !
    !*************** ABBREVIATIONS ***************************************
    ! kvec    : stores the number of the G-vectors, that have been used to
    !           construct the local orbitals
    ! nkvec   : stores the number of G-vectors that have been found and
    !           accepted during the construction of the local orbitals.
    ! enough  : enough is set to .true. when enough G-vectors have been
    !           accepted.
    ! linindq : if the norm of that part of a local orbital (contructed 
    !           with a trial G-vector) that is orthogonal to the previous
    !           ones is larger than linindq, then this G-vector is 
    !           accepted.
    !*********************************************************************
    !
    USE m_constants
    USE m_types
    IMPLICIT NONE

    TYPE(t_atoms),INTENT(IN)   :: atoms
36
    TYPE(t_sym),INTENT(IN)     :: sym
37 38 39 40 41 42 43 44 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
    !     .. 
    !     .. Scalar Arguments ..
    REAL,    INTENT (IN) :: con1,cph ,rph
    INTEGER, INTENT (IN) :: k,na,ntyp,nv
    LOGICAL, INTENT (IN) :: l_lo1
    LOGICAL, INTENT (OUT):: enough
    !     ..
    !     .. Array Arguments ..
    INTEGER, INTENT (IN)::  kvec(2* (2*atoms%llod+1),atoms%nlod )
    REAL,    INTENT (IN) :: alo1(atoms%nlod),blo1(atoms%nlod),clo1(atoms%nlod)
    COMPLEX, INTENT (IN) :: ylm( (atoms%lmaxd+1)**2 )
    COMPLEX, INTENT (OUT):: alo(-atoms%llod:atoms%llod,2* (2*atoms%llod+1),atoms%nlod)  
    COMPLEX, INTENT (OUT):: blo(-atoms%llod:atoms%llod,2* (2*atoms%llod+1),atoms%nlod)  
    COMPLEX, INTENT (OUT):: clo(-atoms%llod:atoms%llod,2* (2*atoms%llod+1),atoms%nlod)  
    INTEGER,INTENT (INOUT):: nkvec(atoms%nlod)
    !     ..
    !     .. Local Scalars ..
    COMPLEX term1
    REAL,PARAMETER:: linindq=1.e-4
    INTEGER l,lo ,mind,ll1,lm,m
    LOGICAL linind
    !     ..
    !
    !---> the whole program is in hartree units, therefore 1/wronskian is
    !---> (rmt**2)/2. the factor i**l, which usually appears in the a, b
    !---> and c coefficients, is included in the t-matrices. thus, it does
    !---> not show up in the formula above.
    !
    !-abccoflo1
    IF ( l_lo1) THEN
       DO lo = 1,atoms%nlo(ntyp)
          IF ( (nkvec(lo).EQ.0).AND.(atoms%llo(lo,ntyp).EQ.0) ) THEN
             enough = .FALSE.
             nkvec(lo) = 1
             m = 0
             clo(m,nkvec(lo),lo) = con1* ((atoms%rmt(ntyp)**2)/2) / SQRT(fpi_const)
             alo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*alo1(lo)
             blo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*blo1(lo)
             clo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*clo1(lo)
             IF (kvec(nkvec(lo),lo)/=k)  CALL juDFT_error("abccoflo:1"&
                  &           ,calledby ="abccoflo")

          ENDIF
       ENDDO
    ELSE
       enough = .TRUE.
       term1 = con1* ((atoms%rmt(ntyp)**2)/2)*CMPLX(rph,cph)
       DO lo = 1,atoms%nlo(ntyp)
85
          IF (sym%invsat(na).EQ.0) THEN
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
             IF ((nkvec(lo)).LT. (2*atoms%llo(lo,ntyp)+1)) THEN
                enough = .FALSE.
                nkvec(lo) = nkvec(lo) + 1
                l = atoms%llo(lo,ntyp)
                ll1 = l*(l+1) + 1
                DO m = -l,l
                   lm = ll1 + m
                   clo(m,nkvec(lo),lo) = term1*ylm(lm)
                END DO
                IF ( kvec(nkvec(lo),lo) == k ) THEN
                   DO m = -l,l
                      alo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*alo1(lo)
                      blo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*blo1(lo)
                      clo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*clo1(lo)
                   END DO
                   !                  WRITE(6,9000) nkvec(lo),k,lo,na,
                   !     +                          (clo(m,nkvec(lo),lo),m=-l,l)
                   ! 9000             format(2i4,2i2,7(' (',e9.3,',',e9.3,')'))
                ELSE
                   nkvec(lo) = nkvec(lo) - 1
                ENDIF
             ENDIF
          ELSE
109
             IF ((sym%invsat(na).EQ.1) .OR. (sym%invsat(na).EQ.2)) THEN
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 141 142 143 144 145 146 147 148 149 150 151 152
                !           only invsat=1 is needed invsat=2 for testing
                IF ((nkvec(lo)).LT. (2* (2*atoms%llo(lo,ntyp)+1))) THEN
                   enough = .FALSE.
                   nkvec(lo) = nkvec(lo) + 1
                   l = atoms%llo(lo,ntyp)
                   ll1 = l*(l+1) + 1
                   DO m = -l,l
                      lm = ll1 + m
                      clo(m,nkvec(lo),lo) = term1*ylm(lm)
                   END DO
                   IF ( kvec(nkvec(lo),lo) == k ) THEN
                      DO m = -l,l
                         !                            if(l.eq.1) then
                         !               WRITE(*,*)'k=',k,' clotmp=',clo(m,nkvec(lo),lo)
                         !               WRITE(*,*)'clo1=',clo1(lo),' term1=',term1
                         !                            endif
                         alo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*alo1(lo)
                         blo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*blo1(lo)
                         clo(m,nkvec(lo),lo) = clo(m,nkvec(lo),lo)*clo1(lo)
                         !                        kvec(nkvec(lo),lo) = k
                      END DO
                   ELSE
                      nkvec(lo) = nkvec(lo) - 1
                   END IF
                END IF
             END IF
          END IF
       END DO
       IF ((k.EQ.nv) .AND. (.NOT.enough)) THEN
          WRITE (6,FMT=*)&
               &     'abccoflo did not find enough linearly independent'
          WRITE (6,FMT=*)&
               &     'clo coefficient-vectors. the linear independence'
          WRITE (6,FMT=*) 'quality, linindq, is set to: ',linindq,'.'
          WRITE (6,FMT=*) 'this value might be to large.'
          CALL juDFT_error&
               &        ("abccoflo: did not find enough lin. ind. clo-vectors"&
               &        ,calledby ="abccoflo")
       END IF
    ENDIF  ! abccoflo1

  END SUBROUTINE abccoflo
END MODULE m_abccoflo