subvxc.F90 21.2 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 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 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
MODULE m_subvxc
CONTAINS
  SUBROUTINE subvxc(lapw,bk, DIMENSION,input,jsp,vr0, atoms,usdus, hybrid, el,ello,sym,&
       nlot_d,kveclo, cell, sphhar, stars,xcpot,mpi,oneD,hamovlp,vx)


    USE m_intgr,     ONLY : intgr3
    USE m_constants
    USE m_gaunt,     ONLY : gaunt1
    USE m_wrapper
    USE m_loddop
    USE m_radflo
    USE m_radfun
    USE m_abcof3
    USE m_types
    IMPLICIT NONE
    TYPE(t_xcpot),INTENT(IN)     :: xcpot
    TYPE(t_mpi),INTENT(IN)       :: mpi
    TYPE(t_dimension),INTENT(IN) :: DIMENSION
    TYPE(t_oneD),INTENT(IN)      :: oneD
    TYPE(t_hybrid),INTENT(INOUT) :: hybrid
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_sym),INTENT(IN)       :: sym
    TYPE(t_stars),INTENT(IN)     :: stars
    TYPE(t_cell),INTENT(IN)      :: cell
    TYPE(t_sphhar),INTENT(IN)    :: sphhar
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_lapw),INTENT(IN)      :: lapw
    TYPE(t_usdus),INTENT(INOUT)  :: usdus
    TYPE(t_potden),INTENT(IN)    :: vx
    !     .. Scalar Arguments ..

    INTEGER, INTENT (IN) :: jsp 
    INTEGER, INTENT (IN) :: nlot_d




    !     .. Array Arguments ..
    INTEGER, INTENT (IN) :: kveclo(nlot_d)

    REAL,    INTENT (IN) :: vr0(atoms%jmtd,atoms%ntype,DIMENSION%jspd)               ! just for radial functions
    REAL,    INTENT (IN) :: el(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT (IN) :: ello(atoms%nlod,atoms%ntype,DIMENSION%jspd)
    REAL,    INTENT (IN) :: bk(3)
    TYPE(t_hamovlp),INTENT(INOUT)::hamovlp

    !     .. Local Scalars ..
    INTEGER               ::  ic,indx,m,ig1,ig2
    INTEGER               ::  nlharm,nnbas,typsym,lm
    INTEGER               ::  noded,nodeu
    INTEGER               ::  nbasf0
    INTEGER               ::  i,j,l,ll,l1,l2 ,m1,m2  ,j1,j2
    INTEGER               ::  ok,p1,p2,lh,mh,pp1,pp2
    INTEGER               ::  igrid,itype,ilharm,istar
    INTEGER               ::  ineq,iatom,ilo,ilop,ieq,icentry
    INTEGER               ::  ikvecat,ikvecprevat,invsfct,ikvec,ikvecp
    INTEGER               ::  lp,mp,pp
    REAL                  ::  a_ex
    REAL                  ::  wronk
    COMPLEX               ::  rc,rr

    !     .. Local Arrays ..
    INTEGER               ::  gg(3)
    INTEGER               ::  pointer_lo(atoms%nlod,atoms%ntype)

    REAL                  ::  integ(0:sphhar%nlhd,hybrid%maxindx,0:atoms%lmaxd,hybrid%maxindx,0:atoms%lmaxd)
    REAL                  ::  grid(atoms%jmtd)
    REAL                  ::  vr(atoms%jmtd,0:sphhar%nlhd)
    REAL                  ::  f(atoms%jmtd,2,0:atoms%lmaxd),g(atoms%jmtd,2,0:atoms%lmaxd)
    REAL                  ::  flo(atoms%jmtd,2,atoms%nlod)
    REAL                  ::  uuilon(atoms%nlod,atoms%ntype),duilon(atoms%nlod,atoms%ntype)
    REAL                  ::  ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)


    REAL                  ::  bas1(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype),&
         bas2(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype)

    COMPLEX               ::  vpw(stars%ng3)
    COMPLEX               ::  vxc(hamovlp%matsize)
    COMPLEX               ::  vrmat(hybrid%maxlmindx,hybrid%maxlmindx)
    COMPLEX               ::  carr(hybrid%maxlmindx,DIMENSION%nvd),carr1(DIMENSION%nvd,DIMENSION%nvd)
    COMPLEX ,ALLOCATABLE  ::  ahlp(:,:,:),bhlp(:,:,:)
    COMPLEX, ALLOCATABLE  ::  bascof(:,:,:)
    COMPLEX               ::  bascof_lo(3,-atoms%llod:atoms%llod,4*atoms%llod+2,atoms%nlod, atoms%nat)



    CALL timestart("subvxc")
    vxc=0


    !  calculate radial functions
    hybrid%nindx      = 2
    DO itype = 1,atoms%ntype


       !
       !--->    generate the radial basis-functions for each l
       !

       WRITE(6,'(a,i3,a)') new_LINE('n')//new_LINE('n')//' wavefunction parameters for atom type',itype,':'
       WRITE(6,'(31x,a,32x,a)') 'radial function','energy derivative'
       WRITE(6,'(a)') '  l    energy            value        '//&
            'derivative    nodes          value        derivative    nodes       norm        wronskian'
       DO l = 0,atoms%lmax(itype)
          CALL radfun(l,itype,jsp,el(l,itype,jsp),vr0(1,itype,jsp),atoms,&
               f(1,1,l),g(1,1,l),usdus,&
               nodeu,noded,wronk)
          WRITE (6,FMT=8010) l,el(l,itype,jsp),usdus%us(l,itype,jsp),&
               usdus%dus(l,itype,jsp),nodeu,usdus%uds(l,itype,jsp),usdus%duds(l,itype,jsp),noded,&
               usdus%ddn(l,itype,jsp),wronk
