mixedbasis.F90 58.2 KB
Newer Older
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
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This subroutine generates the mixed basis set used to evaluate the  !
! exchange term in HF/hybrid functional calculations or EXX           !
! calculations. In the latter case a second mixed basis set is setup  !
! for the OEP integral equation.                                      !
! In all cases the mixed basis consists of IR plane waves             !
!                                                                     !
! IR:                                                                 !
!    M_{\vec{k},\vec{G}} = 1/\sqrt{V} \exp{i(\vec{k}+\vec{G})}        !
!                                                                     !
! which are zero in the MT spheres and radial functions times         !
! spherical harmonics in the MT spheres                               !
!                                                                     !
! MT:                                                                 !
!     a                a                                              !
!    M              = U   Y                                           !
!     PL               PL  LM                                         !
!                                                                     !
!            where     a    a  a                                      !
!                     U  = u  u                                       !
!                      PL   pl  p'l'                                  !
!                                                                     !
!               and    L \in {|l-l'|,...,l+l'}                        !
!                                                                     !
!               and    P counts the different combinations of         !
!                      pl,p'l' which contribute to L                  !
!                                                                     !
!                                               M.Betzinger (09/07)   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Daniel Wortmann's avatar
Daniel Wortmann committed
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
MODULE m_mixedbasis
  use m_judft
CONTAINS

  SUBROUTINE mixedbasis( atoms,kpts,DIMENSION,input, &
       cell,sym,xcpot,hybrid, eig_id, mpi,v,l_restart)

    USE m_radfun,   ONLY : radfun
    USE m_radflo,   ONLY : radflo
    USE m_loddop,   ONLY : loddop
    USE m_util,     ONLY : intgrf_init,intgrf,rorderpf
    USE m_read_core
    USE m_wrapper
    USE m_icorrkeys
    USE m_eig66_io
    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_hybrid),INTENT(INOUT) :: hybrid
    TYPE(t_input),INTENT(IN)     :: input
    TYPE(t_sym),INTENT(IN)       :: sym
    TYPE(t_cell),INTENT(IN)      :: cell
Daniel Wortmann's avatar
Daniel Wortmann committed
55
    TYPE(t_kpts),INTENT(IN)      :: kpts
Daniel Wortmann's avatar
Daniel Wortmann committed
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
    TYPE(t_atoms),INTENT(IN)     :: atoms
    TYPE(t_potden),INTENT(IN)    :: v

    ! - scalars -
    INTEGER,INTENT(IN)              :: eig_id
    LOGICAL,INTENT(INOUT)           :: l_restart



    ! - local scalars -
    INTEGER                         ::  ilo
    INTEGER                         ::  ikpt
    INTEGER                         ::  ispin,itype,l1,l2,l,n,igpt,n1,n2,nn,i,j,ic ,ng      
    INTEGER                         ::  jsp
    INTEGER                         ::  nodem,noded
    INTEGER                         ::  m
    INTEGER                         ::  x,y,z
    INTEGER                         ::  maxindxc,lmaxcd
    INTEGER                         ::  divconq ! use Divide & Conquer algorithm for array sorting (>0: yes, =0: no)
    REAL                            ::  gcutm
    REAL                            ::  wronk
    REAL                            ::  rdum,rdum1,rdum2

    LOGICAL                         ::  ldum,ldum1
    LOGICAL                         ::  exx


    ! - local arrays -
    INTEGER                         ::  g(3)
    INTEGER                         ::  lmaxc(atoms%ntype)
    INTEGER,ALLOCATABLE             ::  nindxc(:,:)
    INTEGER,ALLOCATABLE             ::  ihelp(:) 
    INTEGER,ALLOCATABLE             ::  ptr(:)           ! pointer for array sorting
    INTEGER,ALLOCATABLE             ::  unsrt_pgptm(:,:) ! unsorted pointers to g vectors

    REAL                            ::  kvec(3)
    REAL                            ::  flo(atoms%jmtd,2,atoms%nlod)
    TYPE(t_usdus)                   ::  usdus
    REAL                            ::  uuilon(atoms%nlod,atoms%ntype), duilon(atoms%nlod,atoms%ntype)
    REAL                            ::  ulouilopn(atoms%nlod,atoms%nlod,atoms%ntype)
    REAL                            ::  potatom(atoms%jmtd,atoms%ntype)
    REAL                            ::  el(0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd)
    REAL                            ::  ello(atoms%nlod,atoms%ntype,DIMENSION%jspd)
    REAL                            ::  bashlp(atoms%jmtd)

    REAL   ,ALLOCATABLE             ::  f(:,:,:),df(:,:,:)
    REAL   ,ALLOCATABLE             ::  olap(:,:),work(:),eig(:), eigv(:,:)      
    REAL   ,ALLOCATABLE             ::  bas1(:,:,:,:,:), bas2(:,:,:,:,:)
    REAL   ,ALLOCATABLE             ::  basmhlp(:,:,:,:)
    REAL   ,ALLOCATABLE             ::  gridf(:,:),vr0(:,:,:)
    REAL   ,ALLOCATABLE             ::  core1(:,:,:,:,:), core2(:,:,:,:,:)
    REAL   ,ALLOCATABLE             ::  eig_c(:,:,:,:)
    REAL   ,ALLOCATABLE             ::  length_kg(:,:) ! length of the vectors k + G
109

Daniel Wortmann's avatar
Daniel Wortmann committed
110 111 112
 
    LOGICAL,ALLOCATABLE             ::  selecmat(:,:,:,:)
    LOGICAL,ALLOCATABLE             ::  seleco(:,:),selecu(:,:)
113

Daniel Wortmann's avatar
Daniel Wortmann committed
114 115 116 117 118 119 120 121 122 123 124 125
    CHARACTER, PARAMETER            :: lchar(0:38) =&
         &      (/'s','p','d','f','g','h','i','j','k','l','m','n','o',&
         &        'x','x','x','x','x','x','x','x','x','x','x','x','x',&
         &        'x','x','x','x','x','x','x','x','x','x','x','x','x' /)
    
    CHARACTER(len=2)                ::  nchar
    CHARACTER(len=2)                ::  noel(atoms%ntype)
    CHARACTER(len=10)               ::  fname(atoms%ntype)
    ! writing to a file
    INTEGER, PARAMETER              ::  iounit = 125
    CHARACTER(10), PARAMETER        ::  ioname = 'mixbas'
    LOGICAL                         ::  l_found
126

Daniel Wortmann's avatar
Daniel Wortmann committed
127
    IF ( mpi%irank == 0 ) WRITE(6,'(//A,I2,A)') '### subroutine: mixedbasis ###'
128 129


Daniel Wortmann's avatar
Daniel Wortmann committed
130
    exx = ( xcpot%icorr .EQ. icorr_exx )
131

Daniel Wortmann's avatar
Daniel Wortmann committed
132 133 134 135 136 137 138 139 140 141 142 143
    ! Deallocate arrays which might have been allocated in a previous run of this subroutine
    IF(ALLOCATED(hybrid%ngptm))    DEALLOCATE(hybrid%ngptm)
    IF(ALLOCATED(hybrid%ngptm1))   DEALLOCATE(hybrid%ngptm1)
    IF(ALLOCATED(hybrid%ngptm2))   DEALLOCATE(hybrid%ngptm2)
    IF(ALLOCATED(hybrid%nindxm1))  DEALLOCATE(hybrid%nindxm1)
    IF(ALLOCATED(hybrid%nindxm2))  DEALLOCATE(hybrid%nindxm2)
    IF(ALLOCATED(hybrid%pgptm))    DEALLOCATE(hybrid%pgptm)
    IF(ALLOCATED(hybrid%pgptm1))   DEALLOCATE(hybrid%pgptm1)
    IF(ALLOCATED(hybrid%pgptm2))   DEALLOCATE(hybrid%pgptm2)
    IF(ALLOCATED(hybrid%gptm))     DEALLOCATE(hybrid%gptm)
    IF(ALLOCATED(hybrid%basm1) )   DEALLOCATE(hybrid%basm1)
    IF(ALLOCATED(hybrid%basm2) )   DEALLOCATE(hybrid%basm2)
144

Daniel Wortmann's avatar
Daniel Wortmann committed
145 146
    call usdus%init(atoms,dimension%jspd)

Daniel Wortmann's avatar
Daniel Wortmann committed
147 148 149
    ! If restart is specified read file if it already exists
    ! create it otherwise
    IF ( l_restart ) THEN
150

Daniel Wortmann's avatar
Daniel Wortmann committed
151 152
       ! Test if file exists
       INQUIRE(FILE=ioname,EXIST=l_found)
153

Daniel Wortmann's avatar
Daniel Wortmann committed
154
       IF ( l_found .and..false.) THEN !reading not working yet
155 156 157 158 159 160 161 162 163 164 165 166 167 168 169

          ! Open file
          OPEN(UNIT=iounit,FILE=ioname,FORM='unformatted',STATUS='old')

          ! Read array sizes
          !READ(iounit) kpts%nkptf,hybrid%gptmd
          READ(iounit) hybrid%maxgptm,hybrid%maxindx
          ! Allocate necessary array size and read arrays
          ALLOCATE ( hybrid%ngptm(kpts%nkptf),hybrid%gptm(3,hybrid%gptmd),hybrid%pgptm(hybrid%maxgptm,kpts%nkptf) )
          READ(iounit) hybrid%ngptm,hybrid%gptm,hybrid%pgptm,hybrid%nindx

          ! Read array sizes
          READ(iounit) hybrid%maxgptm1,hybrid%maxlcutm1,hybrid%maxindxm1,hybrid%maxindxp1
          ! Allocate necessary array size and read arrays
          ALLOCATE ( hybrid%ngptm1(kpts%nkptf),hybrid%pgptm1(hybrid%maxgptm1,kpts%nkptf),&
Daniel Wortmann's avatar
Daniel Wortmann committed
170
               &               hybrid%nindxm1(0:hybrid%maxlcutm1,atoms%ntype) )
