spmvec.F90 20.6 KB
Newer Older
1 2 3
      MODULE m_spmvec
      
      CONTAINS
Daniel Wortmann's avatar
Daniel Wortmann committed
4 5
        !Note this module contains a real/complex version of spmvec 
        
6
       
Daniel Wortmann's avatar
Daniel Wortmann committed
7
      SUBROUTINE spmvec_invs(&
Daniel Wortmann's avatar
Daniel Wortmann committed
8 9
     &           atoms,hybrid,&
     &           hybdat,ikpt,kpts,cell,&
10 11 12 13 14 15 16 17
     &           coulomb_mt1,coulomb_mt2,coulomb_mt3,&
     &           coulomb_mtir,vecin,&
     &           vecout)
      
      USE m_wrapper
      USE m_constants
      USE m_types
      IMPLICIT NONE
Daniel Wortmann's avatar
Daniel Wortmann committed
18
      TYPE(t_hybdat),INTENT(IN)   :: hybdat
19
      TYPE(t_hybrid),INTENT(IN)   :: hybrid
Daniel Wortmann's avatar
Daniel Wortmann committed
20 21 22
      TYPE(t_cell),INTENT(IN)     :: cell
      TYPE(t_kpts),INTENT(IN)     :: kpts
      TYPE(t_atoms),INTENT(IN)    :: atoms
23 24 25
      
      ! - scalars -
      INTEGER, INTENT(IN) ::  ikpt 
Daniel Wortmann's avatar
Daniel Wortmann committed
26
          
27
      ! - arrays -