113
       END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
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
       !  8000    FORMAT (1x,/,/,' wavefunction parameters for atom type',i3,':',
       !      +          /,t32,'radial function',t79,'energy derivative',/,t3,
       !      +          'l',t8,'energy',t26,'value',t39,'derivative',t53,
       !      +          'nodes',t68,'value',t81,'derivative',t95,'nodes',t107,
       !      +          'norm',t119,'wronskian')
8010   FORMAT (i3,f10.5,2 (5x,1p,2e16.7,i5),1p,2e16.7)

       bas1(:,1,:,itype)=f(:,1,:)
       bas1(:,2,:,itype)=g(:,1,:)
       bas2(:,1,:,itype)=f(:,2,:)
       bas2(:,2,:,itype)=g(:,2,:)

       !
       !--->   generate the extra radial basis-functions for the local orbitals,
       !--->   if there are any.
       !        
       IF (atoms%nlo(itype).GE.1) THEN

          CALL radflo(&
               atoms,itype,jsp,ello(1,1,jsp),vr0(1,itype,jsp),&
               f,g,mpi,&
               usdus,&
               uuilon,duilon,ulouilopn,flo,.TRUE.)

          DO i=1,atoms%nlo(itype)
             hybrid%nindx(atoms%llo(i,itype),itype) = hybrid%nindx(atoms%llo(i,itype),itype) + 1
             pointer_lo(i,itype)       = hybrid%nindx(atoms%llo(i,itype),itype)
             bas1(:,hybrid%nindx(atoms%llo(i,itype),itype),atoms%llo(i,itype),itype)=&
                  flo(:,1,i)
             bas2(:,hybrid%nindx(atoms%llo(i,itype),itype),atoms%llo(i,itype),itype)=&
                  flo(:,2,i)
          END DO
       END IF
    END DO
148 149


Daniel Wortmann's avatar
Daniel Wortmann committed
150
    ! compute APW coefficients
151