171 172 173 174 175
          ALLOCATE ( hybrid%basm1(atoms%jmtd,hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype) )
          READ(iounit) hybrid%ngptm1,hybrid%pgptm1,hybrid%nindxm1
          READ(iounit) hybrid%basm1

          IF ( exx ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
176 177 178 179 180 181 182 183
             ! Read array sizes
             READ(iounit) hybrid%maxgptm2,hybrid%maxlcutm2,hybrid%maxindxm2,hybrid%maxindxp2
             ! Allocate necessary array size and read arrays
             ALLOCATE ( hybrid%ngptm2(kpts%nkptf),hybrid%pgptm2(hybrid%maxgptm2,kpts%nkptf),&
                  &                 hybrid%nindxm2(0:hybrid%maxlcutm2,atoms%ntype) )
             ALLOCATE ( hybrid%basm2(atoms%jmtd,hybrid%maxindxm2,0:hybrid%maxlcutm2,atoms%ntype) )
             READ(iounit) hybrid%ngptm2,hybrid%pgptm2,hybrid%nindxm2
             READ(iounit) hybrid%basm2
184 185 186 187 188
          END IF

          CLOSE(iounit)

          RETURN
Daniel Wortmann's avatar
Daniel Wortmann committed
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
       END IF
    END IF

    hybrid%maxindx = 0
    DO itype = 1,atoms%ntype
       DO l = 0,atoms%lmax(itype)
          hybrid%maxindx = MAX(hybrid%maxindx,2+COUNT( atoms%llo(:atoms%nlo(itype),itype) .EQ. l))
       END DO
    END DO
    !       maxindx   = maxval( nlo   ) + 2



    ! initialize gridf for radial integration
    CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,gridf)

    !
    ! read in energy parameters from file eig
    ! to avoid meaningless energy parameters which occur in the case
    ! that the energy parameter is set to the atomic prinicipal 
    ! quantum number
    ! (el0 and ello0 are just the values in the enpara file)
    ! 


    DO jsp=1,DIMENSION%jspd
       CALL read_eig(eig_id,1,jsp,el=el(:,:,jsp),ello=ello(:,:,jsp))
    ENDDO

    ALLOCATE ( vr0(atoms%jmtd,atoms%ntype,DIMENSION%jspd) )

    vr0(:,:,:) = v%mt(:,0,:,:) 
   
    ! calculate radial basisfunctions u and u' with
    ! the spherical part of the potential vr0 and store them in
    ! bas1 = large component ,bas2 = small component

    ALLOCATE( f(atoms%jmtd,2,0:atoms%lmaxd), df(atoms%jmtd,2,0:atoms%lmaxd) )
    ALLOCATE( bas1(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd) )
    ALLOCATE( bas2(atoms%jmtd,hybrid%maxindx,0:atoms%lmaxd,atoms%ntype,DIMENSION%jspd) )

    DO itype=1,atoms%ntype
       ng = atoms%jri(itype)
       DO ispin=1,DIMENSION%jspd
233
          DO l=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
234 235
             CALL radfun(l,itype,ispin,el(l,itype,ispin),vr0(:,itype,ispin),atoms,f(:,:,l),df(:,:,l),&
                  &                  usdus,nodem,noded,wronk)
