abclocdn_pulay.F90 9.76 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
MODULE m_abclocdnpulay
  USE m_juDFT
CONTAINS
  SUBROUTINE abclocdn_pulay(&
       &                          atoms,sym,&
       &                          noco,ccchi,kspin,iintsp,&
       &                          con1,phase,ylm,ntyp,na,k,fgp,&
14
       &                          s,nv,ne,nbasf0,alo1,blo1,clo1,&
15
       &                          kvec,nkvec,enough,acof,bcof,ccof,&
16
       &                          acoflo,bcoflo,aveccof,bveccof,cveccof,zMat)
17 18 19 20 21 22 23 24 25 26
    !
    !*********************************************************************
    ! for details see abclocdn; calles by to_pulay
    !*********************************************************************
    !
    USE m_types
    IMPLICIT NONE
    TYPE(t_noco),INTENT(IN)   :: noco
    TYPE(t_sym),INTENT(IN)    :: sym
    TYPE(t_atoms),INTENT(IN)  :: atoms
27
    TYPE(t_zMat),INTENT(IN)   :: zMat
28 29 30 31 32 33 34 35
    !     ..
    !     .. Scalar Arguments ..
    INTEGER, INTENT (IN) :: iintsp
    INTEGER, INTENT (IN) :: k,na,ne,ntyp,nv,kspin
    REAL,    INTENT (IN) :: con1 ,s
    COMPLEX, INTENT (IN) :: phase
    !     ..
    !     .. Array Arguments ..
Daniel Wortmann's avatar
Daniel Wortmann committed
36
    INTEGER, INTENT (IN) :: nbasf0(atoms%nlod,atoms%nat) 
37
    INTEGER, INTENT (IN) :: kvec(2*(2*atoms%llod+1),atoms%nlod)
Daniel Wortmann's avatar
Daniel Wortmann committed
38 39
    REAL,    INTENT (IN) :: alo1(atoms%nlod,atoms%ntype),blo1(atoms%nlod,atoms%ntype)
    REAL,    INTENT (IN) :: clo1(atoms%nlod,atoms%ntype)
40 41
    REAL,    INTENT (IN) :: fgp(3)
    COMPLEX, INTENT (IN) :: ylm( (atoms%lmaxd+1)**2 ),ccchi(2)
