dvanqq2.f90 9.42 KB
Newer Older
sponce's avatar
sponce committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
  !
  ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino  
  !
  ! Copyright (C) 2001 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/dvanqq (QE)
  ! 
  !----------------------------------------------------------------------
  SUBROUTINE dvanqq2
  !----------------------------------------------------------------------
15 16 17 18 19 20 21 22 23 24 25
  !!
  !! New
  !! This routine calculates two integrals of the Q functions and
  !! its derivatives with c V_loc and V_eff which are used
  !! to compute term dV_bare/dtau * psi  in addusdvqpsi.
  !! The result is stored in int1,int2. The routine is called
  !! for each q in nqc. 
  !! 
  !! RM - Nov/Dec 2014 
  !! Imported the noncolinear case implemented by xlzhang
  !!
sponce's avatar
sponce committed
26 27 28
  !
  USE ions_base,        ONLY : nat, ityp, ntyp => nsp
  USE io_global,        ONLY : stdout
29 30
  USE pwcom,            ONLY : lspinorb, nspin, domag
  USE cell_base,        ONLY : tpiba2, omega, tpiba
31
  USE gvect,            ONLY : ngm, gg, nl, g, eigts1, eigts2, eigts3, mill
sponce's avatar
sponce committed
32 33 34
  USE scf,              ONLY : v, vltot
  USE noncollin_module, ONLY : noncolin
  USE kinds,            ONLY : DP
sponce's avatar
sponce committed
35
  USE phcom,            ONLY : int1, int2, int4, int5, &
sponce's avatar
sponce committed
36 37 38
                               int5_so, vlocq, int4_nc
  USE qpoint,           ONLY : xq, eigqts
  USE uspp_param,       ONLY : upf, lmaxq, nh
sponce's avatar
sponce committed
39
  USE mp_global,        ONLY : intra_pool_comm
sponce's avatar
sponce committed
40 41 42 43 44 45 46 47 48
  USE mp,               ONLY : mp_sum
  USE uspp,             ONLY : okvan
  USE fft_base,         ONLY : dfftp
  USE fft_interfaces,   ONLY : fwfft

  implicit none
  !
  !   And the local variables
  !
49
  INTEGER :: nt, na, nb, ig, nta, ntb, ir, ih, jh, ijh, ipol, jpol, is, nspin0
sponce's avatar
sponce committed
50
  !
51
  real(kind=DP), ALLOCATABLE :: qmod (:), qmodg (:), qpg (:,:), &
sponce's avatar
sponce committed
52 53 54 55 56 57
       ylmkq (:,:), ylmk0 (:,:)
  ! the modulus of q+G
  ! the modulus of G
  ! the  q+G vectors
  ! the spherical harmonics

58 59
  COMPLEX(kind=DP) :: fact, fact1, ZDOTC
  COMPLEX(kind=DP), ALLOCATABLE :: aux1 (:), aux2 (:),&
sponce's avatar
sponce committed
60 61
       aux3 (:), aux5 (:), veff (:,:), sk(:)
  ! work space
62
  complex(kind=DP), ALLOCATABLE, TARGET :: qgm(:)
sponce's avatar
sponce committed
63
  ! the augmentation function at G
64
  complex(kind=DP), POINTER :: qgmq (:)
sponce's avatar
sponce committed
65 66
  ! the augmentation function at q+G
  logical :: exst
sponce's avatar
sponce committed
67
  ! 
sponce's avatar
sponce committed
68
  IF (.not.okvan) RETURN
sponce's avatar
sponce committed
69
  ! 
sponce's avatar
sponce committed
70 71
  nspin0=nspin
  IF (nspin==4.and..not.domag) nspin0=1
sponce's avatar
sponce committed
72
  ! 
