incdrhoscf_nc.f90 7.25 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
!
! Copyright (C) 2001-2016 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
subroutine incdrhoscf_nc (drhoscf, weight, ik, dbecsum, dpsi)
  !-----------------------------------------------------------------------
  !
  !     This routine computes the change of the charge density due to the
  !     perturbation. It is called at the end of the computation of the
  !     change of the wavefunction for a given k point.
  !
  USE kinds,                ONLY : DP
  USE ions_base,            ONLY : nat
  USE cell_base,            ONLY : omega
19
  USE fft_base,             ONLY : dffts, dfftp
20 21 22 23 24
  USE fft_interfaces,       ONLY : invfft
  USE lsda_mod,             ONLY : nspin
  USE spin_orb,             ONLY : domag
  USE noncollin_module,     ONLY : npol, nspin_mag
  USE uspp_param,           ONLY : nhm
25
  USE wvfct,                ONLY : npwx, nbnd
26
  USE wavefunctions, ONLY : evc
27 28
  USE klist,                ONLY : ngk,igk_k
  USE qpoint,               ONLY : ikks, ikqs
29 30 31
  USE control_lr,           ONLY : nbnd_occ
  USE mp_bands,             ONLY : me_bgrp, inter_bgrp_comm, ntask_groups
  USE mp,                   ONLY : mp_sum
32
  USE fft_helper_subroutines
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56

  IMPLICIT NONE
  !
  ! I/O variables
  INTEGER, INTENT(IN) :: ik
  ! input: the k point
  REAL(DP), INTENT(IN) :: weight
  ! input: the weight of the k point
  COMPLEX(DP), INTENT(IN) :: dpsi(npwx*npol,nbnd)
  ! input: the perturbed wfcs at the given k point
  COMPLEX(DP), INTENT(INOUT) :: drhoscf (dfftp%nnr,nspin_mag), dbecsum (nhm,nhm,nat,nspin)
  ! input/output: the accumulated change of the charge density and dbecsum
  !
  !   here the local variable
  !
  REAL(DP) :: wgt
  ! the effective weight of the k point
  !
  COMPLEX(DP), ALLOCATABLE :: psi (:,:), dpsic (:,:)
  ! the wavefunctions in real space
  ! the change of wavefunctions in real space
  !
  COMPLEX(DP), ALLOCATABLE :: tg_psi (:,:), tg_dpsi (:,:), tg_drho(:,:)
  !
57
  INTEGER :: npw, npwq, ikk, ikq
58
  INTEGER :: ibnd, jbnd, ir, ir3, ig, incr, v_siz, idx, ioff, ioff_tg, nxyp
59
  INTEGER :: ntgrp, right_inc
60 61 62 63 64 65 66 67 68
  ! counters
  !
  CALL start_clock ('incdrhoscf')
  !
  ALLOCATE (dpsic(dffts%nnr, npol))
  ALLOCATE (psi  (dffts%nnr, npol))
  !
  wgt = 2.d0 * weight / omega
  ikk = ikks(ik)
69 70 71
  ikq = ikqs(ik)
  npw = ngk(ikk)
  npwq= ngk(ikq)
72 73
  incr = 1
  !
74
  IF (dffts%has_task_groups) THEN
75
     !
76
     v_siz = dffts%nnr_tg
77 78 79 80 81
     !
     ALLOCATE( tg_psi( v_siz, npol ) )
     ALLOCATE( tg_dpsi( v_siz, npol ) )
     ALLOCATE( tg_drho( v_siz, nspin_mag ) )
     !
82
     incr  = fftx_ntgrp(dffts)
83 84 85 86 87 88 89 90
     !
  ENDIF
  !
  ! dpsi contains the   perturbed wavefunctions of this k point
  ! evc  contains the unperturbed wavefunctions of this k point
  !
  do ibnd = 1, nbnd_occ(ikk), incr

91
     IF (dffts%has_task_groups) THEN
92 93 94 95 96 97
        !
        tg_drho=(0.0_DP, 0.0_DP)
        tg_psi=(0.0_DP, 0.0_DP)
        tg_dpsi=(0.0_DP, 0.0_DP)
        !
        ioff   = 0
98 99
        CALL tg_get_recip_inc( dffts, right_inc )
        ntgrp = fftx_ntgrp( dffts )
100
        !
101
        DO idx = 1, ntgrp
102
           !
103
           ! ... dtgs%nogrp ffts at the same time. We prepare both
104 105 106 107 108
           ! evc (at k) and dpsi (at k+q)
           !
           IF( idx + ibnd - 1 <= nbnd_occ(ikk) ) THEN
              !
              DO ig = 1, npw
109 110
                 tg_psi( dffts%nl( igk_k( ig,ikk ) ) + ioff, 1 ) = evc( ig, idx+ibnd-1 )
                 tg_psi( dffts%nl( igk_k( ig,ikk ) ) + ioff, 2 ) = evc( npwx+ig, idx+ibnd-1 )
111 112
              END DO
              DO ig = 1, npwq
