mixedbasis.F90 34.5 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 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)   !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
36

Daniel Wortmann's avatar
Daniel Wortmann committed
37 38
MODULE m_mixedbasis

39
CONTAINS
40

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102
SUBROUTINE mixedbasis(atoms,kpts,DIMENSION,input,cell,sym,xcpot,hybrid,enpara,mpi,v,l_restart)

   USE m_judft
   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_inbuild),  INTENT(IN)    :: xcpot
   TYPE(t_mpi),            INTENT(IN)    :: mpi
   TYPE(t_dimension),      INTENT(IN)    :: dimension
   TYPE(t_hybrid),         INTENT(INOUT) :: hybrid
   TYPE(t_enpara),         INTENT(IN)    :: enpara
   TYPE(t_input),          INTENT(IN)    :: input
   TYPE(t_sym),            INTENT(IN)    :: sym
   TYPE(t_cell),           INTENT(IN)    :: cell
   TYPE(t_kpts),           INTENT(IN)    :: kpts
   TYPE(t_atoms),          INTENT(IN)    :: atoms
   TYPE(t_potden),         INTENT(IN)    :: v

   LOGICAL,                INTENT(INOUT) :: l_restart

   ! local type variables
   TYPE(t_usdus)                   ::  usdus

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

   LOGICAL                         ::  ldum,ldum1

   ! - 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)
   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                            ::  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
Daniel Wortmann's avatar
Daniel Wortmann committed
103
 
104 105
   LOGICAL, ALLOCATABLE            ::  selecmat(:,:,:,:)
   LOGICAL, ALLOCATABLE            ::  seleco(:,:),selecu(:,:)
106

107 108 109
   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' /)
Daniel Wortmann's avatar
Daniel Wortmann committed
110
    
111 112 113 114 115 116 117
   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
118

119
   IF (mpi%irank == 0) WRITE(6,'(//A,I2,A)') '### subroutine: mixedbasis ###'
120

121
   IF (xcpot%is_name("exx")) CALL judft_error("EXX is not implemented in this version",calledby='mixedbasis')
Daniel Wortmann's avatar
Daniel Wortmann committed
122
    
123 124 125 126 127 128 129 130
   ! 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
131
   
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
   CALL usdus%init(atoms,dimension%jspd)

   ! If restart is specified read file if it already exists. If not create it.
   IF (l_restart) THEN

      ! Test if file exists
      INQUIRE(FILE=ioname,EXIST=l_found)

      IF (l_found.AND..FALSE.) THEN !reading not working yet
         ! 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))
         ALLOCATE (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))
         ALLOCATE (hybrid%nindxm1(0:hybrid%maxlcutm1,atoms%ntype))
         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
160
       
161
         CLOSE(iounit)
Daniel Wortmann's avatar
Daniel Wortmann committed
162

163 164 165
         RETURN
      END IF
   END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
166

167 168 169 170 171 172 173
   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
Daniel Wortmann's avatar
Daniel Wortmann committed
174

175 176
   ! initialize gridf for radial integration
   CALL intgrf_init(atoms%ntype,atoms%jmtd,atoms%jri,atoms%dx,atoms%rmsh,gridf)
Daniel Wortmann's avatar
Daniel Wortmann committed
177

178
   ALLOCATE (vr0(atoms%jmtd,atoms%ntype,DIMENSION%jspd))
Daniel Wortmann's avatar
Daniel Wortmann committed
179

180
   vr0(:,:,:) = v%mt(:,0,:,:)
Daniel Wortmann's avatar
Daniel Wortmann committed
181
   
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 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 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 350 351 352
   ! 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
         DO l = 0, atoms%lmax(itype)
            CALL radfun(l,itype,ispin,enpara%el0(l,itype,ispin),vr0(:,itype,ispin),atoms,&
                        f(:,:,l),df(:,:,l),usdus,nodem,noded,wronk)
         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
            CALL radflo(atoms,itype,ispin,enpara%ello0(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
         END IF
      END DO
   END DO

   DEALLOCATE(f,df)

   ! the radial functions are normalized
   DO ispin = 1, DIMENSION%jspd
      DO itype = 1, atoms%ntype
         DO l = 0, atoms%lmax(itype)
            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
         END DO
      END DO
   END DO

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

   ! construct G-vectors with cutoff smaller than gcutm
   gcutm = hybrid%gcutm1
   ALLOCATE (hybrid%ngptm(kpts%nkptf))

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

   rdum1 = MAXVAL((/ (SQRT(SUM(MATMUL(kpts%bkf(:,ikpt),cell%bmat)**2)),ikpt=1,kpts%nkptf) /))

   ! 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)
         DO y = -n1, n1
            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
                  kvec = kpts%bkf(:,ikpt)
                  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
         END DO
      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)
         DO y = -n1, n1
            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
                  kvec = kpts%bkf(:,ikpt)
                  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
         END DO
      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
