hsfock.F90 11.4 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
CONTAINS

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

   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
   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)    :: 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
91
   INTEGER                 ::  nbasfcn
92 93 94 95 96 97 98 99
   INTEGER                 ::  nsymop
   INTEGER                 ::  nkpt_EIBZ
   INTEGER                 ::  ncstd
   INTEGER                 ::  ok
   REAL                    ::  a_ex

   ! local arrays
   INTEGER                 ::  nsest(hybrid%nbands(nk)),indx_sest(hybrid%nbands(nk),hybrid%nbands(nk))
100 101
   INTEGER                 ::  rrot(3,3,sym%nsym)
   INTEGER                 ::  psym(sym%nsym) ! Note: psym is only filled up to index nsymop
102

103
   INTEGER,ALLOCATABLE     ::  parent(:)
104 105 106 107
   INTEGER,ALLOCATABLE     ::  pointer_EIBZ(:)
   INTEGER,ALLOCATABLE     ::  n_q(:)

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

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

120 121
   ! initialize weighting factor for HF exchange part
   a_ex=xcpot%get_exchange_weight()
122

123
   ! read in lower triangle part of overlap matrix from direct acces file olap
124 125
   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)
126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
   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
141 142 143 144 145 146 147

   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

148 149 150
      ALLOCATE(parent(kpts%nkptf), stat=ok)
      IF(ok.NE.0) STOP 'mhsfock: failure allocation parent'
      parent = 0
151

152
      CALL timestart("symm_hf")
153
      CALL symm_hf_init(sym,kpts,nk,nsymop,rrot,psym)
154

155
      CALL symm_hf(kpts,nk,sym,dimension,hybdat,eig_irr,atoms,hybrid,cell,lapw,jsp,mpi,&
156
                   rrot,nsymop,psym,nkpt_EIBZ,n_q,parent,pointer_EIBZ,nsest,indx_sest)
157
      CALL timestop("symm_hf")
158

159 160 161
      ! remove weights(wtkpt) in w_iks
      DO ikpt=1,kpts%nkptf
         DO iband=1,dimension%neigd
162
            ikpt0 = kpts%bkp(ikpt)
163 164 165 166 167 168 169 170 171 172
            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,&
173
                               mpi,ex)
174
      CALL timestop("valence exchange calculation")
175

176 177 178 179 180 181 182 183 184 185 186
      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

187 188 189 190
      CALL timestart("core exchange calculation")

      ! calculate contribution from the core states to the HF exchange
      IF (xcpot%is_name("hse").OR.xcpot%is_name("vhse")) THEN
191
#ifdef CPP_NEVER           
192 193
         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)
194
#endif
195 196 197 198 199 200 201 202 203 204 205
         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
206
      IF(dimension%neigd.LT.hybrid%nbands(nk)) STOP " mhsfock: neigd  < nbands(nk) ;trafo from wavefunctions to APW requires at least nbands(nk)"
207

208 209 210 211
      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

212
      call olap%multiply(z,trafo)
213

214
      CALL invtrafo%alloc(olap%l_real,hybrid%nbands(nk),nbasfcn)
215
      CALL trafo%TRANSPOSE(invtrafo)
216 217
      IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)

218 219 220 221 222 223 224 225 226 227
      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

228 229 230 231 232 233 234 235 236 237 238
      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

239 240
      CALL ex%multiply(invtrafo,tmp)
      CALL trafo%multiply(tmp,v_x)
Daniel Wortmann's avatar
Daniel Wortmann committed
241
        
242
      CALL timestop("time for performing T^-1*mat_ex*T^-1*")
243

244 245 246 247 248 249 250 251 252 253 254
      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

255
      CALL symmetrizeh(atoms,kpts%bkf(:,nk),dimension,jsp,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,v_x)
256 257 258 259 260 261 262 263 264 265 266

      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
267

268 269
      CALL write_v_x(v_x,kpts%nkpt*(jsp-1) + nk)
   END IF ! hybrid%l_calhf
270

271
   CALL timestop("total time hsfock")
272

273
END SUBROUTINE hsfock
274

275
END MODULE m_hsfock