Daniel Wortmann's avatar
Daniel Wortmann committed
42 43 44
    COMPLEX, INTENT (INOUT) :: acof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, INTENT (INOUT) :: bcof(:,0:,:)!(nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, INTENT (INOUT) :: ccof(-atoms%llod:,:,:,:)!(-llod:llod,nobd,atoms%nlod,atoms%nat)
45 46
    COMPLEX, INTENT (INOUT) :: acoflo(-atoms%llod:,:,:,:)
    COMPLEX, INTENT (INOUT) :: bcoflo(-atoms%llod:,:,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
47 48 49 50 51
    COMPLEX, INTENT (INOUT) :: aveccof(:,:,0:,:)!(3,nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, INTENT (INOUT) :: bveccof(:,:,0:,:)!(3,nobd,0:dimension%lmd,atoms%nat)
    COMPLEX, INTENT (INOUT) :: cveccof(:,-atoms%llod:,:,:,:)!(3,-atoms%llod:llod,nobd,atoms%nlod,atoms%nat)
    LOGICAL, INTENT (OUT) :: enough(atoms%nat)
    INTEGER, INTENT (INOUT) :: nkvec(atoms%nlod,atoms%nat)
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
    !     ..
    !     .. Local Scalars ..
    COMPLEX ctmp,term1
    REAL,PARAMETER:: linindq=1.0e-4,eps=1.e-30
    INTEGER i,ie,l,ll1,lm,lo ,mind,nbasf,na2,lmp,m
    LOGICAL linind
    !     ..
    !     .. Local Arrays ..
    COMPLEX clotmp(-atoms%llod:atoms%llod)
    !     ..
    enough(na) = .TRUE.
    term1 = con1* ((atoms%rmt(ntyp)**2)/2)*phase
    !
    !---> 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.
    !
    DO lo = 1,atoms%nlo(ntyp)
       l = atoms%llo(lo,ntyp)
       IF (.NOT.((s.LE.eps).AND.(l.GE.1))) THEN
          IF (atoms%invsat(na).EQ.0) THEN
             IF ((nkvec(lo,na)).LT. (2*atoms%llo(lo,ntyp)+1)) THEN
                enough(na) = .FALSE.
                nkvec(lo,na) = nkvec(lo,na) + 1
                nbasf = nbasf0(lo,na) + nkvec(lo,na)
                l = atoms%llo(lo,ntyp)
                ll1 = l* (l+1)
                DO m = -l,l
                   clotmp(m) = term1*CONJG(ylm(ll1+m+1))
                END DO
                !
                IF ( kvec(nkvec(lo,na),lo) == k ) THEN
                   DO ie = 1,ne
                      DO m = -l,l
                         lm = ll1 + m
                         IF (noco%l_noco) THEN
                            IF (noco%l_ss) THEN
90
                               ctmp = clotmp(m)* ccchi(iintsp)*zMat%z_c(kspin+nbasf,ie)
91
                            ELSE
92
                               ctmp = clotmp(m)*( ccchi(1)*zMat%z_c(nbasf,ie)+ccchi(2)*zMat%z_c(kspin+nbasf,ie) )
93 94
                            ENDIF
                         ELSE
95
                            IF (zmat%l_real) THEN
96
                               ctmp = zMat%z_r(nbasf,ie)*clotmp(m)
97
                            ELSE
98
                               ctmp = zMat%z_c(nbasf,ie)*clotmp(m)
99
                            ENDIF
100 101 102 103 104 105 106 107 108 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
                         ENDIF
                         acof(ie,lm,na)     = acof(ie,lm,na) +ctmp*alo1(lo,ntyp)
                         bcof(ie,lm,na)     = bcof(ie,lm,na) +ctmp*blo1(lo,ntyp)
                         ccof(m,ie,lo,na)   = ccof(m,ie,lo,na) +ctmp*clo1(lo,ntyp)
                         acoflo(m,ie,lo,na) = acoflo(m,ie,lo,na) +ctmp*alo1(lo,ntyp)
                         bcoflo(m,ie,lo,na) = bcoflo(m,ie,lo,na) +ctmp*blo1(lo,ntyp)
                         DO i = 1,3
                            aveccof(i,ie,lm,na)=aveccof(i,ie,lm,na) +fgp(i)*ctmp*alo1(lo,ntyp)
                            bveccof(i,ie,lm,na)=bveccof(i,ie,lm,na) +fgp(i)*ctmp*blo1(lo,ntyp)
                            cveccof(i,m,ie,lo,na) =cveccof(i,m,ie,lo,na) +fgp(i)*ctmp*clo1(lo,ntyp)
                         ENDDO
                      END DO
                   END DO
                   !                    write(6,9000) nbasf,k,lo,na,
                   !     +                          (clo1(lo,ntyp)*clotmp(m),m=-l,l)
                   ! 9000               format(2i4,2i2,7(' (',e9.3,',',e9.3,')'))
                ELSE
                   nkvec(lo,na) = nkvec(lo,na) - 1
                END IF
             END IF
          ELSEIF (atoms%invsat(na).EQ.1) THEN
             IF ((nkvec(lo,na)).LT. (2* (2*atoms%llo(lo,ntyp)+1))) THEN
                enough(na) = .FALSE.
                nkvec(lo,na) = nkvec(lo,na) + 1
                nbasf = nbasf0(lo,na) + nkvec(lo,na)
                l = atoms%llo(lo,ntyp)
                ll1 = l* (l+1)
                DO m = -l,l
                   clotmp(m) = term1*CONJG(ylm(ll1+m+1))
                END DO
                !
                IF ( kvec(nkvec(lo,na),lo) == k ) THEN
                   !                     write(*,*)'k vector nr ',k,' has been accepted'
                   DO ie = 1,ne
                      DO m = -l,l
                         lm = ll1 + m
                         IF (noco%l_noco) THEN
                            IF (noco%l_ss) THEN
138
                               ctmp = clotmp(m)*ccchi(iintsp)*zMat%z_c(kspin+nbasf,ie)
139
                            ELSE
140
                               ctmp = clotmp(m)*( ccchi(1)*zMat%z_c(nbasf,ie)+ccchi(2)*zMat%z_c(kspin+nbasf,ie) )
141 142
                            ENDIF
                         ELSE
143
                            IF (zmat%l_real) THEN
144
                               ctmp = zMat%z_r(nbasf,ie)*clotmp(m)
145
                            ELSE
146
                               ctmp = zMat%z_c(nbasf,ie)*clotmp(m)
147
                            END IF
148 149 150 151 152 153 154 155 156 157 158
                         ENDIF
                         acof(ie,lm,na) = acof(ie,lm,na) +ctmp*alo1(lo,ntyp)
                         bcof(ie,lm,na) = bcof(ie,lm,na) +ctmp*blo1(lo,ntyp)
                         ccof(m,ie,lo,na) = ccof(m,ie,lo,na) +ctmp*clo1(lo,ntyp)
                         acoflo(m,ie,lo,na) = acoflo(m,ie,lo,na) +ctmp*alo1(lo,ntyp)
                         bcoflo(m,ie,lo,na) = bcoflo(m,ie,lo,na) +ctmp*blo1(lo,ntyp)
                         DO i = 1,3
                            aveccof(i,ie,lm,na)=aveccof(i,ie,lm,na) +fgp(i)*ctmp*alo1(lo,ntyp)
                            bveccof(i,ie,lm,na)=bveccof(i,ie,lm,na) +fgp(i)*ctmp*blo1(lo,ntyp)
                            cveccof(i,m,ie,lo,na)=cveccof(i,m,ie,lo,na)+fgp(i)*ctmp*clo1(lo,ntyp)
                         ENDDO
159
                         IF (noco%l_soc.AND.sym%invs) THEN
160
                            ctmp = zMat%z_c(nbasf,ie) * CONJG(clotmp(m))*(-1)**(l-m)
161 162 163 164 165 166 167 168 169 170 171 172 173
                            na2 = sym%invsatnr(na)
                            lmp = ll1 - m
                            acof(ie,lmp,na2) = acof(ie,lmp,na2) +ctmp*alo1(lo,ntyp)
                            bcof(ie,lmp,na2) = bcof(ie,lmp,na2) +ctmp*blo1(lo,ntyp)
                            ccof(-m,ie,lo,na2) = ccof(-m,ie,lo,na2) + ctmp*clo1(lo,ntyp)
                            acoflo(-m,ie,lo,na2) = acoflo(-m,ie,lo,na2) +ctmp*alo1(lo,ntyp)
                            bcoflo(-m,ie,lo,na2) = bcoflo(-m,ie,lo,na2) +ctmp*blo1(lo,ntyp)
                            DO i = 1,3
                               aveccof(i,ie,lmp,na2)=aveccof(i,ie,lmp,na2)-fgp(i)*ctmp*alo1(lo,ntyp)
                               bveccof(i,ie,lmp,na2)=bveccof(i,ie,lmp,na2)-fgp(i)*ctmp*blo1(lo,ntyp)
                               cveccof(i,-m,ie,lo,na2) =cveccof(i,-m,ie,lo,na2) -fgp(i)*ctmp*clo1(lo,ntyp)
                            ENDDO
                         ENDIF
174 175 176 177 178 179 180
                      ENDDO ! loop over m
                   ENDDO    ! loop over eigenstates (ie)
                ELSE
                   nkvec(lo,na) = nkvec(lo,na) - 1
                END IF   ! linind
             END IF      ! nkvec(lo,na) < 2*(2*atoms%llo + 1)
          ELSE
181
             CALL juDFT_error("invsat =/= 0 or 1",calledby ="abclocdn_pulay")
182 183 184 185 186 187
          ENDIF
       ELSE
          enough(na) = .FALSE.
       ENDIF ! s > eps  & l >= 1  
    END DO
    IF ((k.EQ.nv) .AND. (.NOT.enough(na))) THEN
188 189
       WRITE (6,FMT=*) 'abclocdn did not find enough linearly independent'
       WRITE (6,FMT=*) 'ccof coefficient-vectors. the linear independence'
190 191
       WRITE (6,FMT=*) 'quality, linindq, is set to: ',linindq,'.'
       WRITE (6,FMT=*) 'this value might be to large.'
192
       CALL juDFT_error("did not find enough lin. ind. ccof-vectors" ,calledby ="abclocdn_pulay")
193 194 195 196
    END IF

  END SUBROUTINE abclocdn_pulay
END MODULE m_abclocdnpulay