mixedbasis.F90 36 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
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_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
54
    TYPE(t_kpts),INTENT(IN)      :: kpts
Daniel Wortmann's avatar
Daniel Wortmann committed
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
    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
Daniel Wortmann's avatar
Daniel Wortmann committed
70
    INTEGER                         ::  m,nk,ok
Daniel Wortmann's avatar
Daniel Wortmann committed
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
    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
108

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

Daniel Wortmann's avatar
Daniel Wortmann committed
113 114 115 116 117 118 119 120 121 122 123 124
    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
125

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


129
    exx = xcpot%is_name("exx")
Daniel Wortmann's avatar
Daniel Wortmann committed
130 131
    if (exx) call judft_error("EXX is not implemented in this version",calledby='mixedbasis.F90')
    
Daniel Wortmann's avatar
Daniel Wortmann committed
132 133 134 135 136 137 138 139
    ! 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%nindxm1))  DEALLOCATE(hybrid%nindxm1)
    IF(ALLOCATED(hybrid%pgptm))    DEALLOCATE(hybrid%pgptm)
    IF(ALLOCATED(hybrid%pgptm1))   DEALLOCATE(hybrid%pgptm1)
    IF(ALLOCATED(hybrid%gptm))     DEALLOCATE(hybrid%gptm)
    IF(ALLOCATED(hybrid%basm1) )   DEALLOCATE(hybrid%basm1)
Daniel Wortmann's avatar
Daniel Wortmann committed
140
   
Daniel Wortmann's avatar
Daniel Wortmann committed
141 142
    call usdus%init(atoms,dimension%jspd)

Daniel Wortmann's avatar
Daniel Wortmann committed
143 144 145
    ! If restart is specified read file if it already exists
    ! create it otherwise
    IF ( l_restart ) THEN
146

Daniel Wortmann's avatar
Daniel Wortmann committed
147 148
       ! Test if file exists
       INQUIRE(FILE=ioname,EXIST=l_found)
149

Daniel Wortmann's avatar
Daniel Wortmann committed
150
       IF ( l_found .and..false.) THEN !reading not working yet
151 152 153 154 155 156 157 158 159 160 161 162 163 164 165

          ! 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
166
               &               hybrid%nindxm1(0:hybrid%maxlcutm1,atoms%ntype) )
167 168 169 170
          ALLOCATE ( hybrid%basm1(atoms%jmtd,hybrid%maxindxm1,0:hybrid%maxlcutm1,atoms%ntype) )
          READ(iounit) hybrid%ngptm1,hybrid%pgptm1,hybrid%nindxm1
          READ(iounit) hybrid%basm1

Daniel Wortmann's avatar
Daniel Wortmann committed
171
       
172 173 174
          CLOSE(iounit)

          RETURN
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
       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
219
          DO l=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
220 221
             CALL radfun(l,itype,ispin,el(l,itype,ispin),vr0(:,itype,ispin),atoms,f(:,:,l),df(:,:,l),&
                  &                  usdus,nodem,noded,wronk)
222 223 224 225 226 227 228 229 230
          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
231 232 233 234 235 236 237 238 239 240 241 242 243 244
             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
245 246 247

          END IF

Daniel Wortmann's avatar
Daniel Wortmann committed
248 249
       END DO
    END DO
250

Daniel Wortmann's avatar
Daniel Wortmann committed
251 252 253 254 255
    DEALLOCATE(f,df)

    ! the radial functions are normalized
    DO ispin=1,DIMENSION%jspd
       DO itype=1,atoms%ntype
256
          DO l=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
257 258 259 260 261 262 263 264 265 266
             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
267
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
268 269 270
       END DO
    END DO

Daniel Wortmann's avatar
Daniel Wortmann committed
271
  
Daniel Wortmann's avatar
Daniel Wortmann committed
272 273 274 275 276 277 278 279
    !
    ! - - - - - - SETUP OF THE MIXED BASIS IN THE IR - - - - - - -
    !


    !
    ! construct G-vectors with cutoff smaller than gcutm
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
280
    gcutm = hybrid%gcutm1
Daniel Wortmann's avatar
Daniel Wortmann committed
281 282 283 284 285 286
    ALLOCATE ( hybrid%ngptm(kpts%nkptf) )

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

Daniel Wortmann's avatar
Daniel Wortmann committed
287
    rdum1 = MAXVAL( (/ (SQRT(SUM(MATMUL(kpts%bkf(:,ikpt),cell%bmat)**2)),ikpt=1,kpts%nkptf ) /) )
Daniel Wortmann's avatar
Daniel Wortmann committed
288 289 290 291 292 293 294 295 296

    ! 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)
297
          DO y = -n1,n1
Daniel Wortmann's avatar
Daniel Wortmann committed
298 299 300 301 302 303 304
             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
