1! 2! Copyright (C) 2001-2016 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!----------------------------------------------------------------------- 9FUNCTION lr_dot(x,y) 10 !--------------------------------------------------------------------- 11 ! 12 ! This subroutine calculates a dot product of the conjugate 13 ! of a complex vector x and a complex vector y 14 ! (sums over the bands and k-points). 15 ! 16 ! Brent Walker, ICTP, 2004 17 ! Modified by Osman Baris Malcioglu, SISSA, 2009 18 ! Modified by Iurii Timrov, SISSA, 2013 (extension to EELS) 19 ! Modified by PG, 2020: replacement of zdotc with dot_product 20 ! 21 USE kinds, ONLY : dp 22 USE io_global, ONLY : stdout 23 USE klist, ONLY : nks, xk, wk, ngk 24 USE lsda_mod, ONLY : nspin 25 USE wvfct, ONLY : npwx,nbnd,wg 26 USE gvecw, ONLY : gcutw 27 USE control_flags, ONLY : gamma_only 28 USE gvect, ONLY : gstart, ngm, g 29 USE mp, ONLY : mp_sum 30 USE mp_global, ONLY : inter_pool_comm, intra_bgrp_comm 31 USE noncollin_module, ONLY : noncolin, npol 32 USE control_lr, ONLY : nbnd_occ 33 USE qpoint, ONLY : nksq 34 ! 35 IMPLICIT NONE 36 ! 37 COMPLEX(kind=dp) :: x(npwx*npol,nbnd,nksq), & 38 y(npwx*npol,nbnd,nksq) 39 COMPLEX(kind=dp) :: lr_dot 40 REAL(kind=dp) :: temp_gamma, degspin 41 INTEGER :: ibnd, ik 42 REAL(kind=dp), EXTERNAL :: DDOT 43 ! 44 CALL start_clock ('lr_dot') 45 ! 46 lr_dot = (0.0d0,0.0d0) 47 temp_gamma = 0.0d0 48 ! 49 IF (nspin==2) THEN 50 degspin = 1.0d0 51 ELSE 52 degspin = 2.0d0 53 ENDIF 54 ! 55 IF (gamma_only) THEN 56 ! 57 CALL lr_dot_gamma() 58 lr_dot = cmplx(temp_gamma,0.0d0,dp) 59 ! 60 ELSEIF (noncolin) THEN 61 ! 62 degspin = 1.0d0 63 CALL lr_dot_k_nc() 64 ! 65 ELSE 66 ! 67 CALL lr_dot_k() 68 ! 69 ENDIF 70 ! 71 lr_dot = lr_dot/degspin 72 ! 73 CALL stop_clock ('lr_dot') 74 ! 75 RETURN 76 ! 77CONTAINS 78 ! 79 SUBROUTINE lr_dot_gamma 80 ! 81 ! Optical case: gamma_only 82 ! Noncollinear case is not implemented 83 ! 84 DO ibnd=1,nbnd 85 ! 86 temp_gamma = temp_gamma + 2.D0*wg(ibnd,1)*DDOT(2*ngk(1),x(:,ibnd,1),1,y(:,ibnd,1),1) 87 ! 88 ! G=0 has been accounted twice, so we subtract one contribution. 89 ! 90 IF (gstart==2) temp_gamma = temp_gamma - wg(ibnd,1)*dble(x(1,ibnd,1))*dble(y(1,ibnd,1)) 91 ! 92 ENDDO 93 ! 94#if defined(__MPI) 95 CALL mp_sum(temp_gamma, intra_bgrp_comm) 96#endif 97 ! 98 RETURN 99 ! 100 END SUBROUTINE lr_dot_gamma 101 ! 102 SUBROUTINE lr_dot_k_nc 103 ! 104 ! Noncollinear case 105 ! 106 USE qpoint, ONLY : ikks, ikqs 107 ! 108 IMPLICIT NONE 109 INTEGER :: ios 110 INTEGER :: ik, & 111 ikk, & ! index of the point k 112 ikq, & ! index of the point k+q 113 npwq ! number of the plane-waves at point k+q 114 ! 115 DO ik = 1, nksq 116 ! 117 ikk = ikks(ik) 118 ikq = ikqs(ik) 119 npwq = ngk(ikq) 120 ! 121 DO ibnd = 1, nbnd_occ(ikk) 122 ! 123 lr_dot = lr_dot + wk(ikk) * dot_product(x(:,ibnd,ik),y(:,ibnd,ik)) 124 ! 125 ENDDO 126 ! 127 ENDDO 128 ! 129#if defined(__MPI) 130 CALL mp_sum(lr_dot, inter_pool_comm) 131 CALL mp_sum(lr_dot, intra_bgrp_comm) 132#endif 133 ! 134 RETURN 135 ! 136 END SUBROUTINE lr_dot_k_nc 137 ! 138 SUBROUTINE lr_dot_k 139 ! 140 ! collinear k point case 141 ! 142 USE qpoint, ONLY : ikks, ikqs 143 ! 144 IMPLICIT NONE 145 INTEGER :: ios 146 INTEGER :: ik, & 147 ikk, & ! index of the point k 148 ikq, & ! index of the point k+q 149 npwq ! number of the plane-waves at point k+q 150 ! 151 DO ik = 1, nksq 152 ! 153 ikk = ikks(ik) 154 ikq = ikqs(ik) 155 npwq = ngk(ikq) 156 ! 157 DO ibnd = 1, nbnd_occ(ikk) 158 ! 159 lr_dot = lr_dot + wk(ikk) * & 160 dot_product( x(1:npwq,ibnd,ik), y(1:npwq,ibnd,ik) ) 161 ! 162 ENDDO 163 ! 164 ENDDO 165 ! 166#if defined(__MPI) 167 CALL mp_sum(lr_dot, inter_pool_comm) 168 CALL mp_sum(lr_dot, intra_bgrp_comm) 169#endif 170 ! 171 RETURN 172 ! 173 END SUBROUTINE lr_dot_k 174 ! 175END FUNCTION lr_dot 176!----------------------------------------------------------------------- 177 178!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 179!! Debugging subroutines 180!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 181 182SUBROUTINE check_vector_gamma (x) 183 !---------------------------------------------------------------------------- 184 ! Checks the inner product for a given vector, and its imaginary and real component 185 ! input: evc 186 ! output : screen output 187 ! 188 USE kinds, ONLY : dp 189 USE mp_global, ONLY : inter_pool_comm, intra_bgrp_comm 190 USE mp, ONLY : mp_sum 191 USE klist , ONLY : ngk 192 USE gvect, ONLY : gstart 193 USE io_global, ONLY : stdout 194 ! 195 IMPLICIT NONE 196 !input/output 197 COMPLEX(kind=dp),INTENT(in) :: x(:) 198 ! 199 ! local variables 200 ! 201 REAL(kind=dp) :: temp_gamma 202 REAL(kind=dp), EXTERNAL :: DDOT 203 ! 204 temp_gamma = 2.D0*DDOT(2*ngk(1),x(:),1,x(:),1) 205 ! 206 IF (gstart==2) temp_gamma = temp_gamma - dble(x(1))*dble(x(1)) 207 ! 208#if defined(__MPI) 209 CALL mp_sum(temp_gamma, intra_bgrp_comm) 210#endif 211 ! 212 WRITE(stdout,'("<x> = ",E15.8)') temp_gamma 213 ! 214 RETURN 215 ! 216END SUBROUTINE check_vector_gamma 217 218SUBROUTINE check_vector_f (x) 219 !----------------------------------------------------------------------- 220 ! 221 ! Checks the inner product for a given vector, and its imaginary and real component 222 ! input: evc 223 ! output: screen output 224 ! 225 USE kinds, ONLY : dp 226 USE mp_global, ONLY : inter_pool_comm, intra_bgrp_comm 227 USE mp, ONLY : mp_sum 228 USE klist , ONLY : ngk 229 USE gvect, ONLY : gstart 230 USE io_global, ONLY : stdout 231 ! 232 IMPLICIT NONE 233 !input/output 234 COMPLEX(kind=dp),INTENT(in) :: x(:) 235 ! 236 ! local variables 237 ! 238 COMPLEX(kind=dp) :: temp_f 239 ! 240 temp_f = dot_product( x(1:ngk(1)), x(1:ngk(1)) ) 241 ! 242#if defined(__MPI) 243 CALL mp_sum(temp_f, intra_bgrp_comm) 244#endif 245 ! 246 WRITE(stdout,'("<x> = ",2E15.8,1X)') temp_f 247 ! 248 RETURN 249 ! 250END SUBROUTINE check_vector_f 251 252SUBROUTINE check_all_bands_gamma (x,sx,nbnd1,nbnd2) 253 !---------------------------------------------------------------------- 254 ! 255 ! Checks all bands of given KS states for orthogonality 256 ! input: evc and sevc 257 ! output : screen output 258 ! 259 USE kinds, ONLY : dp 260 USE mp_global, ONLY : inter_pool_comm, intra_bgrp_comm 261 USE mp, ONLY : mp_sum 262 USE klist , ONLY : ngk 263 USE io_global, ONLY : stdout 264 USE gvect, ONLY : gstart 265 ! 266 IMPLICIT NONE 267 !input/output 268 INTEGER, INTENT(in) :: nbnd1,nbnd2 !Total number of bands for x and sx 269 COMPLEX(kind=dp),INTENT(in) :: x(:,:), sx(:,:) 270 ! 271 ! local variables 272 ! 273 INTEGER :: ibnd, jbnd 274 REAL(kind=dp) :: temp_gamma 275 REAL(kind=dp), EXTERNAL :: DDOT 276 ! 277 DO ibnd=1,nbnd1 278 DO jbnd=ibnd,nbnd2 279 ! 280 temp_gamma = 2.D0*DDOT(2*ngk(1),x(:,ibnd),1,sx(:,jbnd),1) 281 ! 282 IF (gstart==2) temp_gamma = temp_gamma - dble(x(1,ibnd))*dble(sx(1,jbnd)) 283 ! 284#if defined(__MPI) 285 CALL mp_sum(temp_gamma, intra_bgrp_comm) 286#endif 287 ! 288 WRITE(stdout,'("<x,",I02,"|S|x,",I02,"> =",E15.8)') ibnd,jbnd,temp_gamma 289 ENDDO 290 ENDDO 291 ! 292 RETURN 293 ! 294END SUBROUTINE check_all_bands_gamma 295 296SUBROUTINE check_density_gamma (rx,nbnd) 297 !--------------------------------------------------------------------------------- 298 ! 299 ! Checks the contirbution of a given function transformed into real space 300 ! input: revc 301 ! output : screen output 302 ! 303 USE kinds, ONLY : dp 304 USE mp_global, ONLY : inter_pool_comm, intra_bgrp_comm 305 USE mp, ONLY : mp_sum 306 USE wvfct, ONLY : wg 307 USE fft_base, ONLY : dfftp 308 USE io_global, ONLY : stdout 309 USE cell_base, ONLY : omega 310 ! 311 IMPLICIT NONE 312 !input/output 313 INTEGER, INTENT(in) :: nbnd !Total number of bands for x and sx 314 COMPLEX(kind=dp),INTENT(in) :: rx(:,:) 315 ! 316 ! local variables 317 ! 318 INTEGER :: ibnd 319 REAL(kind=dp) :: temp_gamma,w1,w2 320 ! 321 DO ibnd=1,nbnd,2 322 w1 = wg(ibnd,1)/omega 323 ! 324 IF (ibnd<nbnd) THEN 325 w2 = wg(ibnd+1,1)/omega 326 ELSE 327 w2 = w1 328 ENDIF 329 temp_gamma = sum(w1*dble(rx(1:dfftp%nnr,ibnd))*dble(rx(1:dfftp%nnr,ibnd))& 330 + w2*aimag(rx(1:dfftp%nnr,ibnd))*aimag(rx(1:dfftp%nnr,ibnd))) 331#if defined(__MPI) 332 CALL mp_sum(temp_gamma, intra_bgrp_comm) 333#endif 334 WRITE(stdout,'("Contribution of bands ",I02," and ",I02," to total density",E15.8)') ibnd,ibnd+1,temp_gamma 335 ! 336 ENDDO 337 ! 338 RETURN 339 ! 340END SUBROUTINE check_density_gamma 341