Daniel Wortmann's avatar
Daniel Wortmann committed
28 29
      REAL   , INTENT(IN) ::  coulomb_mt1(hybrid%maxindxm1-1,hybrid%maxindxm1-1,&
     &                                    0:hybrid%maxlcutm1,atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
30
      REAL   , INTENT(IN) ::  coulomb_mt2(hybrid%maxindxm1-1,-hybrid%maxlcutm1:hybrid%maxlcutm1,&
Daniel Wortmann's avatar
Daniel Wortmann committed
31
     &                                    0:hybrid%maxlcutm1+1,atoms%nat)
Daniel Wortmann's avatar
Daniel Wortmann committed
32
      REAL   , INTENT(IN) ::  coulomb_mt3(hybrid%maxindxm1-1,atoms%nat,atoms%nat)
33 34 35 36 37
#ifdef CPP_IRCOULOMBAPPROX
      REAL   , INTENT(IN) ::  coulomb_mtir(:,:)
#else
      REAL   , INTENT(IN) ::  coulomb_mtir(:)
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
38 39
      REAL   , INTENT(IN) ::  vecin (:)!(hybrid%nbasm)
      REAL   , INTENT(OUT)::  vecout(:)!(hybrid%nbasm)
Daniel Wortmann's avatar
Daniel Wortmann committed
40 41 42 43 44 45 46 47 48 49 50 51 52 53
      
      ! - local scalars -
      INTEGER             ::  itype,ieq,iatom,ishift
      INTEGER             ::  itype1,ieq1,iatom1,ishift1
      INTEGER             ::  indx0,indx1,indx2,indx3,indx4
      INTEGER             ::  i,ibasm,igptm
      INTEGER             ::  l
      INTEGER             ::  n,m
      
      REAL                ::  gnorm
      
      COMPLEX, PARAMETER  ::  img = (0d0,1d0)
      ! - local arrays -
      
Daniel Wortmann's avatar
Daniel Wortmann committed
54
      REAL                ::  vecinhlp (hybrid%nbasm(ikpt))
Daniel Wortmann's avatar
Daniel Wortmann committed
55 56 57 58 59 60
      REAL   ,ALLOCATABLE ::  coulhlp(:,:)
 
  
      
      vecinhlp = vecin
      
Daniel Wortmann's avatar
Daniel Wortmann committed
61
      CALL reorder(hybrid%nbasm(ikpt),hybrid%nbasp,atoms,hybrid%lcutm1, hybrid%maxlcutm1,hybrid%nindxm1,1, vecinhlp)
Daniel Wortmann's avatar
Daniel Wortmann committed
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 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 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
  

      ibasm = 0
      iatom = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          DO l = 0,hybrid%lcutm1(itype)
            DO m = -l,l
              ibasm = ibasm + hybrid%nindxm1(l,itype) - 1
            END DO
          END DO
        END DO
      END DO
      
      
      ! compute vecout for the indices from 0:ibasm
      iatom = 0
      indx1 = 0; indx2 = 0; indx3 = ibasm
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          DO l = 0,hybrid%lcutm1(itype)
            DO m = -l,l
              indx1 = indx1 + 1
              indx2 = indx2 + hybrid%nindxm1(l,itype) - 1
              indx3 = indx3 + 1

              
              vecout(indx1:indx2) = matmul(coulomb_mt1(:hybrid%nindxm1(l,itype)-1,:hybrid%nindxm1(l,itype)-1, l,itype),&
     &                vecinhlp(indx1:indx2))

              
              
              vecout(indx1:indx2) = vecout(indx1:indx2) + coulomb_mt2(:hybrid%nindxm1(l,itype)-1,m,l,iatom) * vecinhlp(indx3) 
            
              indx1 = indx2
            END DO

            
          END DO
        END DO
      END DO
     
      IF( indx2 .ne. ibasm ) STOP 'spmvec: error counting basis functions'
      
      IF( ikpt .eq. 1) THEN       
        iatom = 0
        indx0 = 0
        DO itype = 1,atoms%ntype
          ishift = sum( (/ ((2*l+1)*(hybrid%nindxm1(l,itype)-1), l=0,hybrid%lcutm1(itype) )/) )
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            l = 0
            m = 0
            
            indx1 = indx0 + 1 
            indx2 = indx1 + hybrid%nindxm1(l,itype) - 2
            
            iatom1 = 0
            indx3  = ibasm
            DO itype1 = 1,atoms%ntype
              ishift1 = (hybrid%lcutm1(itype1)+1)**2
              DO ieq1 = 1,atoms%neq(itype1)
                iatom1 = iatom1 + 1
                indx4  = indx3 + (ieq1 - 1)*ishift1 + 1
                IF( iatom .eq. iatom1 ) CYCLE
                
                vecout(indx1:indx2) = vecout(indx1:indx2) + coulomb_mt3(:hybrid%nindxm1(l,itype)-1, iatom1,iatom) * vecinhlp(indx4)

              END DO
              indx3 = indx3 + atoms%neq(itype1)*ishift1
            END DO
            
            IF( indx3 .ne. hybrid%nbasp ) STOP 'spmvec: error counting index indx3'
     
            vecout(indx1:indx2) = vecout(indx1:indx2) + coulomb_mt2(:hybrid%nindxm1(l,itype)-1,0, hybrid%maxlcutm1+1,iatom) * vecinhlp(indx3+1)
     
            indx0 = indx0 + ishift
          END DO

        END DO
      END IF

      ! compute vecout for the index-range from ibasm+1:nbasm
      
#ifdef CPP_IRCOULOMBAPPROX
        
      indx0 = sum( (/ ( ((2*l+1)*atoms%neq(itype),l=0,hybrid%lcutm1(itype)), itype=1,atoms%ntype ) /) )+ hybrid%ngptm
      indx1 = sum( (/ ( ((2*l+1)*atoms%neq(itype),l=0,hybrid%lcutm1(itype)), itype=1,atoms%ntype ) /) )

      CALL DGEMV('N',indx1,indx0,1d0,coulomb_mtir,(hybrid%maxlcutm1+1)**2*atoms,&
     &          vecinhlp(ibasm+1:),1,0d0,vecout(ibasm+1:),1)
      
      CALL DGEMV('T',indx1,hybrid,1d0,coulomb_mtir(:indx1,indx1+1:),&
     &          indx1,vecinhlp(ibasm+1:),1,0d0,vecout(ibasm+indx1+1:),1)

!       vecout(ibasm+1:ibasm+indx1) = matmul( coulomb_mtir(:indx1,:indx0),vecinhlp(ibasm+1:ibasm+indx0) )
!       vecout(ibasm+indx1+1:ibasm+indx0) = matmul( conjg(transpose(coulomb_mtir(:indx1,indx1+1:indx0))),
!      &                                            vecinhlp(ibasm+1:ibasm+indx1) )

      
      indx0 = ibasm + indx1
      IF( indx0 .ne. hybrid%nbasp ) STOP 'spmvec: error indx0'
      DO i = 1,hybrid%ngptm 
        indx0 = indx0 + 1
        igptm = hybrid%pgptm(i)
        gnorm = sqrt(sum(matmul(kpts%bk(:) + hybrid%gptm(:,igptm),cell%bmat)**2))
        IF( gnorm .eq. 0 ) CYCLE
        vecout(indx0) = vecout(indx0) + fpi*vecinhlp(indx0)/gnorm
      END DO
      
174
#else
Daniel Wortmann's avatar
Daniel Wortmann committed
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246

      indx1 = sum( (/ ( ((2*l+1)*atoms%neq(itype),l=0,hybrid%lcutm1(itype)),&
     &                                      itype=1,atoms%ntype ) /) )+ hybrid%ngptm(ikpt)
      CALL dspmv('U',indx1,1d0,coulomb_mtir,vecinhlp(ibasm+1),1,0d0, vecout(ibasm+1),1)

#endif
      iatom = 0
      indx1 = ibasm; indx2 = 0; indx3 = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
          DO l = 0,hybrid%lcutm1(itype)
            n = hybrid%nindxm1(l,itype)
            DO m = -l,l
              indx1 = indx1 + 1
              indx2 = indx2 + 1
              indx3 = indx3 + n - 1
              
              vecout(indx1) = vecout(indx1) + dotprod(coulomb_mt2(:n-1,m,l,iatom), vecinhlp(indx2:indx3))
              indx2 = indx3
            END DO
             
          END DO
        END DO
      END DO
      
      
      IF( ikpt .eq. 1 ) THEN
        iatom = 0
        indx0 = 0
        DO itype = 1,atoms%ntype
          ishift = sum( (/ ((2*l+1)*(hybrid%nindxm1(l,itype)-1), l=0,hybrid%lcutm1(itype) )/) )
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            indx1 = indx0 + 1
            indx2 = indx1 + hybrid%nindxm1(0,itype) - 2
            vecout(hybrid%nbasp+1) = vecout(hybrid%nbasp+1) + dotprod(coulomb_mt2(:hybrid%nindxm1(0,itype)-1,0, hybrid%maxlcutm1+1,iatom), vecinhlp(indx1:indx2))
        
           indx0 = indx0 + ishift
          END DO
        END DO
         
        iatom = 0
        indx0 = ibasm
        DO itype = 1,atoms%ntype
          ishift = (hybrid%lcutm1(itype)+1)**2
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            indx1 = indx0 + 1
            
            iatom1 = 0
            indx2  = 0
            DO itype1 = 1,atoms%ntype
              ishift1 = sum( (/ ( (2*l+1)*(hybrid%nindxm1(l,itype1)-1), l=0,hybrid%lcutm1(itype1)) /))
              DO ieq1 = 1,atoms%neq(itype1)
                iatom1 = iatom1 + 1
                IF( iatom1 .eq. iatom ) CYCLE
                
                indx3  = indx2 + (ieq1 - 1)*ishift1 + 1
                indx4  = indx3 + hybrid%nindxm1(0,itype1) - 2

                vecout(indx1) = vecout(indx1) + dotprod(coulomb_mt3(:hybrid%nindxm1(0,itype1)-1,iatom,iatom1), vecinhlp(indx3:indx4))

              END DO
              indx2 = indx2 + atoms%neq(itype1)*ishift1
            END DO
            indx0 = indx0 + ishift
          END DO
        END DO
        IF( indx0 .ne. hybrid%nbasp ) STOP 'spmvec: error index counting (indx0)'
      END IF
      
Daniel Wortmann's avatar
Daniel Wortmann committed
247
      CALL reorder(hybrid%nbasm(ikpt),hybrid%nbasp,atoms,hybrid%lcutm1,hybrid%maxlcutm1,hybrid%nindxm1,&
Daniel Wortmann's avatar
Daniel Wortmann committed
248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277
     &             2,&
     &             vecout)
      

    END SUBROUTINE spmvec_invs

    SUBROUTINE spmvec_noinvs(&
     &           atoms,hybrid,&
     &           hybdat,ikpt,kpts,cell,&
     &           coulomb_mt1,coulomb_mt2,coulomb_mt3,&
     &           coulomb_mtir,vecin,&
     &           vecout)
      
      USE m_wrapper
      USE m_constants
      USE m_types
      IMPLICIT NONE
      TYPE(t_hybdat),INTENT(IN)   :: hybdat
      TYPE(t_hybrid),INTENT(IN)   :: hybrid
      TYPE(t_cell),INTENT(IN)     :: cell
      TYPE(t_kpts),INTENT(IN)     :: kpts
      TYPE(t_atoms),INTENT(IN)    :: atoms
      
      ! - scalars -
      INTEGER, INTENT(IN) ::  ikpt 
          
      ! - arrays -
      REAL   , INTENT(IN) ::  coulomb_mt1(hybrid%maxindxm1-1,hybrid%maxindxm1-1,&
     &                                    0:hybrid%maxlcutm1,atoms%ntype)
      COMPLEX, INTENT(IN) ::  coulomb_mt2(hybrid%maxindxm1-1,-hybrid%maxlcutm1:hybrid%maxlcutm1,&
Daniel Wortmann's avatar
Daniel Wortmann committed
278 279
     &                                    0:hybrid%maxlcutm1+1,atoms%nat)
      COMPLEX, INTENT(IN) ::  coulomb_mt3(hybrid%maxindxm1-1,atoms%nat,atoms%nat)
280 281 282 283 284
#ifdef CPP_IRCOULOMBAPPROX
      COMPLEX, INTENT(IN) ::  coulomb_mtir(:,:)
#else
      COMPLEX, INTENT(IN) ::  coulomb_mtir(:)
#endif
Daniel Wortmann's avatar
Daniel Wortmann committed
285 286
      COMPLEX, INTENT(IN) ::  vecin (:)!(hybrid%nbasm)
      COMPLEX, INTENT(OUT)::  vecout(:)!(hybrid%nbasm)
287 288 289 290 291 292 293 294 295 296 297 298 299 300
      
      ! - local scalars -
      INTEGER             ::  itype,ieq,iatom,ishift
      INTEGER             ::  itype1,ieq1,iatom1,ishift1
      INTEGER             ::  indx0,indx1,indx2,indx3,indx4
      INTEGER             ::  i,ibasm,igptm
      INTEGER             ::  l
      INTEGER             ::  n,m
      
      REAL                ::  gnorm
      
      COMPLEX, PARAMETER  ::  img = (0d0,1d0)
      ! - local arrays -
      
Daniel Wortmann's avatar
Daniel Wortmann committed
301
      REAL                ::  vecr(hybrid%maxindxm1-1),veci(hybrid%maxindxm1-1)
Daniel Wortmann's avatar
Daniel Wortmann committed
302
      COMPLEX             ::  vecinhlp (hybrid%nbasm(ikpt))
303 304 305 306 307 308
      COMPLEX,ALLOCATABLE ::  coulhlp(:,:)
 
  
      
      vecinhlp = vecin
      
Daniel Wortmann's avatar
Daniel Wortmann committed
309
      CALL reorder(hybrid%nbasm(ikpt),hybrid%nbasp,atoms,hybrid%lcutm1, hybrid%maxlcutm1,hybrid%nindxm1,1, vec_c=vecinhlp)
310 311 312 313 314 315 316
  

      ibasm = 0
      iatom = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
317
          DO l = 0,hybrid%lcutm1(itype)
318
            DO m = -l,l
Daniel Wortmann's avatar
Daniel Wortmann committed
319
              ibasm = ibasm + hybrid%nindxm1(l,itype) - 1
320 321 322 323 324 325 326 327 328 329 330 331
            END DO
          END DO
        END DO
      END DO
      
      
      ! compute vecout for the indices from 0:ibasm
      iatom = 0
      indx1 = 0; indx2 = 0; indx3 = ibasm
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
332
          DO l = 0,hybrid%lcutm1(itype)
333 334
            DO m = -l,l
              indx1 = indx1 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
335
              indx2 = indx2 + hybrid%nindxm1(l,itype) - 1
336 337 338
              indx3 = indx3 + 1

              
Daniel Wortmann's avatar
Daniel Wortmann committed
339
              vecout(indx1:indx2) = matmul(coulomb_mt1(:hybrid%nindxm1(l,itype)-1,:hybrid%nindxm1(l,itype)-1, l,itype),&
340 341 342
     &                vecinhlp(indx1:indx2))
              
              
Daniel Wortmann's avatar
Daniel Wortmann committed
343
              vecout(indx1:indx2) = vecout(indx1:indx2) + coulomb_mt2(:hybrid%nindxm1(l,itype)-1,m,l,iatom) * vecinhlp(indx3) 
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359
            
              indx1 = indx2
            END DO

            
          END DO
        END DO
      END DO
     
      IF( indx2 .ne. ibasm ) &
     &  STOP 'spmvec: error counting basis functions'
      
      IF( ikpt .eq. 1) THEN       
        iatom = 0
        indx0 = 0
        DO itype = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
360
          ishift = sum( (/ ((2*l+1)*(hybrid%nindxm1(l,itype)-1), l=0,hybrid%lcutm1(itype) )/) )
361 362 363 364 365 366
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            l = 0
            m = 0
            
            indx1 = indx0 + 1 
Daniel Wortmann's avatar
Daniel Wortmann committed
367
            indx2 = indx1 + hybrid%nindxm1(l,itype) - 2
368 369 370 371
            
            iatom1 = 0
            indx3  = ibasm
            DO itype1 = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
372
              ishift1 = (hybrid%lcutm1(itype1)+1)**2
373 374 375 376 377 378
              DO ieq1 = 1,atoms%neq(itype1)
                iatom1 = iatom1 + 1
                indx4  = indx3 + (ieq1 - 1)*ishift1 + 1
                IF( iatom .eq. iatom1 ) CYCLE
                
                vecout(indx1:indx2) = vecout(indx1:indx2)&
Daniel Wortmann's avatar
Daniel Wortmann committed
379
     &                              + coulomb_mt3(:hybrid%nindxm1(l,itype)-1, iatom1,iatom) * vecinhlp(indx4)
380 381 382 383 384

              END DO
              indx3 = indx3 + atoms%neq(itype1)*ishift1
            END DO
            
Daniel Wortmann's avatar
Daniel Wortmann committed
385
            IF( indx3 .ne. hybrid%nbasp ) STOP 'spmvec: error counting index indx3'
386
     
Daniel Wortmann's avatar
Daniel Wortmann committed
387
            vecout(indx1:indx2) = vecout(indx1:indx2) + coulomb_mt2(:hybrid%nindxm1(l,itype)-1,0, hybrid%maxlcutm1+1,iatom) * vecinhlp(indx3+1)
388 389 390 391 392 393 394 395 396 397 398
     
            indx0 = indx0 + ishift
          END DO

        END DO
      END IF

      ! compute vecout for the index-range from ibasm+1:nbasm
      
#ifdef CPP_IRCOULOMBAPPROX
        
Daniel Wortmann's avatar
Daniel Wortmann committed
399 400
      indx0 = sum( (/ ( ((2*l+1)*atoms%neq(itype),l=0,hybrid%lcutm1(itype)), itype=1,atoms%ntype ) /) )+ hybrid%ngptm
      indx1 = sum( (/ ( ((2*l+1)*atoms%neq(itype),l=0,hybrid%lcutm1(itype)), itype=1,atoms%ntype ) /) )
401 402

      CALL ZGEMV('N',indx1,indx0,(1d0,0d0),coulomb_mtir,&
Daniel Wortmann's avatar
Daniel Wortmann committed
403
     &          (hybrid%maxlcutm1+1)**2*atoms,vecinhlp(ibasm+1:),&
404 405 406 407 408 409 410 411 412 413 414 415
     &          1,(0d0,0d0),vecout(ibasm+1:),1)
      
      CALL ZGEMV('C',indx1,hybrid,(1d0,0d0),coulomb_mtir(:indx1,indx1+1:)&
     &          ,indx1,vecinhlp(ibasm+1:),1,(0d0,0d0),&
     &          vecout(ibasm+indx1+1:),1)

!       vecout(ibasm+1:ibasm+indx1) = matmul( coulomb_mtir(:indx1,:indx0),vecinhlp(ibasm+1:ibasm+indx0) )
!       vecout(ibasm+indx1+1:ibasm+indx0) = matmul( conjg(transpose(coulomb_mtir(:indx1,indx1+1:indx0))),
!      &                                            vecinhlp(ibasm+1:ibasm+indx1) )

      
      indx0 = ibasm + indx1
Daniel Wortmann's avatar
Daniel Wortmann committed
416
      IF( indx0 .ne. hybrid%nbasp ) STOP 'spmvec: error indx0'
417 418 419 420 421 422 423 424 425 426
      DO i = 1,hybrid%ngptm 
        indx0 = indx0 + 1
        igptm = hybrid%pgptm(i)
        gnorm = sqrt(sum(matmul(kpts%bk(:) + hybrid%gptm(:,igptm),cell%bmat)**2))
        IF( gnorm .eq. 0 ) CYCLE
        vecout(indx0) = vecout(indx0) + fpi*vecinhlp(indx0)/gnorm
      END DO
      
#else

Daniel Wortmann's avatar
Daniel Wortmann committed
427
      indx1 = sum( (/ ( ((2*l+1)*atoms%neq(itype),l=0,hybrid%lcutm1(itype)),&
428 429 430 431 432 433 434 435 436
     &                                      itype=1,atoms%ntype ) /) )+ hybrid%ngptm(ikpt)
      call zhpmv('U',indx1,(1d0,0d0),coulomb_mtir,vecinhlp(ibasm+1), 1,(0d0,0d0),vecout(ibasm+1),1)

#endif
      iatom = 0
      indx1 = ibasm; indx2 = 0; indx3 = 0
      DO itype = 1,atoms%ntype
        DO ieq = 1,atoms%neq(itype)
          iatom = iatom + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
437 438
          DO l = 0,hybrid%lcutm1(itype)
            n = hybrid%nindxm1(l,itype)
439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
            DO m = -l,l
              indx1 = indx1 + 1
              indx2 = indx2 + 1
              indx3 = indx3 + n - 1
              
              vecout(indx1) = vecout(indx1) + dotprod(coulomb_mt2(:n-1,m,l,iatom), vecinhlp(indx2:indx3))
              indx2 = indx3
            END DO
             
          END DO
        END DO
      END DO
      
      
      IF( ikpt .eq. 1 ) THEN
        iatom = 0
        indx0 = 0
        DO itype = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
457
          ishift = sum( (/ ((2*l+1)*(hybrid%nindxm1(l,itype)-1), l=0,hybrid%lcutm1(itype) )/) )
458 459 460
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            indx1 = indx0 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
461
            indx2 = indx1 + hybrid%nindxm1(0,itype) - 2
Daniel Wortmann's avatar
Daniel Wortmann committed
462
            vecout(hybrid%nbasp+1) = vecout(hybrid%nbasp+1) + dotprod(coulomb_mt2(:hybrid%nindxm1(0,itype)-1,0, hybrid%maxlcutm1+1,iatom), vecinhlp(indx1:indx2))
463 464 465 466 467 468 469 470
        
           indx0 = indx0 + ishift
          END DO
        END DO
         
        iatom = 0
        indx0 = ibasm
        DO itype = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
471
          ishift = (hybrid%lcutm1(itype)+1)**2
472 473 474 475 476 477 478
          DO ieq = 1,atoms%neq(itype)
            iatom = iatom + 1
            indx1 = indx0 + 1
            
            iatom1 = 0
            indx2  = 0
            DO itype1 = 1,atoms%ntype
Daniel Wortmann's avatar
Daniel Wortmann committed
479
              ishift1 = sum( (/ ( (2*l+1)*(hybrid%nindxm1(l,itype1)-1), l=0,hybrid%lcutm1(itype1)) /))
480 481 482 483 484
              DO ieq1 = 1,atoms%neq(itype1)
                iatom1 = iatom1 + 1
                IF( iatom1 .eq. iatom ) CYCLE
                
                indx3  = indx2 + (ieq1 - 1)*ishift1 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
485
                indx4  = indx3 + hybrid%nindxm1(0,itype1) - 2
486

Daniel Wortmann's avatar
Daniel Wortmann committed
487
                vecout(indx1) = vecout(indx1) + dotprod(coulomb_mt3(:hybrid%nindxm1(0,itype1)-1,iatom,iatom1), vecinhlp(indx3:indx4))
488 489 490 491 492 493 494

              END DO
              indx2 = indx2 + atoms%neq(itype1)*ishift1
            END DO
            indx0 = indx0 + ishift
          END DO
        END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
495
        IF( indx0 .ne. hybrid%nbasp ) STOP 'spmvec: error index counting (indx0)'
496 497
      END IF
      
Daniel Wortmann's avatar
Daniel Wortmann committed
498
      CALL reorder(hybrid%nbasm(ikpt),hybrid%nbasp,atoms,hybrid%lcutm1,hybrid%maxlcutm1,hybrid%nindxm1,&
Daniel Wortmann's avatar
Daniel Wortmann committed
499
     &             2,&
Daniel Wortmann's avatar
Daniel Wortmann committed
500
     &             vec_c=vecout)
501 502
      

Daniel Wortmann's avatar
Daniel Wortmann committed
503 504
    END SUBROUTINE spmvec_noinvs

505
      
Daniel Wortmann's avatar
Daniel Wortmann committed
506 507
      SUBROUTINE reorder(nbasm,nbasp,atoms,lcutm, maxlcutm,nindxm,imode, vec_r,vec_c)
      USE m_types
508 509 510 511 512 513 514 515 516 517 518 519
      IMPLICIT NONE
      TYPE(t_atoms),INTENT(IN)   :: atoms
      
      ! - scalars -
      INTEGER, INTENT(IN)   ::  maxlcutm
      INTEGER, INTENT(IN)   ::  nbasm,nbasp
      INTEGER, INTENT(IN)   ::  imode
      
        
      ! - arrays -
      INTEGER, INTENT(IN)   ::  lcutm(atoms%ntype)
      INTEGER, INTENT(IN)   ::  nindxm(0:maxlcutm,atoms%ntype)
Daniel Wortmann's avatar
Daniel Wortmann committed
520 521
      REAL   , INTENT(INOUT),OPTIONAL::  vec_r(nbasm)
      COMPLEX, INTENT(INOUT),OPTIONAL::  vec_c(nbasm)
522 523 524 525 526
      ! - local scalars -
      INTEGER               ::  itype,ieq
      INTEGER               ::  indx1,indx2
      INTEGER               ::  l
      INTEGER               ::  n,m
Daniel Wortmann's avatar
Daniel Wortmann committed
527
      LOGICAL               :: l_real
528
      ! - local arrays -
Daniel Wortmann's avatar
Daniel Wortmann committed
529 530
      REAL                  ::  vechlp_r(nbasm)
      COMPLEX               ::  vechlp_c(nbasm)
531

Daniel Wortmann's avatar
Daniel Wortmann committed
532
      l_real=PRESENT(vec_r)
533 534 535

      IF( imode .ne. 1 .and. imode .ne. 2 ) STOP 'reorder: imode equals neither 1 nor 2'

Daniel Wortmann's avatar
Daniel Wortmann committed
536 537 538 539 540
      if (l_real) THEN
         vechlp_r = vec_r
      else
         vechlp_c = vec_c
      end if
541 542 543 544 545 546 547 548 549 550 551
      
      IF( imode .eq. 1 ) THEN
        indx1 = 0
        indx2 = 0
        DO itype = 1,atoms%ntype
          DO ieq = 1,atoms%neq(itype)
            DO l = 0,lcutm(itype)
              DO m = -l,l
                DO n = 1,nindxm(l,itype)-1
                  indx1 = indx1 + 1
                  indx2 = indx2 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
552 553 554 555 556
                  if (l_real) THEN
                     vec_r(indx1) = vechlp_r(indx2)
                  else
                     vec_c(indx1) = vechlp_c(indx2)
                  endif
557 558 559 560 561 562 563 564 565 566 567 568 569 570
                END DO
                indx2 = indx2 + 1
              END DO
            END DO
          END DO
        END DO
      
        indx2 = 0
        DO itype = 1,atoms%ntype
          DO ieq = 1,atoms%neq(itype)
            DO l = 0,lcutm(itype)
              DO m = -l,l
                indx1 = indx1 + 1
                indx2 = indx2 + nindxm(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
571 572 573 574 575 576
                if (l_real) THEN
                   vec_r(indx1) = vechlp_r(indx2)
                else
                   vec_c(indx1) = vechlp_c(indx2)
                endif
                
577 578 579 580 581 582 583 584 585 586 587 588 589 590
              END DO
            END DO
          END DO
        END DO        
      ELSE IF (imode .eq. 2) THEN
        indx1 = 0
        indx2 = 0
        DO itype = 1,atoms%ntype
          DO ieq = 1,atoms%neq(itype)
            DO l = 0,lcutm(itype)
              DO m = -l,l
                DO n = 1,nindxm(l,itype)-1
                  indx1 = indx1 + 1
                  indx2 = indx2 + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
591 592 593 594 595
                  if (l_real) THEN
                   vec_r(indx2) = vechlp_r(indx1)
                else
                   vec_c(indx2) = vechlp_c(indx1)
                endif
596 597 598 599 600 601 602 603 604 605 606 607 608 609
                END DO
                indx2 = indx2 + 1
              END DO
            END DO
          END DO
        END DO
        
        indx2 = 0
        DO itype = 1,atoms%ntype
          DO ieq = 1,atoms%neq(itype)
            DO l = 0,lcutm(itype)
              DO m = -l,l
                indx1 = indx1 + 1
                indx2 = indx2 + nindxm(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
610 611 612 613 614
                if (l_real) THEN
                   vec_r(indx2) = vechlp_r(indx1)
                else
                   vec_c(indx2) = vechlp_c(indx1)
                endif
615 616 617 618 619 620 621 622 623 624
              END DO
            END DO
          END DO
        END DO        
      END IF
      !IR must not be rearranged

      END SUBROUTINE reorder
      
      END MODULE