LCOV - code coverage report
Current view: top level - hybrid - mixedbasis.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 408 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 2 0.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       8             : ! This subroutine generates the mixed basis set used to evaluate the  !
       9             : ! exchange term in HF/hybrid functional calculations or EXX           !
      10             : ! calculations. In the latter case a second mixed basis set is setup  !
      11             : ! for the OEP integral equation.                                      !
      12             : ! In all cases the mixed basis consists of IR plane waves             !
      13             : !                                                                     !
      14             : ! IR:                                                                 !
      15             : !    M_{\vec{k},\vec{G}} = 1/\sqrt{V} \exp{i(\vec{k}+\vec{G})}        !
      16             : !                                                                     !
      17             : ! which are zero in the MT spheres and radial functions times         !
      18             : ! spherical harmonics in the MT spheres                               !
      19             : !                                                                     !
      20             : ! MT:                                                                 !
      21             : !     a                a                                              !
      22             : !    M              = U   Y                                           !
      23             : !     PL               PL  LM                                         !
      24             : !                                                                     !
      25             : !            where     a    a  a                                      !
      26             : !                     U  = u  u                                       !
      27             : !                      PL   pl  p'l'                                  !
      28             : !                                                                     !
      29             : !               and    L \in {|l-l'|,...,l+l'}                        !
      30             : !                                                                     !
      31             : !               and    P counts the different combinations of         !
      32             : !                      pl,p'l' which contribute to L                  !
      33             : !                                                                     !
      34             : !                                               M.Betzinger (09/07)   !
      35             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      36             : 
      37             : MODULE m_mixedbasis
      38             : 
      39             : CONTAINS
      40             : 
      41           0 :    SUBROUTINE mixedbasis(atoms, kpts, DIMENSION, input, cell, sym, xcpot, hybrid, enpara, mpi, v, l_restart)
      42             : 
      43             :       USE m_judft
      44             :       USE m_radfun, ONLY: radfun
      45             :       USE m_radflo, ONLY: radflo
      46             :       USE m_loddop, ONLY: loddop
      47             :       USE m_util, ONLY: intgrf_init, intgrf, rorderpf
      48             :       USE m_read_core
      49             :       USE m_wrapper
      50             :       USE m_eig66_io
      51             :       USE m_types
      52             : 
      53             :       IMPLICIT NONE
      54             : 
      55             :       TYPE(t_xcpot_inbuild), INTENT(IN)    :: xcpot
      56             :       TYPE(t_mpi), INTENT(IN)    :: mpi
      57             :       TYPE(t_dimension), INTENT(IN)    :: dimension
      58             :       TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      59             :       TYPE(t_enpara), INTENT(IN)    :: enpara
      60             :       TYPE(t_input), INTENT(IN)    :: input
      61             :       TYPE(t_sym), INTENT(IN)    :: sym
      62             :       TYPE(t_cell), INTENT(IN)    :: cell
      63             :       TYPE(t_kpts), INTENT(IN)    :: kpts
      64             :       TYPE(t_atoms), INTENT(IN)    :: atoms
      65             :       TYPE(t_potden), INTENT(IN)    :: v
      66             : 
      67             :       LOGICAL, INTENT(INOUT) :: l_restart
      68             : 
      69             :       ! local type variables
      70           0 :       TYPE(t_usdus)                   ::  usdus
      71             : 
      72             :       ! local scalars
      73             :       INTEGER                         ::  ilo, ikpt, ispin, itype, l1, l2, l, n, igpt, n1, n2, nn, i, j, ic, ng
      74             :       INTEGER                         ::  jsp, nodem, noded, m, nk, ok, x, y, z, maxindxc, lmaxcd
      75             :       INTEGER                         ::  divconq ! use Divide & Conquer algorithm for array sorting (>0: yes, =0: no)
      76             :       REAL                            ::  gcutm, wronk, rdum, rdum1, rdum2
      77             : 
      78             :       LOGICAL                         ::  ldum, ldum1
      79             : 
      80             :       ! - local arrays -
      81             :       INTEGER                         ::  g(3)
      82             :       INTEGER                         ::  lmaxc(atoms%ntype)
      83             :       INTEGER, ALLOCATABLE             ::  nindxc(:, :)
      84           0 :       INTEGER, ALLOCATABLE             ::  ihelp(:)
      85           0 :       INTEGER, ALLOCATABLE             ::  ptr(:)           ! pointer for array sorting
      86           0 :       INTEGER, ALLOCATABLE             ::  unsrt_pgptm(:, :) ! unsorted pointers to g vectors
      87             : 
      88             :       REAL                            ::  kvec(3)
      89           0 :       REAL                            ::  flo(atoms%jmtd, 2, atoms%nlod)
      90           0 :       REAL                            ::  uuilon(atoms%nlod, atoms%ntype), duilon(atoms%nlod, atoms%ntype)
      91           0 :       REAL                            ::  ulouilopn(atoms%nlod, atoms%nlod, atoms%ntype)
      92             :       REAL                            ::  potatom(atoms%jmtd, atoms%ntype)
      93           0 :       REAL                            ::  bashlp(atoms%jmtd)
      94             : 
      95           0 :       REAL, ALLOCATABLE               ::  f(:, :, :), df(:, :, :)
      96           0 :       REAL, ALLOCATABLE               ::  olap(:, :), work(:), eig(:), eigv(:, :)
      97           0 :       REAL, ALLOCATABLE               ::  bas1(:, :, :, :, :), bas2(:, :, :, :, :)
      98           0 :       REAL, ALLOCATABLE               ::  basmhlp(:, :, :, :)
      99           0 :       REAL, ALLOCATABLE               ::  gridf(:, :), vr0(:, :, :)
     100             :       REAL, ALLOCATABLE               ::  core1(:, :, :, :, :), core2(:, :, :, :, :)
     101             :       REAL, ALLOCATABLE               ::  eig_c(:, :, :, :)
     102           0 :       REAL, ALLOCATABLE               ::  length_kg(:, :) ! length of the vectors k + G
     103             : 
     104           0 :       LOGICAL, ALLOCATABLE            ::  selecmat(:, :, :, :)
     105           0 :       LOGICAL, ALLOCATABLE            ::  seleco(:, :), selecu(:, :)
     106             : 
     107             :       CHARACTER, PARAMETER            :: lchar(0:38) = (/'s', 'p', 'd', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o', &
     108             :                                                          'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', &
     109             :                                                          'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x', 'x'/)
     110             : 
     111             :       CHARACTER(len=2)                ::  nchar
     112             :       CHARACTER(len=2)                ::  noel(atoms%ntype)
     113             :       CHARACTER(len=10)               ::  fname(atoms%ntype)
     114             :       ! writing to a file
     115             :       INTEGER, PARAMETER              ::  iounit = 125
     116             :       CHARACTER(10), PARAMETER        ::  ioname = 'mixbas'
     117             :       LOGICAL                         ::  l_found
     118             : 
     119           0 :       IF (mpi%irank == 0) WRITE (6, '(//A,I2,A)') '### subroutine: mixedbasis ###'
     120             : 
     121           0 :       IF (xcpot%is_name("exx")) CALL judft_error("EXX is not implemented in this version", calledby='mixedbasis')
     122             : 
     123             :       ! Deallocate arrays which might have been allocated in a previous run of this subroutine
     124           0 :       IF (ALLOCATED(hybrid%ngptm)) DEALLOCATE (hybrid%ngptm)
     125           0 :       IF (ALLOCATED(hybrid%ngptm1)) DEALLOCATE (hybrid%ngptm1)
     126           0 :       IF (ALLOCATED(hybrid%nindxm1)) DEALLOCATE (hybrid%nindxm1)
     127           0 :       IF (ALLOCATED(hybrid%pgptm)) DEALLOCATE (hybrid%pgptm)
     128           0 :       IF (ALLOCATED(hybrid%pgptm1)) DEALLOCATE (hybrid%pgptm1)
     129           0 :       IF (ALLOCATED(hybrid%gptm)) DEALLOCATE (hybrid%gptm)
     130           0 :       IF (ALLOCATED(hybrid%basm1)) DEALLOCATE (hybrid%basm1)
     131             : 
     132           0 :       CALL usdus%init(atoms, input%jspins)
     133             : 
     134             :       ! If restart is specified read file if it already exists. If not create it.
     135           0 :       IF (l_restart) THEN
     136             : 
     137             :          ! Test if file exists
     138           0 :          INQUIRE (FILE=ioname, EXIST=l_found)
     139             : 
     140             :          IF (l_found .AND. .FALSE.) THEN !reading not working yet
     141             :             ! Open file
     142             :             OPEN (UNIT=iounit, FILE=ioname, FORM='unformatted', STATUS='old')
     143             : 
     144             :             ! Read array sizes
     145             :             !READ(iounit) kpts%nkptf,hybrid%gptmd
     146             :             READ (iounit) hybrid%maxgptm, hybrid%maxindx
     147             :             ! Allocate necessary array size and read arrays
     148             :             ALLOCATE (hybrid%ngptm(kpts%nkptf), hybrid%gptm(3, hybrid%gptmd))
     149             :             ALLOCATE (hybrid%pgptm(hybrid%maxgptm, kpts%nkptf))
     150             :             READ (iounit) hybrid%ngptm, hybrid%gptm, hybrid%pgptm, hybrid%nindx
     151             : 
     152             :             ! Read array sizes
     153             :             READ (iounit) hybrid%maxgptm1, hybrid%maxlcutm1, hybrid%maxindxm1, hybrid%maxindxp1
     154             :             ! Allocate necessary array size and read arrays
     155             :             ALLOCATE (hybrid%ngptm1(kpts%nkptf), hybrid%pgptm1(hybrid%maxgptm1, kpts%nkptf))
     156             :             ALLOCATE (hybrid%nindxm1(0:hybrid%maxlcutm1, atoms%ntype))
     157             :             ALLOCATE (hybrid%basm1(atoms%jmtd, hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype))
     158             :             READ (iounit) hybrid%ngptm1, hybrid%pgptm1, hybrid%nindxm1
     159             :             READ (iounit) hybrid%basm1
     160             : 
     161             :             CLOSE (iounit)
     162             : 
     163             :             RETURN
     164             :          END IF
     165             :       END IF
     166             : 
     167           0 :       hybrid%maxindx = 0
     168           0 :       DO itype = 1, atoms%ntype
     169           0 :          DO l = 0, atoms%lmax(itype)
     170           0 :             hybrid%maxindx = MAX(hybrid%maxindx, 2 + COUNT(atoms%llo(:atoms%nlo(itype), itype) == l))
     171             :          END DO
     172             :       END DO
     173             :       ! maxindx = maxval(nlo) + 2
     174             : 
     175             :       ! initialize gridf for radial integration
     176           0 :       CALL intgrf_init(atoms%ntype, atoms%jmtd, atoms%jri, atoms%dx, atoms%rmsh, gridf)
     177             : 
     178           0 :       ALLOCATE (vr0(atoms%jmtd, atoms%ntype, input%jspins))
     179             : 
     180           0 :       vr0(:, :, :) = v%mt(:, 0, :, :)
     181             : 
     182             :       ! calculate radial basisfunctions u and u' with
     183             :       ! the spherical part of the potential vr0 and store them in
     184             :       ! bas1 = large component ,bas2 = small component
     185             : 
     186           0 :       ALLOCATE (f(atoms%jmtd, 2, 0:atoms%lmaxd), df(atoms%jmtd, 2, 0:atoms%lmaxd))
     187           0 :       ALLOCATE (bas1(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype, input%jspins))
     188           0 :       ALLOCATE (bas2(atoms%jmtd, hybrid%maxindx, 0:atoms%lmaxd, atoms%ntype, input%jspins))
     189             : 
     190           0 :       DO itype = 1, atoms%ntype
     191           0 :          ng = atoms%jri(itype)
     192           0 :          DO ispin = 1, input%jspins
     193           0 :             DO l = 0, atoms%lmax(itype)
     194             :                CALL radfun(l, itype, ispin, enpara%el0(l, itype, ispin), vr0(:, itype, ispin), atoms, &
     195           0 :                            f(:, :, l), df(:, :, l), usdus, nodem, noded, wronk)
     196             :             END DO
     197           0 :             bas1(1:ng, 1, 0:atoms%lmaxd, itype, ispin) = f(1:ng, 1, 0:atoms%lmaxd)
     198           0 :             bas2(1:ng, 1, 0:atoms%lmaxd, itype, ispin) = f(1:ng, 2, 0:atoms%lmaxd)
     199           0 :             bas1(1:ng, 2, 0:atoms%lmaxd, itype, ispin) = df(1:ng, 1, 0:atoms%lmaxd)
     200           0 :             bas2(1:ng, 2, 0:atoms%lmaxd, itype, ispin) = df(1:ng, 2, 0:atoms%lmaxd)
     201             : 
     202           0 :             hybrid%nindx(:, itype) = 2
     203             :             ! generate radial functions for local orbitals
     204           0 :             IF (atoms%nlo(itype) >= 1) THEN
     205             :                CALL radflo(atoms, itype, ispin, enpara%ello0(1, 1, ispin), vr0(:, itype, ispin), &
     206           0 :                            f, df, mpi, usdus, uuilon, duilon, ulouilopn, flo)
     207             : 
     208           0 :                DO ilo = 1, atoms%nlo(itype)
     209           0 :                   hybrid%nindx(atoms%llo(ilo, itype), itype) = hybrid%nindx(atoms%llo(ilo, itype), itype) + 1
     210           0 :                   bas1(1:ng, hybrid%nindx(atoms%llo(ilo, itype), itype), atoms%llo(ilo, itype), itype, ispin) = flo(1:ng, 1, ilo)
     211           0 :                   bas2(1:ng, hybrid%nindx(atoms%llo(ilo, itype), itype), atoms%llo(ilo, itype), itype, ispin) = flo(1:ng, 2, ilo)
     212             :                END DO
     213             :             END IF
     214             :          END DO
     215             :       END DO
     216             : 
     217           0 :       DEALLOCATE (f, df)
     218             : 
     219             :       ! the radial functions are normalized
     220           0 :       DO ispin = 1, input%jspins
     221           0 :          DO itype = 1, atoms%ntype
     222           0 :             DO l = 0, atoms%lmax(itype)
     223           0 :                DO i = 1, hybrid%nindx(l, itype)
     224             :                   rdum = intgrf(bas1(:, i, l, itype, ispin)**2 + bas2(:, i, l, itype, ispin)**2, &
     225           0 :                                 atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     226           0 :                   bas1(:atoms%jri(itype), i, l, itype, ispin) = bas1(:atoms%jri(itype), i, l, itype, ispin)/SQRT(rdum)
     227           0 :                   bas2(:atoms%jri(itype), i, l, itype, ispin) = bas2(:atoms%jri(itype), i, l, itype, ispin)/SQRT(rdum)
     228             :                END DO
     229             :             END DO
     230             :          END DO
     231             :       END DO
     232             : 
     233             :       ! - - - - - - SETUP OF THE MIXED BASIS IN THE IR - - - - - - -
     234             : 
     235             :       ! construct G-vectors with cutoff smaller than gcutm
     236           0 :       gcutm = hybrid%gcutm1
     237           0 :       ALLOCATE (hybrid%ngptm(kpts%nkptf))
     238             : 
     239           0 :       hybrid%ngptm = 0
     240           0 :       i = 0
     241           0 :       n = -1
     242             : 
     243           0 :       rdum1 = MAXVAL((/(SQRT(SUM(MATMUL(kpts%bkf(:, ikpt), cell%bmat)**2)), ikpt=1, kpts%nkptf)/))
     244             : 
     245             :       ! a first run for the determination of the dimensions of the fields gptm,pgptm
     246             : 
     247             :       DO
     248           0 :          n = n + 1
     249           0 :          ldum = .FALSE.
     250           0 :          DO x = -n, n
     251           0 :             n1 = n - ABS(x)
     252           0 :             DO y = -n1, n1
     253           0 :                n2 = n1 - ABS(y)
     254           0 :                DO z = -n2, n2, MAX(2*n2, 1)
     255           0 :                   g = (/x, y, z/)
     256           0 :                   rdum = SQRT(SUM(MATMUL(g, cell%bmat)**2)) - rdum1
     257           0 :                   IF (rdum > gcutm) CYCLE
     258             :                   ldum1 = .FALSE.
     259           0 :                   DO ikpt = 1, kpts%nkptf
     260           0 :                      kvec = kpts%bkf(:, ikpt)
     261           0 :                      rdum = SUM(MATMUL(kvec + g, cell%bmat)**2)
     262             : 
     263           0 :                      IF (rdum <= gcutm**2) THEN
     264           0 :                         IF (.NOT. ldum1) THEN
     265           0 :                            i = i + 1
     266           0 :                            ldum1 = .TRUE.
     267             :                         END IF
     268             : 
     269           0 :                         hybrid%ngptm(ikpt) = hybrid%ngptm(ikpt) + 1
     270           0 :                         ldum = .TRUE.
     271             :                      END IF
     272             :                   END DO
     273             :                END DO
     274             :             END DO
     275             :          END DO
     276           0 :          IF (.NOT. ldum) EXIT
     277             :       END DO
     278             : 
     279           0 :       hybrid%gptmd = i
     280           0 :       hybrid%maxgptm = MAXVAL(hybrid%ngptm)
     281             : 
     282           0 :       ALLOCATE (hybrid%gptm(3, hybrid%gptmd))
     283           0 :       ALLOCATE (hybrid%pgptm(hybrid%maxgptm, kpts%nkptf))
     284             : 
     285           0 :       hybrid%gptm = 0
     286           0 :       hybrid%pgptm = 0
     287           0 :       hybrid%ngptm = 0
     288             : 
     289           0 :       i = 0
     290           0 :       n = -1
     291             : 
     292             :       ! Allocate and initialize arrays needed for G vector ordering
     293           0 :       ALLOCATE (unsrt_pgptm(hybrid%maxgptm, kpts%nkptf))
     294           0 :       ALLOCATE (length_kG(hybrid%maxgptm, kpts%nkptf))
     295           0 :       length_kG = 0
     296           0 :       unsrt_pgptm = 0
     297             : 
     298             :       DO
     299           0 :          n = n + 1
     300           0 :          ldum = .FALSE.
     301           0 :          DO x = -n, n
     302           0 :             n1 = n - ABS(x)
     303           0 :             DO y = -n1, n1
     304           0 :                n2 = n1 - ABS(y)
     305           0 :                DO z = -n2, n2, MAX(2*n2, 1)
     306           0 :                   g = (/x, y, z/)
     307           0 :                   rdum = SQRT(SUM(MATMUL(g, cell%bmat)**2)) - rdum1
     308           0 :                   IF (rdum > gcutm) CYCLE
     309             :                   ldum1 = .FALSE.
     310           0 :                   DO ikpt = 1, kpts%nkptf
     311           0 :                      kvec = kpts%bkf(:, ikpt)
     312           0 :                      rdum = SUM(MATMUL(kvec + g, cell%bmat)**2)
     313             : 
     314           0 :                      IF (rdum <= (gcutm)**2) THEN
     315           0 :                         IF (.NOT. ldum1) THEN
     316           0 :                            i = i + 1
     317           0 :                            hybrid%gptm(:, i) = g
     318             :                            ldum1 = .TRUE.
     319             :                         END IF
     320             : 
     321           0 :                         hybrid%ngptm(ikpt) = hybrid%ngptm(ikpt) + 1
     322           0 :                         ldum = .TRUE.
     323             : 
     324             :                         ! Position of the vector is saved as pointer
     325           0 :                         unsrt_pgptm(hybrid%ngptm(ikpt), ikpt) = i
     326             :                         ! Save length of vector k + G for array sorting
     327           0 :                         length_kG(hybrid%ngptm(ikpt), ikpt) = rdum
     328             :                      END IF
     329             :                   END DO
     330             :                END DO
     331             :             END DO
     332             :          END DO
     333           0 :          IF (.NOT. ldum) EXIT
     334             :       END DO
     335             : 
     336             :       ! Sort pointers in array, so that shortest |k+G| comes first
     337           0 :       DO ikpt = 1, kpts%nkptf
     338           0 :          ALLOCATE (ptr(hybrid%ngptm(ikpt)))
     339             :          ! Divide and conquer algorithm for arrays > 1000 entries
     340           0 :          divconq = MAX(0, INT(1.443*LOG(0.001*hybrid%ngptm(ikpt))))
     341             :          ! create pointers which correspond to a sorted array
     342           0 :          CALL rorderpf(ptr, length_kG(1:hybrid%ngptm(ikpt), ikpt), hybrid%ngptm(ikpt), divconq)
     343             :          ! rearrange old pointers
     344           0 :          DO igpt = 1, hybrid%ngptm(ikpt)
     345           0 :             hybrid%pgptm(igpt, ikpt) = unsrt_pgptm(ptr(igpt), ikpt)
     346             :          END DO
     347           0 :          DEALLOCATE (ptr)
     348             :       END DO
     349           0 :       DEALLOCATE (unsrt_pgptm)
     350           0 :       DEALLOCATE (length_kG)
     351             : 
     352             :       ! construct IR mixed basis set for the representation of the non local exchange elements with cutoff gcutm1
     353             : 
     354             :       ! first run to determine dimension of pgptm1
     355           0 :       ALLOCATE (hybrid%ngptm1(kpts%nkptf))
     356           0 :       hybrid%ngptm1 = 0
     357           0 :       DO igpt = 1, hybrid%gptmd
     358           0 :          g = hybrid%gptm(:, igpt)
     359           0 :          DO ikpt = 1, kpts%nkptf
     360           0 :             kvec = kpts%bkf(:, ikpt)
     361           0 :             rdum = SUM(MATMUL(kvec + g, cell%bmat)**2)
     362           0 :             IF (rdum <= hybrid%gcutm1**2) THEN
     363           0 :                hybrid%ngptm1(ikpt) = hybrid%ngptm1(ikpt) + 1
     364             :             END IF
     365             :          END DO
     366             :       END DO
     367           0 :       hybrid%maxgptm1 = MAXVAL(hybrid%ngptm1)
     368             : 
     369           0 :       ALLOCATE (hybrid%pgptm1(hybrid%maxgptm1, kpts%nkptf))
     370           0 :       hybrid%pgptm1 = 0
     371           0 :       hybrid%ngptm1 = 0
     372             : 
     373             :       ! Allocate and initialize arrays needed for G vector ordering
     374           0 :       ALLOCATE (unsrt_pgptm(hybrid%maxgptm1, kpts%nkptf))
     375           0 :       ALLOCATE (length_kG(hybrid%maxgptm1, kpts%nkptf))
     376           0 :       length_kG = 0
     377           0 :       unsrt_pgptm = 0
     378           0 :       DO igpt = 1, hybrid%gptmd
     379           0 :          g = hybrid%gptm(:, igpt)
     380           0 :          DO ikpt = 1, kpts%nkptf
     381           0 :             kvec = kpts%bkf(:, ikpt)
     382           0 :             rdum = SUM(MATMUL(kvec + g, cell%bmat)**2)
     383           0 :             IF (rdum <= hybrid%gcutm1**2) THEN
     384           0 :                hybrid%ngptm1(ikpt) = hybrid%ngptm1(ikpt) + 1
     385           0 :                unsrt_pgptm(hybrid%ngptm1(ikpt), ikpt) = igpt
     386           0 :                length_kG(hybrid%ngptm1(ikpt), ikpt) = rdum
     387             :             END IF
     388             :          END DO
     389             :       END DO
     390             : 
     391             :       ! Sort pointers in array, so that shortest |k+G| comes first
     392           0 :       DO ikpt = 1, kpts%nkptf
     393           0 :          ALLOCATE (ptr(hybrid%ngptm1(ikpt)))
     394             :          ! Divide and conquer algorithm for arrays > 1000 entries
     395           0 :          divconq = MAX(0, INT(1.443*LOG(0.001*hybrid%ngptm1(ikpt))))
     396             :          ! create pointers which correspond to a sorted array
     397           0 :          CALL rorderpf(ptr, length_kG(1:hybrid%ngptm1(ikpt), ikpt), hybrid%ngptm1(ikpt), divconq)
     398             :          ! rearrange old pointers
     399           0 :          DO igpt = 1, hybrid%ngptm1(ikpt)
     400           0 :             hybrid%pgptm1(igpt, ikpt) = unsrt_pgptm(ptr(igpt), ikpt)
     401             :          END DO
     402           0 :          DEALLOCATE (ptr)
     403             :       END DO
     404           0 :       DEALLOCATE (unsrt_pgptm)
     405           0 :       DEALLOCATE (length_kG)
     406             : 
     407           0 :       IF (mpi%irank == 0) THEN
     408           0 :          WRITE (6, '(/A)') 'Mixed basis'
     409           0 :          WRITE (6, '(A,I5)') 'Number of unique G-vectors: ', hybrid%gptmd
     410           0 :          WRITE (6, *)
     411           0 :          WRITE (6, '(3x,A)') 'IR Plane-wave basis with cutoff of gcutm (hybrid%gcutm1/2*input%rkmax):'
     412           0 :          WRITE (6, '(5x,A,I5)') 'Maximal number of G-vectors:', hybrid%maxgptm
     413           0 :          WRITE (6, *)
     414           0 :          WRITE (6, *)
     415           0 :          WRITE (6, '(3x,A)') 'IR Plane-wave basis for non-local exchange potential:'
     416           0 :          WRITE (6, '(5x,A,I5)') 'Maximal number of G-vectors:', hybrid%maxgptm1
     417           0 :          WRITE (6, *)
     418             :       END IF
     419             : 
     420             :       ! - - - - - - - - Set up MT product basis for the non-local exchange potential  - - - - - - - - - -
     421             : 
     422           0 :       IF (mpi%irank == 0) THEN
     423           0 :          WRITE (6, '(A)') 'MT product basis for non-local exchange potential:'
     424           0 :          WRITE (6, '(A)') 'Reduction due to overlap (quality of orthonormality, should be < 1.0E-06)'
     425             :       END IF
     426             : 
     427           0 :       hybrid%maxlcutm1 = MAXVAL(hybrid%lcutm1)
     428             : 
     429           0 :       ALLOCATE (hybrid%nindxm1(0:hybrid%maxlcutm1, atoms%ntype))
     430           0 :       ALLOCATE (seleco(hybrid%maxindx, 0:atoms%lmaxd), selecu(hybrid%maxindx, 0:atoms%lmaxd))
     431           0 :       ALLOCATE (selecmat(hybrid%maxindx, 0:atoms%lmaxd, hybrid%maxindx, 0:atoms%lmaxd))
     432           0 :       hybrid%nindxm1 = 0    !!! 01/12/10 jij%M.b.
     433             : 
     434             :       ! determine maximal indices of (radial) mixed-basis functions (->nindxm1)
     435             :       ! (will be reduced later-on due to overlap)
     436           0 :       hybrid%maxindxp1 = 0
     437           0 :       DO itype = 1, atoms%ntype
     438           0 :          seleco = .FALSE.
     439           0 :          selecu = .FALSE.
     440           0 :          seleco(1, 0:hybrid%select1(1, itype)) = .TRUE.
     441           0 :          selecu(1, 0:hybrid%select1(3, itype)) = .TRUE.
     442           0 :          seleco(2, 0:hybrid%select1(2, itype)) = .TRUE.
     443           0 :          selecu(2, 0:hybrid%select1(4, itype)) = .TRUE.
     444             : 
     445             :          ! include local orbitals
     446           0 :          IF (hybrid%maxindx >= 3) THEN
     447           0 :             seleco(3:, :) = .TRUE.
     448           0 :             selecu(3:, :) = .TRUE.
     449             :          END IF
     450             : 
     451           0 :          DO l = 0, hybrid%lcutm1(itype)
     452           0 :             n = 0
     453           0 :             M = 0
     454             : 
     455             :             !
     456             :             ! valence * valence
     457             :             !
     458             : 
     459             :             ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
     460             :             selecmat = RESHAPE((/((((seleco(n1, l1) .AND. selecu(n2, l2), &
     461             :                                      n1=1, hybrid%maxindx), l1=0, atoms%lmaxd), n2=1, hybrid%maxindx), l2=0, atoms%lmaxd)/), &
     462           0 :                                (/hybrid%maxindx, atoms%lmaxd + 1, hybrid%maxindx, atoms%lmaxd + 1/))
     463             : 
     464           0 :             DO l1 = 0, atoms%lmax(itype)
     465           0 :                DO l2 = 0, atoms%lmax(itype)
     466           0 :                   IF (l >= ABS(l1 - l2) .AND. l <= l1 + l2) THEN
     467           0 :                      DO n1 = 1, hybrid%nindx(l1, itype)
     468           0 :                         DO n2 = 1, hybrid%nindx(l2, itype)
     469           0 :                            M = M + 1
     470           0 :                            IF (selecmat(n1, l1, n2, l2)) THEN
     471           0 :                               n = n + 1
     472           0 :                               selecmat(n2, l2, n1, l1) = .FALSE. ! prevent double counting of products (a*b = b*a)
     473             :                            END IF
     474             :                         END DO
     475             :                      END DO
     476             :                   END IF
     477             :                END DO
     478             :             END DO
     479           0 :             IF (n == 0 .AND. mpi%irank == 0) &
     480             :                WRITE (6, '(A)') 'mixedbasis: Warning!  No basis-function product of '//lchar(l)// &
     481           0 :                '-angular momentum defined.'
     482           0 :             hybrid%maxindxp1 = MAX(hybrid%maxindxp1, M)
     483           0 :             hybrid%nindxm1(l, itype) = n*input%jspins
     484             :          END DO
     485             :       END DO
     486           0 :       hybrid%maxindxm1 = MAXVAL(hybrid%nindxm1)
     487             : 
     488           0 :       ALLOCATE (hybrid%basm1(atoms%jmtd, hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype))
     489           0 :       hybrid%basm1 = 0
     490             : 
     491             :       ! Define product bases and reduce them according to overlap
     492             : 
     493           0 :       DO itype = 1, atoms%ntype
     494           0 :          seleco = .FALSE.
     495           0 :          selecu = .FALSE.
     496           0 :          seleco(1, 0:hybrid%select1(1, itype)) = .TRUE.
     497           0 :          selecu(1, 0:hybrid%select1(3, itype)) = .TRUE.
     498           0 :          seleco(2, 0:hybrid%select1(2, itype)) = .TRUE.
     499           0 :          selecu(2, 0:hybrid%select1(4, itype)) = .TRUE.
     500             :          ! include lo's
     501           0 :          IF (hybrid%maxindx >= 3) THEN
     502           0 :             seleco(3:, :) = .TRUE.
     503           0 :             selecu(3:, :) = .TRUE.
     504             :          END IF
     505           0 :          IF (atoms%ntype > 1 .AND. mpi%irank == 0)&
     506           0 :               &    WRITE (6, '(6X,A,I3)') 'Atom type', itype
     507           0 :          ng = atoms%jri(itype)
     508           0 :          DO l = 0, hybrid%lcutm1(itype)
     509           0 :             n = hybrid%nindxm1(l, itype)
     510             :             ! allow for zero product-basis functions for
     511             :             ! current l-quantum number
     512           0 :             IF (n == 0) THEN
     513           0 :                IF (mpi%irank == 0) WRITE (6, '(6X,A,'':   0 ->   0'')') lchar(l)
     514             :                CYCLE
     515             :             END IF
     516             : 
     517             :             ! set up the overlap matrix
     518           0 :             ALLOCATE (olap(n, n), eigv(n, n), work(3*n), eig(n), ihelp(n))
     519           0 :             ihelp = 1 ! initialize to avoid a segfault
     520           0 :             i = 0
     521             : 
     522             :             ! valence*valence
     523             : 
     524             :             ! Condense seleco and seleco into selecmat (each product corresponds to a matrix element)
     525             :             selecmat = RESHAPE((/((((seleco(n1, l1) .AND. selecu(n2, l2), &
     526             :                                      n1=1, hybrid%maxindx), l1=0, atoms%lmaxd), n2=1, hybrid%maxindx), l2=0, atoms%lmaxd)/), &
     527           0 :                                (/hybrid%maxindx, atoms%lmaxd + 1, hybrid%maxindx, atoms%lmaxd + 1/))
     528             : 
     529           0 :             DO l1 = 0, atoms%lmax(itype)
     530           0 :                DO l2 = 0, atoms%lmax(itype)
     531           0 :                   IF (l < ABS(l1 - l2) .OR. l > l1 + l2) CYCLE
     532             : 
     533           0 :                   DO n1 = 1, hybrid%nindx(l1, itype)
     534           0 :                      DO n2 = 1, hybrid%nindx(l2, itype)
     535             : 
     536           0 :                         IF (selecmat(n1, l1, n2, l2)) THEN
     537           0 :                            DO ispin = 1, input%jspins
     538           0 :                               i = i + 1
     539           0 :                               IF (i > n) call judft_error('got too many product functions', hint='This is a BUG, please report', calledby='mixedbasis')
     540             : 
     541             :                               hybrid%basm1(:ng, i, l, itype) &
     542             :                                  = (bas1(:ng, n1, l1, itype, ispin) &
     543             :                                     *bas1(:ng, n2, l2, itype, ispin) &
     544             :                                     + bas2(:ng, n1, l1, itype, ispin) &
     545           0 :                                     *bas2(:ng, n2, l2, itype, ispin))/atoms%rmsh(:ng, itype)
     546             : 
     547             :                               !normalize basm1
     548             :                               rdum = SQRT(intgrf(hybrid%basm1(:, i, l, itype)**2, &
     549           0 :                                                  atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf))
     550             : 
     551           0 :                               hybrid%basm1(:ng, i, l, itype) = hybrid%basm1(:ng, i, l, itype)/rdum
     552             : 
     553             :                            END DO !ispin
     554             :                            ! prevent double counting of products (a*b = b*a)
     555           0 :                            selecmat(n2, l2, n1, l1) = .FALSE.
     556             :                         END IF
     557             : 
     558             :                      END DO !n2
     559             :                   END DO !n1
     560             : 
     561             :                END DO !l2
     562             :             END DO  !l1
     563             : 
     564           0 :             IF (i /= n) call judft_error('counting error for product functions', hint='This is a BUG, please report', calledby='mixedbasis')
     565             : 
     566             :             ! In order to get ride of the linear dependencies in the
     567             :             ! radial functions basm1 belonging to fixed l and itype
     568             :             ! the overlap matrix is diagonalized and those eigenvectors
     569             :             ! with a eigenvalue greater then hybrid%tolerance1 are retained
     570             : 
     571             :             ! Calculate overlap
     572           0 :             olap = 0
     573           0 :             DO n2 = 1, n
     574           0 :                DO n1 = 1, n2
     575             :                   olap(n1, n2) = intgrf(hybrid%basm1(:, n1, l, itype)*hybrid%basm1(:, n2, l, itype), &
     576           0 :                                         atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     577           0 :                   olap(n2, n1) = olap(n1, n2)
     578             :                END DO
     579             :             END DO
     580             : 
     581             :             ! Diagonalize
     582           0 :             CALL diagonalize(eigv, eig, olap)
     583             : 
     584             :             ! Get rid of linear dependencies (eigenvalue <= hybrid%tolerance1)
     585           0 :             nn = 0
     586           0 :             DO i = 1, hybrid%nindxm1(l, itype)
     587           0 :                IF (eig(i) > hybrid%tolerance1) THEN
     588           0 :                   nn = nn + 1
     589           0 :                   ihelp(nn) = i
     590             :                END IF
     591             :             END DO
     592           0 :             hybrid%nindxm1(l, itype) = nn
     593           0 :             eig = eig(ihelp)
     594           0 :             eigv(:, :) = eigv(:, ihelp)
     595             : 
     596           0 :             DO i = 1, ng
     597           0 :                hybrid%basm1(i, 1:nn, l, itype) = MATMUL(hybrid%basm1(i, 1:n, l, itype), eigv(:, 1:nn))/SQRT(eig(:nn))
     598             :             END DO
     599             : 
     600             :             ! Add constant function to l=0 basis and then do a Gram-Schmidt orthonormalization
     601           0 :             IF (l == 0) THEN
     602             : 
     603             :                ! Check if basm1 must be reallocated
     604           0 :                IF (nn + 1 > SIZE(hybrid%basm1, 2)) THEN
     605           0 :                   ALLOCATE (basmhlp(atoms%jmtd, nn + 1, 0:hybrid%maxlcutm1, atoms%ntype))
     606           0 :                   basmhlp(:, 1:nn, :, :) = hybrid%basm1
     607           0 :                   DEALLOCATE (hybrid%basm1)
     608           0 :                   ALLOCATE (hybrid%basm1(atoms%jmtd, nn + 1, 0:hybrid%maxlcutm1, atoms%ntype))
     609           0 :                   hybrid%basm1(:, 1:nn, :, :) = basmhlp(:, 1:nn, :, :)
     610           0 :                   DEALLOCATE (basmhlp)
     611             :                END IF
     612             : 
     613           0 :                hybrid%basm1(:ng, nn + 1, 0, itype) = atoms%rmsh(:ng, itype)/SQRT(atoms%rmsh(ng, itype)**3/3)
     614           0 :                DO i = nn, 1, -1
     615           0 :                   DO j = i + 1, nn + 1
     616             :                      hybrid%basm1(:ng, i, 0, itype) = hybrid%basm1(:ng, i, 0, itype) &
     617             :                                                       - intgrf(hybrid%basm1(:ng, i, 0, itype)*hybrid%basm1(:ng, j, 0, itype), &
     618             :                                                                atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf) &
     619           0 :                                                       *hybrid%basm1(:ng, j, 0, itype)
     620             : 
     621             :                   END DO
     622             :                   hybrid%basm1(:ng, i, 0, itype) = hybrid%basm1(:ng, i, 0, itype) &
     623             :                                                    /SQRT(intgrf(hybrid%basm1(:ng, i, 0, itype)**2, &
     624           0 :                                                                 atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf))
     625             :                END DO
     626           0 :                nn = nn + 1
     627           0 :                DEALLOCATE (olap)
     628           0 :                ALLOCATE (olap(nn, nn))
     629           0 :                hybrid%nindxm1(l, itype) = nn
     630             :             END IF
     631             : 
     632             :             ! Check orthonormality of product basis
     633           0 :             rdum = 0
     634           0 :             DO i = 1, nn
     635           0 :                DO j = 1, i
     636             :                   rdum1 = intgrf(hybrid%basm1(:, i, l, itype)*hybrid%basm1(:, j, l, itype), &
     637           0 :                                  atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     638             : 
     639           0 :                   IF (i == j) THEN
     640           0 :                      rdum1 = ABS(1 - rdum1)
     641           0 :                      rdum = rdum + rdum1**2
     642             :                   ELSE
     643           0 :                      rdum1 = ABS(rdum1)
     644           0 :                      rdum = rdum + 2*rdum1**2
     645             :                   END IF
     646             : 
     647           0 :                   IF (rdum1 > 1e-6) THEN
     648           0 :                      IF (mpi%irank == 0) THEN
     649             :                         WRITE (6, '(A)') 'mixedbasis: Bad orthonormality of ' &
     650           0 :                            //lchar(l)//'-product basis. Increase tolerance.'
     651           0 :                         WRITE (6, '(12X,A,F9.6,A,2(I3.3,A))') 'Deviation of', &
     652           0 :                            rdum1, ' occurred for (', i, ',', j, ')-overlap.'
     653             :                      END IF
     654           0 :                      CALL judft_error("Bad orthonormality of product basis", hint='Increase tolerance', calledby='mixedbasis')
     655             :                   END IF
     656             : 
     657             :                END DO
     658             :             END DO
     659             : 
     660           0 :             IF (mpi%irank == 0) THEN
     661           0 :                WRITE (6, '(6X,A,I4,'' ->'',I4,''   ('',ES8.1,'' )'')') lchar(l)//':', n, nn, SQRT(rdum)/nn
     662             :             END IF
     663             : 
     664           0 :             DEALLOCATE (olap, eigv, work, eig, ihelp)
     665             : 
     666             :          END DO !l
     667           0 :          IF (mpi%irank == 0) WRITE (6, '(6X,A,I7)') 'Total:', SUM(hybrid%nindxm1(0:hybrid%lcutm1(itype), itype))
     668             :       END DO ! itype
     669             : 
     670           0 :       hybrid%maxindxm1 = MAXVAL(hybrid%nindxm1)
     671             : 
     672           0 :       ALLOCATE (basmhlp(atoms%jmtd, hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype))
     673             :       basmhlp(1:atoms%jmtd, 1:hybrid%maxindxm1, 0:hybrid%maxlcutm1, 1:atoms%ntype) &
     674           0 :          = hybrid%basm1(1:atoms%jmtd, 1:hybrid%maxindxm1, 0:hybrid%maxlcutm1, 1:atoms%ntype)
     675           0 :       DEALLOCATE (hybrid%basm1)
     676           0 :       ALLOCATE (hybrid%basm1(atoms%jmtd, hybrid%maxindxm1, 0:hybrid%maxlcutm1, atoms%ntype))
     677           0 :       hybrid%basm1 = basmhlp
     678             : 
     679           0 :       DEALLOCATE (basmhlp, seleco, selecu, selecmat)
     680             : 
     681             :       !
     682             :       ! now we build linear combinations of the radial functions
     683             :       ! such that they possess no moment except one radial function in each l-channel
     684             :       !
     685           0 :       IF (mpi%irank == 0) THEN
     686             :          WRITE (6, '(/,A,/,A)') 'Build linear combinations of radial '// &
     687           0 :             'functions in each l-channel,', &
     688             :             'such that they possess no multipolmoment'// &
     689           0 :             ' except the last function:'
     690             : 
     691           0 :          WRITE (6, '(/,17x,A)') 'moment  (quality of orthonormality)'
     692             :       END IF
     693           0 :       DO itype = 1, atoms%ntype
     694           0 :          ng = atoms%jri(itype)
     695             : 
     696           0 :          IF (atoms%ntype > 1 .AND. mpi%irank == 0) WRITE (6, '(6X,A,I3)') 'Atom type', itype
     697             : 
     698           0 :          DO l = 0, hybrid%lcutm1(itype)
     699             :             ! determine radial function with the largest moment
     700             :             ! this function is used to build the linear combinations
     701           0 :             rdum = 0
     702           0 :             DO i = 1, hybrid%nindxm1(l, itype)
     703             :                rdum1 = intgrf(atoms%rmsh(:ng, itype)**(l + 1)*hybrid%basm1(:ng, i, l, itype), &
     704           0 :                               atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     705           0 :                IF (ABS(rdum1) > rdum) THEN
     706           0 :                   n = i
     707           0 :                   rdum = rdum1
     708             :                END IF
     709             :             END DO
     710             : 
     711             :             ! rearrange order of radial functions such that the last function possesses the largest moment
     712           0 :             j = 0
     713           0 :             bashlp(:ng) = hybrid%basm1(:ng, n, l, itype)
     714           0 :             DO i = 1, hybrid%nindxm1(l, itype)
     715           0 :                IF (i == n) CYCLE
     716           0 :                j = j + 1
     717           0 :                hybrid%basm1(:ng, j, l, itype) = hybrid%basm1(:ng, i, l, itype)
     718             :             END DO
     719           0 :             hybrid%basm1(:ng, hybrid%nindxm1(l, itype), l, itype) = bashlp(:ng)
     720             : 
     721             :          END DO
     722             : 
     723           0 :          DO l = 0, hybrid%lcutm1(itype)
     724           0 :             IF (mpi%irank == 0) WRITE (6, '(6X,A)') lchar(l)//':'
     725             : 
     726           0 :             IF (hybrid%nindxm1(l, itype) == 0) THEN
     727           0 :                IF (mpi%irank == 0) WRITE (6, '(6X,A,'':   0 ->    '')') lchar(l)
     728             :                CYCLE
     729             :             END IF
     730             : 
     731           0 :             n = hybrid%nindxm1(l, itype)
     732           0 :             DO i = 1, hybrid%nindxm1(l, itype)
     733           0 :                IF (i == n) CYCLE
     734             :                ! calculate moment of radial function i
     735             :                rdum1 = intgrf(atoms%rmsh(:ng, itype)**(l + 1)*hybrid%basm1(:ng, i, l, itype), &
     736           0 :                               atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     737             : 
     738             :                rdum = intgrf(atoms%rmsh(:ng, itype)**(l + 1)*hybrid%basm1(:ng, n, l, itype), &
     739           0 :                              atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     740             : 
     741           0 :                bashlp(:ng) = hybrid%basm1(:ng, n, l, itype)
     742             : 
     743           0 :                IF (SQRT(rdum**2 + rdum1**2) <= 1E-06 .AND. mpi%irank == 0) &
     744           0 :                   WRITE (6, *) 'Warning: Norm is smaller thann 1E-06!'
     745             : 
     746             :                ! change function n such that n is orthogonal to i
     747             :                ! since the functions basm1 have been orthogonal on input
     748             :                ! the linear combination does not destroy the orthogonality to the residual functions
     749             :                hybrid%basm1(:ng, n, l, itype) = rdum/SQRT(rdum**2 + rdum1**2)*bashlp(:ng) &
     750           0 :                                                 + rdum1/SQRT(rdum**2 + rdum1**2)*hybrid%basm1(:ng, i, l, itype)
     751             : 
     752             :                ! combine basis function i and n so that they possess no momemt
     753             :                hybrid%basm1(:ng, i, l, itype) = rdum1/SQRT(rdum**2 + rdum1**2)*bashlp(:ng) &
     754           0 :                                                 - rdum/SQRT(rdum**2 + rdum1**2)*hybrid%basm1(:ng, i, l, itype)
     755             : 
     756             :                rdum1 = intgrf(atoms%rmsh(:ng, itype)**(l + 1)*hybrid%basm1(:ng, i, l, itype), &
     757           0 :                               atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     758             : 
     759           0 :                IF (rdum1 > 1E-10) call judft_error('moment of radial function does not vanish', calledby='mixedbasis')
     760             : 
     761           0 :                IF (mpi%irank == 0) WRITE (6, '(6x,I4,'' ->  '',ES8.1)') i, rdum1
     762             :             END DO
     763             : 
     764             :             ! test orthogonality
     765           0 :             rdum = 0
     766           0 :             DO i = 1, hybrid%nindxm1(l, itype)
     767           0 :                DO j = 1, hybrid%nindxm1(l, itype)
     768             :                   rdum1 = intgrf(hybrid%basm1(:ng, i, l, itype)*hybrid%basm1(:ng, j, l, itype), &
     769           0 :                                  atoms%jri, atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf)
     770           0 :                   IF (i /= j) THEN
     771           0 :                      rdum = rdum + rdum1
     772             :                   ELSE
     773           0 :                      rdum = rdum + ABS(1 - rdum1)
     774             :                   END IF
     775             :                END DO
     776             :             END DO
     777           0 :             IF (mpi%irank == 0) &
     778           0 :                WRITE (6, '(6x,I4,'' ->'',f10.5,''   ('',ES8.1,'' )'')') n, &
     779             :                intgrf(atoms%rmsh(:ng, itype)**(l + 1)*hybrid%basm1(:ng, n, l, itype), atoms%jri, &
     780           0 :                       atoms%jmtd, atoms%rmsh, atoms%dx, atoms%ntype, itype, gridf), rdum
     781             :          END DO
     782             :       END DO
     783             : 
     784           0 :       DO itype = 1, atoms%ntype
     785           0 :          IF (ANY(hybrid%nindxm1(0:hybrid%lcutm1(itype), itype) == 0)) call judft_error('any hybrid%nindxm1 eq 0', calledby='mixedbasis')
     786             :       END DO
     787             : 
     788             :       !count basis functions
     789           0 :       hybrid%nbasp = 0
     790           0 :       DO itype = 1, atoms%ntype
     791           0 :          DO i = 1, atoms%neq(itype)
     792           0 :             DO l = 0, hybrid%lcutm1(itype)
     793           0 :                DO M = -l, l
     794           0 :                   DO j = 1, hybrid%nindxm1(l, itype)
     795           0 :                      hybrid%nbasp = hybrid%nbasp + 1
     796             :                   END DO
     797             :                END DO
     798             :             END DO
     799             :          END DO
     800             :       END DO
     801           0 :       hybrid%maxbasm1 = hybrid%nbasp + hybrid%maxgptm
     802           0 :       DO nk = 1, kpts%nkptf
     803           0 :          hybrid%nbasm(nk) = hybrid%nbasp + hybrid%ngptm(nk)
     804             :       END DO
     805             : 
     806           0 :       hybrid%maxlmindx = MAXVAL((/(SUM((/(hybrid%nindx(l, itype)*(2*l + 1), l=0, atoms%lmax(itype))/)), itype=1, atoms%ntype)/))
     807             : 
     808             :       !
     809             :       !
     810             :       ! Write all information to a file to enable restarting
     811             :       !
     812           0 :       IF (l_restart .AND. mpi%irank == 0) THEN
     813             : 
     814           0 :          OPEN (UNIT=iounit, FILE=ioname, FORM='unformatted', STATUS='replace')
     815             : 
     816           0 :          WRITE (iounit) kpts%nkptf, hybrid%gptmd
     817           0 :          WRITE (iounit) hybrid%maxgptm, hybrid%maxindx
     818           0 :          WRITE (iounit) hybrid%ngptm, hybrid%gptm, hybrid%pgptm, hybrid%nindx
     819             : 
     820           0 :          WRITE (iounit) hybrid%maxgptm1, hybrid%maxlcutm1, hybrid%maxindxm1, hybrid%maxindxp1
     821           0 :          WRITE (iounit) hybrid%ngptm1, hybrid%pgptm1, hybrid%nindxm1
     822           0 :          WRITE (iounit) hybrid%basm1
     823             : 
     824           0 :          CLOSE (iounit)
     825             : 
     826           0 :          l_restart = .FALSE.
     827             : 
     828             :       END IF
     829             : 
     830           0 :    END SUBROUTINE mixedbasis
     831             : 
     832           0 :    SUBROUTINE read_atompotential(atoms, fname, potential)
     833             : 
     834             :       USE m_types
     835             :       IMPLICIT NONE
     836             :       TYPE(t_atoms), INTENT(IN)   :: atoms
     837             : 
     838             :       ! - scalars -
     839             :       CHARACTER(10), INTENT(IN) ::   fname(atoms%ntype)
     840             : 
     841             :       ! - arrays -
     842             :       REAL, INTENT(OUT)::   potential(atoms%jmtd, atoms%ntype)
     843             :       ! - local scalars -
     844             :       INTEGER                 ::   i, j, itype, ic, m
     845             :       REAL                    ::   rdum, r, n
     846             :       ! - local arrays -
     847           0 :       REAL, ALLOCATABLE    ::   v_atom(:), r_atom(:)
     848             : 
     849           0 :       potential = 0
     850           0 :       DO itype = 1, atoms%ntype
     851           0 :          OPEN (331, file=fname(itype), form='formatted')
     852             : 
     853           0 :          ic = 0
     854             :          DO
     855           0 :             READ (331, *) rdum
     856           0 :             ic = ic + 1
     857           0 :             IF (rdum > (atoms%rmt(itype) + 0.2)) EXIT
     858             :          END DO
     859             : 
     860           0 :          ALLOCATE (v_atom(ic), r_atom(ic))
     861           0 :          REWIND (331)
     862             : 
     863           0 :          DO i = 1, ic
     864           0 :             READ (331, *) r_atom(i), v_atom(i)
     865             :          END DO
     866             : 
     867             :          ! transfer potential from grid used in the atom program
     868             :          ! to those used here
     869           0 :          DO i = 1, atoms%jri(itype)
     870           0 :             r = atoms%rmsh(i, itype)
     871           0 :             DO j = 1, ic - 1
     872           0 :                IF (r_atom(j) > r) THEN
     873           0 :                   M = (v_atom(j) - v_atom(j + 1))/(r_atom(j) - r_atom(j + 1))
     874           0 :                   n = v_atom(j) - M*r_atom(j)
     875           0 :                   potential(i, itype) = n + M*r
     876             :                   !               WRITE(333,'(4f25.10)') n + m*r_atom(j),v_atom(j),n + m*r_atom(j+1),v_atom(j+1)
     877           0 :                   EXIT
     878             :                END IF
     879             :             END DO
     880             :          END DO
     881             : 
     882           0 :          DEALLOCATE (v_atom, r_atom)
     883           0 :          CLOSE (331)
     884             :       END DO
     885             : 
     886           0 :    END SUBROUTINE read_atompotential
     887             : 
     888             : END MODULE m_mixedbasis
     889             : 

Generated by: LCOV version 1.13