305
                   kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
306 307 308 309 310 311 312 313 314 315 316 317 318
                   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
319
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349
       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)
350
          DO y = -n1,n1
Daniel Wortmann's avatar
Daniel Wortmann committed
351 352 353 354 355 356 357
             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
358
                   kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377
                   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
378
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
       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
412
          kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
          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
434
          kvec = kpts%bkf(:,ikpt)
Daniel Wortmann's avatar
Daniel Wortmann committed
435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461
          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 )

Daniel Wortmann's avatar
Daniel Wortmann committed
462
    
Daniel Wortmann's avatar
Daniel Wortmann committed
463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
    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,*)
    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) 
514 515
          n = 0
          M = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
516
        
517 518 519 520 521
          !
          ! valence * valence
          !

          ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Daniel Wortmann's avatar
Daniel Wortmann committed
522 523 524
          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/) )
525 526

          DO l1=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
527 528 529 530 531 532 533 534 535 536 537 538 539 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
             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)
571 572 573
          n = hybrid%nindxm1(l,itype)
          ! allow for zero product-basis functions for
          ! current l-quantum number
Daniel Wortmann's avatar
Daniel Wortmann committed
574 575 576
          IF(n.EQ.0) THEN
             IF ( mpi%irank==0 ) WRITE(6,'(6X,A,'':   0 ->   0'')') lchar(l)
             CYCLE 
577 578 579 580 581 582 583
          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

Daniel Wortmann's avatar
Daniel Wortmann committed
584
       
585 586 587
          ! valence*valence

          ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Daniel Wortmann's avatar
Daniel Wortmann committed
588 589 590
          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/) )
591 592

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

Daniel Wortmann's avatar
Daniel Wortmann committed
596 597
                DO n1=1,hybrid%nindx(l1,itype)
                   DO n2=1,hybrid%nindx(l2,itype)
598

Daniel Wortmann's avatar
Daniel Wortmann committed
599 600 601 602
                      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')
603

Daniel Wortmann's avatar
Daniel Wortmann committed
604 605 606 607 608
                            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)
609

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

Daniel Wortmann's avatar
Daniel Wortmann committed
614 615 616 617 618 619
                            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
620

Daniel Wortmann's avatar
Daniel Wortmann committed
621 622
                   END DO !n2
                END DO !n1
623 624


Daniel Wortmann's avatar
Daniel Wortmann committed
625
             END DO !l2
626 627
          END DO  !l1

Daniel Wortmann's avatar
Daniel Wortmann committed
628
          IF(i .NE. n) call judft_error('counting error for product functions',hint='This is a BUG, please report',calledby='mixedbasis')
629 630 631 632 633 634 635 636 637 638

          ! 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
639 640 641 642 643
             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
644 645 646 647 648 649 650 651
          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
652 653 654 655
             IF(eig(i).GT.hybrid%tolerance1) THEN 
                nn        = nn + 1
                ihelp(nn) = i
             END IF
656 657 658 659 660 661
          END DO
          hybrid%nindxm1(l,itype) = nn
          eig              = eig(ihelp)
          eigv(:,:)        = eigv(:,ihelp)

          DO i=1,ng
Daniel Wortmann's avatar
Daniel Wortmann committed
662
             hybrid%basm1(i,1:nn,l,itype) = MATMUL(hybrid%basm1(i,1:n,l,itype), eigv(:,1:nn))/SQRT(eig(:nn))
663 664 665
          END DO

          ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
Daniel Wortmann's avatar
Daniel Wortmann committed
666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686
          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
687
                hybrid%basm1(:ng,i,0,itype) = hybrid%basm1(:ng,i,0,itype)&
Daniel Wortmann's avatar
Daniel Wortmann committed
688 689 690 691 692 693 694
                     / 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
695 696 697 698 699
          END IF

          ! Check orthonormality of product basis
          rdum = 0
          DO i=1,nn
Daniel Wortmann's avatar
Daniel Wortmann committed
700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722
             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
723 724 725
          END DO

          IF ( mpi%irank == 0 ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
726
             WRITE(6,'(6X,A,I4,'' ->'',I4,''   ('',ES8.1,'' )'')') lchar(l)//':',n,nn,SQRT(rdum)/nn
727 728 729 730
          END IF

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

Daniel Wortmann's avatar
Daniel Wortmann committed
731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764
       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)
765 766 767 768
          ! 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
769 770 771 772 773 774
             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
775
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
776

777 778 779 780
          ! 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
781 782 783
             IF( i .EQ. n ) CYCLE
             j = j + 1
             hybrid%basm1(:ng,j,l,itype) = hybrid%basm1(:ng,i,l,itype)
784 785
          END DO
          hybrid%basm1(:ng,hybrid%nindxm1(l,itype),l,itype) = bashlp(:ng)