Daniel Wortmann's avatar
Daniel Wortmann committed
152 153 154
    !  calculate bascof
    ALLOCATE( ahlp(DIMENSION%nvd,0:DIMENSION%lmd,atoms%nat),bhlp(DIMENSION%nvd,0:DIMENSION%lmd,atoms%nat),stat=ok)
    IF( ok .NE. 0 ) STOP 'subvxc: error in allocation of ahlp/bhlp'
155

Daniel Wortmann's avatar
Daniel Wortmann committed
156 157
    CALL abcof3( input,atoms,sym,jsp,cell, bk,lapw,&
         usdus, kveclo,oneD,ahlp,bhlp,bascof_lo)
158

Daniel Wortmann's avatar
Daniel Wortmann committed
159 160 161 162
    ALLOCATE( bascof(DIMENSION%nvd,2*(DIMENSION%lmd+1),atoms%nat), stat=ok )
    IF( ok .NE. 0 ) STOP 'subvxc: error in allocation of bascof'
    bascof = 0
    ic     = 0
163

Daniel Wortmann's avatar
Daniel Wortmann committed
164 165
    DO itype=1,atoms%ntype
       DO ieq=1,atoms%neq(itype)
166 167 168
          ic   = ic + 1
          indx = 0
          DO l=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
169 170 171 172 173 174 175 176 177 178 179 180
             ll = l*(l+1)
             DO M=-l,l
                lm=ll+M 
                DO i=1,2
                   indx = indx + 1
                   IF( i .EQ. 1) THEN
                      bascof(:,indx,ic) = ahlp(:,lm,ic)
                   ELSE IF( i .EQ. 2 ) THEN
                      bascof(:,indx,ic) = bhlp(:,lm,ic)
                   END IF
                END DO
             END DO
181
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
182 183 184 185
       END DO
    END DO

    DEALLOCATE( ahlp,bhlp )
186

Daniel Wortmann's avatar
Daniel Wortmann committed
187 188 189
    ! Loop over atom types
    iatom = 0
    DO itype = 1,atoms%ntype
190

Daniel Wortmann's avatar
Daniel Wortmann committed
191 192
       typsym = atoms%ntypsy( SUM(atoms%neq(:itype-1))+1 )
       nlharm = sphhar%nlh(typsym)
193

Daniel Wortmann's avatar
Daniel Wortmann committed
194 195
       ! Calculate vxc = vtot - vcoul
       DO l=0,nlharm
196
          DO i=1,atoms%jri(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
197 198 199 200 201 202 203
             IF(l.EQ.0) THEN
                !               vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) -  vrcou(i,0,itype,jsp)   
                vr(i,0)=  vx%mt(i,0,itype,jsp)*sfp_const/atoms%rmsh(i,itype)  !
             ELSE                                              ! vxc = vtot - vcoul
                !               vr(i,l)= vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
                vr(i,l)=  vx%mt(i,l,itype,jsp)      
             END IF
204
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
205 206
       END DO

207

Daniel Wortmann's avatar
Daniel Wortmann committed
208 209 210
       ! Calculate MT contribution to vxc matrix elements
       ! Precompute auxiliary radial integrals
       DO ilharm = 0,nlharm
211 212
          i = 0
          DO l1 = 0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228
             DO p1 = 1,2
                i = i + 1
                j = 0
                DO l2 = 0,atoms%lmax(itype)     
                   DO p2 = 1,2
                      j = j + 1
                      IF( j .LE. i) THEN
                         DO igrid = 1,atoms%jri(itype)
                            grid(igrid)=vr(igrid,ilharm)*(bas1(igrid,p1,l1,itype)*bas1(igrid,p2,l2,itype)+ bas2(igrid,p1,l1,itype)*bas2(igrid,p2,l2,itype) )
                         END DO

                         CALL intgr3(grid,atoms%rmsh(:,itype),atoms%dx(itype),atoms%jri(itype),integ(ilharm,p1,l1,p2,l2) ) ! numerical integration

                         integ(ilharm,p2,l2,p1,l1)=integ(ilharm,p1,l1,p2,l2)
                      END IF
                   END DO
