exchange_core.F90 28.9 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
!     Calculates the HF exchange term
!
!                                          s          s*          s            s*
!                                       phi    (r) phi     (r) phi     (r') phi    (r')
!                         occ.             n_1k       n'k+q       n'k+q        n_2k
!     exchange(n,q)  =  - SUM  INT INT  ------------------------------------------- dr dr'
!                         k,n'                           | r - r' |
!
!                         occ                  s          s    ~        ~       s         s
!                    =  - SUM  SUM  v     < phi      | phi     M    > < M    phi     | phi      >
!                         k,n' I,J   k,IJ      n'k+q      n_1k  q,I      q,J    n_2k      n'k+q
!
!     for the different combinations of n_1 and n_2 and  n' runs over the core states.
!     ( n_1,n_2:  valence-valence, core-core,core-valence )
!
!     It is done directly without employing the mixed basis set.


Daniel Wortmann's avatar
Daniel Wortmann committed
19 20 21 22 23
MODULE m_exchange_core

CONTAINS

  SUBROUTINE exchange_vccv(nk,atoms, hybrid,hybdat, DIMENSION,jsp,lapw,&
24
       maxbands,mnobd,mpi,degenerat,symequivalent,results,&
Daniel Wortmann's avatar
Daniel Wortmann committed
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
       ex_vv_r,ex_vv_c,l_real)


    USE m_constants   
    USE m_util
    USE m_wrapper
    USE m_types
    USE m_io_hybrid
    IMPLICIT NONE

    TYPE(t_hybdat),INTENT(IN)   :: hybdat
    TYPE(t_results),INTENT(INOUT)   :: results
    TYPE(t_mpi),INTENT(IN)   :: mpi
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_hybrid),INTENT(IN)   :: hybrid
    TYPE(t_atoms),INTENT(IN)   :: atoms
    TYPE(t_lapw),INTENT(IN)   :: lapw

    !     -scalars -
    INTEGER,INTENT(IN)      :: jsp 
    INTEGER,INTENT(IN)      ::nk  ,maxbands, mnobd
    !     - arays -
