hsfock.F90 12.1 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
MODULE m_hsfock
8 9

   USE m_judft
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29
! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c
!     This module is the driver routine for the calculation of the Hartree    c
!     Fock exchange term by using the mixed basis set.                        c
!                                                                             c
!     hsfock                                                                  c
!         |                                                                   c
!         |- symm.F:                                                          c
!         |  calculates the irreducible representation                        c
!         |                                                                   c
!         |- wavefproducts.F:                 s      s*                       c
!         |  computes the repsentation of phi    phi       in the mixed basis c
!         |                                  n,k    n',k+q                    c
!         |                                                                   c
!         |- exchange.F:                                                      c
!         |  calculates valence-valence part of the exchange matrix (mat_ex), c
!         |                                                                   c
!         |- exchange_core.F                                                  c
!         |  calculate valence-core contribution                              c
!                                                                             c
!     variables:                                                              c
Daniel Wortmann's avatar
Daniel Wortmann committed
30 31
!         kpts%nkptf   :=   number of kpoints                                      c
!         kpts%nkpt   :=   number of irreducible kpoints                          c
32 33 34 35 36 37 38 39 40 41 42 43 44
!         nbands  :=   number of bands for which the exchange matrix (mat_ex) c
!                      in the space of the wavefunctions is calculated        c
!         te_hfex :=   hf exchange contribution to the total energy           c
!         mnobd   :=   maximum number of occupied bands                       c
!         parent  :=   parent(ikpt) points to the symmetry equivalent point   c
!                      under the little group of kpoint nk                    c
!         symop   :=   symop(ikpt) points to the symmetry operation, which    c
!                      maps parent(ikpt) on ikpt                              c
!                                                                             c
!                                                                             c
!                                               M.Betzinger (09/07)           c
! c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c c

45 46 47 48 49 50
CONTAINS

SUBROUTINE hsfock(nk,atoms,hybrid,lapw,dimension,kpts,jsp,input,hybdat,eig_irr,sym,cell,noco,&
                  results,it,mnobd,xcpot,mpi,irank2,isize2,comm)

   USE m_types
51
   USE m_symm_hf
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
   USE m_util          ,ONLY: intgrf,intgrf_init
   USE m_exchange_valence_hf
   USE m_exchange_core
   USE m_symmetrizeh
   USE m_wrapper
   USE m_hsefunctional ,ONLY: exchange_vccvHSE,exchange_ccccHSE
   USE m_io_hybrid

   IMPLICIT NONE

   TYPE(t_hybdat),        INTENT(IN)    :: hybdat
   TYPE(t_results),       INTENT(INOUT) :: results
   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_input),         INTENT(IN)    :: input
   TYPE(t_noco),          INTENT(IN)    :: noco
   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_lapw),          INTENT(IN)    :: lapw

   ! scalars
   INTEGER,               INTENT(IN)    :: jsp 
   INTEGER,               INTENT(IN)    :: it
   INTEGER,               INTENT(IN)    :: irank2 ,isize2,comm
   INTEGER,               INTENT(IN)    :: nk
   INTEGER,               INTENT(IN)    :: mnobd

   ! arrays
   REAL,                  INTENT(IN)    :: eig_irr(dimension%neigd,kpts%nkpt)

   ! local scalars
   INTEGER                 ::  i,j,ic,ic1,l,itype,n,nn
   INTEGER                 ::  iband,iband1,iband2
   INTEGER                 ::  ikpt,ikpt0
   INTEGER                 ::  irec
   INTEGER                 ::  irecl_olap,irecl_z,irecl_vx
92
   INTEGER                 ::  maxndb, nbasfcn
93 94 95 96 97 98 99 100 101 102
   INTEGER                 ::  nddb
   INTEGER                 ::  nsymop
   INTEGER                 ::  nkpt_EIBZ
   INTEGER                 ::  ncstd
   INTEGER                 ::  ok
   REAL                    ::  a_ex

   ! local arrays
   INTEGER                 ::  degenerat(hybrid%ne_eig(nk))
   INTEGER                 ::  nsest(hybrid%nbands(nk)),indx_sest(hybrid%nbands(nk),hybrid%nbands(nk))
103 104
   INTEGER                 ::  rrot(3,3,sym%nsym)
   INTEGER                 ::  psym(sym%nsym) ! Note: psym is only filled up to index nsymop
105 106 107 108 109 110

   INTEGER,ALLOCATABLE     ::  parent(:),symop(:)
   INTEGER,ALLOCATABLE     ::  pointer_EIBZ(:)
   INTEGER,ALLOCATABLE     ::  n_q(:)

   REAL                    ::  wl_iks(dimension%neigd,kpts%nkptf)
111
    