229
                END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
230 231

             END DO
232
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
233 234 235 236 237 238 239
       END DO

       ! Calculate muffin tin contribution to vxc matrix
       vrmat=0

       j1=0
       DO l1 = 0,atoms%lmax(itype) ! loop: left basis function
240
          DO m1 = -l1,l1
Daniel Wortmann's avatar
Daniel Wortmann committed
241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
             DO p1 = 1,2
                j1 = j1+1
                j2 = 0
                DO l2 = 0,atoms%lmax(itype) ! loop: right basis function
                   DO m2 = -l2,l2
                      DO p2 = 1,2
                         j2 = j2+1
                         rr = 0
                         DO ilharm = 0,nlharm ! loop: lattice harmonics of vxc
                            l = sphhar%llh(ilharm,typsym)
                            DO i = 1,sphhar%nmem(ilharm,typsym)
                               M  = sphhar%mlh(i,ilharm,typsym)
                               rc = sphhar%clnu(i,ilharm,typsym)* gaunt1(l1,l,l2,m1,M,m2,atoms%lmaxd)
                               rr = rr+integ(ilharm,p1,l1,p2,l2)*rc
                            END DO
                         END DO

                         rc           = CMPLX(0,1)**(l2-l1) ! adjusts to a/b/ccof-scaling
                         vrmat(j1,j2) = rr*rc
260

Daniel Wortmann's avatar
Daniel Wortmann committed
261 262
                      END DO
                   END DO
263 264
                END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
265
             END DO
266
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
267 268
       END DO
       nnbas = j1
269

Daniel Wortmann's avatar
Daniel Wortmann committed
270 271
       !        ! Project on bascof
       DO ineq = 1,atoms%neq(itype)
272 273
          iatom = iatom+1

Daniel Wortmann's avatar
Daniel Wortmann committed
274
          carr (:nnbas,:lapw%nv(jsp)) = CONJG(MATMUL(vrmat(:nnbas,:nnbas), TRANSPOSE(bascof(:lapw%nv(jsp),:nnbas,iatom)) ))
275

Daniel Wortmann's avatar
Daniel Wortmann committed
276
          carr1(:lapw%nv(jsp),:lapw%nv(jsp)) = MATMUL(bascof(:lapw%nv(jsp),:nnbas,iatom),carr(:nnbas,:lapw%nv(jsp)) )
277 278
          ic    = 0
          DO j = 1,lapw%nv(jsp)
Daniel Wortmann's avatar
Daniel Wortmann committed
279 280 281 282 283 284 285 286
             !            carr(:nnbas) =  matmul(vrmat(:nnbas,:nnbas),
             !     +                             bascof(j,:nnbas,iatom) )
             DO i = 1,j
                ic = ic + 1
                vxc(ic) = vxc(ic) + carr1(i,j)
                !             vxc(ic) = vxc(ic) + conjg(dotprod ( bascof(i,:nnbas,iatom),
                !     +                                           carr(:nnbas) ))
             END DO
287
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303
       END DO

    END DO ! End loop over atom types

    ! ---------------------------------------------------------------
    ! Calculate plane wave contribution
    DO i=1,stars%ng3
       vpw(i)= vx%pw(i,jsp) 
       !         vpw(i)=vpwtot(i)-vpwcou(i,jsp)      
    END DO

    ! Calculate vxc-matrix,  left basis function (ig1)
    !                        right basis function (ig2)
    ic = 0
    DO ig1=1,lapw%nv(jsp)
       DO ig2=1,ig1
304 305 306 307 308
          ic = ic + 1
          gg(1)=lapw%k1(ig1,jsp)-lapw%k1(ig2,jsp)
          gg(2)=lapw%k2(ig1,jsp)-lapw%k2(ig2,jsp)
          gg(3)=lapw%k3(ig1,jsp)-lapw%k3(ig2,jsp)
          istar=stars%ig(gg(1),gg(2),gg(3))
