dvqpsi_us3.f90 6.97 KB
Newer Older
sponce's avatar
sponce committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
  !
  ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino 
  ! 
  ! Copyright (C) 2001-2003 PWSCF 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 .
  !
  ! adapted from PH/dvqpsi_us (QE)
  !
  !----------------------------------------------------------------------
  SUBROUTINE dvqpsi_us3 (ik, uact, addnlcc, xxk, xq0)
  !----------------------------------------------------------------------
  !
  ! This routine calculates dV_bare/dtau * psi for one perturbation
  ! with a given q. The displacements are described by a vector u.
  ! The result is stored in dvpsi. The routine is called for each k point
  ! and for each pattern u. It computes simultaneously all the bands.
  !
  ! RM - Nov/Dec 2014 
  ! Imported the noncolinear case implemented by xlzhang
  !
  USE kinds,                 ONLY : DP
  USE ions_base,             ONLY : nat, ityp
  USE cell_base,             ONLY : tpiba
  USE fft_base,              ONLY: dfftp, dffts
  USE fft_interfaces,        ONLY: fwfft, invfft
  USE gvect,                 ONLY : eigts1, eigts2, eigts3, mill, g, nl, &
                                    ngm
  USE gvecs,                 ONLY : ngms, doublegrid, nls
  USE lsda_mod,              ONLY : lsda, isk
  USE noncollin_module,      ONLY : npol
  use uspp_param,            ONLY : upf
35
  USE wvfct,                 ONLY : nbnd, npwx
sponce's avatar
sponce committed
36 37 38 39
  USE wavefunctions_module,  ONLY : evc
  USE nlcc_ph,               ONLY : drc
  USE uspp,                  ONLY : nlcc_any
  USE eqv,                   ONLY : dvpsi, dmuxc, vlocq
40 41 42
  USE qpoint,                ONLY : eigqts, npwq !, ikks
  USE klist,                 ONLY : ngk
  USE elph2,                 ONLY : igkq, igk
sponce's avatar
sponce committed
43 44 45 46 47

  implicit none
  !
  !   The dummy variables
  !
48
  integer :: ik, npw
sponce's avatar
sponce committed
49
  ! input: the k point
50
  real(kind=DP) :: xq0(3), xxk(3)
sponce's avatar
sponce committed
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
  complex(DP) :: uact (3 * nat)
  ! input: the pattern of displacements
  logical :: addnlcc
  !
  !   And the local variables
  !

  integer :: na, mu, ig, nt, ibnd, ir, is, ip
  ! counter on atoms
  ! counter on modes
  ! the point k
  ! counter on G vectors
  ! the type of atom
  ! counter on bands
  ! counter on real mesh

  complex(kind=DP) :: gtau, gu, fact, u1, u2, u3, gu0
  complex(kind=DP) , ALLOCATABLE, TARGET :: aux (:)
  complex(kind=DP) , ALLOCATABLE :: aux1 (:), aux2 (:)
  complex(kind=DP) , POINTER :: auxs (:)
  ! work space
  ! SP : seems to be useless ?
  !logical :: htg


  CALL start_clock ('dvqpsi_us')
  !IF (nlcc_any) THEN
  ! SP - changed according to QE5/PH/dvqpsi_us
  !htg = dffts%have_task_groups
80 81 82
  npw  = ngk(ik)


sponce's avatar
sponce committed
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 133 134 135 136 137 138 139 140 141 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 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230
  dffts%have_task_groups=.FALSE.
  IF (nlcc_any.AND.addnlcc) THEN
     ALLOCATE (aux( dfftp%nnr))
     IF (doublegrid) then
        ALLOCATE (auxs(dffts%nnr))
     ELSE
        auxs => aux
     ENDIF
  ENDIF
  ALLOCATE (aux1(dffts%nnr))
  ALLOCATE (aux2(dffts%nnr))
  !
  !    We start by computing the contribution of the local potential.
  !    The computation of the derivative of the local potential is done in
  !    reciprocal space while the product with the wavefunction is done in
  !    real space
  !