112 113 114 115
   TYPE(t_mat)             :: olap,trafo,invtrafo,ex,tmp,v_x,z
   COMPLEX                 ::  exch(dimension%neigd,dimension%neigd)
   COMPLEX,ALLOCATABLE     ::  carr(:)
   COMPLEX,ALLOCATABLE     ::  rep_c(:,:,:,:,:)
116
      
117
   CALL timestart("total time hsfock")
118
    
119
   ! preparations
120

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

124 125
   ! initialize weighting factor for HF exchange part
   a_ex=xcpot%get_exchange_weight()
126

127
   ! read in lower triangle part of overlap matrix from direct acces file olap
128 129
   nbasfcn = MERGE(lapw%nv(1)+lapw%nv(2)+2*atoms%nlotot,lapw%nv(1)+atoms%nlotot,noco%l_noco)
   call olap%alloc(sym%invs,nbasfcn)
130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
   call read_olap(olap, kpts%nkpt*(jsp-1)+nk)
   IF (olap%l_real) THEN
      DO i=1,nbasfcn
        DO j=1,i
          olap%data_r(i,j) = olap%data_r(j,i)
        END DO
      END DO
   ELSE
      DO i=1,nbasfcn
        DO j=1,i
          olap%data_c(i,j) = CONJG(olap%data_c(j,i))
        END DO
      END DO
      olap%data_c=conjg(olap%data_c)
   END IF
145 146 147 148 149 150 151 152 153 154 155

   IF(hybrid%l_calhf) THEN
      ncstd = sum( (/ ( (hybdat%nindxc(l,itype)*(2*l+1)*atoms%neq(itype),l=0,hybdat%lmaxc(itype)), itype = 1,atoms%ntype) /) )
      IF( nk .eq. 1 .and. mpi%irank == 0 ) WRITE(*,*) 'calculate new HF matrix'
      IF( nk .eq. 1 .and. jsp .eq. 1 .and. input%imix .gt. 10) CALL system('rm -f broyd*')
      ! calculate all symmetrie operations, which yield k invariant

      ALLOCATE( parent(kpts%nkptf),symop(kpts%nkptf) ,stat=ok)
      IF( ok .ne. 0 ) STOP 'mhsfock: failure allocation parent/symop'
      parent = 0 ; symop = 0

156 157 158 159 160
      CALL timestart("symm_hf")
      CALL symm_hf_init(sym,kpts,nk,irank2,nsymop,rrot,psym)

      ALLOCATE(rep_c(-hybdat%lmaxcd:hybdat%lmaxcd,-hybdat%lmaxcd:hybdat%lmaxcd,0:hybdat%lmaxcd,nsymop,atoms%nat), stat=ok)
      IF(ok.NE.0) STOP 'hsfock: failure allocation rep_c'
161

162 163
      CALL symm_hf(kpts,nk,sym,dimension,hybdat,eig_irr,atoms,hybrid,cell,lapw,jsp,mpi,irank2,&
                   rrot,nsymop,psym,nkpt_EIBZ,n_q,parent,symop,degenerat,pointer_EIBZ,maxndb,nddb,nsest,indx_sest,rep_c)
164
      CALL timestop("symm_hf")
165

166 167 168
      ! remove weights(wtkpt) in w_iks
      DO ikpt=1,kpts%nkptf
         DO iband=1,dimension%neigd
169
            ikpt0 = kpts%bkp(ikpt)
170 171 172 173 174 175 176 177 178 179 180 181 182
            wl_iks(iband,ikpt) = results%w_iks(iband,ikpt0,jsp) / (kpts%wtkpt(ikpt0)*kpts%nkptf)
         END DO
      END DO

      ! calculate contribution from valence electrons to the
      ! HF exchange
      CALL timestart("valence exchange calculation")
      ex%l_real=sym%invs
      CALL exchange_valence_hf(nk,kpts,nkpt_EIBZ, sym,atoms,hybrid,cell,dimension,input,jsp,hybdat,mnobd,lapw,&
                               eig_irr,results,parent,pointer_EIBZ,n_q,wl_iks,it,xcpot,noco,nsest,indx_sest,&
                               mpi,irank2,isize2,comm,ex)
      DEALLOCATE (rep_c)
      CALL timestop("valence exchange calculation")
183

184 185 186 187 188 189 190 191 192 193 194
      WRITE(1224,'(a,i7)') 'kpoint: ', nk
      DO i = 1, ex%matsize1
         DO j = 1, i
            IF (ex%l_real) THEN
               WRITE(1224,'(2i7,2f15.8)') i, j, ex%data_r(j,i) !ex%data_r(i,j), ex%data_r(j,i)
            ELSE
               WRITE(1224,'(2i7,4f15.8)') i, j, ex%data_c(j,i) !ex%data_c(i,j), ex%data_c(j,i)
            ENDIF
         END DO
      END DO