Daniel Wortmann's avatar
Daniel Wortmann committed
309 310
          IF(istar.NE.0) THEN
             vxc(ic)= vxc(ic) + stars%rgphs(gg(1),gg(2),gg(3))*vpw(istar)
311
          ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
312 313 314
             IF ( mpi%irank == 0 ) WRITE(6,'(A,/6I5)') 'Warning: Gi-Gj not in any star:',&
                  lapw%k1(ig1,jsp),lapw%k2(ig1,jsp),lapw%k3(ig1,jsp),&
                  lapw%k1(ig2,jsp),lapw%k2(ig2,jsp),lapw%k3(ig2,jsp)
315
          ENDIF
Daniel Wortmann's avatar
Daniel Wortmann committed
316 317 318 319 320 321
       ENDDO
    ENDDO

    !    
    ! -------------------------------------------------------------------
    ! Calculate local orbital contribution
322

Daniel Wortmann's avatar
Daniel Wortmann committed
323
    IF( ANY( atoms%nlo .NE. 0) ) THEN 
324

Daniel Wortmann's avatar
Daniel Wortmann committed
325 326 327 328 329
       nbasf0      = lapw%nv(jsp)*(lapw%nv(jsp)+1)/2    ! number of pure APW contributions
       icentry     = nbasf0                   ! icentry counts the entry in the matrix vxc
       iatom       = 0
       ikvecat     = 0
       ikvecprevat = 0
330

Daniel Wortmann's avatar
Daniel Wortmann committed
331
       DO itype=1,atoms%ntype
332

Daniel Wortmann's avatar
Daniel Wortmann committed
333
          typsym = atoms%ntypsy(SUM(atoms%neq(:itype-1))+1)
334 335 336 337
          nlharm = sphhar%nlh(typsym)

          ! Calculate vxc = vtot - vcoul
          DO l=0,nlharm
Daniel Wortmann's avatar
Daniel Wortmann committed
338 339 340 341 342 343 344 345 346
             DO i=1,atoms%jri(itype)
                IF(l.EQ.0) THEN
                   !                 vr(i,0)= vrtot(i,0,itype)*sfp/rmsh(i,itype) -  vrcou(i,0,itype,jsp)
                   vr(i,0)=  vx%mt(i,0,itype,jsp)*sfp_const/atoms%rmsh(i,itype)  !
                ELSE                                              ! vxc = vtot - vcoul
                   vr(i,l)=  vx%mt(i,l,itype,jsp)                    !
                   !                 vr(i,l)=  vrtot(i,l,itype)-vrcou(i,l,itype,jsp)
                END IF
             END DO
347 348 349 350
          END DO

          ! Precompute auxiliary radial integrals
          DO ilharm=0,nlharm
Daniel Wortmann's avatar
Daniel Wortmann committed
351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367
             i = 0
             DO l1=0,atoms%lmax(itype)
                DO p1=1,hybrid%nindx(l1,itype)
                   i = i + 1
                   j = 0
                   DO l2=0,atoms%lmax(itype)     
                      DO p2=1,hybrid%nindx(l2,itype)
                         j = j + 1
                         IF( j .LE. i) THEN
                            DO igrid=1,atoms%jri(itype)
                               grid(igrid)=vr(igrid,ilharm)*(bas1(igrid,p1,l1,itype)*bas1(igrid,p2,l2,itype)+ bas2(igrid,p1,l1,itype)*bas2(igrid,p2,l2,itype))
                            END DO

                            CALL intgr3(grid,atoms%rmsh(:,itype),atoms%dx(itype),atoms%jri(itype),integ(ilharm,p1,l1,p2,l2) ) ! numerical integration

                            integ(ilharm,p2,l2,p1,l1) = integ(ilharm,p1,l1,p2,l2)
                         END IF
368
                      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
369
                   END DO
370 371

                END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
372
             END DO