Daniel Wortmann's avatar
Daniel Wortmann committed
353

354 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 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
   ! 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
         kvec = kpts%bkf(:,ikpt)
         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
         kvec = kpts%bkf(:,ikpt)
         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
406

Daniel Wortmann's avatar
Daniel Wortmann committed
407
    
408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433
   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.
Daniel Wortmann's avatar
Daniel Wortmann committed
434

435 436
   ! determine maximal indices of (radial) mixed-basis functions (->nindxm1) 
   ! (will be reduced later-on due to overlap)
Daniel Wortmann's avatar
Daniel Wortmann committed
437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452
    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) 
453 454
          n = 0
          M = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
455
        
456 457 458 459 460
          !
          ! valence * valence
          !

          ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Daniel Wortmann's avatar
Daniel Wortmann committed
461 462 463
          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/) )
464 465

          DO l1=0,atoms%lmax(itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
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
             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)
510 511 512
          n = hybrid%nindxm1(l,itype)
          ! allow for zero product-basis functions for
          ! current l-quantum number
Daniel Wortmann's avatar
Daniel Wortmann committed
513 514 515
          IF(n.EQ.0) THEN
             IF ( mpi%irank==0 ) WRITE(6,'(6X,A,'':   0 ->   0'')') lchar(l)
             CYCLE 
516 517 518 519 520 521 522
          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
523
       
524 525 526
          ! valence*valence

          ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
Daniel Wortmann's avatar
Daniel Wortmann committed
527 528 529
          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/) )
530 531

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

Daniel Wortmann's avatar
Daniel Wortmann committed
535 536
                DO n1=1,hybrid%nindx(l1,itype)
                   DO n2=1,hybrid%nindx(l2,itype)
537

Daniel Wortmann's avatar
Daniel Wortmann committed
538 539 540 541
                      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')
542

Daniel Wortmann's avatar
Daniel Wortmann committed
543 544 545 546 547
                            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)
548

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

Daniel Wortmann's avatar
Daniel Wortmann committed
553 554 555 556 557 558
                            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
559

Daniel Wortmann's avatar
Daniel Wortmann committed
560 561
                   END DO !n2
                END DO !n1
562 563


Daniel Wortmann's avatar
Daniel Wortmann committed
564
             END DO !l2
565 566
          END DO  !l1

Daniel Wortmann's avatar
Daniel Wortmann committed
567
          IF(i .NE. n) call judft_error('counting error for product functions',hint='This is a BUG, please report',calledby='mixedbasis')
568 569 570 571 572 573 574 575 576 577

          ! 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
578 579 580 581 582
             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
583 584 585 586 587 588 589 590
          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
591 592 593 594
             IF(eig(i).GT.hybrid%tolerance1) THEN 
                nn        = nn + 1
                ihelp(nn) = i
             END IF
595 596 597 598 599 600
          END DO
          hybrid%nindxm1(l,itype) = nn
          eig              = eig(ihelp)
          eigv(:,:)        = eigv(:,ihelp)

          DO i=1,ng
Daniel Wortmann's avatar
Daniel Wortmann committed
601
             hybrid%basm1(i,1:nn,l,itype) = MATMUL(hybrid%basm1(i,1:n,l,itype), eigv(:,1:nn))/SQRT(eig(:nn))
602 603 604
          END DO

          ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
Daniel Wortmann's avatar
Daniel Wortmann committed
605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625
          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
626
                hybrid%basm1(:ng,i,0,itype) = hybrid%basm1(:ng,i,0,itype)&
Daniel Wortmann's avatar
Daniel Wortmann committed
627 628 629 630 631 632 633
                     / 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
634 635 636 637 638
          END IF

          ! Check orthonormality of product basis
          rdum = 0
          DO i=1,nn
Daniel Wortmann's avatar
Daniel Wortmann committed
639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
             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
662 663 664
          END DO

          IF ( mpi%irank == 0 ) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
665
             WRITE(6,'(6X,A,I4,'' ->'',I4,''   ('',ES8.1,'' )'')') lchar(l)//':',n,nn,SQRT(rdum)/nn
666 667 668 669
          END IF

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

Daniel Wortmann's avatar
Daniel Wortmann committed
670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703
       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)
704 705 706 707
          ! 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