Daniel Wortmann's avatar
Daniel Wortmann committed
47
    INTEGER,INTENT(IN)      ::  degenerat(hybrid%ne_eig(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
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
    LOGICAL,INTENT(IN)      :: l_real
    REAL    ,INTENT(INOUT)  ::  ex_vv_r(:,:,:)!(maxbands,mnobd,nkpti)
    COMPLEX ,INTENT(INOUT)  ::  ex_vv_c(:,:,:)!(maxbands,mnobd,nkpti)
    LOGICAL                 ::  symequivalent(COUNT(degenerat .GE. 1), COUNT(degenerat .GE. 1))

    !     - local scalars -
    INTEGER                 ::  iatom,ieq,itype,ic,l,l1,l2, ll,lm ,m1,m2,p1,p2,n,n1,n2,i,j
    INTEGER                 ::  iband1,iband2,ndb1,ndb2,ic1,ic2
    INTEGER                 ::  m

    REAL                    ::  time1,time2
    REAL                    ::  rdum
    REAL                    ::  sum_offdia

    COMPLEX                 ::  cdum

    !     - local arrays -
    INTEGER,ALLOCATABLE     ::  larr(:),larr2(:)
    INTEGER,ALLOCATABLE     ::  parr(:),parr2(:)

    REAL                    ::  integrand(atoms%jmtd)
    REAL                    ::  primf1(atoms%jmtd),primf2(atoms%jmtd)
    REAL,ALLOCATABLE        ::  fprod(:,:),fprod2(:,:)
    REAL,ALLOCATABLE        ::  integral(:,:)

    COMPLEX                 ::  cmt(DIMENSION%neigd,hybrid%maxlmindx,atoms%nat)
Daniel Wortmann's avatar
Daniel Wortmann committed
74
    COMPLEX                 ::  exchange(hybrid%nbands(nk),hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
75 76
    COMPLEX,ALLOCATABLE     ::  carr(:,:),carr2(:,:),carr3(:,:)

Daniel Wortmann's avatar
Daniel Wortmann committed
77
    LOGICAL                 ::  ldum(hybrid%nbands(nk),hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
78

79 80
    WRITE(6,'(A)') new_LINE('n') // new_LINE('n') // '### valence-core-core-valence exchange ###'
    WRITE(6,'(A)') new_LINE('n') // '        k-point       band    exchange (core contribution)'
Daniel Wortmann's avatar
Daniel Wortmann committed
81 82 83 84 85 86 87 88

    ! read in mt wavefunction coefficients from file cmt
    CALL read_cmt(cmt,nk)
    ALLOCATE ( fprod(atoms%jmtd,5),larr(5),parr(5) )

    ! generate ldum(nbands(nk),nbands(nk)), which is true if the corresponding matrix entry is non-zero
    ic1  = 0
    ldum = .FALSE.
Daniel Wortmann's avatar
Daniel Wortmann committed
89
    DO iband1 = 1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
90 91
       ndb1 = degenerat(iband1)
       IF( ndb1 .GE. 1 ) THEN
92 93
          ic1 = ic1 + 1
          ic2 = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
94
          DO iband2 = 1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
95 96 97 98 99 100 101 102 103 104
             ndb2 = degenerat(iband2)
             IF( ndb2 .GE. 1 ) THEN
                ic2 = ic2 + 1
                IF( symequivalent(ic2,ic1) ) THEN
                   IF( ndb1 .NE. ndb2 ) STOP 'exchange: failure symequivalent'
                   DO i = 0,ndb1-1
                      DO j = 0,ndb2 - 1
                         ldum(iband1+i,iband2+j) = .TRUE.
                      END DO
                   END DO
105

Daniel Wortmann's avatar
Daniel Wortmann committed
106 107
                END IF
             END IF
108
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
109 110 111 112 113 114 115 116
       END IF
    END DO

    exchange = 0
    iatom = 0
    rdum  = 0
    DO itype = 1,atoms%ntype
       DO ieq = 1,atoms%neq(itype)
117
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
118
          DO l1 = 0,hybdat%lmaxc(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
             DO p1 = 1,hybdat%nindxc(l1,itype)

                DO l = 0,hybrid%lcutm1(itype)

                   ! Define core-valence product functions

                   n = 0
                   DO l2 = 0,atoms%lmax(itype)
                      IF(l.LT.ABS(l1-l2).OR.l.GT.l1+l2) CYCLE

                      DO p2 = 1,hybrid%nindx(l2,itype)
                         n = n + 1
                         M = SIZE(fprod,2)
                         IF(n.GT.M) THEN
                            ALLOCATE ( fprod2(atoms%jmtd,M),larr2(M),parr2(M) )
                            fprod2 = fprod ; larr2 = larr ; parr2 = parr
                            DEALLOCATE ( fprod,larr,parr )
                            ALLOCATE ( fprod(atoms%jmtd,M+5),larr(M+5),parr(M+5) )
                            fprod(:,:M) = fprod2
                            larr(:M)    = larr2
                            parr(:M)    = parr2
                            DEALLOCATE ( fprod2,larr2,parr2 )
                         END IF
142 143
                         fprod(:atoms%jri(itype),n) = ( hybdat%core1(:atoms%jri(itype),p1,l1,itype) *hybdat%bas1 (:atoms%jri(itype),p2,l2,itype)&
                              +hybdat%core2(:atoms%jri(itype),p1,l1,itype) *hybdat%bas2 (:atoms%jri(itype),p2,l2,itype) )/ atoms%rmsh(:atoms%jri(itype),itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
144 145 146 147
                         larr(n)    = l2
                         parr(n)    = p2
                      END DO
                   END DO
148

Daniel Wortmann's avatar
Daniel Wortmann committed
149
                   ! Evaluate radial integrals (special part of Coulomb matrix : contribution from single MT)
150

Daniel Wortmann's avatar
Daniel Wortmann committed
151
                   ALLOCATE ( integral(n,n),carr(n,hybrid%nbands(nk)), carr2(n,lapw%nv(jsp)),carr3(n,lapw%nv(jsp)) )
152

Daniel Wortmann's avatar
Daniel Wortmann committed
153 154
                   DO i = 1,n
                      CALL primitivef(primf1,fprod(:,i)*atoms%rmsh(:,itype)**(l+1) ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd, itype,atoms%ntype)
155
                      CALL primitivef(primf2,fprod(:atoms%jri(itype),i)/atoms%rmsh(:atoms%jri(itype),itype)**l ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd,-itype,atoms%ntype)  ! -itype is to enforce inward integration
156

157 158
                      primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype)) / atoms%rmsh(:atoms%jri(itype),itype)**l
                      primf2(:atoms%jri(itype)) = primf2(:atoms%jri(itype)) * atoms%rmsh(:atoms%jri(itype),itype)**(l+1)
Daniel Wortmann's avatar
Daniel Wortmann committed
159 160 161 162 163 164
                      DO j = 1,n
                         integrand     = fprod(:,j) * (primf1 + primf2)
                         integral(i,j) = fpi_const/(2*l+1) * intgrf(integrand,atoms%jri,atoms%jmtd,atoms%rmsh,&
                              atoms%dx,atoms%ntype,itype,hybdat%gridf)
                      END DO
                   END DO
165

Daniel Wortmann's avatar
Daniel Wortmann committed
166
                   ! Add everything up
167

Daniel Wortmann's avatar
Daniel Wortmann committed
168 169 170
                   DO m1 = -l1,l1
                      DO M = -l,l
                         m2 = m1 + M
171

Daniel Wortmann's avatar
Daniel Wortmann committed
172 173
                         carr = 0
                         ic   = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
174
                         DO n1=1,hybrid%nbands(nk)
175

Daniel Wortmann's avatar
Daniel Wortmann committed
176 177 178
                            DO i = 1,n
                               ll      = larr(i)
                               IF(ABS(m2).GT.ll) CYCLE
179

Daniel Wortmann's avatar
Daniel Wortmann committed
180
                               lm = SUM((/ ((2*l2+1)*hybrid%nindx(l2,itype),l2=0,ll-1) /)) + (m2+ll)*hybrid%nindx(ll,itype) + parr(i)
181

Daniel Wortmann's avatar
Daniel Wortmann committed
182
                               carr(i,n1) = cmt(n1,lm,iatom) * gaunt(l1,ll,l,m1,m2,M,hybdat%maxfac,hybdat%fac,hybdat%sfac)
183

Daniel Wortmann's avatar
Daniel Wortmann committed
184 185 186 187 188 189 190 191
                            END DO
                            DO n2=1,n1
                               IF( ldum(n2,n1) ) THEN
                                  ic =ic + 1
                                  exchange(n2,n1) = exchange(n2,n1) + dot_PRODUCT( carr(:,n1), MATMUL(integral,carr(:,n2)) )
                               END IF
                            END DO
                         END DO
192
                      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
193 194 195 196
                   END DO

                   DEALLOCATE ( integral,carr,carr2,carr3 )

197
                END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
198 199 200 201
             END DO
          END DO
       END DO
    END DO
202 203


Daniel Wortmann's avatar
Daniel Wortmann committed
204 205 206 207 208 209 210 211 212 213 214 215 216 217

    IF (l_real) THEN
       IF( ANY(ABS(AIMAG(exchange)).GT.1d-10) ) THEN
          IF ( mpi%irank == 0 ) WRITE(6,'(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
          WRITE(*,*) MAXVAL(ABS(AIMAG(exchange)))
          STOP 'exchangeCore: Unusually large imaginary component.'
       END IF
    ENDIF

    ! add the core-valence contribution to the exchange matrix ex_vv

    !      ic         = 0
    sum_offdia = 0
    IF (l_real) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
218 219
       DO n1=1,hybrid%nobd(nk)
          DO n2=1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
220 221 222 223 224
             ex_vv_r(n2,n1,nk) = ex_vv_r(n2,n1,nk) - exchange(n1,n2)
             IF( n1 /= n2) sum_offdia = sum_offdia + 2*ABS(exchange(n1,n2))
          END DO
       END DO
    ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
225 226
       DO n1=1,hybrid%nobd(nk)
          DO n2=1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
227 228
             ex_vv_c(n2,n1,nk) = ex_vv_c(n2,n1,nk) - exchange(n1,n2)
             IF( n1 /= n2) sum_offdia = sum_offdia + 2*ABS(exchange(n1,n2))
229
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
230 231 232
       END DO
    END IF

Daniel Wortmann's avatar
Daniel Wortmann committed
233
    DO n1=1,hybrid%nobd(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
234 235 236
       results%te_hfex%core = results%te_hfex%core - results%w_iks(n1,nk,jsp)*exchange(n1,n1)
    END DO

237
    WRITE(6,'(A,F20.15)') 'sum of the absolut real part of the non diagonal elements',sum_offdia
Daniel Wortmann's avatar
Daniel Wortmann committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265


  END SUBROUTINE exchange_vccv

  SUBROUTINE exchange_vccv1(nk,atoms, hybrid,hybdat, DIMENSION,jsp, lapw,&
       nsymop,nsest,indx_sest,mpi, a_ex,results, mat_ex)

    USE m_constants  
    USE m_util
    USE m_wrapper
    USE m_types
    USE m_io_hybrid
    IMPLICIT NONE

    TYPE(t_hybdat),INTENT(IN)   :: hybdat
    TYPE(t_results),INTENT(INOUT)   :: results
    TYPE(t_mpi),INTENT(IN)   :: mpi
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_hybrid),INTENT(IN)   :: hybrid
    TYPE(t_atoms),INTENT(IN)   :: atoms
    TYPE(t_lapw),INTENT(IN)   :: lapw

    !     -scalars -
    INTEGER,INTENT(IN)      :: jsp 
    INTEGER,INTENT(IN)      :: nk  
    INTEGER,INTENT(IN)      ::  nsymop
    REAL,INTENT(IN)         ::  a_ex
    !     - arays -
Daniel Wortmann's avatar
Daniel Wortmann committed
266
    INTEGER,INTENT(IN)      ::  nsest(hybrid%nbands(nk)),indx_sest(hybrid%nbands(nk),hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289

    TYPE(t_mat),INTENT(INOUT):: mat_ex
    !     - local scalars - 
    INTEGER                 ::  iatom,ieq,itype,ic,l,l1,l2
    INTEGER                 ::  ll,lm ,m1,m2,p1,p2,n,n1,n2,nn2,i,j
    INTEGER                 ::  iband1,iband2,ndb1,ndb2,ic1,ic2
    INTEGER                 ::  m

    REAL                    ::  time1,time2
    REAL                    ::  rdum
    REAL                    ::  sum_offdia

    COMPLEX                 ::  cdum
    !     - local arrays -
    INTEGER,ALLOCATABLE     ::  larr(:),larr2(:)
    INTEGER,ALLOCATABLE     ::  parr(:),parr2(:)

    REAL                    ::  integrand(atoms%jmtd)
    REAL                    ::  primf1(atoms%jmtd),primf2(atoms%jmtd)
    REAL,ALLOCATABLE        ::  fprod(:,:),fprod2(:,:)
    REAL,ALLOCATABLE        ::  integral(:,:)

    COMPLEX                 ::  cmt(DIMENSION%neigd,hybrid%maxlmindx,atoms%nat)
Daniel Wortmann's avatar
Daniel Wortmann committed
290
    COMPLEX                 ::  exchange(hybrid%nbands(nk),hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
291 292
    COMPLEX,ALLOCATABLE     ::  carr(:,:),carr2(:,:),carr3(:,:)

Daniel Wortmann's avatar
Daniel Wortmann committed
293
    LOGICAL                 ::  ldum(hybrid%nbands(nk),hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
294 295 296 297 298 299 300 301 302 303 304 305 306 307



    ! read in mt wavefunction coefficients from file cmt

    CALL read_cmt(cmt,nk)

    ALLOCATE ( fprod(atoms%jmtd,5),larr(5),parr(5) )

    exchange = 0
    iatom    = 0
    rdum     = 0
    DO itype = 1,atoms%ntype
       DO ieq = 1,atoms%neq(itype)
308
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
309
          DO l1 = 0,hybdat%lmaxc(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
             DO p1 = 1,hybdat%nindxc(l1,itype)

                DO l = 0,hybrid%lcutm1(itype) 

                   ! Define core-valence product functions

                   n = 0
                   DO l2 = 0,atoms%lmax(itype)
                      IF(l.LT.ABS(l1-l2).OR.l.GT.l1+l2) CYCLE

                      DO p2 = 1,hybrid%nindx(l2,itype)
                         n = n + 1
                         M = SIZE(fprod,2)
                         IF(n.GT.M) THEN
                            ALLOCATE ( fprod2(atoms%jmtd,M),larr2(M),parr2(M) )
                            fprod2 = fprod ; larr2 = larr ; parr2 = parr
                            DEALLOCATE ( fprod,larr,parr )
                            ALLOCATE ( fprod(atoms%jmtd,M+5),larr(M+5),parr(M+5) )
                            fprod(:,:M) = fprod2
                            larr(:M)    = larr2
                            parr(:M)    = parr2
                            DEALLOCATE ( fprod2,larr2,parr2 )
                         END IF
333 334
                         fprod(:atoms%jri(itype),n) = ( hybdat%core1(:atoms%jri(itype),p1,l1,itype) *hybdat%bas1 (:atoms%jri(itype),p2,l2,itype) &
                              +hybdat%core2(:atoms%jri(itype),p1,l1,itype) *hybdat%bas2 (:atoms%jri(itype),p2,l2,itype) )/ atoms%rmsh(:atoms%jri(itype),itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
335 336 337 338
                         larr(n)    = l2
                         parr(n)    = p2
                      END DO
                   END DO
339

Daniel Wortmann's avatar
Daniel Wortmann committed
340
                   ! Evaluate radial integrals (special part of Coulomb matrix : contribution from single MT)
341

Daniel Wortmann's avatar
Daniel Wortmann committed
342
                   ALLOCATE ( integral(n,n),carr(n,hybrid%nbands(nk)), carr2(n,lapw%nv(jsp)),carr3(n,lapw%nv(jsp)) )
343

Daniel Wortmann's avatar
Daniel Wortmann committed
344
                   DO i = 1,n
345 346
                      CALL primitivef(primf1,fprod(:atoms%jri(itype),i)*atoms%rmsh(:atoms%jri(itype),itype)**(l+1) ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd, itype,atoms%ntype)
                      CALL primitivef(primf2,fprod(:atoms%jri(itype),i)/atoms%rmsh(:atoms%jri(itype),itype)**l ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd,-itype,atoms%ntype)  ! -itype is to enforce inward integration
347

348 349
                      primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype)) / atoms%rmsh(:atoms%jri(itype),itype)**l
                      primf2(:atoms%jri(itype)) = primf2(:atoms%jri(itype)) * atoms%rmsh(:atoms%jri(itype),itype)**(l+1)
Daniel Wortmann's avatar
Daniel Wortmann committed
350 351 352 353
                      DO j = 1,n
                         integrand     = fprod(:,j) * (primf1 + primf2)
                         integral(i,j) = fpi_const/(2*l+1) * intgrf(integrand,atoms%jri,atoms%jmtd,atoms%rmsh,&
                              atoms%dx,atoms%ntype,itype,hybdat%gridf)
354
                      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
355 356 357 358 359 360 361 362 363
                   END DO

                   ! Add everything up

                   DO m1 = -l1,l1
                      DO M = -l,l
                         m2 = m1 + M

                         carr = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
364
                         DO n1=1,hybrid%nbands(nk)
365

Daniel Wortmann's avatar
Daniel Wortmann committed
366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381
                            DO i = 1,n
                               ll      = larr(i)
                               IF(ABS(m2).GT.ll) CYCLE

                               lm = SUM((/ ((2*l2+1)*hybrid%nindx(l2,itype),l2=0,ll-1) /))&
                                    + (m2+ll)*hybrid%nindx(ll,itype) + parr(i)

                               carr(i,n1) = cmt(n1,lm,iatom) * gaunt(l1,ll,l,m1,m2,M,hybdat%maxfac,hybdat%fac,hybdat%sfac)

                            END DO
                            DO n2=1,nsest(n1)!n1
                               nn2 = indx_sest(n2,n1)
                               exchange(nn2,n1) = exchange(nn2,n1) + dot_PRODUCT( carr(:,n1), MATMUL(integral,carr(:,nn2)) )

                            END DO
                         END DO
382
                      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
383
                   END DO
384

Daniel Wortmann's avatar
Daniel Wortmann committed
385
                   DEALLOCATE ( integral,carr,carr2,carr3 )
386

Daniel Wortmann's avatar
Daniel Wortmann committed
387 388
                END DO
             END DO
389
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
390 391 392 393 394 395 396 397 398 399 400 401 402 403
       END DO
    END DO




    IF (mat_ex%l_real) THEN
       IF( ANY(ABS(AIMAG(exchange)).GT.1d-10) ) THEN
          IF (mpi%irank == 0) WRITE(6,'(A)') 'exchangeCore: Warning! Unusually large imaginary component.'
          WRITE(*,*) MAXVAL(ABS(AIMAG(exchange)))
          STOP 'exchangeCore: Unusually large imaginary component.'
       END IF
    ENDIF

Daniel Wortmann's avatar
Daniel Wortmann committed
404
    DO n1=1,hybrid%nobd(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469
       results%te_hfex%core = results%te_hfex%Core - a_ex * results%w_iks(n1,nk,jsp)*exchange(n1,n1)
    END DO

    ! add the core-valence contribution to the exchange matrix mat_ex
    ! factor 1/nsymop is needed due to the symmetrization in symmetrizeh 

    ic         = 0
    sum_offdia = 0
    IF (mat_ex%l_real) THEN
       mat_ex%data_r=mat_ex%data_r+exchange/nsymop
    ELSE
       mat_ex%data_c=mat_ex%data_c+CONJG(exchange)/nsymop
    END IF

  END SUBROUTINE exchange_vccv1


  SUBROUTINE exchange_cccc(nk,atoms,hybdat, ncstd, sym,kpts,a_ex,mpi, results)


    USE m_constants   
    USE m_util
    USE m_wrapper
    USE m_gaunt
    USE m_trafo
    USE m_types
    USE m_io_hybrid
    IMPLICIT NONE
    TYPE(t_hybdat),INTENT(IN)   :: hybdat
    TYPE(t_results),INTENT(INOUT)   :: results
    TYPE(t_mpi),INTENT(IN)   :: mpi
    TYPE(t_sym),INTENT(IN)   :: sym
    TYPE(t_kpts),INTENT(IN)   :: kpts
    TYPE(t_atoms),INTENT(IN)   :: atoms

    ! - scalars - 
    INTEGER,INTENT(IN)    ::  nk,  ncstd

    REAL   ,INTENT(IN)    ::  a_ex

    ! - arays -

    ! - local scalars - 
    INTEGER               ::  itype,ieq,icst,icst1,icst2,iatom,iatom0
    INTEGER               ::  l1,l2,l,ll,llmax
    INTEGER               ::  m1,m2 ,mm,m
    INTEGER               ::  n1,n2,n

    REAL                  ::  rdum,rdum1
    ! - local arrays -
    INTEGER               ::  point(hybdat%maxindxc,-hybdat%lmaxcd:hybdat%lmaxcd, 0:hybdat%lmaxcd,atoms%nat)
    REAL                  ::  rprod(atoms%jmtd),primf1(atoms%jmtd),primf2(atoms%jmtd), integrand(atoms%jmtd)
    COMPLEX               ::  exch(ncstd,ncstd)

    !       IF ( irank == 0 ) THEN
    !         WRITE(6,'(//A)') '### core-core-core-core exchange ###'
    !         WRITE(6,'(/A)') '        k-point       band    exchange'
    !       END IF


    ! set up point
    icst = 0
    iatom= 0
    DO itype = 1,atoms%ntype
       DO ieq = 1,atoms%neq(itype)
470
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
471
          DO l = 0,hybdat%lmaxc(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
472 473 474 475 476 477
             DO M = -l,l
                DO n = 1,hybdat%nindxc(l,itype)
                   icst = icst + 1
                   point(n,M,l,iatom) = icst
                END DO
             END DO
478
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
479 480 481 482 483 484 485 486 487
       END DO
    END DO

    llmax  = 2*hybdat%lmaxcd 
    exch   = 0
    iatom0 = 0
    DO itype = 1,atoms%ntype

       DO l1 = 0,hybdat%lmaxc(itype)  ! left core state
Daniel Wortmann's avatar
Daniel Wortmann committed
488
          DO l2 = 0,hybdat%lmaxc(itype)  ! right core state
Daniel Wortmann's avatar
Daniel Wortmann committed
489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505
             DO l = 0,hybdat%lmaxc(itype)   ! occupied core state

                DO ll = ABS(l1-l),l1+l
                   IF( ll .LT. ABS(l-l2) .OR. ll .GT. l+l2 ) CYCLE
                   IF( MOD(l+l1+ll,2) .NE. 0 ) CYCLE
                   IF( MOD(l+l2+ll,2) .NE. 0 ) CYCLE

                   DO m1 = -l1,l1
                      m2 = m1
                      IF( ABS(m2) .GT. l2 ) CYCLE
                      DO M = -l,l
                         mm = M - m1
                         IF( ABS(mm) .GT. ll ) CYCLE
                         rdum = fpi_const/(2*ll+1)*gaunt1(l,ll,l1,M,mm,m1,llmax) *gaunt1(l,ll,l2,M,mm,m2,llmax)

                         DO n = 1,hybdat%nindxc(l,itype)
                            DO n2 = 1,hybdat%nindxc(l2,itype)
506 507
                               rprod(:atoms%jri(itype)) = ( hybdat%core1(:atoms%jri(itype),n,l,itype)*hybdat%core1(:atoms%jri(itype),n2,l2,itype)&
                                    +hybdat%core2(:atoms%jri(itype),n,l,itype)*hybdat%core2(:atoms%jri(itype),n2,l2,itype) ) / atoms%rmsh(:atoms%jri(itype),itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
508 509

                               CALL primitivef(primf1,rprod(:)*atoms%rmsh(:,itype)**(ll+1) ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd, itype,atoms%ntype)
510
                               CALL primitivef(primf2,rprod(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype),itype)**ll ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd,-itype,atoms%ntype)  ! -itype is to enforce inward integration
Daniel Wortmann's avatar
Daniel Wortmann committed
511 512 513 514 515 516 517


                               primf1 = primf1 / atoms%rmsh(:,itype)**ll
                               primf2 = primf2 * atoms%rmsh(:,itype)**(ll+1)

                               DO n1 = 1,hybdat%nindxc(l1,itype)

518 519
                                  rprod(:atoms%jri(itype)) = ( hybdat%core1(:atoms%jri(itype),n,l,itype)*hybdat%core1(:atoms%jri(itype),n1,l1,itype)&
                                       +hybdat%core2(:atoms%jri(itype),n,l,itype)*hybdat%core2(:atoms%jri(itype),n1,l1,itype) ) / atoms%rmsh(:atoms%jri(itype),itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543

                                  integrand     = rprod * (primf1 + primf2)

                                  rdum1 = rdum*intgrf(integrand,atoms%jri,atoms%jmtd,&
                                       atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf)

                                  iatom = iatom0
                                  DO ieq = 1,atoms%neq(itype)
                                     iatom = iatom + 1
                                     icst1    = point(n1,m1,l1,iatom)
                                     icst2    = point(n2,m2,l2,iatom)
                                     exch(icst1,icst2) = exch(icst1,icst2)+rdum1
                                  END DO
                               END DO  !n1

                            END DO  !n2
                         END DO  !n

                      END DO  !M
                   END DO  !m1

                END DO  !ll

             END DO  !l
544
          END DO  !l2
Daniel Wortmann's avatar
Daniel Wortmann committed
545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
       END DO  !l1
       iatom0 = iatom0 + atoms%neq(itype)
    END DO  !itype      

    IF (sym%invs) THEN
       CALL symmetrize(exch,ncstd,ncstd,3,.FALSE., atoms,hybdat%lmaxc,hybdat%lmaxcd,hybdat%nindxc, sym)
       IF( ANY( ABS(AIMAG(exch)) .GT. 1E-6 ) ) STOP 'exchange_cccc: exch possesses significant imaginary part'
    ENDIF
    !       DO icst = 1,ncstd
    !         IF ( irank == 0 )
    !      &    WRITE(6,'(    ''  ('',F5.3,'','',F5.3,'','',F5.3,'')'',I4,1X,F12.5)')bkpt,icst,REAL(exch(icst,icst))*(-27.211608)
    !       END DO


    ! add core exchange contributions to the te_hfex

    DO icst1 = 1,ncstd
       results%te_hfex%core = results%te_hfex%core - a_ex*kpts%wtkpt(nk)*exch(icst1,icst1)
    END DO

  END SUBROUTINE exchange_cccc

  SUBROUTINE exchange_cccv(nk,atoms,hybdat, hybrid,DIMENSION,maxbands,ncstd,&
       bkpt,sym,mpi, exch_cv_r, exch_cv_c,l_real )

    USE m_constants   
    USE m_util
    USE m_wrapper
    USE m_gaunt
    USE m_trafo
    USE m_io_hybrid
    USE m_types
    IMPLICIT NONE
    TYPE(t_hybdat),INTENT(IN)   :: hybdat
    TYPE(t_mpi),INTENT(IN)   :: mpi
    TYPE(t_dimension),INTENT(IN)   :: DIMENSION
    TYPE(t_hybrid),INTENT(IN)   :: hybrid
    TYPE(t_sym),INTENT(IN)   :: sym
    TYPE(t_atoms),INTENT(IN)   :: atoms
    ! - scalars - 
    INTEGER,INTENT(IN)    ::  nk  ,ncstd
    INTEGER,INTENT(IN)    :: maxbands

    ! - arays -
    REAL   ,INTENT(IN)    ::  bkpt(3)
    LOGICAL,INTENT(IN)    :: l_real
    REAL   ,INTENT(INOUT) ::  exch_cv_r(:,:,:)!(maxbands,ncstd,nkpti)
    COMPLEX,INTENT(INOUT) ::  exch_cv_c(:,:,:) !(maxbands,ncstd,nkpti)
    ! - local scalars - 
    INTEGER               ::  itype,ieq,icst,icst1,icst2,iatom,iatom0
    INTEGER               ::    iatom1,iband
    INTEGER               ::  l1,l2,l,ll,llmax
    INTEGER               ::  lm2,lmp2
    INTEGER               ::  m1,m2 ,mm,m
    INTEGER               ::  n1,n2,n,nn

    REAL                  ::  rdum0,rdum1,rdum2,rdum3,rdum4
    COMPLEX               ::  cdum
    COMPLEX,PARAMETER     ::  img=(0d0,1d0)
    ! - local arrays -
    INTEGER               ::  point(hybdat%maxindxc,-hybdat%lmaxcd:hybdat%lmaxcd, 0:hybdat%lmaxcd,atoms%nat)
    INTEGER               ::  lmstart(0:atoms%lmaxd,atoms%ntype)
    REAL                  ::  rprod(atoms%jmtd),primf1(atoms%jmtd),primf2(atoms%jmtd)
    REAL                  ::  integrand(atoms%jmtd)
    COMPLEX               ::  cexp(atoms%nat)
Daniel Wortmann's avatar
Daniel Wortmann committed
610 611
    COMPLEX               ::  exch(hybrid%nbands(nk),ncstd)
    COMPLEX               ::  cmt(DIMENSION%neigd,hybrid%maxlmindx,atoms%nat),carr(hybrid%nbands(nk))
Daniel Wortmann's avatar
Daniel Wortmann committed
612 613 614 615 616 617 618 619 620 621 622

    IF ( mpi%irank == 0 ) THEN
       WRITE(6,'(//A)') '### core-core-core-valence exchange  ###'
       WRITE(6,'(/A)') '        k-point       band    exchange'
    END IF

    ! set up point
    icst = 0
    iatom= 0
    DO itype = 1,atoms%ntype
       DO ieq = 1,atoms%neq(itype)
623
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
624
          DO l = 0,hybdat%lmaxc(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
625 626 627 628 629 630
             DO M = -l,l
                DO n = 1,hybdat%nindxc(l,itype)
                   icst = icst + 1
                   point(n,M,l,iatom) = icst
                END DO
             END DO
631
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647
       END DO
    END DO

    ! lmstart = lm start index for each l-quantum number and atom type (for cmt-coefficients)
    DO itype = 1,atoms%ntype
       DO l = 0,atoms%lmax(itype)
          lmstart(l,itype) = SUM( (/ (hybrid%nindx(ll,itype)*(2*ll+1),ll=0,l-1) /) )
       END DO
    END DO

    ! read in cmt coefficient at k-point nk

    CALL read_cmt(cmt,nk)
    iatom = 0
    DO itype = 1,atoms%ntype
       DO ieq = 1,atoms%neq(itype)
648
          iatom       = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
649 650 651
          cexp(iatom) = EXP ( img*tpi_const* dot_PRODUCT(bkpt(:),atoms%taual(:,iatom)) )
       END DO
    END DO
652

Daniel Wortmann's avatar
Daniel Wortmann committed
653
    cmt = CONJG(cmt)
654

Daniel Wortmann's avatar
Daniel Wortmann committed
655
    llmax = MAX(2*hybdat%lmaxcd,atoms%lmaxd)
656

Daniel Wortmann's avatar
Daniel Wortmann committed
657 658 659 660 661
    exch   = 0
    iatom0 = 0
    DO itype = 1,atoms%ntype

       DO l1 = 0,hybdat%lmaxc(itype)  ! left core state
662
          DO l2 = 0,atoms%lmax(itype)  ! right valence state
Daniel Wortmann's avatar
Daniel Wortmann committed
663
             DO l = 0,hybdat%lmaxc(itype)   ! occupied core state
664

Daniel Wortmann's avatar
Daniel Wortmann committed
665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688
                DO ll = ABS(l1-l),l1+l
                   IF( ll .LT. ABS(l-l2) .OR. ll .GT. l+l2 ) CYCLE
                   IF( MOD(l+l1+ll,2) .NE. 0 ) CYCLE
                   IF( MOD(l+l2+ll,2) .NE. 0 ) CYCLE

                   !                 WRITE(*,*) 'l1,l2,l,ll',l1,l2,l,ll
                   rdum0 = fpi_const/(2*ll+1)

                   DO m1 = -l1,l1
                      m2 = m1
                      IF( ABS(m2) .GT. l2 ) CYCLE
                      lm2 = lmstart(l2,itype) +(m2+l2)*hybrid%nindx(l2,itype)

                      DO M = -l,l
                         mm = M - m1
                         IF( ABS(M-m1) .GT. ll ) CYCLE

                         rdum1 = gaunt1(l,ll,l1,M,mm,m1,llmax)&
                              *gaunt1(l,ll,l2,M,mm,m1,llmax)*rdum0

                         DO n = 1,hybdat%nindxc(l,itype)
                            DO n2 = 1,hybrid%nindx(l2,itype)
                               lmp2 = lm2 + n2

689 690
                               rprod(:atoms%jri(itype)) = ( hybdat%core1(:atoms%jri(itype),n,l,itype)*hybdat%bas1(:atoms%jri(itype),n2,l2,itype)&
                                    +hybdat%core2(:atoms%jri(itype),n,l,itype)*hybdat%bas2(:atoms%jri(itype),n2,l2,itype) ) / atoms%rmsh(:atoms%jri(itype),itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
691

692 693
                               CALL primitivef(primf1,rprod(:atoms%jri(itype))*atoms%rmsh(:atoms%jri(itype),itype)**(ll+1) ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd, itype,atoms%ntype)
                               CALL primitivef(primf2,rprod(:atoms%jri(itype))/atoms%rmsh(:atoms%jri(itype),itype)**ll ,atoms%rmsh,atoms%dx,atoms%jri,atoms%jmtd,-itype,atoms%ntype)  ! -itype is to enforce inward integration
Daniel Wortmann's avatar
Daniel Wortmann committed
694 695 696 697




698
                               primf1(:atoms%jri(itype)) = primf1(:atoms%jri(itype)) / atoms%rmsh(:atoms%jri(itype),itype)**ll
Daniel Wortmann's avatar
Daniel Wortmann committed
699 700 701 702 703
                               primf2 = primf2 * atoms%rmsh(:,itype)**(ll+1)

                               DO n1 = 1,hybdat%nindxc(l1,itype)

                                  rprod(:) = ( hybdat%core1(:,n,l,itype)*hybdat%core1(:,n1,l1,itype)&
704
                                       +hybdat%core2(:,n,l,itype)*hybdat%core2(:,n1,l1,itype) ) / atoms%rmsh(:atoms%jri(itype),itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
705 706 707 708 709 710 711 712 713 714

                                  integrand     = rprod * (primf1 + primf2)

                                  rdum2 = rdum1*intgrf(integrand,atoms%jri,atoms%jmtd, atoms%rmsh,atoms%dx,atoms%ntype,itype,hybdat%gridf)

                                  iatom = iatom0
                                  DO ieq = 1,atoms%neq(itype)
                                     iatom = iatom + 1
                                     icst1 = point(n1,m1,l1,iatom)
                                     cdum  = rdum2*cexp(iatom)
Daniel Wortmann's avatar
Daniel Wortmann committed
715
                                     DO iband = 1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732

                                        exch(iband,icst1) = exch(iband,icst1) + cdum*cmt(iband,lmp2,iatom)

                                     END DO
                                  END DO

                               END DO  !n1

                            END DO  !n2
                         END DO  !n

                      END DO  !M
                   END DO  !m1

                END DO  !ll

             END DO  !l
733
          END DO  !l2
Daniel Wortmann's avatar
Daniel Wortmann committed
734 735 736 737 738 739 740
       END DO  !l1
       iatom0 = iatom0 + atoms%neq(itype)
    END DO  !itype      


    IF (l_real) THEN
       !symmetrize core-wavefunctions such that phi(-r) = phi(r)*
Daniel Wortmann's avatar
Daniel Wortmann committed
741
       CALL symmetrize( exch,hybrid%nbands(nk),ncstd,2,.FALSE., atoms,hybdat%lmaxc,hybdat%lmaxcd,hybdat%nindxc, sym)
Daniel Wortmann's avatar
Daniel Wortmann committed
742 743 744 745 746 747

       IF( ANY( ABS(AIMAG(exch)) .GT. 1E-6 ) ) STOP 'exchange_cccv: exch possesses significant imaginary part'
    ENDIF

    IF (l_real) THEN
       DO icst = 1,ncstd
Daniel Wortmann's avatar
Daniel Wortmann committed
748
          DO iband = 1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
749 750 751 752 753
             exch_cv_r(iband,icst,nk) = exch_cv_r(iband,icst,nk) - exch(iband,icst)
          END DO
       END DO
    ELSE
       DO icst = 1,ncstd
Daniel Wortmann's avatar
Daniel Wortmann committed
754
          DO iband = 1,hybrid%nbands(nk)
Daniel Wortmann's avatar
Daniel Wortmann committed
755 756 757 758 759 760 761 762 763
             exch_cv_c(iband,icst,nk) = exch_cv_c(iband,icst,nk) - exch(iband,icst)
          END DO
       END DO
    END IF

  END SUBROUTINE exchange_cccv


END MODULE m_exchange_core