373 374 375 376
          END DO


          DO ieq = 1,atoms%neq(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
377 378
             iatom = iatom + 1
             IF( (atoms%invsat(iatom).EQ.0) .OR. (atoms%invsat(iatom) .EQ. 1) ) THEN
379

Daniel Wortmann's avatar
Daniel Wortmann committed
380 381
                IF( atoms%invsat(iatom) .EQ. 0 ) invsfct = 1
                IF( atoms%invsat(iatom) .EQ. 1 ) invsfct = 2
382

Daniel Wortmann's avatar
Daniel Wortmann committed
383 384 385 386 387 388 389 390 391

                DO ilo = 1,atoms%nlo(itype)
                   l1 = atoms%llo(ilo,itype)
                   DO ikvec = 1,invsfct*(2*l1+1)

                      DO m1 = -l1,l1
                         DO p1 = 1,3
                            IF( p1 .EQ. 3) THEN
                               pp1 = pointer_lo(ilo,itype)
392
                            ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
393
                               pp1 = p1
394 395
                            END IF

Daniel Wortmann's avatar
Daniel Wortmann committed
396 397 398 399 400 401 402 403 404 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 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525
                            IF( hybrid%nindx(l1,itype) .LE. 2) STOP 'subvxc: error hybrid%nindx'

                            lm = 0

                            !loop over APW
                            DO l2 = 0,atoms%lmax(itype)
                               DO m2 = -l2,l2
                                  DO p2 = 1,2
                                     lm = lm + 1

                                     rr = 0
                                     DO ilharm = 0,nlharm
                                        lh = sphhar%llh(ilharm,typsym)
                                        DO i = 1,sphhar%nmem(ilharm,typsym)
                                           mh = sphhar%mlh(i,ilharm,typsym)
                                           rc = sphhar%clnu(i,ilharm,typsym)* gaunt1(l1,lh,l2,m1,mh,m2,atoms%lmaxd)
                                           rr = rr+integ(ilharm,p2,l2,pp1,l1)*rc
                                        END DO
                                     END DO

                                     rc = CMPLX(0d0,1d0)**(l2-l1) ! adjusts to a/b/ccof-scaling

                                     ! ic counts the entry in vxc
                                     ic = icentry
                                     DO i=1,lapw%nv(jsp)
                                        ic = ic + 1
                                        IF (hamovlp%l_real) THEN
                                           vxc(ic) = vxc(ic) + invsfct * REAL(rr*rc*bascof(i,lm,iatom) * CONJG(bascof_lo(p1,m1,ikvec,ilo, iatom)))
                                        ELSE
                                           vxc(ic) = vxc(ic) + rr*rc*bascof(i,lm,iatom) *CONJG(bascof_lo(p1,m1,ikvec,ilo, iatom))
                                        ENDIF
                                     END DO

                                  END DO  !p2
                               END DO  ! m2
                            END DO ! l2 ->  loop over APW


                            ! calcualte matrix-elements with local orbitals at the same atom
                            IF( ic .NE. icentry + lapw%nv(jsp) ) STOP 'subvxc: error counting ic'

                            ic = ic + ikvecprevat

                            DO ilop = 1,(ilo-1)
                               lp = atoms%llo(ilop,itype)

                               DO ikvecp = 1,invsfct*(2*lp+1)

                                  ic = ic + 1

                                  DO mp = -lp,lp
                                     DO pp = 1,3
                                        IF ( pp .EQ. 3) THEN
                                           pp2 = pointer_lo(ilop,itype)
                                        ELSE
                                           pp2 = pp
                                        END IF

                                        rr = 0
                                        DO ilharm = 0,nlharm
                                           lh = sphhar%llh(ilharm,typsym)
                                           DO i = 1,sphhar%nmem(ilharm,typsym)
                                              mh = sphhar%mlh(i,ilharm,typsym)
                                              rc = sphhar%clnu(i,ilharm,typsym)* gaunt1(l1,lh,lp,m1,mh,mp,atoms%lmaxd)
                                              rr = rr+integ(ilharm,pp2,lp,pp1,l1)*rc
                                           END DO
                                        END DO

                                        rc = CMPLX(0d0,1d0)**(lp-l1) ! adjusts to a/b/ccof-scaling

                                        IF (hamovlp%l_real) THEN

                                           vxc(ic) = vxc(ic) + invsfct * REAL( rr*rc*bascof_lo(pp,mp, ikvecp,ilop,iatom) * CONJG(bascof_lo(p1,m1,ikvec,ilo, iatom)) )
                                        ELSE
                                           vxc(ic) = vxc(ic) + rr*rc*bascof_lo(pp,mp,ikvecp, ilop,iatom) *CONJG(bascof_lo(p1,m1,ikvec,ilo, iatom))
                                        ENDIF
                                     END DO ! pp
                                  END DO ! mp

                               END DO !ikvecp
                            END DO ! ilop

                            ! calculate matrix-elements of one local orbital with itself

                            DO ikvecp = 1,ikvec
                               ic = ic + 1

                               lp   = l1
                               ilop = ilo
                               DO mp = -lp,lp
                                  DO pp = 1,3
                                     IF ( pp .EQ. 3) THEN
                                        pp2 = pointer_lo(ilop,itype)
                                     ELSE
                                        pp2 = pp
                                     END IF

                                     rr = 0
                                     DO ilharm = 0,nlharm
                                        lh = sphhar%llh(ilharm,typsym)
                                        DO i = 1,sphhar%nmem(ilharm,typsym)
                                           mh = sphhar%mlh(i,ilharm,typsym)
                                           rc = sphhar%clnu(i,ilharm,typsym)* gaunt1(l1,lh,lp,m1,mh,mp,atoms%lmaxd)
                                           rr = rr+integ(ilharm,pp2,lp,pp1,l1)*rc
                                        END DO
                                     END DO

                                     rc = CMPLX(0d0,1d0)**(lp-l1) ! adjusts to a/b/ccof-scaling

                                     IF (hamovlp%l_real) THEN
                                        vxc(ic) = vxc(ic) + invsfct*REAL( rr*rc* bascof_lo(pp,mp,ikvecp,ilop,iatom) * CONJG(bascof_lo(p1,m1,ikvec,ilo, iatom)) )
                                     ELSE
                                        vxc(ic) = vxc(ic) + rr*rc*bascof_lo(pp,mp,ikvecp,ilop, iatom) * CONJG(bascof_lo(p1,m1,ikvec,ilo, iatom))
                                     ENDIF
                                  END DO ! pp
                               END DO ! mp

                            END DO ! ikvecp


                         END DO  ! p1
                      END DO  ! m1

                      icentry = ic       
                   END DO !ikvec
                   ikvecat = ikvecat + invsfct*(2*l1+1)
                END DO  ! ilo
                ikvecprevat = ikvecprevat + ikvecat
                ikvecat     = 0
             END IF  ! atoms%invsat(iatom)
526 527

          END DO ! ieq
Daniel Wortmann's avatar
Daniel Wortmann committed
528 529 530 531 532
       END DO !itype

    END IF ! if any atoms%llo

    !initialize weighting factor
533 534
    a_ex=xcpot%get_exchange_weight()
    
Daniel Wortmann's avatar
Daniel Wortmann committed
535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551
    IF (hamovlp%l_real) THEN
       DO i=1,hamovlp%matsize
          hamovlp%a_r(i) = hamovlp%a_r(i) - a_ex*REAL(vxc(i))
       ENDDO
    ELSE
       DO i=1,hamovlp%matsize
          hamovlp%a_c(i) = hamovlp%a_c(i) - a_ex*vxc(i)
       ENDDO
    ENDIF

    CALL timestop("subvxc")

    DEALLOCATE( bascof )

  END SUBROUTINE subvxc

END MODULE m_subvxc