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

Matthias Redies's avatar
Matthias Redies committed
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132
   SUBROUTINE hsfock(nk, atoms, hybrid, lapw, dimension, kpts, jsp, input, hybdat, eig_irr, sym, cell, noco, &
                     results, it, mnobd, xcpot, mpi)

      USE m_types
      USE m_symm_hf
      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
      TYPE(t_hybrid), INTENT(INOUT) :: hybrid
      TYPE(t_hybdat), INTENT(INOUT) :: hybdat
      TYPE(t_results), INTENT(INOUT) :: results

      ! 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
      INTEGER                 ::  nbasfcn
      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))
      INTEGER                 ::  rrot(3, 3, sym%nsym)
      INTEGER                 ::  psym(sym%nsym) ! Note: psym is only filled up to index nsymop

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

      REAL                    ::  wl_iks(dimension%neigd, kpts%nkptf)

      TYPE(t_mat)             :: olap, trafo, invtrafo, ex, tmp, v_x, z
      COMPLEX                 ::  exch(dimension%neigd, dimension%neigd)
      COMPLEX, ALLOCATABLE     ::  carr(:)

      CALL timestart("total time hsfock")

      ! preparations

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

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

      ! read in lower triangle part of overlap matrix from direct acces file olap
      call timestart("read in olap")
      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)
      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
133 134
         END DO
      ELSE
Matthias Redies's avatar
Matthias Redies committed
135 136 137 138 139 140
         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)
141
      END IF
Matthias Redies's avatar
Matthias Redies committed
142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167
      call timestop("read in olap")

      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 == 1 .and. mpi%irank == 0) WRITE (*, *) 'calculate new HF matrix'
         IF (nk == 1 .and. jsp == 1 .and. input%imix > 10) CALL system('rm -f broyd*')
         ! calculate all symmetrie operations, which yield k invariant

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

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

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

         ! remove weights(wtkpt) in w_iks
         DO ikpt = 1, kpts%nkptf
            DO iband = 1, dimension%neigd
               ikpt0 = kpts%bkp(ikpt)
               wl_iks(iband, ikpt) = results%w_iks(iband, ikpt0, jsp)/(kpts%wtkpt(ikpt0)*kpts%nkptf)
            END DO
         END DO
168

Matthias Redies's avatar
Matthias Redies committed
169 170 171 172 173 174
         ! 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, &
                                  mpi, ex)
175

Matthias Redies's avatar
Matthias Redies committed
176
         CALL timestart("core exchange calculation")
177

Matthias Redies's avatar
Matthias Redies committed
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
         ! calculate contribution from the core states to the HF exchange
         IF (xcpot%is_name("hse") .OR. xcpot%is_name("vhse")) THEN
            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 < hybrid%nbands(nk)) STOP " mhsfock: neigd  < nbands(nk) ;trafo from wavefunctions to APW requires at least nbands(nk)"

         call z%init(olap%l_real, nbasfcn, dimension%neigd)
194
         call read_z(z, kpts%nkptf*(jsp - 1) + nk)
Matthias Redies's avatar
Matthias Redies committed
195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210
         z%matsize2 = hybrid%nbands(nk) ! reduce "visible matsize" for the following computations

         call olap%multiply(z, trafo)

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

         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
211 212
         ENDDO

Matthias Redies's avatar
Matthias Redies committed
213 214 215 216
         CALL ex%multiply(invtrafo, tmp)
         CALL trafo%multiply(tmp, v_x)

         CALL timestop("time for performing T^-1*mat_ex*T^-1*")
217

Matthias Redies's avatar
Matthias Redies committed
218 219 220
         call timestart("symmetrizeh")
         CALL symmetrizeh(atoms, kpts%bkf(:, nk), dimension, jsp, lapw, sym, hybdat%kveclo_eig, cell, nsymop, psym, v_x)
         call timestop("symmetrizeh")
221

Matthias Redies's avatar
Matthias Redies committed
222 223
         CALL write_v_x(v_x, kpts%nkpt*(jsp - 1) + nk)
      END IF ! hybrid%l_calhf
224

Matthias Redies's avatar
Matthias Redies committed
225
      CALL timestop("total time hsfock")
226

Matthias Redies's avatar
Matthias Redies committed
227
   END SUBROUTINE hsfock
228

229
END MODULE m_hsfock