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!----------------------------------------------------------------------- 9SUBROUTINE lr_ortho(dvpsi, evq, ikk, ikq, sevc, inverse) 10 !--------------------------------------------------------------------- 11 ! 12 ! Inspired by LR_Modules/orthogonalize.f90 13 ! This routine orthogonalizes dvpsi to the valence states. 14 ! It should be quite general. It works for metals and insulators, with 15 ! NC as well as with US PP, both SR or FR. 16 ! 17 ! If inverse=.true. then apply P_c = 1 - |evq><sevc| = 1 - |psi><psi|S 18 ! If inverse=.false. then apply P_c^+ = 1 - |sevc><evq| = 1 - S|psi><psi| 19 ! 20 ! See definitions of P_c and P_c^+ eq.(29) in 21 ! B. Walker and R. Gebauer, J. Chem. Phys. 127, 164106 (2007) 22 ! 23 ! NB: IN/OUT is dvpsi ; sevc is used as a workspace 24 ! 25 ! evc0 -> evq 26 ! sevc0 -> sevc 27 ! 28 ! Note: S^{-1} P_c^+(k) = P_c(k) S^{-1} 29 ! 30 ! This subroutine is a modified version of the the subroutine 31 ! orthogonalize.f90. This subroutine should be removed in the future, 32 ! and orthogonalize.f90 should be used instead. Currently this 33 ! subroutine is used only in several places in lr_dav_routines.f90. 34 ! TODO: Check if lr_ortho can be replaced by orthogonalize in 35 ! lr_dav_routines. 36 ! 37 ! Modified by Osman Baris Malcioglu (2009) 38 ! Modified by Iurii Timrov (2013) 39 ! 40 USE kinds, ONLY : DP 41 use gvect, only : gstart 42 USE klist, ONLY : ngk, lgauss, degauss, ngauss 43 USE noncollin_module, ONLY : noncolin, npol 44 USE wvfct, ONLY : npwx, nbnd, et 45 USE ener, ONLY : ef 46 USE control_lr, ONLY : alpha_pv, nbnd_occ 47 USE uspp, ONLY : vkb, okvan 48 USE mp_global, ONLY : intra_bgrp_comm 49 USE mp, ONLY : mp_sum 50 use lr_variables, ONLY : lr_verbosity 51 use control_flags, only : gamma_only 52 USE io_global, ONLY : stdout 53 ! 54 IMPLICIT NONE 55 INTEGER, INTENT(in) :: ikk, ikq 56 ! the index of the k and k+q points 57 ! (in the optical case ikq = ikk) 58 COMPLEX(DP), INTENT(in) :: evq(npwx*npol,nbnd) 59 COMPLEX(DP), INTENT(inout) :: dvpsi(npwx*npol,nbnd) 60 COMPLEX(DP), INTENT(in) :: sevc(npwx*npol,nbnd) 61 ! sevc is a work space allocated by the calling routine (was called dpsi) 62 LOGICAL, INTENT(in):: inverse 63 LOGICAL:: inverse_mode 64 ! 65 CALL start_clock ('lr_ortho') 66 ! 67 IF (lr_verbosity > 5) WRITE(stdout,'("<lr_ortho>")') 68 ! 69 !if (.not. present(inverse)) then 70 ! inverse_mode=.false. 71 !else 72 inverse_mode = inverse 73 !endif 74 ! 75 IF (gamma_only) THEN 76 ! 77 CALL lr_ortho_gamma() 78 ! 79 ELSEIF (noncolin) THEN 80 ! 81 CALL lr_ortho_noncolin() 82 ! 83 ELSE 84 ! 85 CALL lr_ortho_k() 86 ! 87 ENDIF 88 ! 89 CALL stop_clock ('lr_ortho') 90 ! 91 RETURN 92 ! 93CONTAINS 94 95SUBROUTINE lr_ortho_k() 96 ! 97 ! General K points algorithm 98 ! 99 IMPLICIT NONE 100 101 COMPLEX(DP), ALLOCATABLE :: ps(:,:) 102 INTEGER :: ibnd, jbnd, nbnd_eff 103 REAL(DP) :: wg1, w0g, wgp, wwg, deltae, theta 104 REAL(DP), EXTERNAL :: w0gauss, wgauss 105 106 ALLOCATE(ps(nbnd,nbnd)) 107 ! 108 IF (lgauss) THEN 109 ! 110 ! metallic case 111 ! 112 ps = (0.d0, 0.d0) 113 IF (inverse_mode) THEN 114 CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ(ikk), ngk(ikk), (1.d0,0.d0), & 115 sevc, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd ) 116 ELSE 117 CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ(ikk), ngk(ikk), (1.d0,0.d0), & 118 evq, npwx, dvpsi, npwx, (0.d0,0.d0), ps, nbnd ) 119 ENDIF 120 ! 121 DO ibnd = 1, nbnd_occ (ikk) 122 wg1 = wgauss ((ef-et(ibnd,ikk)) / degauss, ngauss) 123 w0g = w0gauss((ef-et(ibnd,ikk)) / degauss, ngauss) / degauss 124 DO jbnd = 1, nbnd 125 wgp = wgauss ( (ef - et (jbnd, ikq) ) / degauss, ngauss) 126 deltae = et (jbnd, ikq) - et (ibnd, ikk) 127 theta = wgauss (deltae / degauss, 0) 128 wwg = wg1 * (1.d0 - theta) + wgp * theta 129 IF (jbnd <= nbnd_occ (ikq) ) THEN 130 IF (abs (deltae) > 1.0d-5) THEN 131 wwg = wwg + alpha_pv * theta * (wgp - wg1) / deltae 132 ELSE 133 ! 134 ! if the two energies are too close takes the limit 135 ! of the 0/0 ratio 136 ! 137 wwg = wwg - alpha_pv * theta * w0g 138 ENDIF 139 ENDIF 140 ! 141 ps(jbnd,ibnd) = wwg * ps(jbnd,ibnd) 142 ! 143 ENDDO 144 CALL DSCAL (2*ngk(ikk), wg1, dvpsi(1,ibnd), 1) 145 ENDDO 146 ! 147 nbnd_eff = nbnd 148 ! 149 ELSE 150 ! 151 ! insulators 152 ! 153 ps = (0.d0, 0.d0) 154 ! 155 IF (inverse_mode) THEN 156 ! 157 ! ps = <sevc|dvpsi> 158 ! 159 CALL ZGEMM( 'C', 'N', nbnd_occ(ikq), nbnd_occ (ikk), ngk(ikk), & 160 (1.d0,0.d0), sevc, npwx, dvpsi, npwx, & 161 (0.d0,0.d0), ps, nbnd ) 162 ! 163 ELSE 164 ! 165 ! ps = <evq|dvpsi> 166 ! 167 CALL ZGEMM( 'C', 'N', nbnd_occ(ikq), nbnd_occ (ikk), ngk(ikk), & 168 (1.d0,0.d0), evq, npwx, dvpsi, npwx, & 169 (0.d0,0.d0), ps, nbnd ) 170 ENDIF 171 ! 172 nbnd_eff = nbnd_occ(ikk) 173 ! 174 ENDIF 175 ! 176#if defined(__MPI) 177 CALL mp_sum(ps(:,1:nbnd_eff),intra_bgrp_comm) 178#endif 179 ! 180 IF (lgauss) THEN 181 ! 182 ! Metallic case 183 ! 184 if (inverse_mode) then 185 CALL ZGEMM( 'N', 'N', ngk(ikk), nbnd_occ(ikk), nbnd, & 186 (-1.d0,0.d0), evq, npwx, ps, nbnd, (1.0d0,0.d0), & 187 dvpsi, npwx ) 188 else 189 CALL ZGEMM( 'N', 'N', ngk(ikk), nbnd_occ(ikk), nbnd, & 190 (-1.d0,0.d0), sevc, npwx, ps, nbnd, (1.0d0,0.d0), & 191 dvpsi, npwx ) 192 endif 193 ! 194 ELSE 195 ! 196 ! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) in an insulator 197 ! 198 if (inverse_mode) then 199 ! 200 ! |dvspi> = |dvpsi> - |evq><sevc|dvpsi> 201 ! 202 CALL ZGEMM( 'N', 'N', ngk(ikk), nbnd_occ(ikk), nbnd_occ(ikk), & 203 (-1.d0,0.d0), evq, npwx, ps, nbnd, (1.0d0,0.d0), & 204 dvpsi, npwx ) 205 ! 206 else 207 ! 208 ! |dvspi> = |dvpsi> - |sevc><evq|dvpsi> 209 ! 210 CALL ZGEMM( 'N', 'N', ngk(ikk), nbnd_occ(ikk), nbnd_occ(ikk), & 211 (-1.d0,0.d0), sevc, npwx, ps, nbnd, (1.0d0,0.d0), & 212 dvpsi, npwx ) 213 ! 214 endif 215 ! 216 ENDIF 217 ! 218 DEALLOCATE(ps) 219 ! 220END SUBROUTINE lr_ortho_k 221 222SUBROUTINE lr_ortho_gamma() 223 ! 224 ! gamma_only algorithm 225 ! 226 IMPLICIT NONE 227 228 COMPLEX(DP), ALLOCATABLE :: ps_c(:,:) 229 REAL(DP), ALLOCATABLE :: ps(:,:) 230 INTEGER :: ibnd, jbnd, nbnd_eff 231 REAL(DP) :: wg1, w0g, wgp, wwg, deltae, theta 232 REAL(DP), EXTERNAL :: w0gauss, wgauss 233 234 ALLOCATE(ps(nbnd,nbnd)) 235 ALLOCATE(ps_c(nbnd,nbnd)) 236 ! 237 IF (lgauss) THEN 238 ! 239 ! Metals 240 ! 241 CALL errore ('lr_ortho', "degauss with gamma point algorithm is not allowed",1) 242 ! 243 ELSE 244 ! 245 ! Insulators 246 ! 247 ps = 0.d0 248 ! 249 IF (inverse_mode) THEN 250 ! 251 ! ps = 2 * <sevc|dvpsi> 252 ! 253 CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*ngk(1), & 254 2.d0, sevc, 2*npwx, dvpsi, 2*npwx, 0.d0, ps, nbnd ) 255 ! 256 ELSE 257 ! 258 ! ps = 2 * <evq|dvpsi> 259 ! 260 CALL DGEMM( 'C', 'N', nbnd, nbnd ,2*ngk(1), & 261 2.d0, evq, 2*npwx, dvpsi, 2*npwx, 0.d0, ps, nbnd ) 262 ! 263 ENDIF 264 ! 265 nbnd_eff = nbnd 266 ! 267 ! G=0 has been accounted twice, therefore we subtruct one contribution. 268 ! 269 IF (gstart == 2) THEN 270 ! 271 IF (inverse_mode) THEN 272 ! 273 ! ps = ps - <sevc|dvpsi> 274 ! 275 CALL DGER( nbnd, nbnd, -1.D0, sevc, 2*npwx, dvpsi, 2*npwx, ps, nbnd) 276 ! 277 ELSE 278 ! 279 ! ps = ps - <evq|dvpsi> 280 ! 281 CALL DGER( nbnd, nbnd, -1.D0, evq, 2*npwx, dvpsi, 2*npwx, ps, nbnd ) 282 ! 283 ENDIF 284 ! 285 ENDIF 286 ! 287 ENDIF 288 ! 289#if defined(__MPI) 290 CALL mp_sum(ps(:,:),intra_bgrp_comm) 291#endif 292 ! 293 ! In the original code dpsi was used as a storage for sevc, since in 294 ! tddfpt we have it stored in memory as sevc0 this part is obsolote 295 ! 296 ps_c = cmplx(ps, 0.d0, dp) 297 ! 298 IF (lgauss) THEN 299 ! 300 ! Metals 301 ! 302 CALL errore ('lr_ortho', "degauss with gamma point algorithm is not allowed",1) 303 ! 304 ELSE 305 ! 306 ! Insulators: note that nbnd_occ(ikk) = nbnd_occ(ikq) 307 ! 308 IF (inverse_mode) THEN 309 ! 310 ! |dvpsi> = |dvpsi> - |evq><sevc|dvpsi> 311 ! 312 CALL ZGEMM( 'N', 'N', ngk(1), nbnd, nbnd, & 313 (-1.d0,0.d0), evq, npwx, ps_c, nbnd, (1.0d0,0.d0), dvpsi, npwx) 314 ! 315 ELSE 316 ! 317 ! |dvpsi> = |dvpsi> - |sevc><evq|dvpsi> 318 ! 319 CALL ZGEMM( 'N', 'N', ngk(1), nbnd, nbnd, & 320 (-1.d0,0.d0), sevc, npwx, ps_c, nbnd, (1.0d0,0.d0), dvpsi, npwx ) 321 ! 322 ENDIF 323 ! 324 ENDIF 325 ! 326 DEALLOCATE(ps) 327 DEALLOCATE(ps_c) 328 ! 329 RETURN 330 ! 331END SUBROUTINE lr_ortho_gamma 332 333SUBROUTINE lr_ortho_noncolin() 334 ! 335 ! Noncollinear case 336 ! Warning: the inverse mode is not implemented! 337 ! 338 IMPLICIT NONE 339 340 COMPLEX(DP), ALLOCATABLE :: ps(:,:) 341 INTEGER :: ibnd, jbnd, nbnd_eff 342 REAL(DP) :: wg1, w0g, wgp, wwg, deltae, theta 343 REAL(DP), EXTERNAL :: w0gauss, wgauss 344 345 IF (inverse_mode) CALL errore ('lr_ortho', "The inverse mode is not implemented!",1) 346 ! 347 ALLOCATE(ps(nbnd,nbnd)) 348 ! 349 IF (lgauss) THEN 350 ! 351 ! metallic case 352 ! 353 ps = (0.d0, 0.d0) 354 ! 355 CALL ZGEMM( 'C', 'N', nbnd, nbnd_occ (ikk), npwx*npol, (1.d0,0.d0), & 356 evq, npwx*npol, dvpsi, npwx*npol, (0.d0,0.d0), ps, nbnd ) 357 ! 358 DO ibnd = 1, nbnd_occ (ikk) 359 ! 360 wg1 = wgauss ((ef-et(ibnd,ikk)) / degauss, ngauss) 361 w0g = w0gauss((ef-et(ibnd,ikk)) / degauss, ngauss) / degauss 362 ! 363 DO jbnd = 1, nbnd 364 ! 365 wgp = wgauss ( (ef - et (jbnd, ikq) ) / degauss, ngauss) 366 deltae = et (jbnd, ikq) - et (ibnd, ikk) 367 theta = wgauss (deltae / degauss, 0) 368 wwg = wg1 * (1.d0 - theta) + wgp * theta 369 ! 370 IF (jbnd <= nbnd_occ (ikq) ) THEN 371 IF (abs (deltae) > 1.0d-5) THEN 372 wwg = wwg + alpha_pv * theta * (wgp - wg1) / deltae 373 ELSE 374 ! 375 ! if the two energies are too close takes the limit 376 ! of the 0/0 ratio 377 ! 378 wwg = wwg - alpha_pv * theta * w0g 379 ENDIF 380 ENDIF 381 ! 382 ps(jbnd,ibnd) = wwg * ps(jbnd,ibnd) 383 ! 384 ENDDO 385 CALL DSCAL (2*npwx*npol, wg1, dvpsi(1,ibnd), 1) 386 ENDDO 387 ! 388 nbnd_eff = nbnd 389 ! 390 ELSE 391 ! 392 ! insulators 393 ! 394 ps = (0.d0, 0.d0) 395 ! 396 ! ps = <evq|dvpsi> 397 ! 398 CALL ZGEMM( 'C', 'N',nbnd_occ(ikq), nbnd_occ(ikk), npwx*npol, & 399 (1.d0,0.d0), evq, npwx*npol, dvpsi, npwx*npol, (0.d0,0.d0), ps, nbnd ) 400 ! 401 nbnd_eff = nbnd_occ(ikk) 402 ! 403 ENDIF 404 ! 405#if defined(__MPI) 406 CALL mp_sum(ps(:,1:nbnd_eff),intra_bgrp_comm) 407#endif 408 ! 409 ! In the original dpsi was used as a storage for sevc, since in 410 ! tddfpt we have it stored in memory as sevc0 this part is obsolote. 411 ! 412 IF (lgauss) THEN 413 ! 414 ! metallic case 415 ! 416 ! |dvpsi> = |dvpsi> - |sevc><evq|dvpsi> 417 ! 418 CALL ZGEMM( 'N', 'N', npwx*npol, nbnd_occ(ikk), nbnd, & 419 (-1.d0,0.d0), sevc, npwx*npol, ps, nbnd, (1.0d0,0.d0), dvpsi,npwx*npol ) 420 ! 421 ELSE 422 ! 423 ! Insulators: note that nbnd_occ(ikk)=nbnd_occ(ikq) 424 ! 425 ! |dvpsi> = |dvpsi> - |sevc><evq|dvpsi> 426 ! 427 CALL ZGEMM( 'N', 'N', npwx*npol, nbnd_occ(ikk), nbnd_occ(ikk), & 428 (-1.d0,0.d0),sevc,npwx*npol,ps,nbnd,(1.0d0,0.d0), dvpsi,npwx*npol ) 429 ! 430 ENDIF 431 ! 432 DEALLOCATE(ps) 433 ! 434 RETURN 435 ! 436END SUBROUTINE lr_ortho_noncolin 437 438END SUBROUTINE lr_ortho 439!----------------------------------------------------------------------- 440