sponce's avatar
sponce committed
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
  CALL start_clock ('dvanqq2')
  int1(:,:,:,:,:) = (0.d0, 0.d0)
  int2(:,:,:,:,:) = (0.d0, 0.d0)
  int4(:,:,:,:,:) = (0.d0, 0.d0)
  int5(:,:,:,:,:) = (0.d0, 0.d0)
  ALLOCATE (sk  (  ngm))    
  ALLOCATE (aux1(  ngm))    
  ALLOCATE (aux2(  ngm))    
  ALLOCATE (aux3(  ngm))    
  ALLOCATE (aux5(  ngm))    
  ALLOCATE (qmodg( ngm))    
  ALLOCATE (ylmk0( ngm , lmaxq * lmaxq))    
  ALLOCATE (qgm  ( ngm))    
  ALLOCATE (ylmkq(ngm , lmaxq * lmaxq))    
  ALLOCATE (qmod( ngm))    
  ALLOCATE (qgmq( ngm))    
  !
  !     compute spherical harmonics
  !
  CALL ylmr2 (lmaxq * lmaxq, ngm, g, gg, ylmk0)
  DO ig = 1, ngm
     qmodg (ig) = sqrt (gg (ig) )
  ENDDO
  ALLOCATE (qpg (3, ngm))    
  CALL setqmod (ngm, xq, g, qmod, qpg)
  CALL ylmr2 (lmaxq * lmaxq, ngm, qpg, qmod, ylmkq)
  DEALLOCATE (qpg)
  DO ig = 1, ngm
     qmod (ig) = sqrt (qmod (ig) )
  ENDDO
  !   we start by computing the FT of the effective potential
  !
  ALLOCATE (veff ( dfftp%nnr , nspin))    
  DO is = 1, nspin
     IF (nspin.ne.4.or.is==1) THEN
        DO ir = 1, dfftp%nnr
109
           veff (ir, is) = CMPLX (vltot (ir) + v%of_r (ir, is), 0.d0, kind=DP)
sponce's avatar
sponce committed
110 111 112
        ENDDO
     ELSE
        DO ir = 1, dfftp%nnr
113
           veff (ir, is) = CMPLX (v%of_r (ir, is), 0.d0, kind=DP)
sponce's avatar
sponce committed
114 115 116 117 118 119 120
        ENDDO
     ENDIF
     CALL fwfft ('Dense', veff(:,is), dfftp)
  ENDDO
  !
  !     We compute here two of the three integrals needed in the phonon
  !
121
  fact1 = CMPLX (0.d0, - tpiba * omega, kind=DP)