236 237 238 239 240 241 242 243 244
          END DO
          bas1(1:ng,1,0:atoms%lmaxd,itype,ispin) =  f(1:ng,1,0:atoms%lmaxd)
          bas2(1:ng,1,0:atoms%lmaxd,itype,ispin) =  f(1:ng,2,0:atoms%lmaxd)
          bas1(1:ng,2,0:atoms%lmaxd,itype,ispin) = df(1:ng,1,0:atoms%lmaxd)
          bas2(1:ng,2,0:atoms%lmaxd,itype,ispin) = df(1:ng,2,0:atoms%lmaxd)

          hybrid%nindx(:,itype) = 2
          ! generate radial functions for local orbitals
          IF (atoms%nlo(itype).GE.1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
245 246 247 248 249 250 251 252 253 254 255 256 257 258
             CALL radflo( atoms,itype,ispin,&
                  &                   ello(1,1,ispin),vr0(:,itype,ispin),&
                  &                   f,df,mpi,&
                  &                   usdus,&
                  &                   uuilon,duilon,ulouilopn,flo)

             DO ilo=1,atoms%nlo(itype)
                hybrid%nindx(atoms%llo(ilo,itype),itype) =&
                     &          hybrid%nindx(atoms%llo(ilo,itype),itype) + 1
                bas1(1:ng,hybrid%nindx(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),&
                     &                                    itype,ispin) = flo(1:ng,1,ilo)
                bas2(1:ng,hybrid%nindx(atoms%llo(ilo,itype),itype),atoms%llo(ilo,itype),&
                     &                                    itype,ispin) = flo(1:ng,2,ilo)
             END DO
259 260 261

          END IF

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

Daniel Wortmann's avatar
Daniel Wortmann committed
265 266 267 268 269
    DEALLOCATE(f,df)

    ! the radial functions are normalized
    DO ispin=1,DIMENSION%jspd
       DO itype=1,atoms%ntype
270
          DO l=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
271 272 273 274 275 276 277 278 279 280
             DO i=1,hybrid%nindx(l,itype)
                rdum = intgrf(bas1(:,i,l,itype,ispin)**2&
                     &                     +bas2(:,i,l,itype,ispin)**2,&
                     &                      atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)

                bas1(:atoms%jri(itype),i,l,itype,ispin)&
                     &         = bas1(:atoms%jri(itype),i,l,itype,ispin)/SQRT(rdum)
                bas2(:atoms%jri(itype),i,l,itype,ispin)&
                     &         = bas2(:atoms%jri(itype),i,l,itype,ispin)/SQRT(rdum)
             END DO
281
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
       END DO
    END DO

    IF( exx ) THEN
       ! generate core wave functions (-> core1/2(jmtd,nindxc,0:lmaxc,ntype,jspd) )   

       CALL core_init( DIMENSION,input,atoms, lmaxcd,maxindxc)


       ALLOCATE( nindxc(0:lmaxcd,atoms%ntype) )
       ALLOCATE( core1(atoms%jmtd,maxindxc,0:lmaxcd,atoms%ntype,DIMENSION%jspd) )      
       ALLOCATE( core2(atoms%jmtd,maxindxc,0:lmaxcd,atoms%ntype,DIMENSION%jspd) )
       ALLOCATE( eig_c(maxindxc,0:lmaxcd,atoms%ntype,DIMENSION%jspd)      )

       DO ispin = 1,input%jspins
          CALL corewf( atoms,ispin,input,DIMENSION,vr0, lmaxcd,maxindxc,mpi,&
               lmaxc,nindxc,core1(:,:,:,:,ispin),core2(:,:,:,:,ispin),eig_c(:,:,:,ispin))

       END DO
    END IF ! exx

    !
    ! - - - - - - SETUP OF THE MIXED BASIS IN THE IR - - - - - - -
    !


    !
    ! construct G-vectors with cutoff smaller than gcutm
    !
    IF( exx ) THEN
       gcutm = 2*input%rkmax
    ELSE
       gcutm = hybrid%gcutm1
    END IF
    ALLOCATE ( hybrid%ngptm(kpts%nkptf) )

    hybrid%ngptm = 0
    i     = 0
    n     =-1

Daniel Wortmann's avatar
Daniel Wortmann committed
322
    rdum1 = MAXVAL( (/ (SQRT(SUM(MATMUL(kpts%bkf(:,ikpt),cell%bmat)**2)),ikpt=1,kpts%nkptf ) /) )
Daniel Wortmann's avatar
Daniel Wortmann committed
323 324 325 326 327 328 329 330 331

    ! a first run for the determination of the dimensions of the fields
    ! gptm,pgptm

    DO
       n    = n + 1
       ldum = .FALSE.
       DO x = -n,n
          n1 = n-ABS(x)
332
          DO y = -n1,n1
Daniel Wortmann's avatar
Daniel Wortmann committed
333 334 335 336 337 338 339
             n2 = n1-ABS(y)
             DO z = -n2,n2,MAX(2*n2,1)
                g     = (/x,y,z/) 
                rdum  = SQRT(SUM(MATMUL(g,cell%bmat)**2))-rdum1
                IF(rdum.GT. gcutm) CYCLE
                ldum1 = .FALSE.
                DO ikpt = 1,kpts%nkptf
Daniel Wortmann's avatar
Daniel Wortmann committed
340
                   kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
341 342 343 344 345 346 347 348 349 350 351 352 353
                   rdum = SUM(MATMUL(kvec+g,cell%bmat)**2)

                   IF(rdum.LE.gcutm**2) THEN
                      IF(.NOT.ldum1) THEN
                         i                     = i + 1
                         ldum1                 = .TRUE.
                      END IF

                      hybrid%ngptm(ikpt) = hybrid%ngptm(ikpt) + 1
                      ldum        = .TRUE.
                   END IF
                END DO
             END DO
354
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384
       END DO
       IF(.NOT.ldum) EXIT
    END DO

    hybrid%gptmd   = i
    hybrid%maxgptm = MAXVAL(hybrid%ngptm)

    ALLOCATE ( hybrid%gptm(3,hybrid%gptmd) )
    ALLOCATE ( hybrid%pgptm(hybrid%maxgptm,kpts%nkptf) )

    hybrid%gptm  = 0
    hybrid%pgptm = 0
    hybrid%ngptm = 0

    i     = 0
    n     =-1

    !     
    ! Allocate and initialize arrays needed for G vector ordering
    !
    ALLOCATE ( unsrt_pgptm(hybrid%maxgptm,kpts%nkptf) )
    ALLOCATE ( length_kG(hybrid%maxgptm,kpts%nkptf) )
    length_kG   = 0
    unsrt_pgptm = 0

    DO
       n    = n + 1
       ldum = .FALSE.
       DO x = -n,n
          n1 = n-ABS(x)
385
          DO y = -n1,n1
Daniel Wortmann's avatar
Daniel Wortmann committed
386 387 388 389 390 391 392
             n2 = n1-ABS(y)
             DO z = -n2,n2,MAX(2*n2,1)
                g     = (/x,y,z/) 
                rdum  = SQRT(SUM(MATMUL(g,cell%bmat)**2))-rdum1
                IF(rdum.GT. gcutm) CYCLE
                ldum1 = .FALSE.
                DO ikpt = 1,kpts%nkptf
Daniel Wortmann's avatar
Daniel Wortmann committed
393
                   kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412
                   rdum = SUM(MATMUL(kvec+g,cell%bmat)**2)

                   IF(rdum.LE.(gcutm)**2) THEN
                      IF(.NOT.ldum1) THEN
                         i          = i + 1
                         hybrid%gptm(:,i)  = g
                         ldum1      = .TRUE.
                      END IF

                      hybrid%ngptm(ikpt)              = hybrid%ngptm(ikpt) + 1
                      ldum                     = .TRUE.

                      ! Position of the vector is saved as pointer
                      unsrt_pgptm(hybrid%ngptm(ikpt),ikpt) = i
                      ! Save length of vector k + G for array sorting
                      length_kG(hybrid%ngptm(ikpt),ikpt)   = rdum
                   END IF
                END DO
             END DO
413
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
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
       END DO
       IF(.NOT.ldum) EXIT
    END DO

    !
    ! Sort pointers in array, so that shortest |k+G| comes first
    !
    DO ikpt = 1,kpts%nkptf
       ALLOCATE( ptr(hybrid%ngptm(ikpt)) )
       ! Divide and conquer algorithm for arrays > 1000 entries
       divconq = MAX( 0, INT( 1.443*LOG( 0.001*hybrid%ngptm(ikpt) ) ) )
       ! create pointers which correspond to a sorted array
       CALL rorderpf(ptr, length_kG(1:hybrid%ngptm(ikpt),ikpt),hybrid%ngptm(ikpt), divconq )
       ! rearrange old pointers
       DO igpt = 1,hybrid%ngptm(ikpt)
          hybrid%pgptm(igpt,ikpt) = unsrt_pgptm(ptr(igpt),ikpt)
       END DO
       DEALLOCATE( ptr )
    END DO
    DEALLOCATE( unsrt_pgptm )
    DEALLOCATE( length_kG )

    !
    ! construct IR mixed basis set for the representation of the non local exchange elements
    ! with cutoff gcutm1
    !

    ! first run to determine dimension of pgptm1
    ALLOCATE( hybrid%ngptm1(kpts%nkptf) )
    hybrid%ngptm1 = 0
    DO igpt = 1,hybrid%gptmd
       g = hybrid%gptm(:,igpt)
       DO ikpt = 1,kpts%nkptf
Daniel Wortmann's avatar
Daniel Wortmann committed
447
          kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
          rdum = SUM(MATMUL(kvec+g,cell%bmat)**2)
          IF( rdum .LE. hybrid%gcutm1**2 ) THEN
             hybrid%ngptm1(ikpt) = hybrid%ngptm1(ikpt) + 1
          END IF
       END DO
    END DO
    hybrid%maxgptm1 = MAXVAL( hybrid%ngptm1 )

    ALLOCATE( hybrid%pgptm1(hybrid%maxgptm1,kpts%nkptf) )
    hybrid%pgptm1 = 0
    hybrid%ngptm1 = 0
    !     
    ! Allocate and initialize arrays needed for G vector ordering
    !
    ALLOCATE ( unsrt_pgptm(hybrid%maxgptm1,kpts%nkptf) )
    ALLOCATE ( length_kG(hybrid%maxgptm1,kpts%nkptf) )
    length_kG   = 0
    unsrt_pgptm = 0
    DO igpt = 1,hybrid%gptmd
       g = hybrid%gptm(:,igpt)
       DO ikpt = 1,kpts%nkptf
Daniel Wortmann's avatar
Daniel Wortmann committed
469
          kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
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
          rdum = SUM(MATMUL(kvec+g,cell%bmat)**2)
          IF( rdum .LE. hybrid%gcutm1**2 ) THEN
             hybrid%ngptm1(ikpt)                   = hybrid%ngptm1(ikpt) + 1
             unsrt_pgptm(hybrid%ngptm1(ikpt),ikpt) = igpt
             length_kG(hybrid%ngptm1(ikpt),ikpt)   = rdum
          END IF
       END DO
    END DO

    !
    ! Sort pointers in array, so that shortest |k+G| comes first
    !
    DO ikpt = 1,kpts%nkptf
       ALLOCATE( ptr(hybrid%ngptm1(ikpt)) )
       ! Divide and conquer algorithm for arrays > 1000 entries
       divconq = MAX( 0, INT( 1.443*LOG( 0.001*hybrid%ngptm1(ikpt) ) ) )
       ! create pointers which correspond to a sorted array
       CALL rorderpf(ptr, length_kG(1:hybrid%ngptm1(ikpt),ikpt),hybrid%ngptm1(ikpt), divconq )
       ! rearrange old pointers
       DO igpt = 1,hybrid%ngptm1(ikpt)
          hybrid%pgptm1(igpt,ikpt) = unsrt_pgptm(ptr(igpt),ikpt)
       END DO
       DEALLOCATE( ptr )
    END DO
    DEALLOCATE( unsrt_pgptm )
    DEALLOCATE( length_kG )

    IF( exx ) THEN
       !
       ! construct IR mixed basis set for the representation of the response function
       ! with cutoff gcutm2 
       !


       ALLOCATE( hybrid%ngptm2(kpts%nkptf) )
       hybrid%ngptm2 = 0

       DO igpt = 1,hybrid%gptmd
508 509
          g = hybrid%gptm(:,igpt)
          DO ikpt = 1,kpts%nkptf
Daniel Wortmann's avatar
Daniel Wortmann committed
510
             kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
511 512 513 514
             rdum = SUM(MATMUL(kvec+g,cell%bmat)**2)
             IF( rdum .LE. hybrid%gcutm2**2 ) THEN
                hybrid%ngptm2(ikpt) = hybrid%ngptm2(ikpt) + 1
             END IF
515
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
516 517 518 519 520 521 522 523 524 525 526 527 528 529
       END DO
       hybrid%maxgptm2 = MAXVAL( hybrid%ngptm2 )

       ALLOCATE( hybrid%pgptm2(hybrid%maxgptm2,kpts%nkptf) )
       hybrid%pgptm2 = 0
       hybrid%ngptm2 = 0
       !     
       ! Allocate and initialize arrays needed for G vector ordering
       !
       ALLOCATE ( unsrt_pgptm(hybrid%maxgptm2,kpts%nkptf) )
       ALLOCATE ( length_kG(hybrid%maxgptm2,kpts%nkptf) )
       length_kG   = 0
       unsrt_pgptm = 0
       DO igpt = 1,hybrid%gptmd
530 531
          g = hybrid%gptm(:,igpt)
          DO ikpt = 1,kpts%nkptf
Daniel Wortmann's avatar
Daniel Wortmann committed
532
             kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
533 534 535 536 537 538
             rdum = SUM(MATMUL(kvec+g,cell%bmat)**2)
             IF( rdum .LE. hybrid%gcutm2**2 ) THEN
                hybrid%ngptm2(ikpt)                   = hybrid%ngptm2(ikpt) + 1
                unsrt_pgptm(hybrid%ngptm2(ikpt),ikpt) = igpt
                length_kG(hybrid%ngptm2(ikpt),ikpt)   = rdum
             END IF
539
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
540 541 542 543 544 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
       END DO

       !
       ! Sort pointers in array, so that shortest |k+G| comes first
       !
       DO ikpt = 1,kpts%nkptf
          ALLOCATE( ptr(hybrid%ngptm2(ikpt)) )
          ! Divide and conquer algorithm for arrays > 1000 entries
          divconq = MAX( 0, INT( 1.443*LOG( 0.001*hybrid%ngptm2(ikpt) ) ) )
          ! create pointers which correspond to a sorted array
          CALL rorderpf(ptr, length_kG(1:hybrid%ngptm2(ikpt),ikpt),hybrid%ngptm2(ikpt), divconq )
          ! rearrange old pointers
          DO igpt = 1,hybrid%ngptm2(ikpt)
             hybrid%pgptm2(igpt,ikpt) = unsrt_pgptm(ptr(igpt),ikpt)
          END DO
          DEALLOCATE( ptr )
       END DO
       DEALLOCATE( unsrt_pgptm )
       DEALLOCATE( length_kG )

    END IF

    IF ( mpi%irank == 0 ) THEN
       WRITE(6,'(/A)')       'Mixed basis'
       WRITE(6,'(A,I5)')     'Number of unique G-vectors: ',hybrid%gptmd
       WRITE(6,*)
       WRITE(6,'(3x,A)') 'IR Plane-wave basis with cutoff of gcutm (hybrid%gcutm1/2*input%rkmax):'
       WRITE(6,'(5x,A,I5)')  'Maximal number of G-vectors:',hybrid%maxgptm
       WRITE(6,*)
       WRITE(6,*)
       WRITE(6,'(3x,A)') 'IR Plane-wave basis for non-local exchange potential:'
       WRITE(6,'(5x,A,I5)')  'Maximal number of G-vectors:',hybrid%maxgptm1
       WRITE(6,*)
       IF( exx ) THEN
          WRITE(6,'(3x,A)') 'IR Plane-wave basis for response function:'
575 576
          WRITE(6,'(5x,A,I5)')  'Maximal number of G-vectors:',hybrid%maxgptm2
          WRITE(6,*)
Daniel Wortmann's avatar
Daniel Wortmann committed
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 610 611 612 613 614 615 616 617
       END IF
    END IF

    !
    ! - - - - - - - - Set up MT product basis for the non-local exchange potential  - - - - - - - - - -
    !

    IF ( mpi%irank == 0 ) THEN
       WRITE(6,'(A)') 'MT product basis for non-local exchange potential:'
       WRITE(6,'(A)') 'Reduction due to overlap (quality of orthonormality, should be < 1.0E-06)'
    END IF

    hybrid%maxlcutm1 = MAXVAL( hybrid%lcutm1 )


    ALLOCATE ( hybrid%nindxm1(0:hybrid%maxlcutm1,atoms%ntype) )
    ALLOCATE ( seleco(hybrid%maxindx,0:atoms%lmaxd) , selecu(hybrid%maxindx,0:atoms%lmaxd) )
    ALLOCATE ( selecmat(hybrid%maxindx,0:atoms%lmaxd,hybrid%maxindx,0:atoms%lmaxd) )
    hybrid%nindxm1 = 0    !!! 01/12/10 jij%M.b.

    !
    ! determine maximal indices of (radial) mixed-basis functions (->nindxm1) 
    ! (will be reduced later-on due to overlap)
    !

    hybrid%maxindxp1  = 0
    DO itype = 1,atoms%ntype
       seleco = .FALSE.
       selecu = .FALSE.
       seleco(1,0:hybrid%select1(1,itype)) = .TRUE.
       selecu(1,0:hybrid%select1(3,itype)) = .TRUE.
       seleco(2,0:hybrid%select1(2,itype)) = .TRUE.
       selecu(2,0:hybrid%select1(4,itype)) = .TRUE.

       ! include local orbitals 
       IF(hybrid%maxindx.GE.3) THEN
          seleco(3:,:) = .TRUE. 
          selecu(3:,:) = .TRUE.
       END IF

       DO l=0,hybrid%lcutm1(itype) 
618 619 620
          n = 0
          M = 0
          IF( exx ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
621 622 623 624 625 626 627 628 629 630 631 632 633
             !
             ! core * valence products
             !
             DO l1 = 0,lmaxc(itype)
                DO l2 = 0,atoms%lmax(itype)
                   IF( l .LT. ABS(l1-l2) .OR. l .GT. l1+l2 ) CYCLE
                   DO n1 = 1,nindxc(l1,itype)             
                      DO n2 = 1,hybrid%nindx(l2,itype)
                         M = M + 1
                         IF( .NOT. seleco(n2,l2) ) CYCLE
                         n = n + 1
                      END DO
                   END DO
634
                END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
635 636
             END DO

637 638 639 640 641 642 643
          END IF

          !
          ! valence * valence
          !

          ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Daniel Wortmann's avatar
Daniel Wortmann committed
644 645 646
          selecmat = RESHAPE ( (/ ((((seleco(n1,l1).AND.selecu(n2,l2),&
               n1=1,hybrid%maxindx),l1=0,atoms%lmaxd), n2=1,hybrid%maxindx),l2=0,atoms%lmaxd) /) ,&
               (/hybrid%maxindx,atoms%lmaxd+1,hybrid%maxindx,atoms%lmaxd+1/) )
647 648

          DO l1=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692
             DO l2=0,atoms%lmax(itype)
                IF( l.GE.ABS(l1-l2) .AND. l.LE.l1+l2) THEN
                   DO n1=1,hybrid%nindx(l1,itype)
                      DO n2=1,hybrid%nindx(l2,itype)
                         M = M + 1
                         IF(selecmat(n1,l1,n2,l2)) THEN
                            n = n + 1
                            selecmat(n2,l2,n1,l1) = .FALSE. ! prevent double counting of products (a*b = b*a)
                         END IF
                      END DO
                   END DO
                END IF
             END DO
          END DO
          IF(n.EQ.0 .AND. mpi%irank==0) &
               WRITE(6,'(A)') 'mixedbasis: Warning!  No basis-function product of ' // lchar(l) //&
               '-angular momentum defined.'
          hybrid%maxindxp1        = MAX(hybrid%maxindxp1,M)
          hybrid%nindxm1(l,itype) = n*DIMENSION%jspd
       END DO
    END DO
    hybrid%maxindxm1 = MAXVAL(hybrid%nindxm1)

    ALLOCATE ( hybrid%basm1(atoms%jmtd,hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype) )
    hybrid%basm1 = 0

    ! Define product bases and reduce them according to overlap

    DO itype=1,atoms%ntype
       seleco = .FALSE.
       selecu = .FALSE.
       seleco(1,0:hybrid%select1(1,itype)) = .TRUE.
       selecu(1,0:hybrid%select1(3,itype)) = .TRUE.
       seleco(2,0:hybrid%select1(2,itype)) = .TRUE.
       selecu(2,0:hybrid%select1(4,itype)) = .TRUE. 
       ! include lo's
       IF(hybrid%maxindx.GE.3) THEN
          seleco(3:,:) = .TRUE. 
          selecu(3:,:) = .TRUE.
       END IF
       IF(atoms%ntype.GT.1 .AND. mpi%irank==0)&
            &    WRITE(6,'(6X,A,I3)') 'Atom type',itype
       ng = atoms%jri(itype)
       DO l=0,hybrid%lcutm1(itype)
693 694 695
          n = hybrid%nindxm1(l,itype)
          ! allow for zero product-basis functions for
          ! current l-quantum number
Daniel Wortmann's avatar
Daniel Wortmann committed
696 697 698
          IF(n.EQ.0) THEN
             IF ( mpi%irank==0 ) WRITE(6,'(6X,A,'':   0 ->   0'')') lchar(l)
             CYCLE 
699 700 701 702 703 704 705 706
          END IF

          ! set up the overlap matrix
          ALLOCATE (olap(n,n),eigv(n,n),work(3*n),eig(n),ihelp(n))
          ihelp    = 1 ! initialize to avoid a segfault
          i        = 0

          IF( exx ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733
             ! core*valence
             DO l1 = 0,lmaxc(itype)
                DO l2 = 0,atoms%lmax(itype)
                   IF( l .LT. ABS(l1-l2) .OR. l .GT. l1+l2 ) CYCLE
                   DO n1 = 1,nindxc(l1,itype)
                      DO n2 = 1,hybrid%nindx(l2,itype)
                         IF( .NOT. seleco(n2,l2) ) CYCLE
                         DO ispin = 1,DIMENSION%jspd  
                            i = i + 1
                            IF(i.GT.n) call judft_error('got too many product functions',hint='This is a BUG, please report',calledby='mixedbasis')
                            hybrid%basm1(:ng,i,l,itype)  &
                                 = (  core1(:ng,n1,l1,itype,ispin)&
                                 * bas1(:ng,n2,l2,itype,ispin) &
                                 + core2(:ng,n1,l1,itype,ispin)&
                                 * bas2(:ng,n2,l2,itype,ispin)  )/atoms%rmsh(:ng,itype)

                            !normalize basm1
                            rdum=SQRT(intgrf( hybrid%basm1(:,i,l,itype)**2,&
                                 atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) )

                            hybrid%basm1(:ng,i,l,itype)=hybrid%basm1(:ng,i,l,itype)/rdum

                         END DO  !ispin
                      END DO  !n2
                   END DO  !n1
                END DO  !l2
             END DO  !l1
734 735 736 737 738
          END IF !exx

          ! valence*valence

          ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Daniel Wortmann's avatar
Daniel Wortmann committed
739 740 741
          selecmat = RESHAPE ( (/ ((((seleco(n1,l1).AND.selecu(n2,l2),&
               n1=1,hybrid%maxindx),l1=0,atoms%lmaxd), n2=1,hybrid%maxindx),l2=0,atoms%lmaxd) /) ,  &
               (/hybrid%maxindx,atoms%lmaxd+1,hybrid%maxindx,atoms%lmaxd+1/) )
742 743

          DO l1=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
744 745
             DO l2=0,atoms%lmax(itype)
                IF( l .LT. ABS(l1-l2) .OR. l .GT. l1+l2 ) CYCLE
746

Daniel Wortmann's avatar
Daniel Wortmann committed
747 748
                DO n1=1,hybrid%nindx(l1,itype)
                   DO n2=1,hybrid%nindx(l2,itype)
749

Daniel Wortmann's avatar
Daniel Wortmann committed
750 751 752 753
                      IF(selecmat(n1,l1,n2,l2)) THEN
                         DO ispin=1,DIMENSION%jspd
                            i = i + 1
                            IF(i.GT.n) call judft_error('got too many product functions',hint='This is a BUG, please report',calledby='mixedbasis')
754

Daniel Wortmann's avatar
Daniel Wortmann committed
755 756 757 758 759
                            hybrid%basm1(:ng,i,l,itype)  &
                                 = (  bas1(:ng,n1,l1,itype,ispin)&
                                 * bas1(:ng,n2,l2,itype,ispin) &
                                 + bas2(:ng,n1,l1,itype,ispin)&
                                 * bas2(:ng,n2,l2,itype,ispin) )/atoms%rmsh(:ng,itype)
760

Daniel Wortmann's avatar
Daniel Wortmann committed
761 762 763
                            !normalize basm1
                            rdum=SQRT(intgrf(hybrid%basm1(:,i,l,itype)**2,&
                                 atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) )
764

Daniel Wortmann's avatar
Daniel Wortmann committed
765 766 767 768 769 770
                            hybrid%basm1(:ng,i,l,itype) = hybrid%basm1(:ng,i,l,itype)/rdum

                         END DO !ispin
                         ! prevent double counting of products (a*b = b*a)
                         selecmat(n2,l2,n1,l1) = .FALSE. 
                      END IF
771

Daniel Wortmann's avatar
Daniel Wortmann committed
772 773
                   END DO !n2
                END DO !n1
774 775


Daniel Wortmann's avatar
Daniel Wortmann committed
776
             END DO !l2
777 778
          END DO  !l1

Daniel Wortmann's avatar
Daniel Wortmann committed
779
          IF(i .NE. n) call judft_error('counting error for product functions',hint='This is a BUG, please report',calledby='mixedbasis')
780 781 782 783 784 785 786 787 788 789

          ! In order to get ride of the linear dependencies in the 
          ! radial functions basm1 belonging to fixed l and itype
          ! the overlap matrix is diagonalized and those eigenvectors
          ! with a eigenvalue greater then hybrid%tolerance1 are retained


          ! Calculate overlap
          olap = 0
          DO n2=1,n
Daniel Wortmann's avatar
Daniel Wortmann committed
790 791 792 793 794
             DO n1=1,n2
                olap(n1,n2) = intgrf(hybrid%basm1(:,n1,l,itype) *hybrid%basm1(:,n2,l,itype),&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
                olap(n2,n1) = olap(n1,n2)
             END DO
795 796 797 798 799 800 801 802
          END DO

          ! Diagonalize
          CALL diagonalize(eigv,eig,olap)

          ! Get rid of linear dependencies (eigenvalue <= hybrid%tolerance1)
          nn = 0
          DO i=1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
803 804 805 806
             IF(eig(i).GT.hybrid%tolerance1) THEN 
                nn        = nn + 1
                ihelp(nn) = i
             END IF
807 808 809 810 811 812
          END DO
          hybrid%nindxm1(l,itype) = nn
          eig              = eig(ihelp)
          eigv(:,:)        = eigv(:,ihelp)

          DO i=1,ng
Daniel Wortmann's avatar
Daniel Wortmann committed
813
             hybrid%basm1(i,1:nn,l,itype) = MATMUL(hybrid%basm1(i,1:n,l,itype), eigv(:,1:nn))/SQRT(eig(:nn))
814 815 816
          END DO

          ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
Daniel Wortmann's avatar
Daniel Wortmann committed
817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837
          IF(l.EQ.0) THEN

             ! Check if basm1 must be reallocated
             IF(nn+1.GT.SIZE(hybrid%basm1,2)) THEN
                ALLOCATE ( basmhlp(atoms%jmtd,nn+1,0:hybrid%maxlcutm1,atoms%ntype) )
                basmhlp(:,1:nn,:,:) = hybrid%basm1
                DEALLOCATE (hybrid%basm1)
                ALLOCATE ( hybrid%basm1(atoms%jmtd,nn+1,0:hybrid%maxlcutm1,atoms%ntype) )
                hybrid%basm1(:,1:nn,:,:) = basmhlp(:,1:nn,:,:) 
                DEALLOCATE ( basmhlp )
             END IF

             hybrid%basm1(:ng,nn+1,0,itype) = atoms%rmsh(:ng,itype) / SQRT(atoms%rmsh(ng,itype)**3/3)
             DO i = nn,1,-1
                DO j = i+1,nn+1
                   hybrid%basm1(:ng,i,0,itype) = hybrid%basm1(:ng,i,0,itype)&
                        - intgrf(hybrid%basm1(:ng,i,0,itype)*hybrid%basm1(:ng,j,0,itype),&
                        atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)&
                        * hybrid%basm1(:ng,j,0,itype)

                END DO
838
                hybrid%basm1(:ng,i,0,itype) = hybrid%basm1(:ng,i,0,itype)&
Daniel Wortmann's avatar
Daniel Wortmann committed
839 840 841 842 843 844 845
                     / SQRT(intgrf(hybrid%basm1(:ng,i,0,itype)**2,&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf))
             END DO
             nn              = nn + 1
             DEALLOCATE ( olap )
             ALLOCATE   ( olap(nn,nn) )
             hybrid%nindxm1(l,itype) = nn
846 847 848 849 850
          END IF

          ! Check orthonormality of product basis
          rdum = 0
          DO i=1,nn
Daniel Wortmann's avatar
Daniel Wortmann committed
851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873
             DO j=1,i
                rdum1 = intgrf(hybrid%basm1(:,i,l,itype)*hybrid%basm1(:,j,l,itype),&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)

                IF(i.EQ.j) THEN
                   rdum1 = ABS(1-rdum1)
                   rdum  = rdum + rdum1**2
                ELSE
                   rdum1 = ABS(rdum1)
                   rdum  = rdum + 2*rdum1**2
                END IF

                IF(rdum1.GT.1d-6) THEN
                   IF ( mpi%irank == 0 ) THEN
                      WRITE(6,'(A)') 'mixedbasis: Bad orthonormality of ' &
                           //lchar(l)//'-product basis. Increase tolerance.'
                      WRITE(6,'(12X,A,F9.6,A,2(I3.3,A))') 'Deviation of',&
                           rdum1,' occurred for (',i,',',j,')-overlap.'
                   END IF
                   CALL judft_error("Bad orthonormality of product basis",hint='Increase tolerance',calledby='mixedbasis')
                END IF

             END DO
874 875 876
          END DO

          IF ( mpi%irank == 0 ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
877
             WRITE(6,'(6X,A,I4,'' ->'',I4,''   ('',ES8.1,'' )'')') lchar(l)//':',n,nn,SQRT(rdum)/nn
878 879 880 881
          END IF

          DEALLOCATE(olap,eigv,work,eig,ihelp)

Daniel Wortmann's avatar
Daniel Wortmann committed
882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
       END DO !l
       IF ( mpi%irank == 0 ) WRITE(6,'(6X,A,I7)') 'Total:', SUM(hybrid%nindxm1(0:hybrid%lcutm1(itype),itype))
    END DO ! itype

    hybrid%maxindxm1 = MAXVAL(hybrid%nindxm1)

    ALLOCATE( basmhlp(atoms%jmtd,hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype) )
    basmhlp(1:atoms%jmtd,1:hybrid%maxindxm1,0:hybrid%maxlcutm1,1:atoms%ntype)&
         = hybrid%basm1(1:atoms%jmtd,1:hybrid%maxindxm1,0:hybrid%maxlcutm1,1:atoms%ntype)
    DEALLOCATE(hybrid%basm1)
    ALLOCATE( hybrid%basm1(atoms%jmtd,hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype) )
    hybrid%basm1=basmhlp

    DEALLOCATE(basmhlp,seleco,selecu,selecmat)


    !
    ! now we build linear combinations of the radial functions
    ! such that they possess no moment except one radial function in each l-channel
    !
    IF ( mpi%irank == 0 ) THEN
       WRITE(6,'(/,A,/,A)')'Build linear combinations of radial '//&
            'functions in each l-channel,',&
            'such that they possess no multipolmoment'//&
            ' except the last function:'

       WRITE(6,'(/,17x,A)') 'moment  (quality of orthonormality)'
    END IF
    DO itype = 1,atoms%ntype
       ng = atoms%jri(itype)

       IF(atoms%ntype.GT.1 .AND. mpi%irank==0) WRITE(6,'(6X,A,I3)') 'Atom type',itype

       DO l = 0,hybrid%lcutm1(itype)
916 917 918 919
          ! determine radial function with the largest moment
          ! this function is used to build the linear combinations
          rdum = 0
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
920 921 922 923 924 925
             rdum1 = intgrf(atoms%rmsh(:ng,itype)**(l+1)*hybrid%basm1(:ng,i,l,itype),&
                  atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
             IF( ABS(rdum1) .GT.rdum ) THEN
                n = i
                rdum = rdum1
             END IF
926
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
927

928 929 930 931
          ! rearrange order of radial functions such that the last function possesses the largest moment
          j = 0
          bashlp(:ng) = hybrid%basm1(:ng,n,l,itype)
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
932 933 934
             IF( i .EQ. n ) CYCLE
             j = j + 1
             hybrid%basm1(:ng,j,l,itype) = hybrid%basm1(:ng,i,l,itype)
935 936
          END DO
          hybrid%basm1(:ng,hybrid%nindxm1(l,itype),l,itype) = bashlp(:ng)
Daniel Wortmann's avatar
Daniel Wortmann committed
937 938 939 940 941

       END DO


       DO l = 0,hybrid%lcutm1(itype)
942 943
          IF ( mpi%irank == 0 ) WRITE(6,'(6X,A)') lchar(l)//':'

Daniel Wortmann's avatar
Daniel Wortmann committed
944 945 946
          IF(hybrid%nindxm1(l,itype).EQ.0) THEN
             IF( mpi%irank == 0 ) WRITE(6,'(6X,A,'':   0 ->    '')') lchar(l)
             CYCLE
947
          END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
948

949 950
          n = hybrid%nindxm1(l,itype)
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
951 952 953 954
             IF(i.EQ.n) CYCLE
             ! calculate moment of radial function i
             rdum1 = intgrf(atoms%rmsh(:ng,itype)**(l+1)*hybrid%basm1(:ng,i,l,itype),&
                  atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
955

Daniel Wortmann's avatar
Daniel Wortmann committed
956 957
             rdum  = intgrf(atoms%rmsh(:ng,itype)**(l+1)*hybrid%basm1(:ng,n,l,itype),&
                  atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
958

Daniel Wortmann's avatar
Daniel Wortmann committed
959
             bashlp(:ng) = hybrid%basm1(:ng,n,l,itype)
960

Daniel Wortmann's avatar
Daniel Wortmann committed
961 962
             IF( SQRT(rdum**2 + rdum1**2) .LE. 1E-06 .AND. mpi%irank == 0)&
                  WRITE(6,*) 'Warning: Norm is smaller thann 1E-06!'
963

Daniel Wortmann's avatar
Daniel Wortmann committed
964 965 966 967 968
             ! change function n such that n is orthogonal to i
             ! since the functions basm1 have been orthogonal on input
             ! the linear combination does not destroy the orthogonality to the residual functions
             hybrid%basm1(:ng,n,l,itype) = rdum /SQRT(rdum**2+rdum1**2) * bashlp(:ng)&
                  + rdum1/SQRT(rdum**2+rdum1**2) * hybrid%basm1(:ng,i,l,itype)
969

Daniel Wortmann's avatar
Daniel Wortmann committed
970 971 972
             ! combine basis function i and n so that they possess no momemt
             hybrid%basm1(:ng,i,l,itype) = rdum1/SQRT(rdum**2+rdum1**2) * bashlp(:ng)&
                  - rdum /SQRT(rdum**2+rdum1**2) * hybrid%basm1(:ng,i,l,itype)
973

Daniel Wortmann's avatar
Daniel Wortmann committed
974 975
             rdum1 = intgrf(atoms%rmsh(:ng,itype)**(l+1)*hybrid%basm1(:ng,i,l,itype),&
                  atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
976

Daniel Wortmann's avatar
Daniel Wortmann committed
977
             IF( rdum1 .GT. 1E-10 ) call judft_error('moment of radial function does not vanish',calledby='mixedbasis')
978

Daniel Wortmann's avatar
Daniel Wortmann committed
979
             IF ( mpi%irank == 0 ) WRITE(6,'(6x,I4,'' ->  '',ES8.1)') i,rdum1
980 981 982 983 984
          END DO

          ! test orthogonality
          rdum = 0
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
985 986 987 988 989 990 991 992 993
             DO j = 1,hybrid%nindxm1(l,itype) 
                rdum1 = intgrf(hybrid%basm1(:ng,i,l,itype)*hybrid%basm1(:ng,j,l,itype),&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
                IF( i .NE. j ) THEN
                   rdum = rdum + rdum1
                ELSE
                   rdum = rdum + ABS(1 - rdum1)
                END IF
             END DO
994 995
          END DO
          IF ( mpi%irank == 0 )&
Daniel Wortmann's avatar
Daniel Wortmann committed
996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008
               WRITE(6,'(6x,I4,'' ->'',f10.5,''   ('',ES8.1,'' )'')')  n,&
               intgrf(atoms%rmsh(:ng,itype)**(l+1)*hybrid%basm1(:ng,n,l,itype),atoms%jri,&
               atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf),rdum
       END DO
    END DO

    IF( exx ) THEN
       CALL judft_warn("EXX not implemented in this version of FLEUR")
       !
       ! - - - - - - Set up MT product basis for the representation of the response function  - - - - - -
       !

       IF ( mpi%irank == 0 ) THEN
1009 1010 1011
          WRITE(6,*)
          WRITE(6,'(A)')       'MT product basis for response function:'
          WRITE(6,'(A)')       'Reduction due to overlap (quality of'//&
Daniel Wortmann's avatar
Daniel Wortmann committed
1012 1013
               ' orthonormality, should be < 1.0E-06)'
       END IF
1014

Daniel Wortmann's avatar
Daniel Wortmann committed
1015
       hybrid%maxlcutm2  = MAXVAL(hybrid%lcutm2)
1016

Daniel Wortmann's avatar
Daniel Wortmann committed
1017 1018 1019
       !
       ! read in atomic exchange potential
       !
1020

Daniel Wortmann's avatar
Daniel Wortmann committed
1021 1022 1023 1024 1025
       !
       ! to get the name of each element type read in inp file
       !
       OPEN (123,file='inp',form='formatted',status='old')
       DO
1026
          READ(123,*) nchar
Daniel Wortmann's avatar
Daniel Wortmann committed
1027 1028
          IF(nchar .EQ. '**') EXIT
       END DO
1029

Daniel Wortmann's avatar
Daniel Wortmann committed
1030
       DO itype = 1,atoms%ntype
1031
          READ(123,*) noel(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
1032 1033
          IF( noel(itype)(2:2) .EQ. " " ) THEN
             WRITE(fname(itype),*) 'vx-',noel(itype)(1:1),'.tab'
1034
          ELSE
Daniel Wortmann's avatar
Daniel Wortmann committed
1035
             WRITE(fname(itype),*) 'vx-',noel(itype),'.tab'
1036 1037
          END IF
          DO
Daniel Wortmann's avatar
Daniel Wortmann committed
1038 1039
             READ(123,*) nchar
             IF(nchar .EQ. '**') EXIT
1040
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
1041 1042
       END DO
       CLOSE(123)
1043

Daniel Wortmann's avatar
Daniel Wortmann committed
1044
       CALL read_atompotential(atoms, (/(fname(i),i=1,atoms%ntype)/),potatom)
1045

Daniel Wortmann's avatar
Daniel Wortmann committed
1046 1047 1048
       ALLOCATE ( hybrid%nindxm2(0:hybrid%maxlcutm2,atoms%ntype) )
       ALLOCATE ( seleco(hybrid%maxindx,0:atoms%lmaxd) , selecu(hybrid%maxindx,0:atoms%lmaxd) )
       ALLOCATE ( selecmat(hybrid%maxindx,0:atoms%lmaxd,hybrid%maxindx,0:atoms%lmaxd) )
1049 1050


Daniel Wortmann's avatar
Daniel Wortmann committed
1051 1052
       ! determine maximal indices of (radial) mixed-basis functions (->nindxm) 
       ! (will be reduced later-on due to overlap)
1053

Daniel Wortmann's avatar
Daniel Wortmann committed
1054 1055 1056 1057 1058 1059 1060 1061 1062
       hybrid%nindxm2    = 0
       hybrid%maxindxp2  = 0
       DO itype=1,atoms%ntype
          seleco = .FALSE.
          selecu = .FALSE.
          seleco(1,0:hybrid%select2(1,itype)) = .TRUE. !0!1
          selecu(1,0:hybrid%select2(3,itype)) = .TRUE. !0!1
          seleco(2,0:hybrid%select2(2,itype)) = .TRUE. !1!1
          selecu(2,0:hybrid%select2(4,itype)) = .TRUE. !0!1
1063
          DO l=0,hybrid%lcutm2(itype) 
Daniel Wortmann's avatar
Daniel Wortmann committed
1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107
             n = 0
             M = 0
             !             !
             !             ! core * valence products
             !             !
             !             DO l1 = 0,lmaxc(itype)
             !               DO l2 = 0,lmax(itype)
             !                 IF( l .lt. abs(l1-l2) .or. l .gt. l1+l2 ) CYCLE
             !                 DO n1 = 1,nindxc(l1,itype)             
             !                   DO n2 = 1,nindx(l2,itype)
             !                     m = m + 1
             !                     IF( .not. seleco(n2,l2) ) CYCLE
             !                     n = n + 1
             !                   END DO
             !                 END DO
             !               END DO
             !             END DO

             ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
             selecmat = RESHAPE ( (/ ((((seleco(n1,l1).AND.selecu(n2,l2),&
                  n1=1,hybrid%maxindx),l1=0,atoms%lmaxd), n2=1,hybrid%maxindx),l2=0,atoms%lmaxd) /) ,&
                  (/hybrid%maxindx,atoms%lmaxd+1,hybrid%maxindx,atoms%lmaxd+1/) )

             DO l1=0,atoms%lmax(itype)
                DO l2=0,atoms%lmax(itype)
                   IF(l.GE.ABS(l1-l2).AND.l.LE.l1+l2) THEN
                      DO n1=1,hybrid%nindx(l1,itype)
                         DO n2=1,hybrid%nindx(l2,itype)
                            M = M + 1
                            IF(selecmat(n1,l1,n2,l2)) THEN
                               n = n + 1
                               selecmat(n2,l2,n1,l1) = .FALSE. ! prevent double counting of products (a*b = b*a)
                            END IF
                         END DO !n2
                      END DO !n1
                   END IF
                END DO !l2
             END DO !l1
             IF(n.EQ.0 .AND. mpi%irank==0) &
                  WRITE(6,'(A)') 'mixedbasis: Warning! No basis-function'//&
                  ' product of ' // lchar(l) //&
                  '-angular momentum defined.'
             hybrid%maxindxp2        = MAX(hybrid%maxindxp2,M)
             hybrid%nindxm2(l,itype) = n*DIMENSION%jspd
1108
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
1109
       END DO
1110

Daniel Wortmann's avatar
Daniel Wortmann committed
1111
       hybrid%maxindxm2 = MAXVAL(hybrid%nindxm2)
1112 1113


Daniel Wortmann's avatar
Daniel Wortmann committed
1114 1115
       ALLOCATE (hybrid%basm2(atoms%jmtd,hybrid%maxindxm2,0:hybrid%maxlcutm2,atoms%ntype) )
       hybrid%basm2 = 0
1116

Daniel Wortmann's avatar
Daniel Wortmann committed
1117 1118 1119 1120 1121 1122 1123 1124 1125 1126
       ! Define product bases and reduce them according to overlap

       DO itype=1,atoms%ntype
          seleco = .FALSE.
          selecu = .FALSE.
          seleco(1,0:hybrid%select2(1,itype)) = .TRUE. !0!1
          selecu(1,0:hybrid%select2(3,itype)) = .TRUE. !0!1
          seleco(2,0:hybrid%select2(2,itype)) = .TRUE. !1!1
          selecu(2,0:hybrid%select2(4,itype)) = .TRUE. !0!1
          IF(atoms%ntype.GT.1 .AND. mpi%irank==0) WRITE(6,'(6X,A,I3)') 'Atom type',itype
1127 1128
          ng = atoms%jri(itype)
          DO l=0,hybrid%lcutm2(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259
             n = hybrid%nindxm2(l,itype)
             ! allow for zero product-basis functions for
             ! current l-quantum number
             IF(n.EQ.0) THEN
                IF ( mpi%irank == 0 ) WRITE(6,'(6X,A,'':   0 ->   0'')') lchar(l)
                CYCLE 
             END IF

             ! set up the overlap matrix
             ALLOCATE (olap(n,n),eigv(n,n),work(3*n),eig(n),ihelp(n))
             ihelp    = 1 ! initialize to avoid a segfault
             i        = 0

             !             ! core*valence
             !             DO l1 = 0,lmaxc(itype)
             !               DO l2 = 0,lmax(itype)
             !                 IF( l .lt. abs(l1-l2) .or. l .gt. l1+l2 ) CYCLE
             !                 DO n1 = 1,nindxc(l1,itype)
             !                   DO n2 = 1,nindx(l2,itype)
             !                     IF( .not. seleco(n2,l2) ) CYCLE
             !                     DO ispin = 1,jspd  
             !                       i = i + 1
             !                       IF(i.gt.n) THEN
             !                          STOP 'mixedbasis:got too many product functions (bug).'
             !                       END IF
             !                       basm2(:ng,i,l,itype)  
             !      &              = ( core1(:ng,n1,l1,itype,ispin) * bas1(:ng,n2,l2,itype,ispin) 
             !      &                 +core2(:ng,n1,l1,itype,ispin) * bas2(:ng,n2,l2,itype,ispin) )/rmsh(:ng,itype)
             !                     
             !                       !normalize basm1
             !                       rdum=sqrt(intgrf( basm1(:,i,l,itype)**2,
             !      &                                  jri,jmtd,rmsh,dx,ntype,itype,gridf) ) 
             ! 
             !                       basm2(:ng,i,l,itype)=basm2(:ng,i,l,itype)/rdum
             ! 
             !                     END DO  !ispin
             !                   END DO  !n2
             !                 END DO  !n1
             !               END DO  !l2
             !             END DO  !l1



             ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
             selecmat = RESHAPE ( (/ ((((seleco(n1,l1).AND.selecu(n2,l2),&
                  n1=1,hybrid%maxindx),l1=0,atoms%lmaxd), n2=1,hybrid%maxindx),l2=0,atoms%lmaxd) /) ,  &
                  (/hybrid%maxindx,atoms%lmaxd+1,hybrid%maxindx,atoms%lmaxd+1/) )

             DO l1=0,atoms%lmax(itype)
                DO n1 = 1,hybrid%nindx(l1,itype)
                   DO l2=0,atoms%lmax(itype)
                      DO n2 = 1,hybrid%nindx(l2,itype)

                         IF(l.GE.ABS(l1-l2).AND.l.LE.l1+l2) THEN
                            IF(selecmat(n1,l1,n2,l2)) THEN
                               DO ispin=1,DIMENSION%jspd
                                  i = i + 1
                                  IF(i.GT.n) call judft_error('got too many product functions',hint='This is a BUG, please report',calledby='mixedbasis')

                                  hybrid%basm2(:ng,i,l,itype) = ( bas1(:ng,n1,l1,itype,ispin)&
                                       * bas1(:ng,n2,l2,itype,ispin) +   bas2(:ng,n1,l1,itype,ispin)&
                                       * bas2(:ng,n2,l2,itype,ispin) ) / atoms%rmsh(:ng,itype)

                                  !normalize basm2
                                  rdum=SQRT(intgrf(hybrid%basm2(:,i,l,itype)**2,&
                                       atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) )

                                  hybrid%basm2(:ng,i,l,itype)=hybrid%basm2(:ng,i,l,itype)/rdum

                               END DO !ispin
                               ! prevent double counting of products (a*b = b*a)
                               selecmat(n2,l2,n1,l1) = .FALSE. 
                            END IF
                         END IF

                      END DO  !n2
                   END DO  !l2
                END DO !n1
             END DO  !l1

             IF(i.NE.n) call judft_error('counting error for product functions',hint='This is a BUG, please report',calledby='mixedbasis')


             !Add constant function and atomic EXX potential to l=0 basis

             IF(l.EQ.0) THEN

                ! Check if basm2 must be reallocated
                IF(n+2.GT.SIZE(hybrid%basm2,2)) THEN
                   ALLOCATE ( basmhlp(atoms%jmtd,n,0:hybrid%maxlcutm2,atoms%ntype) )
                   basmhlp(:,1:n,:,:) = hybrid%basm2
                   DEALLOCATE (hybrid%basm2)
                   ALLOCATE ( hybrid%basm2(atoms%jmtd,n+2,0:hybrid%maxlcutm2,atoms%ntype) )
                   hybrid%basm2(:,1:n,:,:) = basmhlp(:,1:n,:,:) 
                   DEALLOCATE ( basmhlp )
                END IF

                ! add normalized atomic exx potential
                hybrid%basm2(:ng,n+2,0,itype)= potatom(:ng,itype)*atoms%rmsh(:ng,itype)

                rdum=intgrf(hybrid%basm2(:ng,n+2,0,itype)*hybrid%basm2(:ng,n+2,0,itype),&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)

                hybrid%basm2(:ng,n+2,0,itype) = hybrid%basm2(:ng,n+2,0,itype)/SQRT(rdum)


                ! add normalized constant function
                hybrid%basm2(:ng,n+1,0,itype) = atoms%rmsh(:ng,itype)&
                     / SQRT(atoms%rmsh(ng,itype)**3/3)

                ! orthogonalize the constant function with respect the atomic potential
                hybrid%basm2(:ng,n+1,0,itype) = hybrid%basm2(:ng,n+1,0,itype)&
                     - intgrf(hybrid%basm2(:ng,n+1,0,itype)* hybrid%basm2(:ng,n+2,0,itype),&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)&
                     * hybrid%basm2(:ng,n+2,0,itype)


                hybrid%basm2(:ng,n+1,0,itype) = hybrid%basm2(:ng,n+1,0,itype)&
                     / SQRT( intgrf(hybrid%basm2(:ng,n+1,0,itype)**2,&
                     atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) )

                n               = n + 2
                DEALLOCATE( olap )
                ALLOCATE  ( olap(n,n) )
                hybrid%nindxm2(l,itype)= n
             END IF


             ! Calculate overlap
             olap = 0
             DO n2=1,n
1260
                DO n1=1,n2
Daniel Wortmann's avatar
Daniel Wortmann committed
1261 1262 1263
                   olap(n1,n2) = intgrf(hybrid%basm2(:,n1,l,itype)*hybrid%basm2(:,n2,l,itype),&
                        atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
                   olap(n2,n1) = olap(n1,n2)
1264
                END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
1265
             END DO
1266

Daniel Wortmann's avatar
Daniel Wortmann committed
1267
             IF( l .EQ. 0 ) THEN
1268 1269


Daniel Wortmann's avatar
Daniel Wortmann committed
1270 1271
                ! orthogonalize all function ( 1:n-2 ) with respect to the constant and atomic function
                DO i = n-2,1,-1
1272

Daniel Wortmann's avatar
Daniel Wortmann committed
1273 1274 1275 1276 1277 1278
                   DO j = n-1,n
                      hybrid%basm2(:ng,i,0,itype) = hybrid%basm2(:ng,i,0,itype) - olap(j,i)*hybrid%basm2(:ng,j,0,itype)
                   END DO
                   hybrid%basm2(:ng,i,0,itype) = hybrid%basm2(:ng,i,0,itype) &
                        / SQRT( intgrf(hybrid%basm2(:ng,i,0,itype)**2, atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf) )
                END DO
1279

Daniel Wortmann's avatar
Daniel Wortmann committed
1280 1281 1282 1283 1284 1285 1286 1287 1288
                ! Calculate overlap
                olap = 0
                DO n2=1,n
                   DO n1=1,n2
                      olap(n1,n2) = intgrf(hybrid%basm2(:,n1,l,itype)*hybrid%basm2(:,n2,l,itype),&
                           atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)
                      olap(n2,n1) = olap(n1,n2)
                   END DO
                END DO
1289

Daniel Wortmann's avatar
Daniel Wortmann committed
1290
                n = n - 2
1291

Daniel Wortmann's avatar
Daniel Wortmann committed
1292
             END IF
1293

Daniel Wortmann's avatar
Daniel Wortmann committed
1294 1295 1296 1297
             ! In order to get ride of the linear dependencies in the
             ! radial functions basm2 belonging to fixed l and itype
             ! the overlap matrix is diagonalized and those eigenvectors
             ! with a eigenvalue greater then tolerance are retained
1298 1299


Daniel Wortmann's avatar
Daniel Wortmann committed
1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317
             ! Diagonalize

             CALL diagonalize(eigv,eig,olap(:n,:n))

             ! Get rid of linear dependencies (eigenvalue <= tolerance)

             IF( .FALSE. ) THEN
                IF( l .EQ. 0 ) rdum2 = hybrid%tolerance2*1E07
                ic = INT(rdum2/10.**(8-l))
                rdum2 = rdum2 - ic*10.**(8-l)
                WRITE(*,*) ic
                nn = 0
                DO i= n,1,-1
                   nn        = nn + 1
                   ihelp(nn) = i
                   WRITE(*,*) nn,i
                   IF( nn .EQ. ic ) EXIT
                END DO
1318

Daniel Wortmann's avatar
Daniel Wortmann committed
1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371
             ELSE
                nn = 0
                DO i=1,n !hybrid%nindxm2(l,itype)
                   IF(eig(i).GT.hybrid%tolerance2) THEN 
                      nn        = nn + 1
                      ihelp(nn) = i
                   END IF
                END DO
             END IF

             hybrid%nindxm2(l,itype) = nn
             eig              = eig(ihelp)
             eigv(:,:)        = eigv(:,ihelp)

             DO i=1,ng
                hybrid%basm2(i,1:nn,l,itype) = MATMUL(hybrid%basm2(i,1:n,l,itype), eigv(:,1:nn)) / SQRT(eig(:nn))
             END DO

             IF( l .EQ. 0 ) THEN
                hybrid%nindxm2(l,itype)         = nn + 2
                hybrid%basm2(:,nn+1,l,itype) = hybrid%basm2(:,n+1,l,itype)
                hybrid%basm2(:,nn+2,l,itype) = hybrid%basm2(:,n+2,l,itype)
                nn = nn + 2
             END IF


             ! Check orthonormality of product basis
             rdum = 0
             DO i=1,nn
                DO j=1,i
                   rdum1 = intgrf(hybrid%basm2(:,i,l,itype)*hybrid%basm2(:,j,l,itype),&
                        atoms%jri,atoms%jmtd,atoms%rmsh,atoms%dx,atoms%ntype,itype,gridf)

                   IF(i.EQ.j) THEN
                      rdum1 = ABS(1-rdum1)
                      rdum  = rdum + rdum1**2
                   ELSE
                      rdum1 = ABS(rdum1)
                      rdum  = rdum + 2*rdum1**2
                   ENDIF

                   IF(rdum1.GT.1d-6) THEN
                      IF ( mpi%irank == 0 ) THEN
                         WRITE(6,'(A)') 'mixedbasis: Bad orthonormality of ' &
                              //lchar(l)//'-product basis. Increase tolerance.'
                         WRITE(6,'(12X,A,F9.6,A,2(I3.3,A))') 'Deviation of',&
                              rdum1,' occurred for (',i,',',j,')-overlap.'
                      END IF
                      call judft_error('Bad orthonormality of product basis',hint='Increase tolerance',calledby='mixedbasis')
                   END IF

                END DO
             END DO
1372

Daniel Wortmann's avatar
Daniel Wortmann committed
1373 1374
             IF ( mpi%irank == 0 ) WRITE(6,'(6X,A,I4,'' ->'',I4,''   ('',ES8.1,'' )'')')&
                  lchar(l)//':',n,nn,SQRT(rdum)/nn
1375

Daniel Wortmann's avatar
Daniel Wortmann committed
1376
             DEALLOCATE(olap,eigv,work,eig,ihelp)
1377 1378

          END DO !l
Daniel Wortmann's avatar
Daniel Wortmann committed
1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397
          IF ( mpi%irank == 0 ) WRITE(6,'(6X,A,I7)') 'Total:',SUM(hybrid%nindxm2(0:hybrid%lcutm2(itype),itype))
       END DO ! itype

       hybrid%maxindxm2 = MAXVAL(hybrid%nindxm2)

       ALLOCATE( basmhlp(atoms%jmtd,hybrid%maxindxm2,0:hybrid%maxlcutm2,atoms%ntype) )
       basmhlp(1:atoms%jmtd,1:hybrid%maxindxm2,0:hybrid%maxlcutm2,1:atoms%ntype)&
            = hybrid%basm2(1:atoms%jmtd,1:hybrid%maxindxm2,0:hybrid%maxlcutm2,1:atoms%ntype)
       DEALLOCATE(hybrid%basm2)
       ALLOCATE( hybrid%basm2(atoms%jmtd,hybrid%maxindxm2,0:hybrid%maxlcutm2,atoms%ntype) )
       hybrid%basm2=basmhlp
       DEALLOCATE(basmhlp)
    END IF


    DO itype = 1,atoms%ntype
       IF ( ANY(hybrid%nindxm1(0:hybrid%lcutm1(itype),itype) == 0) ) call judft_error('any hybrid%nindxm1 eq 0',calledby='mixedbasis')
    END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414
    !count basis functions
    hybrid%nbasp   = 0
    DO itype=1,atoms%ntype
       DO i=1,atoms%neq(itype)
          DO l=0,hybrid%lcutm1(itype)
             DO M=-l,l
                DO j=1,hybrid%nindxm1(l,itype)
                   hybrid%nbasp = hybrid%nbasp + 1
                END DO
             END DO
          END DO
       END DO
    END DO
    hybrid%maxbasm1  = hybrid%nbasp  + hybrid%maxgptm
    
    !
    
Daniel Wortmann's avatar
Daniel Wortmann committed
1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434
    !       ic = 0
    !       DO itype = 1,ntype
    !         DO l = 0,lcutm2(itype)
    !           DO i = 1,nindxm2(l,itype)
    !             ic = ic + 1
    !             DO j = 1,jmtd
    !               WRITE(200+ic,'(2f15.10)') rmsh(j,itype),basm2(j,i,l,itype,1)
    !             END DO
    !           END DO
    !         END DO
    !       END DO
    !       
    !       STOP

    !
    ! Write all information to a file to enable restarting
    !
    IF ( l_restart .AND. mpi%irank == 0 ) THEN

       OPEN(UNIT=iounit,FILE=ioname,FORM='unformatted', STATUS='replace')
Daniel Wortmann's avatar
Daniel Wortmann committed
1435
 
Daniel Wortmann's avatar
Daniel Wortmann committed
1436 1437 1438 1439 1440 1441 1442 1443 1444
       WRITE(iounit) kpts%nkptf,hybrid%gptmd
       WRITE(iounit) hybrid%maxgptm,hybrid%maxindx
       WRITE(iounit) hybrid%ngptm,hybrid%gptm,hybrid%pgptm,hybrid%nindx

       WRITE(iounit) hybrid%maxgptm1,hybrid%maxlcutm1,hybrid%maxindxm1,hybrid%maxindxp1
       WRITE(iounit) hybrid%ngptm1,hybrid%pgptm1,hybrid%nindxm1
       WRITE(iounit) hybrid%basm1

       IF ( exx ) THEN
1445 1446 1447
          WRITE(iounit) hybrid%maxgptm2,hybrid%maxlcutm2,hybrid%maxindxm2,hybrid%maxindxp2
          WRITE(iounit) hybrid%ngptm2,hybrid%pgptm2,hybrid%nindxm2
          WRITE(iounit) hybrid%basm2
Daniel Wortmann's avatar
Daniel Wortmann committed
1448
       END IF
1449

Daniel Wortmann's avatar
Daniel Wortmann committed
1450
       CLOSE(iounit)
1451

Daniel Wortmann's avatar
Daniel Wortmann committed
1452
       l_restart = .FALSE.
1453

Daniel Wortmann's avatar
Daniel Wortmann committed
1454
    END IF
1455

Daniel Wortmann's avatar
Daniel Wortmann committed
1456
  END SUBROUTINE mixedbasis
1457 1458


Daniel Wortmann's avatar
Daniel Wortmann committed
1459
  SUBROUTINE read_atompotential(atoms,fname, potential)
1460 1461

    USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481
    IMPLICIT NONE
    TYPE(t_atoms),INTENT(IN)   :: atoms

    ! - scalars -
    CHARACTER(10),INTENT(IN) ::   fname(atoms%ntype)

    ! - arrays -
    REAL        ,INTENT(OUT)::   potential(atoms%jmtd,atoms%ntype)
    ! - local scalars -
    INTEGER                 ::   i,j,itype,ic,m     
    REAL                    ::   rdum,r ,n
    ! - local arrays -
    REAL    ,ALLOCATABLE    ::   v_atom (:),r_atom(:)

    potential =0
    DO itype = 1,atoms%ntype
       OPEN(331,file=fname(itype),form='formatted')

       ic = 0
       DO 
1482 1483
          READ(331,*) rdum
          ic = ic + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
1484 1485
          IF( rdum .GT. (atoms%rmt(itype)+0.2) ) EXIT
       END DO
1486

Daniel Wortmann's avatar
Daniel Wortmann committed
1487 1488
       ALLOCATE( v_atom(ic),r_atom(ic) )
       REWIND(331)
1489

Daniel Wortmann's avatar
Daniel Wortmann committed
1490
       DO i = 1,ic
1491
          READ(331,*) r_atom(i),v_atom(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
1492
       END DO
1493