708 709 710 711 712 713
             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
714
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
715

716 717 718 719
          ! 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
720 721 722
             IF( i .EQ. n ) CYCLE
             j = j + 1
             hybrid%basm1(:ng,j,l,itype) = hybrid%basm1(:ng,i,l,itype)
723 724
          END DO
          hybrid%basm1(:ng,hybrid%nindxm1(l,itype),l,itype) = bashlp(:ng)
Daniel Wortmann's avatar
Daniel Wortmann committed
725 726 727 728 729

       END DO


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

Daniel Wortmann's avatar
Daniel Wortmann committed
732 733 734
          IF(hybrid%nindxm1(l,itype).EQ.0) THEN
             IF( mpi%irank == 0 ) WRITE(6,'(6X,A,'':   0 ->    '')') lchar(l)
             CYCLE
735
          END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
736

737 738
          n = hybrid%nindxm1(l,itype)
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
739 740 741 742
             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)
743

Daniel Wortmann's avatar
Daniel Wortmann committed
744 745
             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)
746

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

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

Daniel Wortmann's avatar
Daniel Wortmann committed
752 753 754 755 756
             ! 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)
757

Daniel Wortmann's avatar
Daniel Wortmann committed
758 759 760
             ! 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)
761

Daniel Wortmann's avatar
Daniel Wortmann committed
762 763
             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)
764

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

Daniel Wortmann's avatar
Daniel Wortmann committed
767
             IF ( mpi%irank == 0 ) WRITE(6,'(6x,I4,'' ->  '',ES8.1)') i,rdum1
768 769 770 771 772
          END DO

          ! test orthogonality
          rdum = 0
          DO i = 1,hybrid%nindxm1(l,itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
773 774 775 776 777 778 779 780 781
             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
782 783
          END DO
          IF ( mpi%irank == 0 )&
Daniel Wortmann's avatar
Daniel Wortmann committed
784 785 786 787 788 789 790 791 792 793 794 795
               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
796 797 798 799 800 801 802 803 804 805 806 807 808 809
    !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
810 811 812 813 814
    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
815

Daniel Wortmann's avatar
Daniel Wortmann committed
816
    !
Daniel Wortmann's avatar
Daniel Wortmann committed
817 818 819 820 821 822
    !
    ! 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
823
 
Daniel Wortmann's avatar
Daniel Wortmann committed
824 825 826 827 828 829 830 831
       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
832
   
Daniel Wortmann's avatar
Daniel Wortmann committed
833
       CLOSE(iounit)
834

Daniel Wortmann's avatar
Daniel Wortmann committed
835
       l_restart = .FALSE.
836

Daniel Wortmann's avatar
Daniel Wortmann committed
837
    END IF
838

Daniel Wortmann's avatar
Daniel Wortmann committed
839
  END SUBROUTINE mixedbasis
840 841


Daniel Wortmann's avatar
Daniel Wortmann committed
842
  SUBROUTINE read_atompotential(atoms,fname, potential)
843 844

    USE m_types
Daniel Wortmann's avatar
Daniel Wortmann committed
845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864
    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 
865 866
          READ(331,*) rdum
          ic = ic + 1
Daniel Wortmann's avatar
Daniel Wortmann committed
867 868
          IF( rdum .GT. (atoms%rmt(itype)+0.2) ) EXIT
       END DO
869

Daniel Wortmann's avatar
Daniel Wortmann committed
870 871
       ALLOCATE( v_atom(ic),r_atom(ic) )
       REWIND(331)
872

Daniel Wortmann's avatar
Daniel Wortmann committed
873
       DO i = 1,ic
874
          READ(331,*) r_atom(i),v_atom(i)
Daniel Wortmann's avatar
Daniel Wortmann committed
875
       END DO
876

Daniel Wortmann's avatar
Daniel Wortmann committed
877 878 879
       ! transfer potential from grid used in the atom program
       ! to those used here
       DO i = 1,atoms%jri(itype)
880 881
          r = atoms%rmsh(i,itype)
          DO j = 1,ic-1
Daniel Wortmann's avatar
Daniel Wortmann committed
882 883 884 885 886 887 888
             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
889
          END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
890
       END DO
891

Daniel Wortmann's avatar
Daniel Wortmann committed
892 893 894
       DEALLOCATE( v_atom,r_atom )
       CLOSE(331) 
    END DO
895

Daniel Wortmann's avatar
Daniel Wortmann committed
896
  END SUBROUTINE read_atompotential
897

Daniel Wortmann's avatar
Daniel Wortmann committed
898
END MODULE m_mixedbasis
899