113 114
                 tg_dpsi( dffts%nl( igk_k( ig,ikq ) ) + ioff, 1 ) = dpsi( ig, idx+ibnd-1 )
                 tg_dpsi( dffts%nl( igk_k( ig,ikq ) ) + ioff, 2 ) = dpsi( npwx+ig, idx+ibnd-1 )
115 116 117 118
              END DO
              !
           END IF
           !
119
           ioff = ioff + right_inc
120 121
           !
        END DO
122 123 124 125
        CALL invfft ('tgWave', tg_psi(:,1), dffts)
        CALL invfft ('tgWave', tg_psi(:,2), dffts)
        CALL invfft ('tgWave', tg_dpsi(:,1), dffts)
        CALL invfft ('tgWave', tg_dpsi(:,2), dffts)
126

127 128 129
        do ir = 1, dffts%nr1x * dffts%nr2x * dffts%my_nr3p
           tg_drho (ir,1) = tg_drho (ir,1) + wgt * (CONJG(tg_psi (ir,1) )*tg_dpsi (ir,1) &
                                                  + CONJG(tg_psi (ir,2) )*tg_dpsi (ir,2) )
130 131
        enddo
        IF (domag) THEN
132 133 134 135 136 137 138
           do ir = 1, dffts%nr1x * dffts%nr2x * dffts%my_nr3p
              tg_drho(ir,2)= tg_drho(ir,2) + wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,2) &
                                                  + CONJG(tg_psi(ir,2))*tg_dpsi(ir,1) )
              tg_drho(ir,3)= tg_drho(ir,3) + wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,2) &
                                                  - CONJG(tg_psi(ir,2))*tg_dpsi(ir,1) ) * (0.d0,-1.d0)
              tg_drho(ir,4)= tg_drho(ir,4) + wgt * (CONJG(tg_psi(ir,1))*tg_dpsi(ir,1) &
                                                  - CONJG(tg_psi(ir,2))*tg_dpsi(ir,2) )
139 140 141 142 143 144
           enddo
        ENDIF
        !
        ! reduce the group charge (equivalent to sum over the bands of the
        ! orbital group)
        !
145
        CALL tg_reduce_rho( drhoscf, tg_drho, dffts )
146 147 148 149 150 151 152 153 154
        !
     ELSE
        !
        ! Normal case: no task groups
        !
        ! FFT to R-space of the unperturbed wfct's evc
        !
        psi = (0.d0, 0.d0)
        do ig = 1, npw
155 156
           psi (dffts%nl (igk_k(ig,ikk) ), 1) = evc (ig, ibnd)
           psi (dffts%nl (igk_k(ig,ikk) ), 2) = evc (ig+npwx, ibnd)
157 158 159 160 161 162 163 164
        enddo
        CALL invfft ('Wave', psi(:,1), dffts)
        CALL invfft ('Wave', psi(:,2), dffts)
        !
        ! FFT to R-space of the perturbed wfct's dpsi
        !
        dpsic = (0.d0, 0.d0)
        do ig = 1, npwq
165 166
           dpsic (dffts%nl (igk_k(ig,ikq)), 1 ) = dpsi (ig, ibnd)
           dpsic (dffts%nl (igk_k(ig,ikq)), 2 ) = dpsi (ig+npwx, ibnd)
167 168 169 170 171 172 173 174 175 176 177 178
        enddo
        CALL invfft ('Wave', dpsic(:,1), dffts)
        CALL invfft ('Wave', dpsic(:,2), dffts)
        !
        ! Calculation of the response charge density
        !
        do ir = 1, dffts%nnr
           drhoscf(ir,1)=drhoscf(ir,1)+wgt*(CONJG(psi(ir,1))*dpsic(ir,1)  +  &
                                            CONJG(psi(ir,2))*dpsic(ir,2) )
        enddo
        IF (domag) THEN
           do ir = 1, dffts%nnr
179 180 181 182 183 184
              drhoscf(ir,2)=drhoscf (ir,2) + wgt * (CONJG(psi(ir,1))*dpsic(ir,2) &
                                                  + CONJG(psi(ir,2))*dpsic(ir,1) )
              drhoscf(ir,3)=drhoscf (ir,3) + wgt * (CONJG(psi(ir,1))*dpsic(ir,2) &
                                                  - CONJG(psi(ir,2))*dpsic(ir,1) ) * (0.d0,-1.d0)
              drhoscf(ir,4)=drhoscf (ir,4) + wgt * (CONJG(psi(ir,1))*dpsic(ir,1) &
                                                  - CONJG(psi(ir,2))*dpsic(ir,2) )
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
           enddo
        END IF
        !
     END IF
     !
  enddo
  !
  ! Ultrasoft contribution
  ! Calculate dbecsum_nc = <evc|vkb><vkb|dpsi>
  !
  CALL addusdbec_nc (ik, weight, dpsi, dbecsum)
  !
  DEALLOCATE(psi)
  DEALLOCATE(dpsic)
  !
200
  IF (dffts%has_task_groups) THEN
201 202 203 204 205 206 207 208 209 210
     DEALLOCATE(tg_psi)
     DEALLOCATE(tg_dpsi)
     DEALLOCATE(tg_drho)
  END IF
  !
  CALL stop_clock ('incdrhoscf')
  !
  RETURN
  !
end subroutine incdrhoscf_nc