! RM - I don't think we need this
!  IF (lgamma) THEN
!     ikk = ik
!  ELSE
!     ikk = 2 * ik - 1
!  ENDIF
 ! ikk = ikks(ik)
  dvpsi(:,:) = (0.d0, 0.d0)
  aux1(:) = (0.d0, 0.d0)
  DO na = 1, nat
     fact = tpiba * (0.d0, -1.d0) * eigqts (na)
     mu = 3 * (na - 1)
     IF (abs (uact (mu + 1) ) + abs (uact (mu + 2) ) + abs (uact (mu + &
          3) ) .gt.1.0d-12) THEN
        nt = ityp (na)
        u1 = uact (mu + 1)
        u2 = uact (mu + 2)
        u3 = uact (mu + 3)
        gu0 = xq0 (1) * u1 + xq0 (2) * u2 + xq0 (3) * u3
        DO ig = 1, ngms
           gtau = eigts1 (mill(1,ig), na) * eigts2 (mill(2,ig), na) * eigts3 ( &
                mill(3,ig), na)
           gu = gu0 + g (1, ig) * u1 + g (2, ig) * u2 + g (3, ig) * u3
           aux1 (nls (ig) ) = aux1 (nls (ig) ) + vlocq (ig, nt) * gu * &
                fact * gtau
        ENDDO
     ENDIF
  ENDDO
  !
  ! add NLCC when present
  !
   IF (nlcc_any.and.addnlcc) THEN
      aux(:) = (0.d0, 0.d0)
      DO na = 1,nat
         fact = tpiba*(0.d0,-1.d0)*eigqts(na)
         mu = 3*(na-1)
         IF (abs(uact(mu+1))+abs(uact(mu+2))  &
                         +abs(uact(mu+3)).gt.1.0d-12) then
            nt=ityp(na)
            u1 = uact(mu+1)
            u2 = uact(mu+2)
            u3 = uact(mu+3)
            gu0 = xq0(1)*u1 +xq0(2)*u2+xq0(3)*u3
            IF (upf(nt)%nlcc) THEN
               DO ig = 1,ngm
                  gtau = eigts1(mill(1,ig),na)*   &
                         eigts2(mill(2,ig),na)*   &
                         eigts3(mill(3,ig),na)
                  gu = gu0+g(1,ig)*u1+g(2,ig)*u2+g(3,ig)*u3
                  aux(nl(ig))=aux(nl(ig))+drc(ig,nt)*gu*fact*gtau
               ENDDO
            ENDIF
         ENDIF
      ENDDO
      CALL invfft ('Dense', aux, dfftp)
      IF (.not.lsda) THEN
         DO ir=1,dfftp%nnr
            aux(ir) = aux(ir) * dmuxc(ir,1,1)
         ENDDO
      ELSE
         ! RM - I think it should be ik to be consistent with dvqpsi_us_only3
         !is=isk(ikk)
         is=isk(ik)
         DO ir=1,dfftp%nnr
            aux(ir) = aux(ir) * 0.5d0 *  &
                 (dmuxc(ir,is,1)+dmuxc(ir,is,2))
         ENDDO
      ENDIF
      CALL fwfft ('Dense', aux, dfftp)
      IF (doublegrid) THEN
         auxs(:) = (0.d0, 0.d0)
         DO ig=1,ngms
            auxs(nls(ig)) = aux(nl(ig))
         ENDDO
      ENDIF
      aux1(:) = aux1(:) + auxs(:)
   ENDIF
  !
  ! Now we compute dV_loc/dtau in real space
  !
  CALL invfft ('Smooth', aux1, dffts)
  DO ibnd = 1, nbnd
     DO ip = 1, npol
        aux2(:) = (0.d0, 0.d0)
        IF ( ip == 1 ) THEN
           DO ig = 1, npw
              aux2 (nls (igk (ig) ) ) = evc (ig, ibnd)
           ENDDO
        ELSE
           DO ig = 1, npw
              aux2 (nls (igk (ig) ) ) = evc (ig+npwx, ibnd)
           ENDDO
        ENDIF
        !
        !  This wavefunction is computed in real space
        !
        CALL invfft ('Wave', aux2, dffts)
        DO ir = 1, dffts%nnr
           aux2 (ir) = aux2 (ir) * aux1 (ir)
        ENDDO
        !
        ! and finally dV_loc/dtau * psi is transformed in reciprocal space
        !
        CALL fwfft ('Wave', aux2, dffts)
        IF ( ip == 1 ) THEN
           DO ig = 1, npwq
              dvpsi (ig, ibnd) = aux2 (nls (igkq (ig) ) )
           ENDDO
        ELSE
           DO ig = 1, npwq
              dvpsi (ig+npwx, ibnd) = aux2 (nls (igkq (ig) ) )
           ENDDO
        ENDIF
     ENDDO
  ENDDO
  ! 
  DEALLOCATE (aux2)
  DEALLOCATE (aux1)
  IF (nlcc_any.and.addnlcc) THEN
     DEALLOCATE (aux)
     IF (doublegrid) DEALLOCATE (auxs)
  ENDIF
  !
  !   We add the contribution of the nonlocal potential in the US form
  !   First a term similar to the KB case.
  !   Then a term due to the change of the D coefficients in the perturbat
  !
  CALL dvqpsi_us_only3 (ik, uact, xxk)
  CALL stop_clock ('dvqpsi_us')
  RETURN
END SUBROUTINE dvqpsi_us3