195 196 197 198 199 200
      CALL timestart("core exchange calculation")
      ! do the rest of the calculation only on master
      IF (irank2 /= 0) RETURN

      ! calculate contribution from the core states to the HF exchange
      IF (xcpot%is_name("hse").OR.xcpot%is_name("vhse")) THEN
201
#ifdef CPP_NEVER           
202 203
         CALL exchange_vccvHSE(nk,atoms,hybrid,hybdat,dimension,jsp,lapw,nsymop,nsest,indx_sest,mpi,a_ex,results,mat_ex%core)
         CALL exchange_ccccHSE(nk,obsolete,atoms,hybdat,ncstd,kpts(:,nk),sym,a_ex,mpi,results%core)
204
#endif
205 206 207 208 209 210 211 212 213 214 215 216 217 218
         STOP "HSE not implemented in hsfock"
      ELSE
         CALL exchange_vccv1(nk,atoms,hybrid,hybdat,dimension,jsp,lapw,nsymop,nsest,indx_sest,mpi,a_ex,results,ex)
         CALL exchange_cccc(nk,atoms,hybdat,ncstd,sym,kpts,a_ex,mpi,results)
      END IF

      DEALLOCATE(n_q)
      CALL timestop("core exchange calculation")

      CALL timestart("time for performing T^-1*mat_ex*T^-1*")
      !calculate trafo from wavefunctions to APW basis
      IF(dimension%neigd.LT.hybrid%nbands(nk)) STOP 'mhsfock: neigd  < nbands(nk) ; '& 
                                                    'trafo from wavefunctions to APW requires at least nbands(nk) '

219 220 221 222
      call z%init(olap%l_real,nbasfcn,dimension%neigd)
      call read_z(z,kpts%nkpt*(jsp-1)+nk)
      z%matsize2 = hybrid%nbands(nk) ! reduce "visible matsize" for the following computations

223
      call olap%multiply(z,trafo)
224

225
      CALL invtrafo%alloc(olap%l_real,hybrid%nbands(nk),nbasfcn)
226
      CALL trafo%TRANSPOSE(invtrafo)
227 228
      IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)

229 230 231 232 233 234 235 236 237 238
      DO i=1,hybrid%nbands(nk)
         DO j=1,i-1
            IF (ex%l_real) THEN
               ex%data_r(i,j)=ex%data_r(j,i)
            ELSE
               ex%data_c(i,j)=conjg(ex%data_c(j,i))
            END IF
         ENDDO
      ENDDO

239 240 241 242 243 244 245 246 247 248 249
      WRITE(1225,'(a,i7)') 'kpoint: ', nk
      DO i = 1, ex%matsize1
         DO j = 1, i
            IF (ex%l_real) THEN
               WRITE(1225,'(2i7,2f15.8)') i, j, ex%data_r(i,j), ex%data_r(j,i)
            ELSE
               WRITE(1225,'(2i7,4f15.8)') i, j, ex%data_c(i,j), ex%data_c(j,i)
            ENDIF
         END DO
      END DO

250 251
      CALL ex%multiply(invtrafo,tmp)
      CALL trafo%multiply(tmp,v_x)
Daniel Wortmann's avatar
Daniel Wortmann committed
252
        
253
      CALL timestop("time for performing T^-1*mat_ex*T^-1*")
254

255 256 257 258 259 260 261 262 263 264 265
      WRITE(1231,'(a,i7)') 'kpoint: ', nk
      DO i = 1, v_x%matsize1
         DO j = 1, i
            IF (v_x%l_real) THEN
               WRITE(1231,'(2i7,1f15.8)') i, j, v_x%data_r(i,j)
            ELSE
               WRITE(1231,'(2i7,2f15.8)') i, j, v_x%data_c(i,j)
            ENDIF
         END DO
      END DO

266
      CALL symmetrizeh(atoms,kpts%bkf(:,nk),dimension,jsp,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,v_x)
267 268 269 270 271 272 273 274 275 276 277

      WRITE(1232,'(a,i7)') 'kpoint: ', nk
      DO i = 1, v_x%matsize1
         DO j = 1, i
            IF (v_x%l_real) THEN
               WRITE(1232,'(2i7,1f15.8)') i, j, v_x%data_r(j,i) ! Note the different indices in comparison to points above. This is wanted!
            ELSE
               WRITE(1232,'(2i7,2f15.8)') i, j, v_x%data_c(j,i) ! Note the different indices in comparison to points above. This is wanted!
            ENDIF
         END DO
      END DO
278

279 280
      CALL write_v_x(v_x,kpts%nkpt*(jsp-1) + nk)
   END IF ! hybrid%l_calhf
281

282
   CALL timestop("total time hsfock")
283

284
END SUBROUTINE hsfock
285

286
END MODULE m_hsfock