hsfock.F90 9.95 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
   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_xcpot_inbuild), INTENT(IN)    :: xcpot
   TYPE(t_mpi),           INTENT(IN)    :: mpi
   TYPE(t_dimension),     INTENT(IN)    :: dimension
   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
Matthias Redies's avatar
Matthias Redies committed
72 73 74
   TYPE(t_hybrid),        INTENT(INOUT) :: hybrid
   TYPE(t_hybdat),        INTENT(INOUT) :: hybdat
   TYPE(t_results),       INTENT(INOUT) :: results
75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90

   ! 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
   call timestart("read in olap")
125 126
   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)
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
   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
142
   call timestop("read in olap")
143 144 145 146 147 148 149

   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

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

154
      CALL timestart("symm_hf")
155
      CALL symm_hf_init(sym,kpts,nk,nsymop,rrot,psym)
156

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

161 162 163
      ! remove weights(wtkpt) in w_iks
      DO ikpt=1,kpts%nkptf
         DO iband=1,dimension%neigd
164
            ikpt0 = kpts%bkp(ikpt)
165 166 167 168 169 170 171 172 173
            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
      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,&
174
                               mpi,ex)
Matthias Redies's avatar
Matthias Redies committed
175
   
176 177 178 179
      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
180
#ifdef CPP_NEVER           
181 182
         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)
183
#endif
184 185 186 187 188 189 190 191 192 193 194
         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
195
      IF(dimension%neigd.LT.hybrid%nbands(nk)) STOP " mhsfock: neigd  < nbands(nk) ;trafo from wavefunctions to APW requires at least nbands(nk)"
196

197 198 199 200
      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

201
      call olap%multiply(z,trafo)
202

203
      CALL invtrafo%alloc(olap%l_real,hybrid%nbands(nk),nbasfcn)
204
      CALL trafo%TRANSPOSE(invtrafo)
205 206
      IF(.NOT.invtrafo%l_real) invtrafo%data_c = CONJG(invtrafo%data_c)

207 208 209 210 211 212 213 214 215 216 217 218
      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

      CALL ex%multiply(invtrafo,tmp)
      CALL trafo%multiply(tmp,v_x)
Daniel Wortmann's avatar
Daniel Wortmann committed
219
        
220
      CALL timestop("time for performing T^-1*mat_ex*T^-1*")
221

Matthias Redies's avatar
Matthias Redies committed
222
      call timestart("symmetrizeh")
223
      CALL symmetrizeh(atoms,kpts%bkf(:,nk),dimension,jsp,lapw,sym,hybdat%kveclo_eig,cell,nsymop,psym,v_x)
Matthias Redies's avatar
Matthias Redies committed
224
      call timestop("symmetrizeh")
225

226 227
      CALL write_v_x(v_x,kpts%nkpt*(jsp-1) + nk)
   END IF ! hybrid%l_calhf
228

229
   CALL timestop("total time hsfock")
230

231
END SUBROUTINE hsfock
232

233
END MODULE m_hsfock