Daniel Wortmann's avatar
Daniel Wortmann committed
786 787 788 789 790

       END DO


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

Daniel Wortmann's avatar
Daniel Wortmann committed
793 794 795
          IF(hybrid%nindxm1(l,itype).EQ.0) THEN
             IF( mpi%irank == 0 ) WRITE(6,'(6X,A,'':   0 ->    '')') lchar(l)
             CYCLE
796
          END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
797

798 799
          n = hybrid%nindxm1(l,itype)
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
800 801 802 803
             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)
804

Daniel Wortmann's avatar
Daniel Wortmann committed
805 806
             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)
807

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

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

Daniel Wortmann's avatar
Daniel Wortmann committed
813 814 815 816 817
             ! 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)
818

Daniel Wortmann's avatar
Daniel Wortmann committed
819 820 821
             ! 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)
822

Daniel Wortmann's avatar
Daniel Wortmann committed
823 824
             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)
825

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

Daniel Wortmann's avatar
Daniel Wortmann committed
828
             IF ( mpi%irank == 0 ) WRITE(6,'(6x,I4,'' ->  '',ES8.1)') i,rdum1
829 830 831 832 833
          END DO

          ! test orthogonality
          rdum = 0
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
834 835 836 837 838 839 840 841 842
             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
843 844
          END DO
          IF ( mpi%irank == 0 )&
Daniel Wortmann's avatar
Daniel Wortmann committed
845 846 847 848 849 850 851 852 853 854 855 856
               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



    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
857 858 859 860 861 862 863 864 865 866 867 868 869 870
    !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
871 872 873 874 875
    DO nk = 1,kpts%nkptf
       hybrid%nbasm(nk) = hybrid%nbasp + hybrid%ngptm(nk)
    END DO

     hybrid%maxlmindx = MAXVAL((/ ( SUM( (/ (hybrid%nindx(l,itype)*(2*l+1), l=0,atoms%lmax(itype)) /) ),itype=1,atoms%ntype) /) )
Daniel Wortmann's avatar
Daniel Wortmann committed
876

Daniel Wortmann's avatar
Daniel Wortmann committed
877
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
878 879 880 881 882 883
    !
    ! 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
884
 
Daniel Wortmann's avatar
Daniel Wortmann committed
885 886 887 888 889 890 891 892
       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

Daniel Wortmann's avatar
Daniel Wortmann committed
893
   
Daniel Wortmann's avatar
Daniel Wortmann committed
894
       CLOSE(iounit)
895

Daniel Wortmann's avatar
Daniel Wortmann committed
896
       l_restart = .FALSE.
897

Daniel Wortmann's avatar
Daniel Wortmann committed
898
    END IF
899

Daniel Wortmann's avatar
Daniel Wortmann committed
900
  END SUBROUTINE mixedbasis
901 902


Daniel Wortmann's avatar
Daniel Wortmann committed
903
  SUBROUTINE read_atompotential(atoms,fname, potential)
904 905

    USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925
    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 
926 927
          READ(331,*) rdum
          ic = ic + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
928 929
          IF( rdum .GT. (atoms%rmt(itype)+0.2) ) EXIT
       END DO
930

Daniel Wortmann's avatar
Daniel Wortmann committed
931 932
       ALLOCATE( v_atom(ic),r_atom(ic) )
       REWIND(331)
933

Daniel Wortmann's avatar
Daniel Wortmann committed
934
       DO i = 1,ic
935
          READ(331,*) r_atom(i),v_atom(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
936
       END DO
937

Daniel Wortmann's avatar
Daniel Wortmann committed
938 939 940
       ! transfer potential from grid used in the atom program
       ! to those used here
       DO i = 1,atoms%jri(itype)
941 942
          r = atoms%rmsh(i,itype)
          DO j = 1,ic-1
Daniel Wortmann's avatar
Daniel Wortmann committed
943 944 945 946 947 948 949
             IF( r_atom(j) .GT. r  ) THEN
                M = (v_atom(j)-v_atom(j+1))/(r_atom(j)-r_atom(j+1))
                n = v_atom(j) - M*r_atom(j)
                potential(i,itype) =n + M*r
                !               WRITE(333,'(4f25.10)') n + m*r_atom(j),v_atom(j),n + m*r_atom(j+1),v_atom(j+1)
                EXIT
             END IF
950
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
951
       END DO
952

Daniel Wortmann's avatar
Daniel Wortmann committed
953 954 955
       DEALLOCATE( v_atom,r_atom )
       CLOSE(331) 
    END DO
956

Daniel Wortmann's avatar
Daniel Wortmann committed
957
  END SUBROUTINE read_atompotential
958

Daniel Wortmann's avatar
Daniel Wortmann committed
959
END MODULE m_mixedbasis
960