sponce's avatar
sponce committed
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 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
  !
  DO ntb = 1, ntyp
     IF (upf(ntb)%tvanp ) THEN
        ijh = 0
        DO ih = 1, nh (ntb)
           DO jh = ih, nh (ntb)
              ijh = ijh + 1
              !
              !    compute the augmentation function
              !
              CALL qvan2 (ngm, ih, jh, ntb, qmodg, qgm, ylmk0)
              CALL qvan2 (ngm, ih, jh, ntb, qmod, qgmq, ylmkq)
              !
              !     NB: for this integral the moving atom and the atom of Q
              !     do not necessarily coincide
              !
              DO nb = 1, nat
                 IF (ityp (nb) == ntb) THEN
                    DO ig = 1, ngm
                       aux1 (ig) = qgmq (ig) * eigts1 (mill(1,ig), nb) &
                                             * eigts2 (mill(2,ig), nb) &
                                             * eigts3 (mill(3,ig), nb)
                    ENDDO
                    DO na = 1, nat
                       fact = eigqts (na) * CONJG(eigqts (nb) )
                       !
                       !    nb is the atom of the augmentation function
                       !
                       nta = ityp (na)
                       DO ig=1, ngm
                          sk(ig)=vlocq(ig,nta) * eigts1(mill(1,ig), na) &
                                               * eigts2(mill(2,ig), na) &
                                               * eigts3(mill(3,ig), na) 
                       ENDDO
                       DO ipol = 1, 3
                          DO ig=1, ngm
                            aux5(ig)= sk(ig) * (g (ipol, ig) + xq (ipol) )
                          ENDDO
                          int2 (ih, jh, ipol, na, nb) = fact * fact1 * &
                                ZDOTC (ngm, aux1, 1, aux5, 1)
                          DO jpol = 1, 3
                             IF (jpol >= ipol) THEN
                                DO ig = 1, ngm
                                   aux3 (ig) = aux5 (ig) * &
                                               (g (jpol, ig) + xq (jpol) )
                                ENDDO
                                int5 (ijh, ipol, jpol, na, nb) = &
                                     CONJG(fact) * tpiba2 * omega * &
                                     ZDOTC (ngm, aux3, 1, aux1, 1)
                             ELSE
                                int5 (ijh, ipol, jpol, na, nb) = &
                                     int5 (ijh, jpol, ipol, na, nb)
                             ENDIF
                          ENDDO
                       ENDDO
                    ENDDO
                    DO ig = 1, ngm
                       aux1 (ig) = qgm (ig) * eigts1 (mill(1,ig), nb) &
                                            * eigts2 (mill(2,ig), nb) &
                                            * eigts3 (mill(3,ig), nb)
                    ENDDO
                    DO is = 1, nspin0
                       DO ipol = 1, 3
                          DO ig = 1, ngm
                             aux2 (ig) = veff (nl (ig), is) * g (ipol, ig)
                          ENDDO
                          int1 (ih, jh, ipol, nb, is) = - fact1 * &
                               ZDOTC (ngm, aux1, 1, aux2, 1)
                          DO jpol = 1, 3
                             IF (jpol >= ipol) THEN
                                DO ig = 1, ngm
                                   aux3 (ig) = aux2 (ig) * g (jpol, ig)
                                ENDDO
                                int4 (ijh, ipol, jpol, nb, is) = - tpiba2 * &
                                     omega * ZDOTC (ngm, aux3, 1, aux1, 1)
                             ELSE
                                int4 (ijh, ipol, jpol, nb, is) = &
                                     int4 (ijh, jpol, ipol, nb, is)
                             ENDIF
                          ENDDO
                       ENDDO
                    ENDDO
                 ENDIF
              ENDDO
           ENDDO
        ENDDO
        DO ih = 1, nh (ntb)
           DO jh = ih + 1, nh (ntb)
              !
              !    We use the symmetry properties of the integral factor
              !
              DO nb = 1, nat
                 IF (ityp (nb) == ntb) THEN
                    DO ipol = 1, 3
                       DO is = 1, nspin
                          int1(jh,ih,ipol,nb,is) = int1(ih,jh,ipol,nb,is)
                       ENDDO
                       DO na = 1, nat
                          int2(jh,ih,ipol,na,nb) = int2(ih,jh,ipol,na,nb)
                       ENDDO
                    ENDDO
                 ENDIF
              ENDDO
           ENDDO
        ENDDO
     ENDIF
  ENDDO
  CALL mp_sum(  int1, intra_pool_comm )
  CALL mp_sum(  int2, intra_pool_comm )
  call mp_sum(  int4, intra_pool_comm )
  call mp_sum(  int5, intra_pool_comm )
  IF (noncolin) THEN
     CALL set_int12_nc(0)
     int4_nc = (0.d0, 0.d0)
     IF (lspinorb) int5_so = (0.d0, 0.d0)
     DO nt = 1, ntyp
        IF ( upf(nt)%tvanp ) THEN
           DO na = 1, nat
              IF (ityp(na)==nt) THEN
                 IF (upf(nt)%has_so)  THEN
                    CALL transform_int4_so(int4,na)
                    CALL transform_int5_so(int5,na)
                 ELSE
                    CALL transform_int4_nc(int4,na)
                    IF (lspinorb) CALL transform_int5_nc(int5,na)
                 ENDIF
              ENDIF
           ENDDO
        ENDIF
     ENDDO
  ENDIF
  !
  DEALLOCATE (veff)
  DEALLOCATE (qgmq)
  DEALLOCATE (qmod)
  DEALLOCATE (ylmkq)
  DEALLOCATE (qgm)
  DEALLOCATE (ylmk0)
  DEALLOCATE (qmodg)
  DEALLOCATE (aux5)
  DEALLOCATE (aux3)
  DEALLOCATE (aux2)
  DEALLOCATE (aux1)
  DEALLOCATE (sk)

  CALL stop_clock ('dvanqq2')
  RETURN
END SUBROUTINE dvanqq2