1 ! 2 ! Copyright (C) 2016-2019 Samuel Ponce', Roxana Margine, Feliciano Giustino 3 ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino 4 ! Copyright (C) 2007-2009 Jesse Noffsinger, Brad Malone, Feliciano Giustino 5 ! 6 ! This file is distributed under the terms of the GNU General Public 7 ! License. See the file `LICENSE' in the root directory of the 8 ! present distribution, or http://www.gnu.org/copyleft.gpl.txt . 9 ! 10 !---------------------------------------------------------------------- 11 MODULE pw2wan2epw 12 !---------------------------------------------------------------------- 13 !! 14 !! This module contains routine to go from PW results to Wannier90 to EPW 15 !! 16 IMPLICIT NONE 17 ! 18 CONTAINS 19 ! 20 !------------------------------------------------------------------------ 21 SUBROUTINE pw2wan90epw 22 !------------------------------------------------------------------------ 23 !! 24 !! This is the interface to the Wannier90 code: see http://www.wannier.org 25 !! Adapted from the code PP/pw2wannier - Quantum-ESPRESSO group 26 !! 27 !! 10/2008 Parellel computation of Amn and Mmn 28 !! 12/2008 Added phase setting of overlap matrix elements 29 !! 02/2009 works with standard nkc1*nkc2*nkc3 grids 30 !! 12/2009 works with USPP 31 !! 12/2014 RM: Imported the noncolinear case implemented by xlzhang 32 !! 06/2016 SP: Debug of SOC + print/reading of nnkp file 33 !! 11/2019 RM: Imported and adapted the SCDM method from pw2wannier 34 !! 35 !------------------------------------------------------------------------ 36 USE io_global, ONLY : stdout 37 USE klist, ONLY : nkstot 38 USE io_files, ONLY : prefix 39 USE epwcom, ONLY : write_wfn, scdm_proj, scdm_entanglement, & 40 scdm_sigma, wannier_plot 41 USE wannierEPW, ONLY : seedname2, ispinw, ikstart, ikstop, iknum 42 USE constants_epw, ONLY : zero, one 43 USE noncollin_module, ONLY : noncolin 44 ! 45 IMPLICIT NONE 46 ! 47 CHARACTER(LEN = 4) :: spin_component 48 !! Determine the spin case 49 CHARACTER(LEN = 256) :: outdir 50 !! Name of the output directory 51 ! 52 outdir = './' 53 seedname2 = prefix 54 spin_component = 'none' 55 ! 56 IF (scdm_proj) THEN 57 IF ((TRIM(scdm_entanglement) /= 'isolated') .AND. & 58 (TRIM(scdm_entanglement) /= 'erfc') .AND. & 59 (TRIM(scdm_entanglement) /= 'gaussian')) THEN 60 CALL errore('pw2wan90epw', & 61 'Can not recognize the choice for scdm_entanglement. ' & 62 //'Valid options are: isolated, erfc and gaussian') 63 ENDIF 64 ENDIF 65 IF (scdm_sigma <= zero) & 66 CALL errore('pw2wan90epw', 'Sigma in the SCDM method must be positive.') 67 ! 68 WRITE(stdout, *) 69 SELECT CASE(TRIM(spin_component)) 70 CASE ('up') 71 WRITE(stdout, *) ' Spin CASE ( up )' 72 ispinw = 1 73 ikstart = 1 74 ikstop = nkstot / 2 75 iknum = nkstot / 2 76 CASE ('down') 77 WRITE(stdout, *) ' Spin CASE ( down )' 78 ispinw = 2 79 ikstart = nkstot / 2 + 1 80 ikstop = nkstot 81 iknum = nkstot / 2 82 CASE DEFAULT 83 IF (noncolin) THEN 84 WRITE(stdout, *) ' Spin CASE ( non-collinear )' 85 ELSE 86 WRITE(stdout, *) ' Spin CASE ( default = unpolarized )' 87 ENDIF 88 ispinw = 0 89 ikstart = 1 90 ikstop = nkstot 91 iknum = nkstot 92 END SELECT 93 ! 94 WRITE(stdout, *) 95 WRITE(stdout, *) ' Initializing Wannier90' 96 WRITE(stdout, *) 97 ! 98 CALL setup_nnkp() 99 CALL ylm_expansion() 100 IF (scdm_proj) THEN 101 CALL compute_amn_with_scdm() 102 ELSE 103 CALL compute_amn_para() 104 ENDIF 105 CALL compute_mmn_para() 106 ! calling of phases_a_m removed because now we don't call setphases_wrap. 107! CALL phases_a_m() 108 CALL write_band() 109 ! 110 WRITE(stdout, *) 111 WRITE(stdout, *) ' Running Wannier90' 112 ! 113 CALL run_wannier() 114 ! 115 IF (wannier_plot) CALL write_plot() 116 ! 117 CALL lib_dealloc() 118 ! 119 !------------------------------------------------------------------------- 120 END SUBROUTINE pw2wan90epw 121 !------------------------------------------------------------------------- 122 ! 123 !------------------------------------------------------------------------- 124 SUBROUTINE lib_dealloc() 125 !----------------------------------------------------------------------- 126 !! 127 !! Routine to de-allocate Wannier related matrices. 128 !! 129 USE wannierEPW, ONLY : atcart, atsym, kpb, g_kpb, center_w, alpha_w, & 130 l_w, mr_w, r_w, zaxis, xaxis, excluded_band, & 131 m_mat, u_mat, u_mat_opt, a_mat, eigval, & 132 lwindow, gf, ig_, zerophase, wann_centers, & 133 wann_spreads 134 ! 135 IMPLICIT NONE 136 ! 137 ! Local variables 138 INTEGER :: ierr 139 !! Error status 140 ! 141 DEALLOCATE(atcart, STAT = ierr) 142 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating atcart', 1) 143 DEALLOCATE(atsym, STAT = ierr) 144 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating atsym', 1) 145 DEALLOCATE(kpb, STAT = ierr) 146 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating kpb', 1) 147 DEALLOCATE(g_kpb, STAT = ierr) 148 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating g_kpb', 1) 149 DEALLOCATE(center_w, STAT = ierr) 150 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating center_w', 1) 151 DEALLOCATE(alpha_w, STAT = ierr) 152 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating alpha_w', 1) 153 DEALLOCATE(l_w, STAT = ierr) 154 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating l_w', 1) 155 DEALLOCATE(mr_w, STAT = ierr) 156 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating mr_w', 1) 157 DEALLOCATE(r_w, STAT = ierr) 158 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating r_w', 1) 159 DEALLOCATE(zaxis, STAT = ierr) 160 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating zaxis', 1) 161 DEALLOCATE(xaxis, STAT = ierr) 162 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating xaxis', 1) 163 DEALLOCATE(excluded_band, STAT = ierr) 164 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating excluded_band', 1) 165 DEALLOCATE(m_mat, STAT = ierr) 166 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating m_mat', 1) 167 DEALLOCATE(u_mat, STAT = ierr) 168 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating u_mat', 1) 169 DEALLOCATE(u_mat_opt, STAT = ierr) 170 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating u_mat_opt', 1) 171 DEALLOCATE(a_mat, STAT = ierr) 172 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating a_mat', 1) 173 DEALLOCATE(eigval, STAT = ierr) 174 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating eigval', 1) 175 DEALLOCATE(lwindow, STAT = ierr) 176 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating lwindow', 1) 177 DEALLOCATE(gf, STAT = ierr) 178 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating gf', 1) 179 DEALLOCATE(ig_, STAT = ierr) 180 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating ig_', 1) 181 DEALLOCATE(zerophase, STAT = ierr) 182 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating zerophase', 1) 183 DEALLOCATE(wann_centers, STAT = ierr) 184 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating wann_centers', 1) 185 DEALLOCATE(wann_spreads, STAT = ierr) 186 IF (ierr /= 0) CALL errore('lib_dealloc', 'Error deallocating wann_spreads', 1) 187 ! 188 !----------------------------------------------------------------------- 189 END SUBROUTINE lib_dealloc 190 !----------------------------------------------------------------------- 191 ! 192 !------------------------------------------------------------------------- 193 SUBROUTINE setup_nnkp() 194 !----------------------------------------------------------------------- 195 !! 196 !! This routine writes and reads the .nnkp file. 197 !! The file specifies 198 !! 1) The initial projections functions in the format 199 !! num_proj 200 !! proj_site(1,i), proj_site(2,i), proj_site(3,i), proj_l(i), proj_m(i), proj_radial(i) 201 !! proj_z(1,i), proj_z(2,i), proj_z(3,i),proj_x(1,i), proj_x(2,i), proj_x(3,i), proj_zona(i) 202 !! proj_s(i), proj_s_qaxis(1,i), proj_s_qaxis(2,i), proj_s_qaxis(3,i) 203 !! 204 !! 2) begin nnkpts: the nearest neighbours of each k-point, 205 !! and therefore provides the information required to 206 !! calculate the M_mn(k,b) matrix elements 207 !! 208 ! ---------------------------------------------------------------------- 209 USE kinds, ONLY : DP 210 USE io_global, ONLY : meta_ionode, stdout, meta_ionode_id 211 USE mp_world, ONLY : world_comm 212 USE cell_base, ONLY : at, bg, alat 213 USE gvect, ONLY : g, gg 214 USE ions_base, ONLY : nat, tau, ityp, atm 215 USE mp, ONLY : mp_bcast, mp_sum 216 USE wvfct, ONLY : nbnd, npwx 217 USE wannierEPW, ONLY : num_nnmax, mp_grid, atcart, atsym, kpb, g_kpb, & 218 center_w, alpha_w, l_w, mr_w, r_w, zaxis, & 219 xaxis, excluded_band, rlatt, glatt, gf, & 220 csph, ig_, iknum, seedname2, kpt_latt, nnb, & 221 num_bands, n_wannier, nexband, nnbx, n_proj, & 222 spin_eig, spin_qaxis, zerophase 223 USE noncollin_module, ONLY : noncolin 224 USE constants_epw, ONLY : bohr, eps6, twopi 225 USE mp_pools, ONLY : intra_pool_comm 226 USE epwcom, ONLY : scdm_proj 227 USE w90_io, ONLY : post_proc_flag 228 USE io_var, ONLY : iunnkp 229 USE elph2, ONLY : ibndstart, ibndend, nbndep, nbndskip 230 ! 231 IMPLICIT NONE 232 ! 233 LOGICAL :: have_nnkp 234 !! Check if the .nnkp file exists. 235 LOGICAL :: found 236 !! Check if the section in the .nnkp file is found. 237 ! 238 INTEGER :: ik 239 !! Counter on k-points 240 INTEGER :: ib 241 !! Counter on b-vectors 242 INTEGER :: ig 243 !! Index on G_k+b vectors 244 INTEGER :: iw 245 !! Counter on number of projections 246 INTEGER :: ia 247 !! Counter on number of atoms 248 INTEGER :: type 249 !! Type of atom 250 INTEGER :: ibnd 251 !! Counter on band index 252 INTEGER :: indexb 253 !! Index of exclude_bands 254 INTEGER :: ipol 255 !! Counter on polarizations 256 INTEGER :: idum 257 !! Dummy index for reading nnkp file 258 INTEGER :: tmp_auto 259 !! Needed for the selection of projections with SCDM 260 INTEGER :: ierr 261 !! Error status 262 INTEGER, ALLOCATABLE :: ig_check(:, :) 263 !! Temporary index on G_k+b vectors 264 INTEGER :: exclude_bands(nbnd) 265 !! Bands excluded from the calcultion of WFs 266 REAL(KIND = DP) :: xnorm 267 !! Norm of xaxis 268 REAL(KIND = DP) :: znorm 269 !! Norm of zaxis 270 REAL(KIND = DP) :: coseno 271 !! Cosine between xaxis and zaxis 272 REAL(KIND = DP) :: gg_ 273 !! Square of g_(3) 274 REAL(KIND = DP) :: g_(3) 275 !! Temporary vector G_k+b, g_(:) = g_kpb(:,ik,ib) 276 ! 277 INTERFACE 278 !! 279 !! SP: An interface is required because the Wannier routine has optional arguments 280 !! 281 !----------------------------------------------------------------------- 282 SUBROUTINE wannier_setup(seed__name, mp_grid_loc, num_kpts_loc, & 283 real_lattice_loc, recip_lattice_loc, kpt_latt_loc, num_bands_tot, & 284 num_atoms_loc, atom_symbols_loc, atoms_cart_loc, gamma_only_loc, spinors_loc, & 285 nntot_loc, nnlist_loc, nncell_loc, num_bands_loc, num_wann_loc, & 286 proj_site_loc, proj_l_loc, proj_m_loc, proj_radial_loc, proj_z_loc, & 287 proj_x_loc, proj_zona_loc, exclude_bands_loc, proj_s_loc, proj_s_qaxis_loc) 288 !----------------------------------------------------------------------- 289 ! 290 USE kinds, ONLY : DP 291 USE wannierEPW, ONLY : num_nnmax 292 ! 293 IMPLICIT NONE 294 ! 295 CHARACTER(LEN = *), INTENT(in) :: seed__name 296 !! Name 297 CHARACTER(LEN = *), INTENT(in) :: atom_symbols_loc(num_atoms_loc) 298 !! Symbol for atoms 299 LOGICAL, INTENT(in) :: gamma_only_loc 300 !! Gamma point only 301 LOGICAL, INTENT(in) :: spinors_loc 302 !! Local spinor 303 INTEGER, INTENT(in) :: mp_grid_loc(3) 304 !! MP grid 305 INTEGER, INTENT(in) :: num_kpts_loc 306 !! Local number of k-points 307 INTEGER, INTENT(in) :: num_bands_tot 308 !! Number of bands 309 INTEGER, INTENT(in) :: num_atoms_loc 310 !! Number of atoms 311 INTEGER, INTENT(out) :: nntot_loc 312 !! Local nntot 313 INTEGER, INTENT(out) :: nnlist_loc(num_kpts_loc, num_nnmax) 314 !! Local nnlist 315 INTEGER, INTENT(out) :: nncell_loc(3, num_kpts_loc, num_nnmax) 316 !! Local ncell 317 INTEGER, INTENT(out) :: num_bands_loc 318 !! Number of bands 319 INTEGER, INTENT(out) :: num_wann_loc 320 !! Number of Wannier functions 321 INTEGER, INTENT(out) :: proj_l_loc(num_bands_tot) 322 !! Projection of l local momentum 323 INTEGER, INTENT(out) :: proj_m_loc(num_bands_tot) 324 !! Local projection 325 INTEGER, INTENT(out) :: proj_radial_loc(num_bands_tot) 326 !! Radial projection 327 INTEGER, INTENT(out) :: exclude_bands_loc(num_bands_tot) 328 !! Exclude bands 329 INTEGER, OPTIONAL, INTENT(out) :: proj_s_loc(num_bands_tot) 330 !! Projection on s 331 REAL(KIND = DP), INTENT(in) :: REAL_lattice_loc(3, 3) 332 !! Lattice 333 REAL(KIND = DP), INTENT(in) :: recip_lattice_loc(3, 3) 334 !! Reciprocal lattice 335 REAL(KIND = DP), INTENT(in) :: kpt_latt_loc(3, num_kpts_loc) 336 !! kpt grid 337 REAL(KIND = DP), INTENT(in) :: atoms_cart_loc(3, num_atoms_loc) 338 !! Cartesian location of atoms 339 REAL(KIND = DP), INTENT(out) :: proj_site_loc(3, num_bands_tot) 340 !! Site projection 341 REAL(KIND = DP), INTENT(out) :: proj_z_loc(3, num_bands_tot) 342 !! Projection on z 343 REAL(KIND = DP), INTENT(out) :: proj_x_loc(3, num_bands_tot) 344 !! Projection on x 345 REAL(KIND = DP), INTENT(out) :: proj_zona_loc(num_bands_tot) 346 !! Projection on z 347 REAL(KIND = DP), OPTIONAL, INTENT(out) :: proj_s_qaxis_loc(3, num_bands_tot) 348 !! Projection s q-axis 349 ! 350 !----------------------------------------------------------------------- 351 END SUBROUTINE wannier_setup 352 !----------------------------------------------------------------------- 353 END INTERFACE 354 ! 355 num_nnmax = 32 356 ! 357 ! aam: translations between PW2Wannier90 and Wannier90 358 ! pw2wannier90 <==> Wannier90 359 ! nbnd num_bands_tot 360 ! n_wannier num_wann 361 ! num_bands num_bands 362 ! nat num_atoms 363 ! iknum num_kpts 364 ! rlatt TRANSPOSE(real_lattice) 365 ! glatt TRANSPOSE(recip_lattice) 366 ! kpt_latt kpt_latt 367 ! nnb nntot 368 ! kpb nnlist 369 ! g_kpb nncell 370 ! mp_grid mp_grid 371 ! center_w proj_site 372 ! l_w,mr_w,r_w proj_l,proj_m,proj_radial 373 ! xaxis,zaxis proj_x,proj_z 374 ! alpha_w proj_zona 375 ! exclude_bands exclude_bands 376 ! atcart atoms_cart 377 ! atsym atom_symbols 378 ! 379 ALLOCATE(atcart(3, nat), STAT = ierr) 380 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating atcart', 1) 381 ALLOCATE(atsym(nat), STAT = ierr) 382 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating atsym', 1) 383 ALLOCATE(kpb(iknum, num_nnmax), STAT = ierr) 384 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating kpb', 1) 385 ALLOCATE(g_kpb(3, iknum, num_nnmax), STAT = ierr) 386 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating b_kpb', 1) 387 ALLOCATE(center_w(3, nbnd), STAT = ierr) 388 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating center_w', 1) 389 ALLOCATE(alpha_w(nbnd), STAT = ierr) 390 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating alpha_w', 1) 391 ALLOCATE(l_w(nbnd), STAT = ierr) 392 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating l_w', 1) 393 ALLOCATE(mr_w(nbnd), STAT = ierr) 394 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating mr_w', 1) 395 ALLOCATE(r_w(nbnd), STAT = ierr) 396 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating r_w', 1) 397 ALLOCATE(zaxis(3, nbnd), STAT = ierr) 398 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating zaxis', 1) 399 ALLOCATE(xaxis(3, nbnd), STAT = ierr) 400 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating xaxis', 1) 401 ALLOCATE(excluded_band(nbnd), STAT = ierr) 402 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating excluded_band', 1) 403 ! 404 ! real lattice (Cartesians, Angstrom) 405 rlatt(:, :) = TRANSPOSE(at(:, :)) * alat * bohr 406 ! reciprocal lattice (Cartesians, Angstrom) 407 glatt(:, :) = TRANSPOSE(bg(:, :)) * twopi / ( alat * bohr ) 408 ! atom coordinates in Cartesian coords and Angstrom units 409 atcart(:, :) = tau(:, :) * bohr * alat 410 ! atom symbols 411 DO ia = 1, nat 412 type = ityp(ia) 413 atsym(ia) = atm(type) 414 ENDDO 415 ! 416 IF (meta_ionode) THEN 417 !postproc_setup = .TRUE. 418 post_proc_flag = .TRUE. 419 CALL wannier_setup(seedname2, mp_grid, iknum, & ! input 420 rlatt, glatt, kpt_latt, nbnd, & ! input 421 nat, atsym, atcart, .FALSE., noncolin, & ! input 422 nnb, kpb, g_kpb, num_bands, n_wannier, & ! output 423 center_w, l_w, mr_w, r_w, zaxis, & ! output 424 xaxis, alpha_w, exclude_bands) ! output 425 ! SP: In wannier_setup, the .nnkp file is produced. 426 ENDIF 427 ! 428 CALL mp_bcast(nnb, meta_ionode_id, world_comm) 429 CALL mp_bcast(kpb, meta_ionode_id, world_comm) 430 CALL mp_bcast(g_kpb, meta_ionode_id, world_comm) 431 CALL mp_bcast(num_bands, meta_ionode_id, world_comm) 432 CALL mp_bcast(n_wannier, meta_ionode_id, world_comm) 433 CALL mp_bcast(center_w, meta_ionode_id, world_comm) 434 CALL mp_bcast(l_w, meta_ionode_id, world_comm) 435 CALL mp_bcast(mr_w, meta_ionode_id, world_comm) 436 CALL mp_bcast(r_w, meta_ionode_id, world_comm) 437 CALL mp_bcast(zaxis, meta_ionode_id, world_comm) 438 CALL mp_bcast(xaxis, meta_ionode_id, world_comm) 439 CALL mp_bcast(alpha_w, meta_ionode_id, world_comm) 440 CALL mp_bcast(exclude_bands, meta_ionode_id, world_comm) 441 CALL mp_bcast(noncolin, meta_ionode_id, world_comm) 442 ! 443 ! SP: Commented because we now write on file the .nnkp file and read from it. 444 ! 445 ! n_proj = nr. of projections (=#WF unless spinors then =#WF/2) 446 !IF (noncolin) THEN 447 ! n_proj = n_wannier/2 448 !ELSE 449 n_proj = n_wannier 450 !ENDIF 451 ! 452 IF (scdm_proj) THEN 453 WRITE(stdout, *) 454 WRITE(stdout, *) ' Initial Wannier auto_projections' 455 WRITE(stdout, *) 456 ELSE 457 WRITE(stdout, *) 458 WRITE(stdout, *) ' Initial Wannier projections' 459 WRITE(stdout, *) 460 ! 461 DO iw = 1, n_proj 462 WRITE(stdout, '(5x, "(", 3f10.5, ") : l = ", i3, " mr = ", i3)') & 463 center_w(:, iw), l_w(iw), mr_w(iw) 464 ENDDO 465 ENDIF 466 ! 467 excluded_band(1:nbnd) = .FALSE. 468 nexband = 0 469 band_loop: DO ibnd = 1, nbnd 470 indexb = exclude_bands(ibnd) 471 IF (indexb > nbnd .OR. indexb < 0) THEN 472 CALL errore('setup_nnkp', 'wrong excluded band index ', 1) 473 ELSEIF (indexb == 0) THEN 474 EXIT band_loop 475 ELSE 476 nexband = nexband + 1 477 excluded_band(indexb) = .TRUE. 478 ENDIF 479 ENDDO band_loop 480 ! 481 ibndstart = 1 482 DO ibnd = 1, nbnd 483 IF (excluded_band(ibnd) .eqv. .FALSE.) THEN 484 ibndstart = ibnd 485 EXIT 486 ENDIF 487 ENDDO 488 ibndend= nbnd 489 DO ibnd = nbnd, 1, -1 490 IF (excluded_band(ibnd) .eqv. .FALSE.) THEN 491 ibndend = ibnd 492 EXIT 493 ENDIF 494 ENDDO 495 nbndep = ibndend - ibndstart + 1 496 nbndskip = ibndstart - 1 497 ! 498 WRITE(stdout, '(/, " - Number of bands is (", i3, ")")') num_bands 499 WRITE(stdout, '(" - Number of total bands is (", i3, ")")') nbnd 500 WRITE(stdout, '(" - Number of excluded bands is (", i3, ")")') nexband 501 WRITE(stdout, '(" - Number of wannier functions is (", i3, ")")') n_wannier 502 ! 503 IF ((nbnd - nexband) /= num_bands ) & 504 CALL errore('setup_nnkp', ' something wrong with num_bands', 1) 505 ! 506 ! Now we read the .nnkp file 507 ! 508 IF (meta_ionode) THEN ! Read nnkp file on ionode only 509 INQUIRE(FILE = TRIM(seedname2)//".nnkp", EXIST = have_nnkp) 510 IF (.NOT. have_nnkp) THEN 511 CALL errore('setup_nnkp', 'Could not find the file ' & 512 //TRIM(seedname2)//'.nnkp', 1 ) 513 ENDIF 514 OPEN(UNIT = iunnkp, FILE = TRIM(seedname2)//".nnkp", FORM = 'formatted') 515 ENDIF 516 ! 517 IF (meta_ionode) THEN ! read from ionode only 518 IF (noncolin) THEN 519 CALL scan_file_to(iunnkp, 'spinor_projections', found) 520 IF (.NOT. found) THEN 521 CALL errore('setup_nnkp', 'Could not find projections block in ' & 522 //TRIM(seedname2)//'.nnkp', 1) 523 ENDIF 524 ELSE 525 CALL scan_file_to(iunnkp, 'projections', found) 526 IF (.NOT. found) THEN 527 CALL errore('setup_nnkp', 'Could not find projections block in ' & 528 //TRIM(seedname2)//'.nnkp', 1) 529 ENDIF 530 ENDIF 531 READ(iunnkp, *) n_proj 532 ENDIF 533 CALL mp_bcast(n_proj, meta_ionode_id, world_comm) 534 ! 535 ALLOCATE(gf(npwx, n_proj), STAT = ierr) 536 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating gf', 1) 537 ALLOCATE(csph(16, n_proj), STAT = ierr) 538 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating csph', 1) 539 IF (noncolin) THEN 540 ALLOCATE(spin_eig(n_proj), STAT = ierr) 541 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating spin_eig', 1) 542 ALLOCATE(spin_qaxis(3, n_proj), STAT = ierr) 543 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating spin_qaxis', 1) 544 ENDIF 545 ! 546 IF (meta_ionode) THEN ! read from ionode only 547 DO iw = 1, n_proj 548 READ(iunnkp, *) (center_w(ipol, iw), ipol = 1, 3), l_w(iw), mr_w(iw), r_w(iw) 549 READ(iunnkp, *) (zaxis(ipol, iw), ipol = 1, 3), (xaxis(ipol, iw), ipol = 1, 3), alpha_w(iw) 550 xnorm = DSQRT(SUM(xaxis(:, iw) * xaxis(:, iw))) 551 IF (xnorm < eps6) CALL errore('setup_nnkp', '|xaxis| < eps', 1) 552 znorm = DSQRT(SUM(zaxis(:, iw) * zaxis(:, iw))) 553 IF (znorm < eps6) CALL errore('setup_nnkp', '|zaxis| < eps', 1) 554 coseno = SUM(xaxis(:, iw) * zaxis(:, iw)) / xnorm / znorm 555 IF (ABS(coseno) > eps6) CALL errore('setup_nnkp', 'xaxis and zaxis are not orthogonal!', 1) 556 IF (alpha_w(iw) < eps6) CALL errore('setup_nnkp', 'zona value must be positive', 1) 557 ! convert wannier center in cartesian coordinates (in unit of alat) 558 CALL cryst_to_cart(1, center_w(:, iw), at, 1) 559 IF (noncolin) THEN 560 READ(iunnkp, *) spin_eig(iw), (spin_qaxis(ipol, iw), ipol = 1, 3) 561 xnorm = DSQRT(SUM(spin_qaxis(:, iw) * spin_qaxis(:, iw))) 562 IF (xnorm < eps6) CALL errore('setup_nnkp', '|xaxis| < eps', 1) 563 spin_qaxis(:, iw) = spin_qaxis(:, iw) / xnorm 564 ENDIF 565 ENDDO 566 ENDIF 567 ! 568 ! automatic projections 569 IF (meta_ionode) THEN 570 CALL scan_file_to(iunnkp, 'auto_projections', found) 571 IF (found) THEN 572 READ(iunnkp, *) n_wannier 573 READ(iunnkp, *) tmp_auto 574 ! 575 IF (scdm_proj) THEN 576 IF (n_proj > 0) THEN 577 WRITE(stdout, '(//, " ****** begin Error message ******",/)') 578 WRITE(stdout, '(/, " Found a projection block, an auto_projections block", /)') 579 WRITE(stdout, '(/, " and scdm_proj = T in the input file. These three options are inconsistent.", /)') 580 WRITE(stdout, '(/, " Please refer to the Wannier90 User guide for correct use of these flags.",/)') 581 WRITE(stdout, '(/, " ****** end Error message ******",//)') 582 CALL errore('setup_nnkp', 'Inconsistent options for projections.', 1) 583 ELSE 584 IF (tmp_auto /= 0) & 585 CALL errore('setup_nnkp', 'Second entry in auto_projections block is not 0. ' // & 586 'See Wannier90 User Guide in the auto_projections section for clarifications.', 1) 587 ENDIF 588 ELSE 589 ! Fire an error whether or not a projections block is found 590 CALL errore('setup_nnkp', 'scdm_proj = F but found an auto_projections block in ' & 591 //TRIM(seedname2)//'.nnkp', 1) 592 ENDIF 593 ELSE 594 IF (scdm_proj) THEN 595 ! Fire an error whether or not a projections block is found 596 CALL errore('setup_nnkp', 'scdm_proj = T but cannot find an auto_projections block in ' & 597 //TRIM(seedname2)//'.nnkp', 1) 598 ENDIF 599 ENDIF 600 ENDIF 601 ! 602 IF (.NOT. scdm_proj) WRITE(stdout, *) ' - All guiding functions are given ' 603 ! 604 ! Broadcast 605 CALL mp_bcast(center_w, meta_ionode_id, world_comm) 606 CALL mp_bcast(l_w, meta_ionode_id, world_comm) 607 CALL mp_bcast(mr_w, meta_ionode_id, world_comm) 608 CALL mp_bcast(r_w, meta_ionode_id, world_comm) 609 CALL mp_bcast(zaxis, meta_ionode_id, world_comm) 610 CALL mp_bcast(xaxis, meta_ionode_id, world_comm) 611 CALL mp_bcast(alpha_w, meta_ionode_id, world_comm) 612 IF (noncolin) THEN 613 CALL mp_bcast(spin_eig, meta_ionode_id, world_comm) 614 CALL mp_bcast(spin_qaxis, meta_ionode_id, world_comm) 615 ENDIF 616 ! 617 IF (meta_ionode) THEN ! read from ionode only 618 CALL scan_file_to(iunnkp, 'nnkpts', found) 619 IF (.NOT. found) THEN 620 CALL errore('setup_nnkp', 'Could not find nnkpts block in ' & 621 //TRIM(seedname2)//'.nnkp', 1) 622 ENDIF 623 READ(iunnkp, *) nnb 624 ENDIF 625 ! 626 ! Broadcast 627 CALL mp_bcast(nnb, meta_ionode_id, world_comm) 628 ! 629 nnbx = 0 630 nnbx = MAX(nnbx, nnb) 631 ALLOCATE(ig_(iknum, nnbx), STAT = ierr) 632 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating ig_', 1) 633 ALLOCATE(ig_check(iknum, nnbx), STAT = ierr) 634 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating ig_check', 1) 635 ALLOCATE(zerophase(iknum, nnb), STAT = ierr) 636 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error allocating zerophase', 1) 637 zerophase = .FALSE. 638 ! 639 ! Read data about neighbours 640 WRITE(stdout, *) 641 WRITE(stdout, *) ' Reading data about k-point neighbours ' 642 WRITE(stdout, *) 643 IF (meta_ionode) THEN 644 DO ik = 1, iknum 645 DO ib = 1, nnb 646 READ(iunnkp, *) idum, kpb(ik, ib), (g_kpb(ipol, ik, ib), ipol = 1, 3) 647 ENDDO 648 ENDDO 649 ENDIF 650 ! 651 ! Broadcast 652 CALL mp_bcast(kpb, meta_ionode_id, world_comm) 653 CALL mp_bcast(g_kpb, meta_ionode_id, world_comm) 654 ! 655 DO ik = 1, iknum 656 DO ib = 1, nnb 657 IF ((g_kpb(1, ik, ib) == 0) .AND. & 658 (g_kpb(2, ik, ib) == 0) .AND. & 659 (g_kpb(3, ik, ib) == 0) ) zerophase(ik, ib) = .TRUE. 660 g_(:) = REAL(g_kpb(:, ik, ib)) 661 CALL cryst_to_cart(1, g_, bg, 1) 662 gg_ = g_(1) * g_(1) + g_(2) * g_(2) + g_(3) * g_(3) 663 ig_(ik, ib) = 0 664 ig = 1 665 DO WHILE (gg(ig) <= gg_ + eps6) 666 IF ((ABS(g(1, ig) - g_(1)) < eps6) .AND. & 667 (ABS(g(2, ig) - g_(2)) < eps6) .AND. & 668 (ABS(g(3, ig) - g_(3)) < eps6) ) ig_(ik, ib) = ig 669 ig = ig + 1 670 ENDDO 671 ENDDO 672 ENDDO 673 ! 674 ig_check(:, :) = ig_(:, :) 675 CALL mp_sum(ig_check, intra_pool_comm) 676 DO ik = 1, iknum 677 DO ib = 1, nnb 678 IF (ig_check(ik, ib) == 0) & 679 CALL errore('setup_nnkp', 'g_kpb vector is not in the list of Gs', 100 * ik + ib) 680 ENDDO 681 ENDDO 682 DEALLOCATE(ig_check, STAT = ierr) 683 IF (ierr /= 0) CALL errore('setup_nnkp', 'Error deallocating ig_check', 1) 684 ! 685 WRITE(stdout, *) ' - All neighbours are found ' 686 WRITE(stdout, *) 687 ! 688 IF (meta_ionode) THEN 689 CALL scan_file_to(iunnkp, 'exclude_bands', found) 690 IF (.NOT. found) THEN 691 CALL errore('setup_nnkp', 'Could not find exclude_bands block in ' & 692 //TRIM(seedname2)//'.nnkp', 1) 693 ENDIF 694 READ(iunnkp, *) nexband 695 excluded_band(1:nbnd) = .FALSE. 696 DO ibnd = 1, nexband 697 READ(iunnkp, *) indexb 698 IF (indexb < 1 .OR. indexb > nbnd) & 699 CALL errore('setup_nnkp', ' wrong excluded band index ', 1) 700 excluded_band(indexb) = .TRUE. 701 ENDDO 702 num_bands = nbnd - nexband 703 ENDIF 704 ! 705 ! Broadcast 706 CALL mp_bcast(nexband, meta_ionode_id, world_comm) 707 CALL mp_bcast(excluded_band, meta_ionode_id, world_comm) 708 CALL mp_bcast(num_bands, meta_ionode_id, world_comm) 709 ! 710 IF (meta_ionode) THEN 711 CLOSE(iunnkp) 712 ENDIF 713 ! 714 RETURN 715 ! 716 !-------------------------------------------------------------------------- 717 END SUBROUTINE setup_nnkp 718 !-------------------------------------------------------------------------- 719 ! 720 !-------------------------------------------------------------------------- 721 SUBROUTINE scan_file_to(iunr, keyword, found) 722 !----------------------------------------------------------------------- 723 ! 724 IMPLICIT NONE 725 ! 726 CHARACTER(LEN = *), INTENT(in) :: keyword 727 !! Keyword searched for 728 LOGICAL, INTENT(out) :: found 729 !! Check if the section in the .nnkp file is found. 730 INTEGER, INTENT(in) :: iunr 731 !! Unit number for file 732 ! 733 ! Local variables 734 CHARACTER(LEN = 80) :: line1, line2 735 !! 736 ! 737 ! Uncommenting the following line the file scan restarts every time 738 ! from the beginning thus making the reading independent on the order 739 ! of data-blocks 740 ! REWIND(iunr) 741 ! 742 10 CONTINUE 743 READ(iunr, *, END = 20) line1, line2 744 IF (line1 /= 'begin') GOTO 10 745 IF (line2 /= keyword) GOTO 10 746 found = .TRUE. 747 RETURN 748 20 found = .FALSE. 749 REWIND(iunr) 750 ! 751 !------------------------------------------------------------------------ 752 END SUBROUTINE scan_file_to 753 !------------------------------------------------------------------------ 754 ! 755 !------------------------------------------------------------------------ 756 SUBROUTINE run_wannier() 757 !----------------------------------------------------------------------- 758 ! 759 USE kinds, ONLY : DP 760 USE io_global, ONLY : stdout, meta_ionode_id, meta_ionode 761 USE ions_base, ONLY : nat 762 USE mp, ONLY : mp_bcast 763 USE mp_world, ONLY : world_comm 764 USE cell_base, ONLY : alat 765 USE io_files, ONLY : prefix 766 USE io_var, ONLY : iuqpeig 767 USE pwcom, ONLY : nkstot 768 USE wannierEPW,ONLY : u_mat, lwindow, wann_centers, wann_spreads, eigval, & 769 n_wannier, spreads, nnb, rlatt, glatt, kpt_latt, & 770 iknum, seedname2, num_bands, u_mat_opt, atsym, a_mat,& 771 atcart, m_mat, mp_grid, excluded_band 772 USE epwcom, ONLY : eig_read 773 USE wvfct, ONLY : nbnd 774 USE constants_epw, ONLY : zero, czero, bohr 775 ! 776 IMPLICIT NONE 777 ! 778 CHARACTER(LEN = 256) :: tempfile 779 !! Temporary file 780 CHARACTER (LEN = 80) :: line 781 !! Temporary character 782 INTEGER :: iw 783 !! Counter on wannier functions 784 INTEGER :: ik 785 !! Counter of k-point index 786 INTEGER :: ibnd 787 !! Counter on bands 788 INTEGER :: ibnd1 789 !! Band index 790 INTEGER :: ios 791 !! Integer variable for I/O control 792 INTEGER :: ierr 793 !! Error status 794 REAL(KIND = DP), ALLOCATABLE :: eigvaltmp(:, :) 795 !! Temporary array containing the eigenvalues (KS or GW) when read from files 796 ! 797 ALLOCATE(u_mat(n_wannier, n_wannier, iknum), STAT = ierr) 798 IF (ierr /= 0) CALL errore('run_wannier', 'Error allocating u_mat', 1) 799 ALLOCATE(u_mat_opt(num_bands, n_wannier, iknum), STAT = ierr) 800 IF (ierr /= 0) CALL errore('run_wannier', 'Error allocating u_mat_opt', 1) 801 ALLOCATE(lwindow(num_bands, iknum), STAT = ierr) 802 IF (ierr /= 0) CALL errore('run_wannier', 'Error allocating lwindow', 1) 803 ALLOCATE(wann_centers(3, n_wannier), STAT = ierr) 804 IF (ierr /= 0) CALL errore('run_wannier', 'Error allocating wann_centers', 1) 805 ALLOCATE(wann_spreads(n_wannier), STAT = ierr) 806 IF (ierr /= 0) CALL errore('run_wannier', 'Error allocating wann_spreads', 1) 807 ! 808 u_mat(:, :, :) = czero 809 u_mat_opt(:, :, :) = czero 810 wann_centers(:, :) = zero 811 wann_spreads(:) = zero 812 ! 813 IF (meta_ionode) THEN 814 ! read in external eigenvalues, e.g. GW 815 IF (eig_read) THEN 816 ALLOCATE(eigvaltmp(nbnd, iknum), STAT = ierr) 817 IF (ierr /= 0) CALL errore('run_wannier', 'Error allocating eigvaltmp', 1) 818 eigvaltmp(:, :) = zero 819 eigval(:, :) = zero 820 WRITE (stdout, '(5x, a, i5, a, i5, a)') "Reading external electronic eigenvalues (", & 821 nbnd, ",", nkstot,")" 822 tempfile = TRIM(prefix)//'.eig' 823 OPEN(iuqpeig, FILE = tempfile, FORM = 'formatted', ACTION = 'read', IOSTAT = ios) 824 IF (ios /= 0) CALL errore('run_wannier', 'error opening'//tempfile, 1) 825 READ(iuqpeig, '(a)') line 826 DO ik = 1, nkstot 827 ! We do not save the k-point for the moment ==> should be read and 828 ! tested against the current one 829 READ(iuqpeig, '(a)') line 830 READ(iuqpeig, *) eigvaltmp(:, ik) 831 ENDDO 832 DO ik = 1, nkstot 833 ibnd1 = 0 834 DO ibnd = 1, nbnd 835 IF (excluded_band(ibnd)) CYCLE 836 ibnd1 = ibnd1 + 1 837 eigval(ibnd1, ik) = eigvaltmp(ibnd, ik) 838 ENDDO 839 ENDDO 840 CLOSE(iuqpeig) 841 DEALLOCATE(eigvaltmp, STAT = ierr) 842 IF (ierr /= 0) CALL errore('run_wannier', 'Error deallocating eigvaltmp', 1) 843 ENDIF 844 845 ! SP : This file is not used for now. Only required to build the UNK file 846 ! tempfile = TRIM(prefix)//'.mmn' 847 ! OPEN(iummn, FILE = tempfile, IOSTAT = ios, FORM = 'unformatted') 848 ! WRITE(iummn) m_mat 849 ! CLOSE(iummn) 850 851 CALL wannier_run(seedname2, mp_grid, iknum, & ! input 852 rlatt, glatt, kpt_latt, num_bands, & ! input 853 n_wannier, nnb, nat, atsym, & ! input 854 atcart, .FALSE., m_mat, a_mat, eigval, & ! input 855 u_mat, u_mat_opt, lwindow, wann_centers, & ! output 856 wann_spreads, spreads) ! output 857 ENDIF 858 ! 859 CALL mp_bcast(u_mat, meta_ionode_id, world_comm) 860 CALL mp_bcast(u_mat_opt, meta_ionode_id, world_comm) 861 CALL mp_bcast(lwindow, meta_ionode_id, world_comm) 862 CALL mp_bcast(wann_centers, meta_ionode_id, world_comm) 863 CALL mp_bcast(wann_spreads, meta_ionode_id, world_comm) 864 CALL mp_bcast(spreads, meta_ionode_id, world_comm) 865 ! 866 ! 867 ! output the results of the wannierization 868 ! 869 WRITE(stdout, *) 870 WRITE(stdout, *) ' Wannier Function centers (cartesian, alat) and spreads (ang):' 871 WRITE(stdout, *) 872 DO iw = 1, n_wannier 873 WRITE(stdout, '(5x, "(", 3f10.5, ") : ",f8.5)') & 874 wann_centers(:, iw) / alat / bohr, wann_spreads(iw) 875 ENDDO 876 WRITE(stdout, *) 877 ! 878 ! store the final minimisation matrix on disk for later use 879 ! 880 CALL write_filukk 881 ! 882 RETURN 883 ! 884 !----------------------------------------------------------------------- 885 END SUBROUTINE run_wannier 886 !----------------------------------------------------------------------- 887 ! 888 !----------------------------------------------------------------------- 889 SUBROUTINE compute_amn_para() 890 !----------------------------------------------------------------------- 891 !! adapted from compute_amn in pw2wannier90.f90 892 !! parallelization on k-points has been added 893 !! 10/2008 Jesse Noffsinger UC Berkeley 894 !! 01/2018 Roxana Margine updated 895 !! 896 USE kinds, ONLY : DP 897 USE io_global, ONLY : stdout, meta_ionode 898 USE klist, ONLY : nks, igk_k 899 USE klist_epw, ONLY : xk_loc 900 USE wvfct, ONLY : nbnd, npw, npwx, g2kin 901 USE wavefunctions, ONLY : evc 902 USE gvect, ONLY : g, ngm 903 USE gvecw, ONLY : gcutw 904 USE uspp, ONLY : nkb, vkb 905 USE becmod, ONLY : becp, calbec, deallocate_bec_type, & 906 allocate_bec_type 907 USE wannierEPW, ONLY : csph, excluded_band, gf, num_bands, & 908 n_wannier, iknum, n_proj, a_mat, & 909 spin_qaxis, spin_eig, gf_spinor, sgf_spinor 910 USE uspp_param, ONLY : upf 911 USE noncollin_module,ONLY : noncolin, npol 912 USE constants_epw, ONLY : czero, cone, zero, eps6 913 USE mp_global, ONLY : my_pool_id, npool, intra_pool_comm, inter_pool_comm 914 USE mp, ONLY : mp_sum 915 USE kfold, ONLY : ktokpmq 916 USE io_epw, ONLY : readwfc 917 ! 918 IMPLICIT NONE 919 ! 920 LOGICAL :: any_uspp 921 !! Check if uspp are used 922 LOGICAL :: spin_z_pos 923 !! Detect if spin quantisation axis is along +z 924 LOGICAL :: spin_z_neg 925 !! Detect if spin quantisation axis is along -z 926 INTEGER :: iw 927 !! Counter on number of projections 928 INTEGER :: ik 929 !! Counter of k-point index 930 INTEGER :: ibnd 931 !! Counter on band index 932 INTEGER :: ibnd1 933 !! Band index 934 INTEGER :: ipool 935 !! Index of current pool 936 INTEGER :: nkq 937 !! Index of k-point in the pool 938 INTEGER :: nkq_abs 939 !! Absolute index of k-point 940 INTEGER :: ik_g 941 !! Temporary index of k-point, ik_g = nkq_abs 942 INTEGER :: ipol 943 !! Index of spin-up/spin-down polarizations 944 INTEGER :: istart 945 !! Starting index on plane waves 946 INTEGER :: ierr 947 !! Error status 948 REAL(KIND = DP) :: zero_vect(3) 949 !! Temporary zero vector 950 COMPLEX(KIND = DP) :: amn 951 !! Element of A_mn matrix 952 COMPLEX(KIND = DP), EXTERNAL :: ZDOTC 953 !! <psi_mk|g_n> 954 COMPLEX(KIND = DP) :: amn_tmp 955 !! Element of A_mn matrix 956 COMPLEX(KIND = DP) :: fac(2) 957 !! Factor for spin quantization axis 958 COMPLEX(KIND = DP), ALLOCATABLE :: sgf(:, :) 959 !! Guiding functions 960 ! 961 !nocolin: we have half as many projections g(r) defined as wannier 962 ! functions. We project onto (1,0) (ie up spin) and then onto 963 ! (0,1) to obtain num_wann projections. jry 964 ! 965 any_uspp = ANY(upf(:)%tvanp) 966 ! 967 ALLOCATE(a_mat(num_bands, n_wannier, iknum), STAT = ierr) 968 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error allocating a_mat', 1) 969 ALLOCATE(sgf(npwx, n_proj), STAT = ierr) 970 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error allocating sgf', 1) 971 a_mat(:, :, :) = czero 972 sgf(:, :) = czero 973 zero_vect(:) = zero 974 ! 975 IF (noncolin) THEN 976 ALLOCATE( gf_spinor(2 * npwx, n_proj), STAT = ierr) 977 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error allocating gf_spinor', 1) 978 ALLOCATE(sgf_spinor(2 * npwx, n_proj), STAT = ierr) 979 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error allocating sgf_spinor', 1) 980 gf_spinor(:, :) = czero 981 sgf_spinor(:, :) = czero 982 ENDIF 983 ! 984 WRITE(stdout, '(5x, a)') 'AMN' 985 ! 986 IF (any_uspp) THEN 987 CALL allocate_bec_type(nkb, n_wannier, becp) 988 ENDIF 989 ! 990#if defined(__MPI) 991 WRITE(stdout, '(6x, a, i5, a, i4, a)') 'k points = ', iknum, ' in ', npool, ' pools' 992#endif 993 ! 994 DO ik = 1, nks 995 ! 996 ! returns in-pool index nkq and absolute index nkq_abs of xk 997 CALL ktokpmq(xk_loc(:, ik), zero_vect, +1, ipool, nkq, nkq_abs) 998 ik_g = nkq_abs 999 ! 1000 WRITE(stdout, '(5x, i8, " of ", i4, a)') ik , nks, ' on ionode' 1001 FLUSH(stdout) 1002 ! 1003 ! read wfc at k 1004 CALL readwfc(my_pool_id + 1, ik, evc) 1005 ! 1006 ! sorts k+G vectors in order of increasing magnitude, up to ecut 1007 CALL gk_sort(xk_loc(1, ik), ngm, g, gcutw, npw, igk_k(1, ik), g2kin) 1008 ! 1009 CALL generate_guiding_functions(ik) ! they are called gf(npw, n_proj) 1010 ! 1011 IF (noncolin) THEN 1012 CALL orient_gf_spinor(npw) 1013 ENDIF 1014 ! 1015 ! USPP 1016 ! 1017 IF (any_uspp) THEN 1018 CALL init_us_2(npw, igk_k(1,ik), xk_loc(1,ik), vkb) 1019 ! below we compute the product of beta functions with trial func. 1020 IF (noncolin) THEN 1021 CALL calbec(npw, vkb, gf_spinor, becp, n_proj) 1022 ELSE 1023 CALL calbec(npw, vkb, gf, becp, n_proj) 1024 ENDIF 1025 ! and we use it for the product S|trial_func> 1026 IF (noncolin) THEN 1027 CALL s_psi(npwx, npw, n_proj, gf_spinor, sgf_spinor) 1028 ELSE 1029 CALL s_psi(npwx, npw, n_proj, gf, sgf) 1030 ENDIF 1031 ELSE 1032 sgf(:, :) = gf(:, :) 1033 ENDIF 1034 ! 1035 IF (noncolin) THEN 1036 ! we do the projection as g(r)*a(r) and g(r)*b(r) 1037 DO iw = 1, n_proj 1038 ! 1039 spin_z_pos = .FALSE. 1040 spin_z_neg = .FALSE. 1041 ! detect if spin quantisation axis is along z 1042 IF ((ABS(spin_qaxis(1, iw) - 0.0d0) < eps6) .AND. & 1043 (ABS(spin_qaxis(2, iw) - 0.0d0) < eps6) .AND. & 1044 (ABS(spin_qaxis(3, iw) - 1.0d0) < eps6)) THEN 1045 spin_z_pos = .TRUE. 1046 ELSEIF ((ABS(spin_qaxis(1, iw) - 0.0d0) < eps6) .AND. & 1047 (ABS(spin_qaxis(2, iw) - 0.0d0) < eps6) .AND. & 1048 (ABS(spin_qaxis(3, iw) + 1.0d0) < eps6)) THEN 1049 spin_z_neg = .TRUE. 1050 ENDIF 1051 ! 1052 IF (spin_z_pos .OR. spin_z_neg) THEN 1053 ! 1054 ibnd1 = 0 1055 DO ibnd = 1, nbnd 1056 IF (excluded_band(ibnd)) CYCLE 1057 ibnd1 = ibnd1 + 1 1058 ! 1059 IF (spin_z_pos) THEN 1060 ipol = (3 - spin_eig(iw)) / 2 1061 ELSE 1062 ipol = (3 + spin_eig(iw)) / 2 1063 ENDIF 1064 istart = (ipol - 1) * npwx + 1 1065 ! 1066 amn = czero 1067 IF (any_uspp) THEN 1068 amn = ZDOTC(npw, evc(1, ibnd), 1, sgf_spinor(1, iw), 1) 1069 amn = amn + ZDOTC(npw, evc(npwx + 1, ibnd), 1, sgf_spinor(npwx + 1, iw), 1) 1070 ELSE 1071 amn = ZDOTC(npw, evc(istart, ibnd), 1, sgf(1, iw), 1) 1072 ENDIF 1073 CALL mp_sum(amn, intra_pool_comm) 1074 ! 1075 ! changed since n_proj is now read from .nnkp file (n_proj=n_wannier) 1076 ! a_mat(ibnd1, iw + n_proj * (ipol - 1), ik_g) = amn 1077 a_mat(ibnd1, iw, ik_g) = amn 1078 ENDDO ! ibnd 1079 ELSE 1080 ! general routine for quantisation axis (a,b,c) 1081 ! 'up' eigenvector is 1/DSQRT(1+c) [c+1,a+ib] 1082 ! 'down' eigenvector is 1/DSQRT(1-c) [c-1,a+ib] 1083 IF (spin_eig(iw) == 1) THEN 1084 fac(1) = (1.0d0 / DSQRT(1 + spin_qaxis(3, iw))) & 1085 * (spin_qaxis(3, iw) + 1) * cone 1086 fac(2) = (1.0d0 / DSQRT(1 + spin_qaxis(3, iw))) & 1087 * CMPLX(spin_qaxis(1, iw), spin_qaxis(2, iw), KIND = DP) 1088 ELSE 1089 fac(1) = (1.0d0 / DSQRT(1 + spin_qaxis(3, iw))) & 1090 * (spin_qaxis(3, iw)) * cone 1091 fac(2) = (1.0d0 / DSQRT(1 - spin_qaxis(3, iw))) & 1092 * CMPLX(spin_qaxis(1, iw), spin_qaxis(2, iw), KIND = DP) 1093 ENDIF 1094 ! 1095 ibnd1 = 0 1096 DO ibnd = 1, nbnd 1097 IF (excluded_band(ibnd)) CYCLE 1098 ibnd1 = ibnd1 + 1 1099 ! 1100 amn = czero 1101 DO ipol = 1, npol 1102 istart = (ipol - 1) * npwx + 1 1103 amn_tmp = czero 1104 IF (any_uspp) THEN 1105 amn_tmp = ZDOTC(npw, evc(istart, ibnd), 1, sgf_spinor(istart, iw), 1) 1106 CALL mp_sum(amn_tmp, intra_pool_comm) 1107 amn = amn + amn_tmp 1108 ELSE 1109 amn_tmp = ZDOTC(npw, evc(istart, ibnd), 1, sgf(1, iw), 1) 1110 CALL mp_sum(amn_tmp, intra_pool_comm) 1111 amn = amn + fac(ipol) * amn_tmp 1112 ENDIF 1113 ENDDO ! ipol 1114 ! 1115 ! changed since n_proj is now read from .nnkp file (n_proj=n_wannier) 1116 !a_mat(ibnd1, iw + n_proj * (ipol - 1), ik_g) = amn 1117 a_mat(ibnd1, iw, ik_g) = amn 1118 ENDDO ! ibnd 1119 ENDIF ! spin_z_pos 1120 ENDDO ! iw (wannier fns) 1121 ELSE ! scalar wavefunction 1122 DO iw = 1, n_proj 1123 ! 1124 ibnd1 = 0 1125 DO ibnd = 1, nbnd 1126 IF (excluded_band(ibnd)) CYCLE 1127 ibnd1 = ibnd1 + 1 1128 ! 1129 amn = czero 1130 amn = ZDOTC(npw, evc(1, ibnd), 1, sgf(1, iw), 1) 1131 CALL mp_sum(amn, intra_pool_comm) 1132 ! 1133 a_mat(ibnd1, iw, ik_g) = amn 1134 ENDDO ! ibnd 1135 ENDDO ! iw (wannier fns) 1136 ENDIF ! scalar wavefunction 1137 ENDDO ! k-points 1138 ! 1139 DEALLOCATE(csph, STAT = ierr) 1140 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error deallocating csph', 1) 1141 DEALLOCATE(sgf, STAT = ierr) 1142 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error deallocating sgf', 1) 1143 IF (noncolin) THEN 1144 DEALLOCATE(gf_spinor, STAT = ierr) 1145 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error deallocating gf_spinor', 1) 1146 DEALLOCATE(sgf_spinor, STAT = ierr) 1147 IF (ierr /= 0) CALL errore('compute_amn_para', 'Error deallocating sgf_spinor', 1) 1148 ENDIF 1149 ! 1150 IF (any_uspp) CALL deallocate_bec_type(becp) 1151 ! 1152 CALL mp_sum(a_mat, inter_pool_comm) 1153 ! 1154 ! RMDB 1155 !IF (meta_ionode) THEN 1156 ! ALLOCATE(a_mat_tmp(nbnd, n_wannier, iknum), STAT = ierr) 1157 ! IF (ierr /= 0) CALL errore('compute_amn_para', 'Error allocating a_mat_tmp', 1) 1158 ! a_mat_tmp(:, :, :) = czero 1159 ! ! 1160 ! DO ik = 1, iknum 1161 ! DO iw = 1, n_proj 1162 ! ! 1163 ! ibnd1 = 0 1164 ! DO ibnd = 1, nbnd 1165 ! IF (excluded_band(ibnd)) CYCLE 1166 ! ibnd1 = ibnd1 + 1 1167 ! ! 1168 ! a_mat_tmp(ibnd, iw, ik) = a_mat(ibnd1, iw, ik) 1169 ! ENDDO ! ibnd 1170 ! ! 1171 ! DO ibnd = 1, nbnd 1172 ! WRITE(501, '(3i5, 2f18.12)') ibnd, iw, ik, a_mat_tmp(ibnd, iw, ik) 1173 ! ENDDO ! ibnd 1174 ! ENDDO ! iw 1175 ! ENDDO ! ik 1176 ! DEALLOCATE(a_mat_tmp, STAT = ierr) 1177 ! IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating a_mat_tmp', 1) 1178 !ENDIF 1179 ! 1180 WRITE(stdout, *) 1181 WRITE(stdout, '(5x, a)') 'AMN calculated' 1182 ! 1183 RETURN 1184 ! 1185 !----------------------------------------------------------------------- 1186 END SUBROUTINE compute_amn_para 1187 !----------------------------------------------------------------------- 1188 ! 1189 !----------------------------------------------------------------------- 1190 SUBROUTINE compute_amn_with_scdm() 1191 !----------------------------------------------------------------------- 1192 ! 1193 !! This subroutine computes automatic Wannier functions using 1194 !! the selected columns of density matrix (SCDM) method. 1195 !! The subroutine is adapted from compute_amn_with_scdm and 1196 !! compute_amn_with_scdm_spinor subroutines in pw2wannier90.f90. 1197 ! 1198 USE kinds, ONLY : DP 1199 USE constants, ONLY : rytoev 1200 USE io_global, ONLY : stdout, meta_ionode, meta_ionode_id 1201 USE wvfct, ONLY : nbnd, npw, npwx, et, g2kin 1202 USE gvecw, ONLY : gcutw 1203 USE wavefunctions, ONLY : evc, psic, psic_nc 1204 USE noncollin_module, ONLY : noncolin 1205 USE wannierEPW, ONLY : n_wannier, iknum, num_bands, nexband, & 1206 excluded_band, kpt_latt, a_mat 1207 USE klist, ONLY : nks, igk_k 1208 USE klist_epw, ONLY : xk_all, xk_loc 1209 USE gvect, ONLY : g, ngm, mill 1210 USE fft_base, ONLY : dffts 1211 USE scatter_mod, ONLY : gather_grid 1212 USE fft_interfaces, ONLY : invfft 1213 USE mp, ONLY : mp_bcast, mp_sum 1214 USE mp_world, ONLY : world_comm 1215 USE mp_global, ONLY : my_pool_id, npool, intra_pool_comm, inter_pool_comm 1216 USE constants_epw, ONLY : zero, czero, one, twopi 1217 USE kfold, ONLY : ktokpmq 1218 USE epwcom, ONLY : scdm_entanglement, scdm_mu, scdm_sigma 1219 USE io_epw, ONLY : readwfc 1220 ! 1221 IMPLICIT NONE 1222 ! 1223 ! Local variables 1224 LOGICAL :: found_gamma 1225 !! find G-point index 1226 INTEGER :: ik 1227 !! Counter over k-points 1228 INTEGER :: ibnd 1229 !! Counter over bands 1230 INTEGER :: ibnd1 1231 !! Band index 1232 INTEGER :: iw 1233 !! Counter over the nr. of projections 1234 INTEGER :: nrtot 1235 !! Size of FFT grid 1236 INTEGER :: ierr 1237 !! Error status 1238 INTEGER :: info 1239 !! if info=0 succesful exit of QR factorization 1240 INTEGER :: lcwork 1241 !! 1242 INTEGER :: ipt, jpt, kpt, lpt 1243 !! Counter over FFT grid 1244 INTEGER :: count_piv_spin 1245 !! 1246 INTEGER :: minmn 1247 !! 1248 INTEGER :: minmn2 1249 !! 1250 INTEGER :: maxmn2 1251 !! 1252 INTEGER :: numbands 1253 !! Nr. of bands 1254 INTEGER :: nbtot 1255 !! Total nr. of bands 1256 INTEGER :: ig 1257 !! Counter on G vectors 1258 INTEGER :: ig_local 1259 !! Counter on G vectors 1260 INTEGER :: ipool 1261 !! Index of current pool 1262 INTEGER :: nkq 1263 !! Index of k-point in the pool 1264 INTEGER :: nkq_abs 1265 !! Absolute index of k-point 1266 INTEGER :: ik_g 1267 !! Temporary index of k-point, ik_g = nkq_abs 1268 INTEGER, ALLOCATABLE :: piv(:) 1269 !! vv: Pivot array in the QR factorization 1270 INTEGER, ALLOCATABLE :: piv_pos(:) 1271 !! jml: Position of piv 1272 INTEGER, ALLOCATABLE :: piv_spin(:) 1273 !! jml: Spin index of piv 1274 REAL(KIND = DP):: pnorm 1275 !! 1276 REAL(KIND = DP):: norm_psi 1277 !! norm of wave function |Psi| 1278 REAL(KIND = DP):: f_gamma 1279 !! generalization of Fermi-Dirac probability for occupied states at G-point 1280 REAL(KIND = DP):: tpi_r_dot_k 1281 !! scalar product i*2*pi*k*r 1282 REAL(KIND = DP):: tpi_r_dot_g 1283 !! scalar prodoct i*2*pi*G*r 1284 REAL(KIND = DP), EXTERNAL :: DDOT 1285 !! Scalar product of two vectors 1286 REAL(KIND = DP) :: zero_vect(3) 1287 !! Temporary zero vector 1288 REAL(KIND = DP) :: g_(3) 1289 !! Temporary vector G 1290 REAL(KIND = DP), ALLOCATABLE :: focc(:) 1291 !! generalization of Fermi-Dirac probability for occupied states 1292 ! 1293 ! vv: Real array for the QR factorization and SVD 1294 REAL(KIND = DP), ALLOCATABLE :: rwork(:) 1295 !! 1296 REAL(KIND = DP), ALLOCATABLE :: rwork2(:) 1297 !! 1298 REAL(KIND = DP), ALLOCATABLE :: singval(:) 1299 !! 1300 REAL(KIND = DP), ALLOCATABLE :: rpos(:, :) 1301 !! real space coords. of FFT grid 1302 REAL(KIND = DP), ALLOCATABLE :: cpos(:, :) 1303 !! real space coords. of FFT grid 1304 ! 1305 COMPLEX(KIND = DP) :: nowfc_tmp 1306 !! 1307 COMPLEX(KIND = DP), ALLOCATABLE :: nowfc(:, :) 1308 !! sum_G (psi(G) * e^(i*G*r)) * focc * phase / norm_psi 1309 COMPLEX(KIND = DP), ALLOCATABLE :: psi_gamma(:, :) 1310 !! psi * f_gamma / norm_psi 1311 COMPLEX(KIND = DP) :: tmp_cwork(2) 1312 !! 1313 COMPLEX(KIND = DP), ALLOCATABLE :: phase(:) 1314 !! exp(i*k*r) 1315 COMPLEX(KIND = DP), ALLOCATABLE :: phase_g(:, :) 1316 !! exp(i*G*r) 1317 ! 1318 ! complex arrays in QR factorization 1319 COMPLEX(KIND = DP), ALLOCATABLE :: qr_tau(:) 1320 !! 1321 COMPLEX(KIND = DP), ALLOCATABLE :: cwork(:) 1322 !! 1323 COMPLEX(KIND = DP), ALLOCATABLE :: umat(:, :) 1324 !! 1325 COMPLEX(KIND = DP), ALLOCATABLE :: vtmat(:, :) 1326 !! 1327 COMPLEX(KIND = DP), ALLOCATABLE :: amat(:, :) 1328 !! Elements of A_mn matrix 1329 !COMPLEX(KIND = DP), ALLOCATABLE :: a_mat_tmp(:, :, :) 1330 !! Temprary variable to store a_mat (used for debugging) 1331 ! 1332 ! vv: Write info about SCDM in output 1333 IF (TRIM(scdm_entanglement) == 'isolated') THEN 1334 WRITE(stdout, '(1x, a, a/)') 'Case : ', TRIM(scdm_entanglement) 1335 ELSEIF ((TRIM(scdm_entanglement) == 'erfc') .OR. & 1336 (TRIM(scdm_entanglement) == 'gaussian')) THEN 1337 WRITE(stdout, '(1x, a, a)') 'Case : ', TRIM(scdm_entanglement) 1338 WRITE(stdout, '(1x, a, f10.3, a/, 1x, a, f10.3, a/)') 'mu = ', scdm_mu, ' eV', 'sigma =', scdm_sigma, ' eV' 1339 ENDIF 1340 ! 1341 ! vv: Allocate all the variables for the SCDM method: 1342 ! 1)For the QR decomposition 1343 ! 2)For the unk's on the real grid 1344 ! 3)For the SVD 1345 ! 1346 IF (TRIM(scdm_entanglement) == 'isolated') THEN 1347 numbands = n_wannier 1348 nbtot = n_wannier + nexband 1349 ELSE 1350 numbands = nbnd - nexband 1351 nbtot = nbnd 1352 ENDIF 1353 ! 1354 nrtot = dffts%nr1 * dffts%nr2 * dffts%nr3 1355 IF (noncolin) THEN 1356 minmn = MIN(numbands, 2 * nrtot) 1357 ALLOCATE(piv(2 * nrtot), STAT = ierr) 1358 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating piv', 1) 1359 ALLOCATE(piv_pos(n_wannier), STAT = ierr) 1360 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating piv_pos', 1) 1361 ALLOCATE(piv_spin(n_wannier), STAT = ierr) 1362 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating piv_spin', 1) 1363 IF (meta_ionode) THEN 1364 ALLOCATE(qr_tau(2 * minmn), STAT = ierr) 1365 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating qr_tau', 1) 1366 ALLOCATE(rwork(2 * 2 * nrtot), STAT = ierr) 1367 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating rwork', 1) 1368 ALLOCATE(psi_gamma(2 * nrtot, numbands), STAT = ierr) 1369 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating psi_gamma', 1) 1370 ENDIF 1371 ELSE 1372 minmn = MIN(numbands, nrtot) 1373 ALLOCATE(piv(nrtot), STAT = ierr) 1374 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating piv', 1) 1375 IF (meta_ionode) THEN 1376 ALLOCATE(qr_tau(2 * minmn), STAT = ierr) 1377 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating qr_tau', 1) 1378 ALLOCATE(rwork(2 * nrtot), STAT = ierr) 1379 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating rwork', 1) 1380 ALLOCATE(psi_gamma(nrtot, numbands), STAT = ierr) 1381 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating psi_gamma', 1) 1382 ENDIF 1383 ENDIF 1384 ! 1385 ALLOCATE(nowfc(n_wannier, numbands), STAT = ierr) 1386 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating nowfc', 1) 1387 ALLOCATE(focc(numbands), STAT = ierr) 1388 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating focc', 1) 1389 minmn2 = MIN(numbands, n_wannier) 1390 maxmn2 = MAX(numbands, n_wannier) 1391 ALLOCATE(rwork2(5 * minmn2), STAT = ierr) 1392 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating rwork2', 1) 1393 ! 1394 ALLOCATE(rpos(nrtot, 3), STAT = ierr) 1395 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating rpos', 1) 1396 ALLOCATE(cpos(n_wannier, 3), STAT = ierr) 1397 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating cpos', 1) 1398 ALLOCATE(phase(n_wannier), STAT = ierr) 1399 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating phase', 1) 1400 ALLOCATE(singval(n_wannier), STAT = ierr) 1401 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating singval', 1) 1402 ALLOCATE(umat(numbands, n_wannier), STAT = ierr) 1403 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating umat', 1) 1404 ALLOCATE(vtmat(n_wannier, n_wannier), STAT = ierr) 1405 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating vtmat', 1) 1406 ALLOCATE(amat(numbands, n_wannier), STAT = ierr) 1407 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating amat', 1) 1408 ! 1409 ALLOCATE(a_mat(numbands, n_wannier, iknum), STAT = ierr) 1410 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating a_mat', 1) 1411 a_mat(:, :, :) = czero 1412 zero_vect(:) = zero 1413 ! 1414 WRITE(stdout, '(5x, a)') 'AMN with SCDM' 1415 ! 1416 ! Check that Gamma-point is first in the list of k-vectors 1417 ! 1418 found_gamma = .FALSE. 1419 found_gamma = kpt_latt(1, 1) == 0.0d0 .AND. & 1420 kpt_latt(2, 1) == 0.0d0 .AND. & 1421 kpt_latt(3, 1) == 0.0d0 1422 IF (.NOT. found_gamma) CALL errore('compute_amn_with_scdm', 'No Gamma point found.', 1) 1423 ! 1424 IF (meta_ionode) THEN 1425 ! read wfc at G-point 1426 ik = 1 1427 CALL readwfc(my_pool_id + 1, ik, evc) 1428 ! 1429 ! sorts k+G vectors in order of increasing magnitude, up to ecut 1430 CALL gk_sort(xk_all(1, ik), ngm, g, gcutw, npw, igk_k(1, ik), g2kin) 1431 ! 1432 ibnd1 = 0 1433 f_gamma = zero 1434 psi_gamma(:, :) = czero 1435 IF (noncolin) THEN 1436 DO ibnd = 1, nbtot 1437 IF (excluded_band(ibnd)) CYCLE 1438 ibnd1 = ibnd1 + 1 1439 ! 1440 ! check ibnd1 <= numbands 1441 IF (ibnd1 > numbands) & 1442 CALL errore('compute_amn_with_scdm', & 1443 'Something wrong with the number of bands. Check exclude_bands.', 1) 1444 IF (TRIM(scdm_entanglement) == 'isolated') THEN 1445 f_gamma = 1.d0 1446 ELSEIF (TRIM(scdm_entanglement) == 'erfc') THEN 1447 f_gamma = 0.5d0 * ERFC((et(ibnd, ik) * rytoev - scdm_mu) / scdm_sigma) 1448 ELSEIF (TRIM(scdm_entanglement) == 'gaussian') THEN 1449 f_gamma = EXP(-1.d0 * ((et(ibnd, ik) * rytoev - scdm_mu)**2) / (scdm_sigma**2)) 1450 ELSE 1451 CALL errore('compute_amn_with_scdm', 'scdm_entanglement value not recognized.', 1) 1452 ENDIF 1453 ! 1454 ! vv: Compute unk's on a real grid (the fft grid) 1455 ! 1456 psic_nc(:, :) = czero 1457 psic_nc(dffts%nl(igk_k(1:npw, ik)), 1) = evc( 1:npw, ibnd) 1458 psic_nc(dffts%nl(igk_k(1:npw, ik)), 2) = evc(1 + npwx:npw + npwx, ibnd) 1459 CALL invfft('Wave', psic_nc(:, 1), dffts) 1460 CALL invfft('Wave', psic_nc(:, 2), dffts) 1461 ! 1462 ! vv: Build Psi_k = Unk * focc at G-point only 1463 pnorm = REAL(SUM(psic_nc(1:nrtot, 1) * CONJG(psic_nc(1:nrtot, 1))), KIND = DP) + & 1464 REAL(SUM(psic_nc(1:nrtot, 2) * CONJG(psic_nc(1:nrtot, 2))), KIND = DP) 1465 norm_psi = DSQRT(pnorm) 1466 psi_gamma( 1:nrtot, ibnd1) = (psic_nc(1:nrtot, 1) / norm_psi) * f_gamma 1467 psi_gamma(1 + nrtot:2 * nrtot, ibnd1) = (psic_nc(1:nrtot, 2) / norm_psi) * f_gamma 1468 ENDDO ! ibnd 1469 ELSE ! scalar 1470 DO ibnd = 1, nbtot 1471 IF (excluded_band(ibnd)) CYCLE 1472 ibnd1 = ibnd1 + 1 1473 ! 1474 ! check ibnd1 <= numbands 1475 IF (ibnd1 > numbands) & 1476 CALL errore('compute_amn_with_scdm', & 1477 'Something wrong with the number of bands. Check exclude_bands.', 1) 1478 IF (TRIM(scdm_entanglement) == 'isolated') THEN 1479 f_gamma = 1.d0 1480 ELSEIF (TRIM(scdm_entanglement) == 'erfc') THEN 1481 f_gamma = 0.5d0 * ERFC((et(ibnd, ik) * rytoev - scdm_mu) / scdm_sigma) 1482 ELSEIF (TRIM(scdm_entanglement) == 'gaussian') THEN 1483 f_gamma = EXP(-1.d0 * ((et(ibnd, ik) * rytoev - scdm_mu)**2) / (scdm_sigma**2)) 1484 ELSE 1485 CALL errore('compute_amn_with_scdm', 'scdm_entanglement value not recognized.', 1) 1486 ENDIF 1487 ! 1488 ! vv: Compute unk's on a real grid (the fft grid) 1489 ! 1490 psic(:) = czero 1491 psic(dffts%nl(igk_k(1:npw, ik))) = evc(1:npw, ibnd) 1492 CALL invfft('Wave', psic, dffts) 1493 ! 1494 ! vv: Build Psi_k = Unk * focc at G-point only 1495 pnorm = REAL(SUM(psic(1:nrtot) * CONJG(psic(1:nrtot))), KIND = DP) 1496 norm_psi = DSQRT(pnorm) 1497 psi_gamma(1:nrtot, ibnd1) = (psic(1:nrtot) / norm_psi) * f_gamma 1498 ENDDO ! ibnd 1499 ENDIF ! noncolin 1500 ! 1501 ! vv: Perform QR factorization with pivoting on Psi_Gamma 1502 ! vv: Preliminary call to define optimal values for lwork and cwork size 1503 ! 1504 piv(:) = 0 1505 qr_tau(:) = czero 1506 tmp_cwork(:) = czero 1507 rwork(:) = zero 1508 IF (noncolin) THEN 1509 CALL ZGEQP3(numbands, 2 * nrtot, TRANSPOSE(CONJG(psi_gamma)), numbands, & 1510 piv, qr_tau, tmp_cwork, -1, rwork, info) 1511 ELSE 1512 CALL ZGEQP3(numbands, nrtot, TRANSPOSE(CONJG(psi_gamma)), numbands, & 1513 piv, qr_tau, tmp_cwork, -1, rwork, info) 1514 ENDIF 1515 IF (info /= 0) CALL errore('compute_amn_with_scdm', 'Error in computing the QR factorization', 1) 1516 ! 1517 lcwork = AINT(REAL(tmp_cwork(1))) 1518 ALLOCATE(cwork(lcwork), STAT = ierr) 1519 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating cwork', 1) 1520 piv(:) = 0 1521 qr_tau(:) = czero 1522 cwork(:) = czero 1523 rwork(:) = zero 1524 ! 1525 IF (noncolin) THEN 1526 CALL ZGEQP3(numbands, 2 * nrtot, TRANSPOSE(CONJG(psi_gamma)), numbands, & 1527 piv, qr_tau, cwork, lcwork, rwork, info) 1528 ELSE 1529 CALL ZGEQP3(numbands, nrtot, TRANSPOSE(CONJG(psi_gamma)), numbands, & 1530 piv, qr_tau, cwork, lcwork, rwork, info) 1531 ENDIF 1532 IF (info /= 0) CALL errore('compute_amn_with_scdm', 'Error in computing the QR factorization', 1) 1533 ! 1534 DEALLOCATE(qr_tau, STAT = ierr) 1535 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating qr_tau', 1) 1536 DEALLOCATE(rwork, STAT = ierr) 1537 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating rwork', 1) 1538 DEALLOCATE(psi_gamma, STAT = ierr) 1539 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating psi_gamma', 1) 1540 DEALLOCATE(cwork, STAT = ierr) 1541 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating cwork', 1) 1542 ENDIF ! meta_ionode 1543 ! 1544 CALL mp_bcast(piv, meta_ionode_id, world_comm) 1545 ! 1546 ! jml: calculate position and spin part of piv 1547 IF (noncolin) THEN 1548 count_piv_spin = 0 1549 piv_pos(:) = 0 1550 piv_spin(:) = 0 1551 DO iw = 1, n_wannier 1552 IF (piv(iw) <= nrtot) THEN 1553 piv_pos(iw) = piv(iw) 1554 piv_spin(iw) = 1 1555 count_piv_spin = count_piv_spin + 1 1556 ELSE 1557 piv_pos(iw) = piv(iw) - nrtot 1558 piv_spin(iw) = 2 1559 ENDIF 1560 ENDDO 1561 WRITE(stdout, '(a, i5)') " Number of pivot points with spin up : ", count_piv_spin 1562 WRITE(stdout, '(a, i5)') " Number of pivot points with spin down: ", n_wannier - count_piv_spin 1563 ENDIF 1564 ! 1565 ! vv: Compute the points 1566 lpt = 0 1567 rpos(:, :) = zero 1568 DO kpt = 0, dffts%nr3 - 1 1569 DO jpt = 0, dffts%nr2 - 1 1570 DO ipt = 0, dffts%nr1 - 1 1571 lpt = lpt + 1 1572 rpos(lpt, 1) = REAL(ipt, KIND = DP) / REAL(dffts%nr1, KIND = DP) 1573 rpos(lpt, 2) = REAL(jpt, KIND = DP) / REAL(dffts%nr2, KIND = DP) 1574 rpos(lpt, 3) = REAL(kpt, KIND = DP) / REAL(dffts%nr3, KIND = DP) 1575 ENDDO 1576 ENDDO 1577 ENDDO 1578 ! 1579 cpos(:, :) = zero 1580 IF (noncolin) THEN 1581 DO iw = 1, n_wannier 1582 cpos(iw, :) = rpos(piv_pos(iw), :) 1583 cpos(iw, :) = cpos(iw, :) - ANINT(cpos(iw, :)) 1584 ENDDO 1585 ELSE 1586 DO iw = 1, n_wannier 1587 cpos(iw, :) = rpos(piv(iw), :) 1588 cpos(iw, :) = cpos(iw, :) - ANINT(cpos(iw, :)) 1589 ENDDO 1590 ENDIF 1591 ! 1592#if defined(__MPI) 1593 WRITE(stdout, '(6x, a, i5, a, i4, a)') 'k points = ', iknum, ' in ', npool, ' pools' 1594#endif 1595 ! 1596 DO ik = 1, nks 1597 ! 1598 ! returns in-pool index nkq and absolute index nkq_abs of xk 1599 CALL ktokpmq(xk_loc(:,ik), zero_vect, +1, ipool, nkq, nkq_abs) 1600 ik_g = nkq_abs 1601 ! 1602 WRITE(stdout, '(5x, i8, " of ", i4, a)') ik , nks, ' on ionode' 1603 FLUSH(stdout) 1604 ! 1605 ! read wfc at k 1606 CALL readwfc(my_pool_id + 1, ik, evc) 1607 ! 1608 ! sorts k+G vectors in order of increasing magnitude, up to ecut 1609 CALL gk_sort(xk_loc(1, ik), ngm, g, gcutw, npw, igk_k(1, ik), g2kin) 1610 ! 1611 ! vv: SCDM method for generating the Amn matrix 1612 ! jml: calculate of psi_nk at pivot points using slow FT 1613 ! This is faster than using invfft because the number of pivot 1614 ! points is much smaller than the number of FFT grid points. 1615 ! 1616 ! jml: calculate phase factors before the loop over bands 1617 ALLOCATE(phase_g(npw, n_wannier), STAT = ierr) 1618 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating phase_g', 1) 1619 phase_g(:, :) = czero 1620 phase(:) = czero 1621 ! 1622 DO iw = 1, n_wannier 1623 !tpi_r_dot_k = twopi * DDOT(3, cpos(iw, :), 1, kpt_latt(:, ik_g), 1) 1624 tpi_r_dot_k = twopi * (cpos(iw, 1) * kpt_latt(1, ik_g) + & 1625 cpos(iw, 2) * kpt_latt(2, ik_g) + & 1626 cpos(iw, 3) * kpt_latt(3, ik_g)) 1627 phase(iw) = CMPLX(COS(tpi_r_dot_k), SIN(tpi_r_dot_k), KIND = DP) 1628 DO ig_local = 1, npw 1629 ig = igk_k(ig_local, ik) 1630 g_(:) = REAL(mill(:, ig), KIND = DP) 1631 !tpi_r_dot_g = twopi * DDOT(3, cpos(iw, :), 1, g_(:), 1) 1632 tpi_r_dot_g = twopi * (cpos(iw, 1) * g_(1) + & 1633 cpos(iw, 2) * g_(2) + & 1634 cpos(iw, 3) * g_(3)) 1635 phase_g(ig_local, iw) = CMPLX(COS(tpi_r_dot_g), SIN(tpi_r_dot_g), KIND = DP) 1636 ENDDO 1637 ENDDO 1638 ! 1639 ! vv: Generate the occupation numbers matrix according to scdm_entanglement 1640 focc(:) = zero 1641 nowfc(:, :) = czero 1642 ibnd1 = 0 1643 IF (noncolin) THEN 1644 DO ibnd = 1, nbtot 1645 IF (excluded_band(ibnd)) CYCLE 1646 ibnd1 = ibnd1 + 1 1647 ! 1648 IF (TRIM(scdm_entanglement) == 'isolated') THEN 1649 focc(ibnd1) = 1.0d0 1650 ELSEIF (TRIM(scdm_entanglement) == 'erfc') THEN 1651 focc(ibnd1) = 0.5d0 * ERFC((et(ibnd, ik_g) * rytoev - scdm_mu) / scdm_sigma) 1652 ELSEIF (TRIM(scdm_entanglement) == 'gaussian') THEN 1653 focc(ibnd1) = EXP(-1.0d0 * ((et(ibnd, ik_g) * rytoev - scdm_mu)**2) / (scdm_sigma**2)) 1654 ELSE 1655 CALL errore('compute_amn_with_scdm', 'scdm_entanglement value not recognized.', 1) 1656 ENDIF 1657 ! 1658 norm_psi = REAL(SUM(evc( 1:npw, ibnd) * CONJG(evc( 1:npw, ibnd)))) + & 1659 REAL(SUM(evc(1 + npwx:npw + npwx, ibnd) * CONJG(evc(1 + npwx:npw + npwx, ibnd)) )) 1660 CALL mp_sum(norm_psi, intra_pool_comm) 1661 norm_psi = DSQRT(norm_psi) 1662 ! 1663 ! jml: nowfc = sum_G (psi(G) * exp(i*G*r)) * focc * phase(iw) / norm_psi 1664 ! 1665 DO iw = 1, n_wannier 1666 IF (piv_spin(iw) == 1) THEN ! spin up 1667 nowfc_tmp = SUM(evc( 1:npw, ibnd) * phase_g(1:npw, iw)) 1668 ELSE 1669 nowfc_tmp = SUM(evc(1 + npwx:npw + npwx, ibnd) * phase_g(1:npw, iw) ) 1670 ENDIF 1671 nowfc(iw, ibnd1) = nowfc_tmp * phase(iw) * focc(ibnd1) / norm_psi 1672 ENDDO ! iw (wannier fns) 1673 ENDDO ! ibnd 1674 ELSE ! scalar wavefunction 1675 DO ibnd = 1, nbtot 1676 IF (excluded_band(ibnd)) CYCLE 1677 ibnd1 = ibnd1 + 1 1678 ! 1679 IF (TRIM(scdm_entanglement) == 'isolated') THEN 1680 focc(ibnd1) = 1.0d0 1681 ELSEIF (TRIM(scdm_entanglement) == 'erfc') THEN 1682 focc(ibnd1) = 0.5d0 * ERFC((et(ibnd, ik_g) * rytoev - scdm_mu) / scdm_sigma) 1683 ELSEIF (TRIM(scdm_entanglement) == 'gaussian') THEN 1684 focc(ibnd1) = EXP(-1.0d0 * ((et(ibnd, ik_g) * rytoev - scdm_mu)**2) / (scdm_sigma**2)) 1685 ELSE 1686 CALL errore('compute_amn_with_scdm', 'scdm_entanglement value not recognized.', 1) 1687 ENDIF 1688 ! 1689 norm_psi = REAL(SUM(evc(1:npw, ibnd) * CONJG(evc(1:npw, ibnd)))) 1690 CALL mp_sum(norm_psi, intra_pool_comm) 1691 norm_psi = DSQRT(norm_psi) 1692 ! 1693 ! jml: nowfc = sum_G (psi(G) * exp(i*G*r)) * focc * phase(iw) / norm_psi 1694 ! 1695 DO iw = 1, n_wannier 1696 nowfc_tmp = SUM(evc(1:npw, ibnd) * phase_g(1:npw, iw)) 1697 nowfc(iw, ibnd1) = nowfc_tmp * phase(iw) * focc(ibnd1) / norm_psi 1698 ENDDO ! iw (wannier fns) 1699 ENDDO ! ibnd 1700 ENDIF ! noncollinear 1701 CALL mp_sum(nowfc, intra_pool_comm) 1702 ! 1703 DEALLOCATE(phase_g, STAT = ierr) 1704 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating phase_g', 1) 1705 ! 1706 singval(:) = zero 1707 umat(:, :) = czero 1708 vtmat(:, :) = czero 1709 tmp_cwork(:) = czero 1710 rwork2(:) = zero 1711 ! 1712 CALL ZGESVD('s', 's', numbands, n_wannier, TRANSPOSE(CONJG(nowfc)), numbands, & 1713 singval, umat, numbands, vtmat, n_wannier, tmp_cwork, -1, rwork2, info) 1714 ! 1715 lcwork = AINT(REAL(tmp_cwork(1))) 1716 ALLOCATE(cwork(lcwork), STAT = ierr) 1717 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating cwork', 1) 1718 singval(:) = zero 1719 umat(:, :) = czero 1720 vtmat(:, :) = czero 1721 cwork(:) = czero 1722 rwork2(:) = zero 1723 ! 1724 ! vv: SVD to generate orthogonal projections 1725 CALL ZGESVD('s', 's', numbands, n_wannier, TRANSPOSE(CONJG(nowfc)), numbands, & 1726 singval, umat, numbands, vtmat, n_wannier, cwork, lcwork, rwork2, info) 1727 IF(info /= 0) CALL errore('compute_amn_with_scdm', 'Error in computing the SVD of the PSI matrix in the SCDM method', 1) 1728 DEALLOCATE(cwork, STAT = ierr) 1729 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating cwork', 1) 1730 ! 1731 amat(:, :) = czero 1732 amat = MATMUL(umat, vtmat) 1733 DO iw = 1, n_wannier 1734 ibnd1 = 0 1735 DO ibnd = 1, nbtot 1736 IF (excluded_band(ibnd)) CYCLE 1737 ibnd1 = ibnd1 + 1 1738 ! 1739 a_mat(ibnd1, iw, ik_g) = amat(ibnd1, iw) 1740 ENDDO ! ibnd 1741 ENDDO ! iw (wannier fns) 1742 ! 1743 ENDDO ! k-points 1744 ! 1745 CALL mp_sum(a_mat, inter_pool_comm) 1746 ! 1747 ! RMDB 1748 !IF (meta_ionode) THEN 1749 ! ALLOCATE(a_mat_tmp(nbtot, n_wannier, iknum), STAT = ierr) 1750 ! IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error allocating a_mat_tmp', 1) 1751 ! a_mat_tmp(:, :, :) = czero 1752 ! ! 1753 ! DO ik = 1, iknum 1754 ! DO iw = 1, n_wannier 1755 ! ! 1756 ! ibnd1 = 0 1757 ! DO ibnd = 1, nbtot 1758 ! IF (excluded_band(ibnd)) CYCLE 1759 ! ibnd1 = ibnd1 + 1 1760 ! ! 1761 ! a_mat_tmp(ibnd, iw, ik) = a_mat(ibnd1, iw, ik) 1762 ! ENDDO ! ibnd 1763 ! ! 1764 ! DO ibnd = 1, nbtot 1765 ! WRITE(501, '(3i5, 2f18.12)') ibnd, iw, ik, a_mat_tmp(ibnd, iw, ik) 1766 ! ENDDO ! ibnd 1767 ! ENDDO ! iw 1768 ! ENDDO ! ik 1769 ! DEALLOCATE(a_mat_tmp, STAT = ierr) 1770 ! IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating a_mat_tmp', 1) 1771 !ENDIF 1772 ! RMDB 1773 ! 1774 DEALLOCATE(piv, STAT = ierr) 1775 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating piv', 1) 1776 DEALLOCATE(nowfc, STAT = ierr) 1777 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating nowfc', 1) 1778 DEALLOCATE(focc, STAT = ierr) 1779 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating focc', 1) 1780 DEALLOCATE(rwork2, STAT = ierr) 1781 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating rwork2', 1) 1782 DEALLOCATE(rpos, STAT = ierr) 1783 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating rpos', 1) 1784 DEALLOCATE(cpos, STAT = ierr) 1785 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating cpos', 1) 1786 DEALLOCATE(phase, STAT = ierr) 1787 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating phase', 1) 1788 DEALLOCATE(singval, STAT = ierr) 1789 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating singval', 1) 1790 DEALLOCATE(umat, STAT = ierr) 1791 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating umat', 1) 1792 DEALLOCATE(vtmat, STAT = ierr) 1793 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating vtmat', 1) 1794 DEALLOCATE(amat, STAT = ierr) 1795 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating amat', 1) 1796 ! 1797 IF (noncolin) THEN 1798 DEALLOCATE(piv_pos, STAT = ierr) 1799 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating piv_pos', 1) 1800 DEALLOCATE(piv_spin, STAT = ierr) 1801 IF (ierr /= 0) CALL errore('compute_amn_with_scdm', 'Error deallocating piv_spin', 1) 1802 ENDIF 1803 ! 1804 WRITE(stdout, *) 1805 WRITE(stdout, '(5x, a)') 'AMN calculated with SCDM' 1806 ! 1807 RETURN 1808 ! 1809 !----------------------------------------------------------------------- 1810 END SUBROUTINE compute_amn_with_scdm 1811 !----------------------------------------------------------------------- 1812 ! 1813 !----------------------------------------------------------------------- 1814 SUBROUTINE orient_gf_spinor(npw) 1815 !----------------------------------------------------------------------- 1816 !! 1817 !! FIXME 1818 !! 1819 USE kinds, ONLY : DP 1820 USE constants_epw, ONLY : czero, cone, eps6 1821 USE wvfct, ONLY : npwx 1822 USE wannierEPW, ONLY : gf, n_proj, spin_qaxis, spin_eig, gf_spinor 1823 ! 1824 IMPLICIT NONE 1825 ! 1826 INTEGER, INTENT(in) :: npw 1827 !! Number of plane waves 1828 ! 1829 ! Local variables 1830 LOGICAL :: spin_z_pos 1831 !! Detect if spin quantisation axis is along +z 1832 LOGICAL :: spin_z_neg 1833 !! Detect if spin quantisation axis is along -z 1834 INTEGER :: iw 1835 !! Counter on number of projections 1836 INTEGER :: ipol 1837 !! Index of spin-up/spin-down polarizations 1838 INTEGER :: istart, iend 1839 !! Starting/ending index on plane waves 1840 COMPLEX(KIND = DP) :: fac(2) 1841 !! Factor 1842 ! 1843 gf_spinor(:, :) = czero 1844 DO iw = 1, n_proj 1845 spin_z_pos = .FALSE. 1846 spin_z_neg = .FALSE. 1847 ! detect if spin quantisation axis is along z 1848 IF ((ABS(spin_qaxis(1, iw) - 0.0d0) < eps6) .AND. & 1849 (ABS(spin_qaxis(2, iw) - 0.0d0) < eps6) .AND. & 1850 (ABS(spin_qaxis(3, iw) - 1.0d0) < eps6)) THEN 1851 spin_z_pos = .TRUE. 1852 ELSEIF ((ABS(spin_qaxis(1, iw) - 0.0d0) < eps6) .AND. & 1853 (ABS(spin_qaxis(2, iw) - 0.0d0) < eps6) .AND. & 1854 (ABS(spin_qaxis(3, iw) + 1.0d0) < eps6)) THEN 1855 spin_z_neg = .TRUE. 1856 ENDIF 1857 IF (spin_z_pos .OR. spin_z_neg) THEN 1858 IF (spin_z_pos) THEN 1859 ipol = (3 - spin_eig(iw)) / 2 1860 ELSE 1861 ipol = (3 + spin_eig(iw)) / 2 1862 ENDIF 1863 istart = (ipol - 1) * npwx + 1 1864 iend = (ipol - 1) * npwx + npw 1865 gf_spinor(istart:iend, iw) = gf(1:npw, iw) 1866 ELSE 1867 IF (spin_eig(iw) == 1) THEN 1868 fac(1) = (1.0d0 / DSQRT(1 + spin_qaxis(3, iw))) & 1869 * (spin_qaxis(3, iw) + 1 ) * cone 1870 fac(2) = (1.0d0 / DSQRT(1 + spin_qaxis(3, iw))) & 1871 * CMPLX(spin_qaxis(1, iw), spin_qaxis(2, iw), KIND = DP) 1872 ELSE 1873 fac(1) = (1.0d0 / DSQRT(1 + spin_qaxis(3, iw))) & 1874 * ( spin_qaxis(3, iw) ) * cone 1875 fac(2) = (1.0d0 / DSQRT(1 - spin_qaxis(3, iw))) & 1876 * CMPLX(spin_qaxis(1, iw), spin_qaxis(2,iw), KIND = DP) 1877 ENDIF 1878 gf_spinor(1:npw, iw) = gf(1:npw, iw) * fac(1) 1879 gf_spinor(1 + npwx:npw + npwx, iw) = gf(1:npw, iw) * fac(2) 1880 ENDIF 1881 ENDDO 1882 ! 1883 RETURN 1884 ! 1885 !----------------------------------------------------------------------- 1886 END SUBROUTINE orient_gf_spinor 1887 !----------------------------------------------------------------------- 1888 ! 1889 !----------------------------------------------------------------------- 1890 SUBROUTINE compute_mmn_para() 1891 !----------------------------------------------------------------------- 1892 !! 1893 !! adapted from compute_mmn in pw2wannier90.f90 1894 !! parallelization on k-points has been added 1895 !! 10/2008 Jesse Noffsinger UC Berkeley 1896 !! 1897 USE kinds, ONLY : DP 1898 USE io_global, ONLY : stdout, meta_ionode 1899 USE io_files, ONLY : diropn, prefix 1900 USE wvfct, ONLY : nbnd, npw, npwx, g2kin 1901 USE wavefunctions, ONLY : evc, psic, psic_nc 1902 USE units_lr, ONLY : lrwfc, iuwfc 1903 USE fft_base, ONLY : dffts 1904 USE fft_interfaces, ONLY : fwfft, invfft 1905 USE klist, ONLY : nkstot, nks, igk_k 1906 USE klist_epw, ONLY : xk_all, xk_loc 1907 USE gvect, ONLY : g, ngm 1908 USE gvecw, ONLY : gcutw 1909 USE cell_base, ONLY : omega, tpiba, bg 1910 USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau 1911 USE uspp, ONLY : nkb, vkb 1912 USE uspp_param, ONLY : upf, lmaxq, nh, nhm 1913 USE becmod, ONLY : becp, calbec, allocate_bec_type, & 1914 deallocate_bec_type 1915 USE noncollin_module,ONLY : noncolin, npol 1916 USE spin_orb, ONLY : lspinorb 1917 USE wannierEPW, ONLY : m_mat, num_bands, nnb, iknum, g_kpb, kpb, ig_, & 1918 excluded_band, write_mmn, zerophase 1919 USE constants_epw, ONLY : czero, cone, twopi, zero 1920 USE io_var, ONLY : iummn 1921#if defined(__NAG) 1922 USE f90_unix_io, ONLY : flush 1923#endif 1924 USE mp, ONLY : mp_sum 1925 USE mp_global, ONLY : my_pool_id, npool, intra_pool_comm, inter_pool_comm 1926 USE kfold, ONLY : ktokpmq 1927 USE io_epw, ONLY : readwfc 1928 USE elph2, ONLY : nbndep 1929 ! 1930 IMPLICIT NONE 1931 ! 1932 CHARACTER (LEN=80) :: filmmn 1933 !! file containing M_mn(k,b) matrix 1934 LOGICAL :: any_uspp 1935 !! Check if uspp are used 1936 LOGICAL :: exst 1937 !! Does the file exist 1938 ! 1939 INTEGER :: ik 1940 !! Counter on k-points 1941 INTEGER :: ikp 1942 !! Index of k-point nearest neighbour 1943 INTEGER :: ib 1944 !! Counter on b-vectors 1945 INTEGER :: npwq 1946 !! Number of planes waves at k+b 1947 INTEGER :: m 1948 !! Counter over bands 1949 INTEGER :: n 1950 !! Counter over bands 1951 INTEGER :: ibnd_m 1952 !! Band index 1953 INTEGER :: ibnd_n 1954 !! Band index 1955 INTEGER :: ih, jh 1956 !! Counters over the beta functions per atomic type 1957 INTEGER :: ikb, jkb, ijkb0 1958 !! Indexes over beta functions 1959 INTEGER :: na 1960 !! Counter over atoms 1961 INTEGER :: nt 1962 !! Number of type of atoms 1963 INTEGER :: ind 1964 !! Counter over nearest neighbours of all k-points 1965 INTEGER :: nbt 1966 !! 1967 INTEGER :: ipool 1968 !! Index of current pool 1969 INTEGER :: nkq 1970 !! Index of k-point in the pool 1971 INTEGER :: nkq_abs 1972 !! Absolute index of k-point 1973 INTEGER :: ipol 1974 !! Index of spin-up/spin-down polarizations 1975 INTEGER :: istart, iend 1976 !! Starting/ending index on plane waves 1977 INTEGER :: ik_g 1978 !! Temporary index of k-point, ik_g = nkq_abs 1979 INTEGER :: ind0 1980 !! Starting index for k-point nearest neighbours in each pool 1981 INTEGER :: ierr 1982 !! Error status 1983 INTEGER, ALLOCATABLE :: igkq(:) 1984 !! 1985 REAL(KIND = DP) :: arg 1986 !! 2*pi*(k+b - k)*R 1987 COMPLEX(KIND = DP) :: mmn 1988 !! Element of M_mn matrix 1989 REAL(KIND = DP), ALLOCATABLE :: dxk(:, :) 1990 !! Temporary vector k+b - k 1991 REAL(KIND = DP), ALLOCATABLE :: qg(:) 1992 !! Square of dxk 1993 REAL(KIND = DP), ALLOCATABLE :: ylm(:, :) 1994 !! Spherical harmonics 1995 REAL(KIND = DP) :: zero_vect(3) 1996 !! Temporary zero vector 1997 REAL(KIND = DP) :: g_(3) 1998 !! Temporary vector G_k+b, g_(:) = g_kpb(:,ik,ib) 1999 COMPLEX(KIND = DP), EXTERNAL :: ZDOTC 2000 !! Scalar product of complex vectors 2001 COMPLEX(KIND = DP) :: phase1 2002 !! e^{i*2*pi*(k+b - k)*R} 2003 COMPLEX(KIND = DP), ALLOCATABLE :: phase(:) 2004 !! Phase 2005 COMPLEX(KIND = DP), ALLOCATABLE :: aux(:) 2006 !! Auxillary variable 2007 COMPLEX(KIND = DP), ALLOCATABLE :: aux_nc(:, :) 2008 !! NC auxillary variable 2009 COMPLEX(KIND = DP), ALLOCATABLE :: evcq(:, :) 2010 !! Wave functions psi_k+b 2011 COMPLEX(KIND = DP), ALLOCATABLE :: Mkb(:, :) 2012 !! Element of M_mn matrix 2013 COMPLEX(KIND = DP), ALLOCATABLE :: becp2(:, :) 2014 !! Beta functions 2015 COMPLEX(KIND = DP), ALLOCATABLE :: becp2_nc(:, :, :) 2016 !! Beta functions, noncolin 2017 COMPLEX(KIND = DP), ALLOCATABLE :: qb(:, :, :, :) 2018 !! Local variables for uspp 2019 COMPLEX(KIND = DP), ALLOCATABLE :: qgm(:) 2020 !! Local variables for uspp 2021 COMPLEX(KIND = DP), ALLOCATABLE :: qq_so(:, :, :, :) 2022 !! Local variables for uspp 2023 ! 2024 any_uspp = ANY(upf(:)%tvanp) 2025 ! 2026 ALLOCATE(phase(dffts%nnr), STAT = ierr) 2027 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating phase', 1) 2028 ALLOCATE(igkq(npwx), STAT = ierr) 2029 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating igkq', 1) 2030 ALLOCATE(evcq(npol * npwx, nbnd), STAT = ierr) 2031 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating evcq', 1) 2032 ! 2033 IF (noncolin) THEN 2034 ALLOCATE(aux_nc(npwx, npol), STAT = ierr) 2035 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating aux_nc', 1) 2036 ELSE 2037 ALLOCATE(aux(npwx), STAT = ierr) 2038 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating aux', 1) 2039 ENDIF 2040 ! 2041 ALLOCATE(m_mat(num_bands, num_bands, nnb, iknum), STAT = ierr) 2042 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating m_mat', 1) 2043 ! 2044 ! close all the wfc files to allow access for each pool to all wfs 2045 CLOSE(UNIT = iuwfc, STATUS = 'keep') 2046 ! 2047 WRITE(stdout, *) 2048 WRITE(stdout, '(5x, a)') 'MMN' 2049 ! 2050 zero_vect(:) = zero 2051 m_mat(:, :, :, :) = czero 2052 ! 2053 ! USPP 2054 ! 2055 IF (any_uspp) THEN 2056 CALL allocate_bec_type(nkb, nbnd, becp) 2057 IF (noncolin) THEN 2058 ALLOCATE(becp2_nc(nkb, 2, nbnd), STAT = ierr) 2059 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating becp2_nc', 1) 2060 ELSE 2061 ALLOCATE(becp2(nkb, nbnd), STAT = ierr) 2062 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating becp2', 1) 2063 ENDIF 2064 ! 2065 ! qb is FT of Q(r) 2066 ! 2067 nbt = nnb * iknum 2068 ! 2069 ALLOCATE(qg(nbt), STAT = ierr) 2070 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating qg', 1) 2071 ALLOCATE(dxk(3, nbt), STAT = ierr) 2072 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating dxk', 1) 2073 ! 2074 ind = 0 2075 DO ik = 1, iknum ! loop over k-points 2076 DO ib = 1, nnb ! loop over nearest neighbours for each k-point 2077 ind = ind + 1 2078 ikp = kpb(ik, ib) 2079 ! 2080 g_(:) = REAL(g_kpb(:, ik, ib)) 2081 ! bring g_ from crystal to cartesian 2082 CALL cryst_to_cart(1, g_, bg, 1) 2083 dxk(:, ind) = xk_all(:, ikp) + g_(:) - xk_all(:, ik) 2084 qg(ind) = SUM(dxk(:, ind) * dxk(:, ind)) 2085 ENDDO 2086 ENDDO 2087 ! 2088 ALLOCATE(ylm(nbt, lmaxq * lmaxq), STAT = ierr) 2089 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating ylm', 1) 2090 ALLOCATE(qgm(nbt), STAT = ierr) 2091 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating qgm', 1) 2092 ALLOCATE(qb(nhm, nhm, ntyp, nbt), STAT = ierr) 2093 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating qb', 1) 2094 ALLOCATE(qq_so(nhm, nhm, 4, ntyp), STAT = ierr) 2095 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating qq_so', 1) 2096 ! 2097 ! get spherical harmonics ylm 2098 CALL ylmr2(lmaxq * lmaxq, nbt, dxk, qg, ylm) 2099 qg(:) = DSQRT(qg(:)) * tpiba 2100 ! 2101 DO nt = 1, ntyp 2102 IF (upf(nt)%tvanp) THEN 2103 DO ih = 1, nh(nt) 2104 DO jh = 1, nh(nt) 2105 CALL qvan2(nbt, ih, jh, nt, qg, qgm, ylm) 2106 qb(ih, jh, nt, 1:nbt) = omega * qgm(1:nbt) 2107 ENDDO 2108 ENDDO 2109 ENDIF 2110 ENDDO 2111 ! 2112 DEALLOCATE(qg, STAT = ierr) 2113 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating qg', 1) 2114 DEALLOCATE(qgm, STAT = ierr) 2115 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating qgm', 1) 2116 DEALLOCATE(ylm, STAT = ierr) 2117 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating ylm', 1) 2118 ! 2119 ENDIF 2120 ! 2121 ALLOCATE(Mkb(nbnd, nbnd), STAT = ierr) 2122 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error allocating Mkb', 1) 2123 ! 2124#if defined(__MPI) 2125 WRITE(stdout, '(6x, a, i5, a, i4, a)') 'k points = ', iknum, ' in ', npool, ' pools' 2126#endif 2127 ! 2128 ! returns in-pool index nkq and absolute index nkq_abs of first k-point in this pool 2129 CALL ktokpmq( xk_loc(:,1), zero_vect, +1, ipool, nkq, nkq_abs ) 2130 ind0 = (nkq_abs - 1) * nnb 2131 ! 2132 ind = ind0 2133 DO ik = 1, nks 2134 ! 2135 ! returns in-pool index nkq and absolute index nkq_abs of xk 2136 CALL ktokpmq(xk_loc(:, ik), zero_vect, +1, ipool, nkq, nkq_abs) 2137 ik_g = nkq_abs 2138 ! 2139 WRITE(stdout, '(5x, i8, " of ", i4, a)') ik , nks, ' on ionode' 2140 FLUSH(stdout) 2141 ! 2142 ! read wfc at k 2143 CALL readwfc(my_pool_id + 1, ik, evc) 2144 ! 2145 ! sorts k+G vectors in order of increasing magnitude, up to ecut 2146 CALL gk_sort(xk_loc(1, ik), ngm, g, gcutw, npw, igk_k(1, ik), g2kin) 2147 ! 2148 ! USPP 2149 ! 2150 IF (any_uspp) THEN 2151 CALL init_us_2(npw, igk_k(1, ik), xk_loc(1, ik), vkb) 2152 ! below we compute the product of beta functions with |psi> 2153 CALL calbec(npw, vkb, evc, becp) 2154 ENDIF 2155 ! 2156 ! 2157 DO ib = 1, nnb ! loop on finite diff b-vectors 2158 ind = ind + 1 2159 ! 2160 ikp = kpb(ik_g, ib) 2161 ! 2162 CALL ktokpmq(xk_all(:, ikp), zero_vect, +1, ipool, nkq, nkq_abs) 2163 ! 2164 ! read wfc at k+b 2165 CALL readwfc(ipool, nkq, evcq) 2166 ! 2167 ! sorts k+G vectors in order of increasing magnitude, up to ecut 2168 CALL gk_sort(xk_all(1, ikp), ngm, g, gcutw, npwq, igkq, g2kin) 2169 ! 2170 ! compute the phase 2171 IF (.NOT. zerophase(ik_g, ib)) THEN 2172 phase(:) = czero 2173 IF (ig_(ik_g, ib) > 0) phase(dffts%nl(ig_(ik_g, ib))) = cone 2174 CALL invfft('Wave', phase, dffts) 2175 ENDIF 2176 ! 2177 ! USPP 2178 ! 2179 IF (any_uspp) THEN 2180 CALL init_us_2(npwq, igkq, xk_all(1, ikp), vkb) 2181 ! below we compute the product of beta functions with |psi> 2182 IF (noncolin) THEN 2183 CALL calbec(npwq, vkb, evcq, becp2_nc) 2184 ! 2185 IF (lspinorb) THEN 2186 qq_so(:, :, :, :) = czero 2187 CALL transform_qq_so(qb(:, :, :, ind), qq_so) 2188 ENDIF 2189 ! 2190 ELSE 2191 CALL calbec(npwq, vkb, evcq, becp2) 2192 ENDIF 2193 ENDIF 2194 ! 2195 Mkb(:, :) = czero 2196 ! 2197 IF (any_uspp) THEN 2198 ijkb0 = 0 2199 DO nt = 1, ntyp 2200 IF (upf(nt)%tvanp) THEN 2201 DO na = 1, nat 2202 ! 2203 arg = twopi * DOT_PRODUCT(dxk(:, ind), tau(:, na)) 2204 phase1 = CMPLX(COS(arg), -SIN(arg), KIND = DP) 2205 ! 2206 IF (ityp(na) == nt) THEN 2207 DO jh = 1, nh(nt) 2208 jkb = ijkb0 + jh 2209 DO ih = 1, nh(nt) 2210 ikb = ijkb0 + ih 2211 ! 2212 DO m = 1, nbnd 2213 IF (excluded_band(m)) CYCLE 2214 IF (noncolin) THEN 2215 DO n = 1, nbnd 2216 IF (excluded_band(n)) CYCLE 2217 IF (lspinorb) THEN 2218 Mkb(m, n) = Mkb(m, n) + phase1 * & 2219 (qq_so(ih, jh, 1, nt) * CONJG(becp%nc(ikb, 1, m)) * becp2_nc(jkb, 1, n) + & 2220 qq_so(ih, jh, 2, nt) * CONJG(becp%nc(ikb, 1, m)) * becp2_nc(jkb, 2, n) + & 2221 qq_so(ih, jh, 3, nt) * CONJG(becp%nc(ikb, 2, m)) * becp2_nc(jkb, 1, n) + & 2222 qq_so(ih, jh, 4, nt) * CONJG(becp%nc(ikb, 2, m)) * becp2_nc(jkb, 2, n)) 2223 ELSE 2224 Mkb(m, n) = Mkb(m, n) + phase1 * qb(ih, jh, nt, ind) * & 2225 (CONJG(becp%nc(ikb, 1, m)) * becp2_nc(jkb, 1, n) + & 2226 CONJG(becp%nc(ikb, 2, m)) * becp2_nc(jkb, 2, n)) 2227 ENDIF 2228 ENDDO 2229 ELSE 2230 DO n = 1, nbnd 2231 IF (excluded_band(n)) CYCLE 2232 Mkb(m, n) = Mkb(m, n) + phase1 * qb(ih, jh, nt, ind) * & 2233 CONJG(becp%k(ikb, m)) * becp2(jkb, n) 2234 ENDDO 2235 ENDIF 2236 ENDDO ! m 2237 ENDDO ! ih 2238 ENDDO ! jh 2239 ijkb0 = ijkb0 + nh(nt) 2240 ENDIF ! ityp 2241 ENDDO ! nat 2242 ELSE ! tvanp 2243 DO na = 1, nat 2244 IF (ityp(na) == nt) ijkb0 = ijkb0 + nh(nt) 2245 ENDDO 2246 ENDIF ! tvanp 2247 ENDDO ! ntyp 2248 ENDIF ! any_uspp 2249 ! 2250 DO m = 1, nbnd ! loop on band m 2251 IF (excluded_band(m)) CYCLE 2252 ! 2253 IF (noncolin) THEN 2254 psic_nc(:, :) = czero 2255 DO ipol = 1, 2 !npol 2256 istart = (ipol - 1) * npwx + 1 2257 iend = (ipol - 1) * npwx + npw 2258 psic_nc(dffts%nl(igk_k(1:npw, ik)), ipol) = evc(istart:iend, m) 2259 IF (.NOT. zerophase(ik_g, ib)) THEN 2260 CALL invfft('Wave', psic_nc(:, ipol), dffts) 2261 psic_nc(1:dffts%nnr, ipol) = psic_nc(1:dffts%nnr, ipol) * & 2262 phase(1:dffts%nnr) 2263 CALL fwfft('Wave', psic_nc(:, ipol), dffts) 2264 ENDIF 2265 aux_nc(1:npwq, ipol) = psic_nc(dffts%nl(igkq(1:npwq)), ipol) 2266 ENDDO 2267 ELSE 2268 psic(:) = czero 2269 psic(dffts%nl(igk_k(1:npw, ik))) = evc(1:npw, m) 2270 IF (.NOT. zerophase(ik_g,ib)) THEN 2271 CALL invfft('Wave', psic, dffts) 2272 psic(1:dffts%nnr) = psic(1:dffts%nnr) * phase(1:dffts%nnr) 2273 CALL fwfft('Wave', psic, dffts) 2274 ENDIF 2275 aux(1:npwq) = psic(dffts%nl(igkq(1:npwq))) 2276 ENDIF 2277 ! aa = 0.d0 2278 ! 2279 ! Mkb(m,n) = Mkb(m,n) + \sum_{ijI} qb_{ij}^I * e^-i(b*tau_I) 2280 ! <psi_m,k1| beta_i,k1 > < beta_j,k2 | psi_n,k2 > 2281 ! 2282 IF (noncolin) THEN 2283 DO n = 1, nbnd 2284 IF (excluded_band(n)) CYCLE 2285 mmn = czero 2286 mmn = mmn + ZDOTC(npwq, aux_nc(1, 1), 1, evcq( 1, n), 1) & 2287 + ZDOTC(npwq, aux_nc(1, 2), 1, evcq(1 + npwx, n), 1) 2288 CALL mp_sum(mmn, intra_pool_comm) 2289 Mkb(m, n) = mmn + Mkb(m, n) 2290 ! aa = aa + ABS(mmn)**2 2291 ENDDO ! n 2292 ELSE ! scalar wavefunction 2293 DO n = 1, nbnd 2294 IF (excluded_band(n)) CYCLE 2295 mmn = ZDOTC(npwq, aux, 1, evcq(1, n), 1) 2296 CALL mp_sum(mmn, intra_pool_comm) 2297 Mkb(m, n) = mmn + Mkb(m, n) 2298 ! aa = aa + ABS(mmn)**2 2299 ENDDO ! n 2300 ENDIF ! scalar wavefunction 2301 ENDDO ! m 2302 ! 2303 ibnd_n = 0 2304 DO n = 1, nbnd 2305 IF (excluded_band(n)) CYCLE 2306 ibnd_n = ibnd_n + 1 2307 ! 2308 ibnd_m = 0 2309 DO m = 1, nbnd 2310 IF (excluded_band(m)) CYCLE 2311 ibnd_m = ibnd_m + 1 2312 ! 2313 m_mat(ibnd_m, ibnd_n, ib, ik_g) = Mkb(m, n) 2314 ENDDO ! m 2315 ENDDO ! n 2316 ! 2317 ENDDO ! ib 2318 ENDDO ! ik 2319 ! 2320 CALL mp_sum(m_mat, inter_pool_comm) 2321 ! 2322 ! RM - write mmn to file (file needed with vme = true) 2323 IF (meta_ionode) THEN 2324 write_mmn = .TRUE. 2325 IF (write_mmn) THEN 2326 ! 2327 filmmn = TRIM(prefix)//'.mmn' 2328 OPEN(UNIT = iummn, FILE = filmmn, FORM = 'formatted') 2329 DO ik = 1, nkstot 2330 DO ib = 1, nnb 2331 ! 2332 DO n = 1, nbndep 2333 DO m = 1, nbndep 2334 WRITE(iummn,*) m_mat(m, n, ib, ik) 2335 ENDDO ! m 2336 ENDDO ! n 2337 ! 2338 ENDDO ! ib 2339 ENDDO ! ik 2340 ! 2341 CLOSE(iummn) 2342 ENDIF !write_mmn 2343 ENDIF 2344 ! 2345 DEALLOCATE(Mkb, STAT = ierr) 2346 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating Mkb', 1) 2347 DEALLOCATE(phase, STAT = ierr) 2348 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating phase', 1) 2349 DEALLOCATE(evcq, STAT = ierr) 2350 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating evcq', 1) 2351 DEALLOCATE(igkq, STAT = ierr) 2352 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating igkq', 1) 2353 IF (noncolin) THEN 2354 DEALLOCATE(aux_nc, STAT = ierr) 2355 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating aux_nc', 1) 2356 ELSE 2357 DEALLOCATE(aux, STAT = ierr) 2358 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating aux', 1) 2359 ENDIF 2360 ! 2361 IF (any_uspp) THEN 2362 DEALLOCATE(dxk, STAT = ierr) 2363 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating dxk', 1) 2364 DEALLOCATE(qb, STAT = ierr) 2365 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating qb', 1) 2366 DEALLOCATE(qq_so, STAT = ierr) 2367 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating qq_so', 1) 2368 CALL deallocate_bec_type(becp) 2369 IF (noncolin) THEN 2370 DEALLOCATE(becp2_nc, STAT = ierr) 2371 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating becp2_nc', 1) 2372 ELSE 2373 DEALLOCATE(becp2, STAT = ierr) 2374 IF (ierr /= 0) CALL errore('compute_mmn_para', 'Error deallocating becp2', 1) 2375 ENDIF 2376 ENDIF 2377 ! 2378 WRITE(stdout, '(5x, a)') 'MMN calculated' 2379 ! 2380 ! reopen wfc here, leaving UNIT = 20 in the same state 2381 iuwfc = 20 2382 CALL diropn(iuwfc, 'wfc', lrwfc, exst) 2383 ! 2384 RETURN 2385 ! 2386 !------------------------------------------------------------------------ 2387 END SUBROUTINE compute_mmn_para 2388 !------------------------------------------------------------------------ 2389 ! 2390 !----------------------------------------------------------------------- 2391 SUBROUTINE compute_pmn_para() 2392 !----------------------------------------------------------------------- 2393 ! 2394 !! 2395 !! Computes dipole matrix elements. 2396 !! This can be used to compute the velocity in the local approximation. 2397 !! The commutator with the non-local pseudopotetial is neglected. 2398 !! 2399 !! 06/2010 Jesse Noffsinger: adapted from compute_amn_para 2400 !! 08/2016 Samuel Ponce: adapted to work with SOC 2401 !! 2402 ! 2403 USE kinds, ONLY : DP 2404 USE io_global, ONLY : stdout 2405 USE mp, ONLY : mp_sum 2406 USE mp_global, ONLY : my_pool_id 2407 USE klist, ONLY : nks, igk_k 2408 USE klist_epw, ONLY : xk_loc 2409 USE wvfct, ONLY : nbnd, npw, npwx, g2kin 2410 USE gvecw, ONLY : gcutw 2411 USE wavefunctions, ONLY : evc 2412 USE gvect, ONLY : g, ngm 2413 USE cell_base, ONLY : tpiba 2414 USE noncollin_module,ONLY : noncolin 2415 USE elph2, ONLY : dmec, ibndstart, ibndend, nbndep 2416 USE constants_epw, ONLY : czero 2417 USE uspp_param, ONLY : upf 2418 USE becmod, ONLY : becp, deallocate_bec_type, allocate_bec_type 2419 USE uspp, ONLY : nkb 2420 USE wannierEPW, ONLY : n_wannier 2421 USE kfold, ONLY : ktokpmq 2422 USE io_epw, ONLY : readwfc 2423 ! 2424 IMPLICIT NONE 2425 ! 2426 LOGICAL :: any_uspp 2427 !! Check if USPP is present 2428 ! 2429 INTEGER :: ik 2430 !! Counter on k-point 2431 INTEGER :: ibnd, jbnd 2432 !! Counter on bands 2433 INTEGER :: ig 2434 !! Counter on G vector 2435 INTEGER :: ierr 2436 !! Error status 2437 ! 2438 COMPLEX(KIND = DP) :: dipole_aux(3, nbnd, nbnd) 2439 !! Auxilary dipole 2440 COMPLEX(KIND = DP) :: caux 2441 !! Wavefunction squared 2442 ! 2443 any_uspp = ANY(upf(:)%tvanp) 2444 ! 2445 IF (any_uspp) CALL errore('pw2wan90epw', & 2446 'dipole matrix calculation not implimented with USPP - set vme=.TRUE.',1) 2447 ! 2448 ALLOCATE(dmec(3, nbndep, nbndep, nks), STAT = ierr) 2449 IF (ierr /= 0) CALL errore('compute_pmn_para', 'Error allocating dmec', 1) 2450 ! 2451 ! initialize 2452 dmec(:, :, :, :) = czero 2453 dipole_aux(:, :, :) = czero 2454 ! 2455 IF (any_uspp) THEN 2456 CALL allocate_bec_type(nkb, n_wannier, becp) 2457 ENDIF 2458 ! 2459 ! computes velocity dmec = v_mn(\alpha,k) in the local approximation 2460 ! acording to [Eqn. 60 of Comp. Phys. Commun. 209, 116 (2016)] 2461 ! v_mn(\alpha,k) = k \delta_mn + \sum_G * CONJG(c_mk(G)) * c_nk(G) * G 2462 ! 2463 DO ik = 1, nks 2464 ! 2465 ! read wfc for the given kpt 2466 CALL readwfc(my_pool_id + 1, ik, evc) 2467 ! 2468 ! sorts k+G vectors in order of increasing magnitude, up to ecut 2469 CALL gk_sort(xk_loc(:, ik), ngm, g, gcutw, npw, igk_k(:, ik), g2kin) 2470 ! 2471 dipole_aux = czero 2472 DO jbnd = ibndstart, ibndend 2473 DO ibnd = ibndstart, ibndend 2474 ! 2475 IF (ibnd == jbnd) CYCLE 2476 ! 2477 ! taken from PP/epsilon.f90 subroutine dipole_calc 2478 DO ig = 1, npw 2479 IF (igk_k(ig, ik) > SIZE(g, 2) .OR. igk_k(ig, ik) < 1) CYCLE 2480 ! 2481 caux = CONJG(evc(ig, ibnd)) * evc(ig, jbnd) 2482 ! 2483 IF (noncolin) THEN 2484 ! 2485 caux = caux + CONJG(evc(ig + npwx, ibnd)) * evc(ig + npwx, jbnd) 2486 ! 2487 ENDIF 2488 ! 2489 dipole_aux(:, ibnd, jbnd) = dipole_aux(:, ibnd, jbnd) + & 2490 (g(:, igk_k(ig, ik))) * caux 2491 ! 2492 ENDDO 2493 ! 2494 ENDDO !bands i 2495 ENDDO ! bands j 2496 ! 2497 ! metal diagonal part 2498 DO ibnd = ibndstart, ibndend 2499 DO ig = 1, npw 2500 IF (igk_k(ig, ik) > SIZE(g, 2) .OR. igk_k(ig, ik) < 1) CYCLE 2501 ! 2502 caux = CONJG(evc(ig, ibnd)) * evc(ig, ibnd) 2503 ! 2504 IF (noncolin) THEN 2505 ! 2506 caux = caux + CONJG(evc(ig + npwx, ibnd)) * evc(ig + npwx, ibnd) 2507 ! 2508 ENDIF 2509 ! 2510 dipole_aux(:, ibnd, ibnd) = dipole_aux(:, ibnd, ibnd) + & 2511 (g(:, igk_k(ig, ik)) + xk_loc(:, ik)) * caux 2512 ! 2513 ENDDO 2514 ENDDO 2515 ! need to divide by 2pi/a to fix the units 2516 DO jbnd = ibndstart, ibndend 2517 DO ibnd = ibndstart, ibndend 2518 dmec(:, ibnd-ibndstart+1, jbnd-ibndstart+1, ik) = dipole_aux(:, ibnd, jbnd) * tpiba 2519 ENDDO 2520 ENDDO 2521 ! 2522 ENDDO ! k-points 2523 ! 2524 IF (any_uspp) CALL deallocate_bec_type(becp) 2525 ! 2526 WRITE(stdout, '(/5x, a)') 'Dipole matrix elements calculated' 2527 WRITE(stdout, *) 2528 !DBSP 2529 !WRITE(stdout, *) 'dmec ', SUM(dmec) 2530 !IF (meta_ionode) THEN 2531 ! WRITE(stdout, *) 'dmec(:, :, :, 1) ',sum(dmec(:, :, :, 1)) 2532 ! WRITE(stdout, *) 'dmec(:, :, :, 2) ',sum(dmec(:, :, :, 2)) 2533 !ENDIF 2534 ! 2535 RETURN 2536 ! 2537 !------------------------------------------------------------------------ 2538 END SUBROUTINE compute_pmn_para 2539 !----------------------------------------------------------------------- 2540 ! 2541 !----------------------------------------------------------------------- 2542 SUBROUTINE write_filukk 2543 !----------------------------------------------------------------------- 2544 ! 2545 ! Here we compute and write out the final ukk matrix which is used by 2546 ! epw.x to localize the electron wavefuctions (and therefore the ep-matrix 2547 ! elements) 2548 ! 10/2008 Jesse Noffsinger UC Berkeley 2549 ! 07/2010 Fixed the rotation for ndimwin when lower bands are not included 2550 ! 2551 USE kinds, ONLY : DP 2552 USE io_var, ONLY : iunukk 2553 USE wvfct, ONLY : nbnd 2554 USE wannierEPW, ONLY : n_wannier, iknum, u_mat, u_mat_opt, lwindow, & 2555 excluded_band, num_bands, wann_centers 2556 USE epwcom, ONLY : filukk 2557 USE constants_epw,ONLY : czero, bohr 2558 USE io_global, ONLY : meta_ionode 2559 USE cell_base, ONLY : alat 2560 USE elph2, ONLY : nbndep, ibndstart, ibndend 2561 ! 2562 IMPLICIT NONE 2563 ! 2564 INTEGER :: ik 2565 !! Counter of k-point index 2566 INTEGER :: ibnd 2567 !! Counter on band index 2568 INTEGER :: ibnd1 2569 !! Band index 2570 INTEGER :: iw 2571 !! Counter on number of Wannier functions 2572 INTEGER :: ierr 2573 !! Error status 2574 INTEGER :: ndimwin(iknum) 2575 !! Number of bands within outer window at each k-point 2576 ! 2577 COMPLEX(KIND = DP), ALLOCATABLE :: u_kc(:, :, :) 2578 !! Rotation matrix 2579 ! 2580 IF (meta_ionode) THEN 2581 ! 2582 ndimwin(:) = 0 2583 DO ik = 1, iknum 2584 DO ibnd = 1, num_bands 2585 IF (lwindow(ibnd, ik)) ndimwin(ik) = ndimwin(ik) + 1 2586 ENDDO 2587 ENDDO 2588 ! 2589 ! get the final rotation matrix, which is the product of the optimal 2590 ! subspace and the rotation among the n_wannier wavefunctions 2591 ! 2592 ALLOCATE(u_kc(nbndep, n_wannier, iknum), STAT = ierr) 2593 IF (ierr /= 0) CALL errore('write_filukk', 'Error allocating u_kc', 1) 2594 u_kc(:, :, :) = czero 2595 ! 2596 DO ik = 1, iknum 2597 ! 2598 u_kc(1:ndimwin(ik), 1:n_wannier, ik) = & 2599 MATMUL(u_mat_opt(1:ndimwin(ik), :, ik), u_mat(:, 1:n_wannier, ik)) 2600 ! 2601 ENDDO 2602 ! 2603 OPEN(UNIT = iunukk, FILE = filukk, FORM = 'formatted') 2604 ! 2605 WRITE(iunukk, *) ibndstart, ibndend 2606 ! 2607 DO ik = 1, iknum 2608 DO ibnd = 1, nbndep 2609 DO iw = 1, n_wannier 2610 WRITE(iunukk, *) u_kc(ibnd, iw, ik) 2611 ENDDO 2612 ENDDO 2613 ENDDO 2614 ! 2615 ! needs also lwindow when disentanglement is used 2616 ! 2617 DO ik = 1, iknum 2618 DO ibnd = 1, nbndep 2619 WRITE(iunukk, *) lwindow(ibnd, ik) 2620 ENDDO 2621 ENDDO 2622 ! 2623 DO ibnd = 1, nbnd 2624 WRITE(iunukk, *) excluded_band(ibnd) 2625 ENDDO 2626 ! 2627 ! Now write the Wannier centers to files 2628 DO iw = 1, n_wannier 2629 ! SP : Need more precision other WS are not determined properly. 2630 !WRITE (iuukk, '(3f12.8)') wann_centers(:, iw) / alat / bohr 2631 WRITE (iunukk, '(3E22.12)') wann_centers(:, iw) / alat / bohr 2632 ENDDO 2633 ! 2634 CLOSE(iunukk) 2635 ! 2636 DEALLOCATE(u_kc, STAT = ierr) 2637 IF (ierr /= 0) CALL errore('write_filukk', 'Error deallocating u_kc', 1) 2638 ! 2639 ENDIF 2640 ! 2641 RETURN 2642 ! 2643 !------------------------------------------------------------------------ 2644 END SUBROUTINE write_filukk 2645 !----------------------------------------------------------------------- 2646 ! 2647 !----------------------------------------------------------------------- 2648 SUBROUTINE phases_a_m 2649 !----------------------------------------------------------------------- 2650 !! 2651 !! We will set phases here on the matrices. It should not affect 2652 !! the spreads and centers found in w90, but it will leave 2653 !! u_mat_opt and u_mat to reflect the known phases. 2654 !! 2655 ! 2656 USE kinds, ONLY : DP 2657 USE mp_global, ONLY : inter_pool_comm 2658 USE mp, ONLY : mp_sum 2659 USE klist, ONLY : nkstot, nks 2660 USE klist_epw, ONLY : xk_loc 2661 USE wvfct, ONLY : nbnd 2662 USE wannierEPW, ONLY : a_mat, m_mat, n_wannier, n_proj, & 2663 nnb, kpb, iknum, excluded_band 2664 USE elph2, ONLY : umat, umat_all 2665 USE constants_epw, ONLY : czero, cone, zero 2666 USE kfold, ONLY : ktokpmq 2667 ! 2668 IMPLICIT NONE 2669 ! 2670 INTEGER :: ik 2671 !! Counter on k-points 2672 INTEGER :: ikb 2673 !! Index of k-point nearest neighbour 2674 INTEGER :: ib 2675 !! Counter on b-vectors 2676 INTEGER :: ipool 2677 !! Index of current pool 2678 INTEGER :: nkq 2679 !! Index of k-point in the pool 2680 INTEGER :: nkq_abs 2681 !! Absolute index of k-point 2682 INTEGER :: ik_g 2683 !! Temporary index of k-point, ik_g = nkq_abs 2684 INTEGER :: m, n 2685 !! Counter over bands 2686 INTEGER :: ibnd_m, ibnd_n 2687 !! Band index 2688 INTEGER :: iw 2689 !! Counter on number of projections 2690 INTEGER :: ierr 2691 !! Error status 2692 ! 2693 REAL(KIND = DP) :: zero_vect(3) 2694 !! Temporary zero vector 2695 ! 2696 COMPLEX(KIND = DP), ALLOCATABLE :: a_mat_tmp(:, :, :) 2697 !! Temporary a_mat matrices 2698 COMPLEX(KIND = DP), ALLOCATABLE :: a_mat_tmp1(:, :, :) 2699 !! Temporary a_mat matrices 2700 COMPLEX(KIND = DP), ALLOCATABLE :: m_mn_tmp1(:, :) 2701 !! Temporary m_mat matrices 2702 COMPLEX(KIND = DP), ALLOCATABLE :: m_mn_tmp2(:, :) 2703 !! Temporary m_mat matrices 2704 COMPLEX(KIND = DP), ALLOCATABLE :: m_mn_tmp3(:, :, :, :) 2705 !! Temporary m_mat matrices 2706 COMPLEX(KIND = DP), ALLOCATABLE :: m_mat_tmp(:, :, :, :) 2707 !! Temporary m_mat matrices 2708 ! 2709 ! RM: Band-dimension of a_mat and m_mat is num_bands while that of 2710 ! umat and umat_all is nbnd. This causes a problem if exclude_bands 2711 ! is used because num_bands = nbnd - nexband 2712 ! a_mat(num_bands,n_wannier,nkstot) in compute_amn_para 2713 ! m_mat(num_bands,n_wannier,nnb,nkstot) in compute_mmn_para 2714 ! umat(nbnd,nbnd,nks) and umat_mat(nbnd,nbnd,nkstot) in setphases_wrap 2715 ! 2716 ALLOCATE(a_mat_tmp(nbnd, n_wannier, iknum), STAT = ierr) 2717 IF (ierr /= 0) CALL errore('phases_a_m', 'Error allocating a_mat_tmp', 1) 2718 ALLOCATE(a_mat_tmp1(nbnd, n_wannier, iknum), STAT = ierr) 2719 IF (ierr /= 0) CALL errore('phases_a_m', 'Error allocating a_mat_tmp1', 1) 2720 ALLOCATE(m_mat_tmp(nbnd, nbnd, nnb, iknum), STAT = ierr) 2721 IF (ierr /= 0) CALL errore('phases_a_m', 'Error allocating m_mat_tmp', 1) 2722 ALLOCATE(m_mn_tmp1(nbnd, nbnd), STAT = ierr) 2723 IF (ierr /= 0) CALL errore('phases_a_m', 'Error allocating m_mn_tmp1', 1) 2724 ALLOCATE(m_mn_tmp2(nbnd, nbnd), STAT = ierr) 2725 IF (ierr /= 0) CALL errore('phases_a_m', 'Error allocating m_mn_tmp2', 1) 2726 ALLOCATE(m_mn_tmp3(nbnd, nbnd, nnb, iknum), STAT = ierr) 2727 IF (ierr /= 0) CALL errore('phases_a_m', 'Error allocating m_mn_tmp3', 1) 2728 ! 2729 ! zero all temporary/work quantities 2730 ! 2731 zero_vect(:) = zero 2732 a_mat_tmp(:, :, :) = czero 2733 a_mat_tmp1(:, :, :) = czero 2734 m_mn_tmp1(:, :) = czero 2735 m_mn_tmp2(:, :) = czero 2736 m_mn_tmp3(:, :, :, :) = czero 2737 m_mat_tmp(:, :, :, :) = czero 2738 ! 2739 ! full size a_mat_tmp1 and m_mat_tmp matrices to nbnd bands 2740 ! 2741 DO ik = 1, nkstot 2742 DO iw = 1, n_proj 2743 ibnd_m = 0 2744 DO m = 1, nbnd 2745 IF (excluded_band(m)) CYCLE 2746 ibnd_m = ibnd_m + 1 2747 a_mat_tmp1(m, iw, ik) = a_mat(ibnd_m, iw, ik) 2748 ENDDO 2749 ENDDO 2750 ! 2751 DO ib = 1, nnb 2752 ibnd_n = 0 2753 DO n = 1, nbnd 2754 IF (excluded_band(n)) CYCLE 2755 ibnd_n = ibnd_n + 1 2756 ibnd_m = 0 2757 DO m = 1, nbnd 2758 IF (excluded_band(m)) CYCLE 2759 ibnd_m = ibnd_m + 1 2760 m_mat_tmp(m, n, ib, ik) = m_mat(ibnd_m, ibnd_n, ib, ik) 2761 ENDDO 2762 ENDDO 2763 ENDDO 2764 ENDDO 2765 ! 2766 DO ik = 1, nks 2767 ! 2768 ! returns in-pool index nkq and absolute index nkq_abs of xk 2769 CALL ktokpmq(xk_loc(:, ik), zero_vect, +1, ipool, nkq, nkq_abs) 2770 ik_g = nkq_abs 2771 ! 2772 ! GF_n are the guiding functions which are our initial guesses 2773 ! Amn(k) = <psi_k,m|GF_n>. 2774 ! We want U(k)^\dagger<psi_k,m|GF_m> 2775 ! 2776 ! CALL ZGEMM('c', 'n', nbnd, n_wannier, nbnd, cone, umat(:, :, ik), & 2777 ! nbnd, a_mat(:, :, ik_g), nbnd, czero, a_mat_tmp(:, :, ik_g), nbnd) 2778 CALL ZGEMM('c', 'n', nbnd, n_wannier, nbnd, cone, umat(:, :, ik), & 2779 nbnd, a_mat_tmp1(:, :, ik_g), nbnd, czero, a_mat_tmp(:, :, ik_g), nbnd) 2780 ! 2781 DO ib = 1,nnb 2782 ikb = kpb(ik_g, ib) 2783 ! 2784 ! Mmn(k,k+b) = <psi_k_m| psi_(k+b)_n> so we need 2785 ! (U(k)^\dagger <psi_k_m| ) * (|psi_k+b_n> U(k+b) 2786 ! = U(k)^\dagger (M_mn) = m_mat_tmp, 2787 ! Mmn(k,k+b)' = m_mat_tmp*U(k+b) 2788 ! 2789 ! CALL ZGEMM('c', 'n', nbnd, nbnd, nbnd, cone, umat(:, :, ik), & 2790 ! nbnd, m_mat(:, :, ib, ik_g), nbnd, czero, m_mn_tmp1(:, :), nbnd) 2791 CALL ZGEMM('c', 'n', nbnd, nbnd, nbnd, cone, umat(:, :, ik), & 2792 nbnd, m_mat_tmp(:, :, ib, ik_g), nbnd, czero, m_mn_tmp1(:, :), nbnd) 2793 CALL ZGEMM('n', 'n', nbnd, nbnd, nbnd, cone, m_mn_tmp1(:, :), & 2794 nbnd, umat_all(:, :, ikb), nbnd, czero, m_mn_tmp2(:, :), nbnd) 2795 ! 2796 ! m_mn_tmp1 = MATMUL(CONJG(TRANSPOSE(umat(:, :, ik))), m_mat(:, :, ib, ik_g)) 2797 ! m_mn_tmp2 = MATMUL(m_mn_tmp1, umat_g(:, :, ikb)) 2798 ! 2799 m_mn_tmp3(:, :, ib, ik_g) = m_mn_tmp2(:, :) 2800 ENDDO 2801 ENDDO 2802 CALL mp_sum(a_mat_tmp, inter_pool_comm) 2803 CALL mp_sum(m_mn_tmp3, inter_pool_comm) 2804 ! 2805 ! a_mat(:, :, :) = a_mat_tmp(:, :, :) 2806 ! m_mat(:, :, :, :) = m_mn_tmp3(:, :, :, :) 2807 ! 2808 ! slim down a_mat and m_mat matrices to num_bands=nbnd-nexband bands 2809 ! 2810 DO ik = 1, nkstot 2811 DO iw = 1, n_proj 2812 ibnd_m = 0 2813 DO m = 1, nbnd 2814 IF (excluded_band(m)) CYCLE 2815 ibnd_m = ibnd_m + 1 2816 a_mat(ibnd_m, iw, ik) = a_mat_tmp(m, iw, ik) 2817 ENDDO 2818 ENDDO 2819 ! 2820 DO ib = 1, nnb 2821 ibnd_n = 0 2822 DO n = 1, nbnd 2823 IF (excluded_band(n)) CYCLE 2824 ibnd_n = ibnd_n + 1 2825 ibnd_m = 0 2826 DO m = 1, nbnd 2827 IF (excluded_band(m)) CYCLE 2828 ibnd_m = ibnd_m + 1 2829 m_mat(ibnd_m, ibnd_n, ib, ik) = m_mn_tmp3(m, n, ib, ik) 2830 ENDDO 2831 ENDDO 2832 ENDDO 2833 ENDDO 2834 ! 2835 DEALLOCATE(a_mat_tmp, STAT = ierr) 2836 IF (ierr /= 0) CALL errore('phases_a_m', 'Error deallocating a_mat_tmp', 1) 2837 DEALLOCATE(a_mat_tmp1, STAT = ierr) 2838 IF (ierr /= 0) CALL errore('phases_a_m', 'Error deallocating a_mat_tmp1', 1) 2839 DEALLOCATE(m_mat_tmp, STAT = ierr) 2840 IF (ierr /= 0) CALL errore('phases_a_m', 'Error deallocating m_mat_tmp', 1) 2841 DEALLOCATE(m_mn_tmp1, STAT = ierr) 2842 IF (ierr /= 0) CALL errore('phases_a_m', 'Error deallocating m_mn_tmp1', 1) 2843 DEALLOCATE(m_mn_tmp2, STAT = ierr) 2844 IF (ierr /= 0) CALL errore('phases_a_m', 'Error deallocating m_mn_tmp2', 1) 2845 DEALLOCATE(m_mn_tmp3, STAT = ierr) 2846 IF (ierr /= 0) CALL errore('phases_a_m', 'Error deallocating m_mn_tmp3', 1) 2847 ! 2848 RETURN 2849 ! 2850 !----------------------------------------------------------------------- 2851 END SUBROUTINE phases_a_m 2852 !----------------------------------------------------------------------- 2853 ! 2854 !----------------------------------------------------------------------- 2855 SUBROUTINE generate_guiding_functions(ik) 2856 !----------------------------------------------------------------------- 2857 !! 2858 !! Guiding functions 2859 !! 2860 ! 2861 USE kinds, ONLY : DP 2862 USE mp_global, ONLY : intra_pool_comm 2863 USE mp, ONLY : mp_sum 2864 USE wvfct, ONLY : npw 2865 USE gvect, ONLY : g 2866 USE cell_base, ONLY : tpiba 2867 USE wannierEPW, ONLY : n_proj, gf, center_w, csph, alpha_w, r_w 2868 USE klist, ONLY : igk_k 2869 USE klist_epw, ONLY : xk_loc 2870 USE constants_epw, ONLY : zero, czero, ci, twopi, eps8 2871 ! 2872 IMPLICIT NONE 2873 ! 2874 INTEGER, INTENT(in) :: ik 2875 !! Index of k-point 2876 ! 2877 ! Local variables 2878 INTEGER, PARAMETER :: lmax = 3 2879 !! Nr of spherical harmonics 2880 INTEGER, PARAMETER :: lmax2 = (lmax + 1)**2 2881 !! Total nr. of spherical harmonics 2882 INTEGER :: iw 2883 !! Counter on wannier projections 2884 INTEGER :: ig 2885 !! Counter on G-vectors 2886 INTEGER :: iig 2887 !! Index of G vector 2888 INTEGER :: lm 2889 !! Counter on spherical harmonics 2890 INTEGER :: l 2891 !! Momentum l 2892 INTEGER :: ierr 2893 !! Error status 2894 REAL(KIND = DP) :: arg 2895 !! 2*pi*(k+G)*r 2896 REAL(KIND = DP) :: anorm 2897 !! Anormal 2898 REAL(KIND = DP), EXTERNAL :: DDOT 2899 !! Scalar product of two vectors 2900 REAL(KIND = DP), ALLOCATABLE :: gk(:, :) 2901 !! k+G vectors 2902 REAL(KIND = DP), ALLOCATABLE :: qg(:) 2903 !! Norm of k+G 2904 REAL(KIND = DP), ALLOCATABLE :: ylm(:, :) 2905 !! Spherical harmonics 2906 REAL(KIND = DP), ALLOCATABLE :: radial(:, :) 2907 !! Radiam 2908 COMPLEX(KIND = DP) :: lphase 2909 !! (-i)^l 2910 COMPLEX(KIND = DP), EXTERNAL :: ZDOTC 2911 !! Scalar product of two complex vectors 2912 COMPLEX(KIND = DP), ALLOCATABLE :: sk(:) 2913 !! e^{-i*2*pi*(k+G)*r} 2914 ! 2915 ALLOCATE(gk(3, npw), STAT = ierr) 2916 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error allocating gk', 1) 2917 ALLOCATE(qg(npw), STAT = ierr) 2918 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error allocating qg', 1) 2919 ALLOCATE(ylm(npw, lmax2), STAT = ierr) 2920 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error allocating ylm', 1) 2921 ALLOCATE(sk(npw), STAT = ierr) 2922 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error allocating sk', 1) 2923 ALLOCATE(radial(npw, 0:lmax), STAT = ierr) 2924 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error allocating radial', 1) 2925 gk(:, :) = zero 2926 qg(:) = zero 2927 ylm(:, :) = zero 2928 sk(:) = czero 2929 radial(:, 0:lmax) = zero 2930 ! 2931 DO ig = 1, npw 2932 gk(:, ig) = xk_loc(:, ik) + g(:, igk_k(ig, ik)) 2933 qg(ig) = SUM(gk(:, ig) * gk(:, ig)) 2934 ENDDO 2935 ! 2936 ! get spherical harmonics ylm up to lmax 2937 CALL ylmr2(lmax2, npw, gk, qg, ylm) 2938 ! define qg as the norm of (k+g) in a.u. 2939 qg(:) = DSQRT(qg(:)) * tpiba 2940 ! 2941 ! RM changed according to QE4.0.3/PP/pw2wannier90 2942 DO iw = 1, n_proj 2943 ! 2944 gf(:, iw) = czero 2945 ! 2946 CALL radialpart(npw, qg, alpha_w(iw), r_w(iw), lmax, radial) 2947 ! 2948 DO lm = 1, lmax2 2949 IF (ABS(csph(lm, iw)) < eps8) CYCLE 2950 l = INT(DSQRT(lm - 1.d0)) 2951 lphase = (- ci)**l 2952 ! 2953 DO ig = 1, npw 2954 gf(ig, iw) = gf(ig, iw) + csph(lm, iw) * ylm(ig, lm) * radial(ig, l) * lphase 2955 ENDDO !ig 2956 ENDDO ! lm 2957 ! 2958 DO ig = 1, npw 2959 iig = igk_k(ig, ik) 2960 arg = twopi * DDOT(3, gk(:, ig), 1, center_w(:, iw), 1) 2961 ! center_w are cartesian coordinates in units of alat 2962 sk(ig) = CMPLX(COS(arg), - SIN(arg), KIND = DP) 2963 gf(ig, iw) = gf(ig, iw) * sk(ig) 2964 ENDDO 2965 anorm = REAL(ZDOTC(npw, gf(1, iw), 1, gf(1, iw), 1)) 2966 CALL mp_sum(anorm, intra_pool_comm) 2967 gf(:, iw) = gf(:, iw) / DSQRT(anorm) 2968 ENDDO 2969 ! 2970 DEALLOCATE(gk, STAT = ierr) 2971 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error deallocating gk', 1) 2972 DEALLOCATE(qg, STAT = ierr) 2973 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error deallocating qg', 1) 2974 DEALLOCATE(ylm, STAT = ierr) 2975 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error deallocating ylm', 1) 2976 DEALLOCATE(sk, STAT = ierr) 2977 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error deallocating sk', 1) 2978 DEALLOCATE(radial, STAT = ierr) 2979 IF (ierr /= 0) CALL errore('generate_guiding_functions', 'Error deallocating radial', 1) 2980 ! 2981 RETURN 2982 ! 2983 !----------------------------------------------------------------------- 2984 END SUBROUTINE generate_guiding_functions 2985 !----------------------------------------------------------------------- 2986 ! 2987 !----------------------------------------------------------------------- 2988 SUBROUTINE write_band() 2989 !----------------------------------------------------------------------- 2990 !! 2991 !! Write bands 2992 !! 2993 USE wvfct, ONLY : nbnd, et 2994 USE constants, ONLY : rytoev 2995 USE wannierEPW, ONLY : ikstart, ikstop, iknum, num_bands, eigval, & 2996 excluded_band 2997 USE constants_epw, ONLY : zero 2998 ! 2999 IMPLICIT NONE 3000 ! 3001 ! Local variables 3002 INTEGER :: ik 3003 !! Counter on k-point 3004 INTEGER ikevc 3005 !! k-point index 3006 INTEGER :: ibnd 3007 !! Counter on bands 3008 INTEGER :: ibnd1 3009 !! Band index 3010 INTEGER :: ierr 3011 !! Error status 3012 ! 3013 ALLOCATE(eigval(num_bands, iknum), STAT = ierr) 3014 IF (ierr /= 0) CALL errore('write_band', 'Error allocating eigval', 1) 3015 eigval(:, :) = zero 3016 ! 3017 DO ik = ikstart, ikstop 3018 ikevc = ik - ikstart + 1 3019 ibnd1 = 0 3020 DO ibnd = 1, nbnd 3021 IF (excluded_band(ibnd)) CYCLE 3022 ibnd1 = ibnd1 + 1 3023 ! 3024 ! RM - same value for rytoev as in wannier90 3025 ! eigval(ibnd1, ikevc) = et(ibnd, ik) * ryd2ev 3026 eigval(ibnd1, ikevc) = et(ibnd, ik) * rytoev 3027 ENDDO 3028 ENDDO 3029 ! 3030 RETURN 3031 ! 3032 !------------------------------------------------------------------------ 3033 END SUBROUTINE write_band 3034 !----------------------------------------------------------------------- 3035 ! 3036 !----------------------------------------------------------------------- 3037 SUBROUTINE write_plot() 3038 !----------------------------------------------------------------------- 3039 !! 3040 !! JN 06/2009: 3041 !! added a couple of calls -- now works with multiple 3042 !! pools/procs (but one proc per pool) 3043 !! 3044 !! HL 03/2020: 3045 !! New implementation based on some parts of 3046 !! QE (subroutine of write_plot in /PP/src/pw2wannier.f90) 3047 !! and Wannier90 (subroutine of plot_wannier in /src/plot.F90) 3048 !! 3049 USE kinds, ONLY : DP 3050 USE io_global, ONLY : stdout, meta_ionode, meta_ionode_id 3051 USE wvfct, ONLY : nbnd, npw, npwx 3052 USE wavefunctions, ONLY : evc, psic, psic_nc 3053 USE wannierEPW, ONLY : excluded_band, u_mat, u_mat_opt, iknum, & 3054 n_wannier, num_bands, lwindow 3055 USE epwcom, ONLY : wannier_plot_supercell, reduce_unk 3056 USE klist, ONLY : nks, igk_k, ngk 3057 USE klist_epw, ONLY : xk_loc 3058 USE fft_base, ONLY : dffts 3059 USE fft_interfaces, ONLY : invfft 3060 USE noncollin_module,ONLY : noncolin, npol 3061 USE constants_epw, ONLY : czero, zero, twopi, ci 3062 USE mp_pools, ONLY : my_pool_id 3063 USE kfold, ONLY : ktokpmq 3064 USE io_epw, ONLY : readwfc 3065 USE mp_world, ONLY : world_comm 3066 USE elph2, ONLY : nbndep, wanplotlist, num_wannier_plot 3067 USE cell_base, ONLY : at, bg, alat, tpiba 3068 USE constants_epw, ONLY : bohr 3069 USE wannierEPW, ONLY : wann_centers 3070 USE epwcom, ONLY : wannier_plot_radius, wannier_plot_scale 3071 USE control_flags, ONLY : iverbosity 3072 USE mp, ONLY : mp_barrier 3073 USE parallel_include 3074 ! 3075 IMPLICIT NONE 3076 ! 3077 ! Local variables 3078 INTEGER :: ik 3079 !! Counter on k-point 3080 INTEGER :: ibnd 3081 !! Counter on bands 3082 INTEGER :: ibnd0 3083 !! Band index 3084 INTEGER :: ibnd1 3085 !! Band index 3086 INTEGER :: ig 3087 !! Counter on G vectors 3088 INTEGER :: ipool 3089 !! Index of current pool 3090 INTEGER :: nkq 3091 !! Index of k-point in the pool 3092 INTEGER :: nkq_abs 3093 !! Absolute index of k-point 3094 INTEGER :: ik_g 3095 !! Temporary index of k-point, ik_g = nkq_abs 3096 INTEGER :: ipol 3097 !! Counter on polarizations 3098 INTEGER :: ispw 3099 !! Starting index of plane waves 3100 INTEGER :: iepw 3101 !! Ending index of plane waves 3102 INTEGER :: ngx 3103 !! Number of soft grids 3104 INTEGER :: ngy 3105 !! Number of soft grids 3106 INTEGER :: ngz 3107 !! Number of soft grids 3108 INTEGER :: n1by2 3109 !! Number of reduced grids 3110 INTEGER :: n2by2 3111 !! Number of reduced grids 3112 INTEGER :: n3by2 3113 !! Number of reduced grids 3114 INTEGER :: nx 3115 !! Counter on primitive-cell grids 3116 INTEGER :: ny 3117 !! Counter on primitive-cell grids 3118 INTEGER :: nz 3119 !! Counter on primitive-cell grids 3120 INTEGER :: nxx 3121 !! Counter on super-cell grids 3122 INTEGER :: nyy 3123 !! Counter on super-cell grids 3124 INTEGER :: nzz 3125 !! Counter on super-cell grids 3126 INTEGER :: loop_w 3127 !! Counter on Wannier functions 3128 INTEGER :: wann_index 3129 !! Index of Wannier functions to plot 3130 INTEGER :: ngridwf 3131 !! Number of grids to describe real-space Wannier function 3132 INTEGER :: ngridwf_max 3133 !! Max. of number of grids to describe real-space Wannier function 3134 INTEGER :: ipoint 3135 !! Real-space grid index 3136 INTEGER :: idx 3137 !! Real-space grid index 3138 INTEGER :: qxx 3139 !! Temporary grid index 3140 INTEGER :: qyy 3141 !! Temporary grid index 3142 INTEGER :: qzz 3143 !! Temporary grid index 3144 INTEGER :: i 3145 !! Counter index 3146 INTEGER :: ierr 3147 !! Error status 3148 INTEGER :: ndimwin(nks) 3149 !! Number of bands within outer window at each k-point 3150 INTEGER :: ngs(3) 3151 !! Size of the supercell for Wannier function plot 3152 INTEGER :: istart(3) 3153 !! First grid to consider in cube files 3154 INTEGER :: iend(3) 3155 !! Last grid to consider in cube files 3156 INTEGER :: ilength(3) 3157 !! Length of grids in cube files 3158 REAL(KIND = DP) :: rstart(3) 3159 !! Start of cube wrt simulation (home) cell origin 3160 REAL(KIND = DP) :: rend(3) 3161 !! End of cube wrt simulation (home) cell origin 3162 REAL(KIND = DP) :: rlength(3) 3163 !! Length of region to consider in cube files 3164 REAL(KIND = DP) :: orig(3) 3165 !! Origin of cube wrt simulation (home) cell in Cart. coords. 3166 REAL(KIND = DP) :: dgrid(3) 3167 !! Grid spacing in each lattice direction 3168 REAL(KIND = DP) :: moda(3) 3169 !! Length of real lattice vectors 3170 REAL(KIND = DP) :: modb(3) 3171 !! Length of reciprocal lattice vectors 3172 REAL(KIND = DP) :: scalfac 3173 !! Scaling factor 3174 REAL(KIND = DP) :: tmax 3175 !! Temporary max. of wavefunction norm 3176 REAL(KIND = DP) :: tmaxx 3177 !! Temporary max. of wavefunction norm 3178 REAL(KIND = DP) :: ratio 3179 !! Im/Re Ratio 3180 REAL(KIND = DP) :: ratmax 3181 !! Max Im/Re Ratio 3182 REAL(KIND = DP) :: zero_vect(3) 3183 !! Temporary zero vector 3184 REAL(KIND = DP) :: xcrys(3) 3185 !! x in crystal coords. 3186 COMPLEX(KIND = DP) :: catmp 3187 !! Temporary complex number 3188 COMPLEX(KIND = DP) :: uptmp 3189 !! Temporary complex number 3190 COMPLEX(KIND = DP) :: dntmp 3191 !! Temporary complex number 3192 COMPLEX(KIND = DP) :: wmod 3193 !! Conversion factor for wavefunction 3194 COMPLEX(KIND = DP), ALLOCATABLE :: u_kc(:, :) 3195 !! Rotation matrix, U^dis*U 3196 COMPLEX(KIND = DP), ALLOCATABLE :: rotwf(:, :, :) 3197 !! Rotated wave function in reciprocal space 3198 COMPLEX(KIND = DP), ALLOCATABLE :: psic_small(:, :) 3199 !! Rotated wave function in reduced real space 3200 COMPLEX(KIND = DP), ALLOCATABLE :: wann_func(:, :, :, :) 3201 !! Real-space Wannier function 3202 ! 3203 CALL start_clock('write_plot') 3204 ! 3205 ngs(:) = wannier_plot_supercell(:) 3206 ! 3207 zero_vect(:) = zero 3208 WRITE(stdout,'(a,/)') ' Writing out Wannier function cube files' 3209 ! 3210 IF (iverbosity == 1) THEN 3211 WRITE(stdout,'(a,f6.3)') 'write_plot: wannier_plot_radius =', & 3212 wannier_plot_radius 3213 WRITE(stdout,'(a,f6.3)') 'write_plot: wannier_plot_scale =', & 3214 wannier_plot_scale 3215 ENDIF 3216 ! 3217 ngx = dffts%nr1 3218 ngy = dffts%nr2 3219 ngz = dffts%nr3 3220 IF (reduce_unk) THEN 3221 ngx = (dffts%nr1 + 1) / 2 3222 ngy = (dffts%nr2 + 1) / 2 3223 ngz = (dffts%nr3 + 1) / 2 3224 ALLOCATE(psic_small(ngx*ngy*ngz, npol), STAT = ierr) 3225 IF (ierr /= 0) CALL errore('write_plot', 'Error allocating psic_small', 1) 3226 psic_small(:, :) = czero 3227 WRITE(stdout, '(a)') 'write_plot: Real-space grids for plotting Wannier functions are reduced' 3228 ENDIF 3229 WRITE(stdout, '(3(a, i5))') 'nr1s =', ngx, ', nr2s =', ngy, ', nr3s =', ngz 3230 WRITE(stdout, '(a, 3i5)') 'write_plot: wannier_plot_supercell =', & 3231 (wannier_plot_supercell(i), i = 1, 3) 3232 ! 3233 ndimwin(:) = 0 3234 ALLOCATE(u_kc(nbndep, n_wannier), STAT = ierr) 3235 IF (ierr /= 0) CALL errore('write_plot', 'Error allocating u_kc', 1) 3236 u_kc(:, :) = czero 3237 ALLOCATE(rotwf(npwx*npol, num_wannier_plot, nks), STAT = ierr) 3238 IF (ierr /= 0) CALL errore('write_plot', 'Error allocating rotwf', 1) 3239 rotwf(:, :, :) = czero 3240 ! 3241 DO ik = 1, nks 3242 ! 3243 npw = ngk(ik) 3244 ! 3245 ! returns in-pool index nkq and absolute index nkq_abs of xk 3246 ! 3247 CALL ktokpmq(xk_loc(:, ik), zero_vect, +1, ipool, nkq, nkq_abs) 3248 ik_g = nkq_abs 3249 ! 3250 DO ibnd = 1, num_bands 3251 IF (lwindow(ibnd, ik_g)) ndimwin(ik) = ndimwin(ik) + 1 3252 ENDDO 3253 ! 3254 ! get the final rotation matrix, which is the product of the optimal 3255 ! subspace and the rotation among the n_wannier wavefunctions 3256 ! 3257 u_kc(1:ndimwin(ik), 1:n_wannier) = & 3258 MATMUL(u_mat_opt(1:ndimwin(ik), :, ik_g), u_mat(:, 1:n_wannier, ik_g)) 3259 ! 3260 ! read wfc at ik 3261 ! 3262 CALL readwfc(my_pool_id + 1, ik, evc) 3263 ! 3264 ibnd0 = 0 3265 ibnd1 = 0 3266 DO ibnd = 1, nbnd 3267 ! 3268 IF (excluded_band(ibnd)) CYCLE 3269 ibnd0 = ibnd0 + 1 3270 IF (lwindow(ibnd0, ik_g)) THEN 3271 ibnd1 = ibnd1 + 1 3272 DO loop_w = 1, num_wannier_plot 3273 DO ig = 1, npw 3274 rotwf(ig, loop_w, ik) = rotwf(ig, loop_w, ik) + & 3275 u_kc(ibnd1, wanplotlist(loop_w)) * evc(ig, ibnd) 3276 IF (noncolin) THEN 3277 rotwf(ig + npwx, loop_w, ik) = rotwf(ig + npwx, loop_w, ik) + & 3278 u_kc(ibnd1, wanplotlist(loop_w)) * evc(ig + npwx, ibnd) 3279 ENDIF 3280 ENDDO 3281 ENDDO 3282 ENDIF 3283 ! 3284 ENDDO 3285 ! 3286 ENDDO 3287 ! 3288 CALL mp_barrier(world_comm) 3289 ! 3290 ! Lengths of real and reciprocal lattice vectors 3291 ! 3292 DO i = 1, 3 3293 moda(i) = DSQRT( at(1, i) * at(1, i) & 3294 + at(2, i) * at(2, i) & 3295 + at(3, i) * at(3, i) ) * alat 3296 modb(i) = DSQRT( bg(1, i) * bg(1, i) & 3297 + bg(2, i) * bg(2, i) & 3298 + bg(3, i) * bg(3, i) ) * tpiba 3299 ENDDO 3300 ! 3301 ! Grid spacing in each lattice direction 3302 ! 3303 dgrid(1) = moda(1) / ngx 3304 dgrid(2) = moda(2) / ngy 3305 dgrid(3) = moda(3) / ngz 3306 ! 3307 DO loop_w = 1, num_wannier_plot 3308 wann_index = wanplotlist(loop_w) 3309 ! 3310 ! Find start and end of cube wrt simulation (home) cell origin 3311 ! 3312 DO i = 1, 3 3313 ! ... in terms of distance along each lattice vector direction i 3314 rstart(i) = (wann_centers(1, wann_index) / bohr / alat * bg(1, i) & 3315 + wann_centers(2, wann_index) / bohr / alat * bg(2, i) & 3316 + wann_centers(3, wann_index) / bohr / alat * bg(3, i)) * moda(i) & 3317 - twopi * wannier_plot_radius / bohr / (moda(i) * modb(i)) 3318 rend(i) = (wann_centers(1, wann_index) / bohr / alat * bg(1, i) & 3319 + wann_centers(2, wann_index) / bohr / alat * bg(2, i) & 3320 + wann_centers(3, wann_index) / bohr / alat * bg(3, i)) * moda(i) & 3321 + twopi * wannier_plot_radius / bohr / (moda(i) * modb(i)) 3322 ENDDO 3323 ! 3324 rlength(:) = rend(:) - rstart(:) 3325 ilength(:) = CEILING(rlength(:) / dgrid(:)) 3326 ! 3327 ! ... in terms of integer gridpoints along each lattice vector direction i 3328 ! 3329 istart(:) = FLOOR(rstart(:) / dgrid(:)) + 1 3330 iend(:) = istart(:) + ilength(:) - 1 3331 ! 3332 ! Origin of cube wrt simulation (home) cell in Cartesian co-ordinates 3333 ! 3334 DO i = 1, 3 3335 orig(i) = & 3336 REAL(istart(1) - 1, KIND = DP) * dgrid(1) * at(i, 1) * alat / moda(1) & 3337 + REAL(istart(2) - 1, KIND = DP) * dgrid(2) * at(i, 2) * alat / moda(2) & 3338 + REAL(istart(3) - 1, KIND = DP) * dgrid(3) * at(i, 3) * alat / moda(3) 3339 ENDDO 3340 ! 3341 DO nzz = 1, ilength(3) 3342 qzz = nzz + istart(3) - 1 3343 IF ((qzz < (-((ngs(3))/2)*ngz)) .OR. & 3344 (qzz > ((ngs(3) + 1)/2)*ngz - 1)) THEN 3345 WRITE(stdout, *) 'Error plotting WF cube. Try one of the following:' 3346 WRITE(stdout, *) ' (1) increase wannier_plot_supercell;' 3347 WRITE(stdout, *) ' (2) decrease wannier_plot_radius;' 3348 CALL errore('write_plot', 'Wrong range in plotting WF cube.', 1) 3349 ENDIF 3350 DO nyy = 1, ilength(2) 3351 qyy = nyy + istart(2) - 1 3352 IF ((qyy < (-((ngs(2))/2)*ngy)) .OR. & 3353 (qyy > ((ngs(2) + 1)/2)*ngy - 1)) THEN 3354 WRITE(stdout, *) 'Error plotting WF cube. Try one of the following:' 3355 WRITE(stdout, *) ' (1) increase wannier_plot_supercell;' 3356 WRITE(stdout, *) ' (2) decrease wannier_plot_radius;' 3357 CALL errore('write_plot', 'Wrong range in plotting WF cube.', 1) 3358 ENDIF 3359 DO nxx = 1, ilength(1) 3360 qxx = nxx + istart(1) - 1 3361 IF ((qxx < (-((ngs(1))/2)*ngx)) .OR. & 3362 (qxx > ((ngs(1) + 1)/2)*ngx - 1)) THEN 3363 WRITE(stdout, *) 'Error plotting WF cube. Try one of the following:' 3364 WRITE(stdout, *) ' (1) increase wannier_plot_supercell;' 3365 WRITE(stdout, *) ' (2) decrease wannier_plot_radius;' 3366 CALL errore('write_plot', 'Wrong range in plotting WF cube.', 1) 3367 ENDIF 3368 ENDDO 3369 ENDDO 3370 ENDDO 3371 ! 3372 ngridwf = ilength(1) * ilength(2) * ilength(3) 3373 IF (loop_w == 1) THEN 3374 ngridwf_max = ngridwf 3375 ALLOCATE(wann_func(ilength(1), ilength(2), ilength(3), npol), STAT = ierr) 3376 IF (ierr /= 0) CALL errore('write_plot', 'Error allocating wann_func', 1) 3377 ELSE 3378 IF (ngridwf > ngridwf_max) THEN 3379 DEALLOCATE(wann_func, STAT = ierr) 3380 IF (ierr /= 0) CALL errore('write_plot', 'Error deallocating wann_func', 1) 3381 ALLOCATE(wann_func(ilength(1), ilength(2), ilength(3), npol), STAT = ierr) 3382 IF (ierr /= 0) CALL errore('write_plot', 'Error allocating wann_func', 1) 3383 ngridwf_max = ngridwf 3384 ENDIF 3385 ENDIF 3386 wann_func(:, :, :, :) = czero 3387 ! 3388 DO ik = 1, nks 3389 ! 3390 npw = ngk(ik) 3391 xcrys = xk_loc(:, ik) 3392 CALL cryst_to_cart(1, xcrys, at, -1) 3393 IF (noncolin) THEN 3394 psic_nc(:, :) = czero 3395 DO ipol = 1, 2 3396 ispw = 1 + npwx * (ipol - 1) 3397 iepw = npw + npwx * (ipol - 1) 3398 psic_nc(dffts%nl(igk_k(1:npw, ik)), ipol) = rotwf(ispw:iepw, loop_w, ik) 3399 CALL invfft('Wave', psic_nc(:,ipol), dffts) 3400 ENDDO 3401 ELSE 3402 psic(:) = czero 3403 psic(dffts%nl(igk_k(1:npw, ik))) = rotwf(1:npw, loop_w, ik) 3404 CALL invfft('Wave', psic, dffts) 3405 ENDIF 3406 ! 3407 IF (reduce_unk) THEN 3408 idx = 0 3409 DO nz = 1, dffts%nr3, 2 3410 DO ny = 1, dffts%nr2, 2 3411 DO nx = 1, dffts%nr1, 2 3412 ipoint = nx + (ny - 1)*dffts%nr1 + (nz - 1)*dffts%nr2*dffts%nr1 3413 idx = idx + 1 3414 IF (noncolin) THEN 3415 psic_small(idx, 1) = psic_nc(ipoint, 1) 3416 psic_small(idx, 2) = psic_nc(ipoint, 2) 3417 ELSE 3418 psic_small(idx, 1) = psic(ipoint) 3419 ENDIF 3420 ENDDO 3421 ENDDO 3422 ENDDO 3423 ENDIF 3424 ! 3425 DO nzz = istart(3), iend(3) 3426 nz = MOD(nzz, ngz) 3427 IF (nz < 1) nz = nz + ngz 3428 qzz = nzz - istart(3) + 1 3429 DO nyy = istart(2), iend(2) 3430 ny = MOD(nyy, ngy) 3431 IF (ny < 1) ny = ny + ngy 3432 qyy = nyy - istart(2) + 1 3433 DO nxx = istart(1), iend(1) 3434 nx = MOD(nxx, ngx) 3435 IF (nx < 1) nx = nx + ngx 3436 qxx = nxx - istart(1) + 1 3437 scalfac = xcrys(1) * DBLE(nxx - 1) / DBLE(ngx) + & 3438 xcrys(2) * DBLE(nyy - 1) / DBLE(ngy) + & 3439 xcrys(3) * DBLE(nzz - 1) / DBLE(ngz) 3440 ipoint = nx + (ny - 1) * ngx + (nz - 1) * ngy * ngx 3441 catmp = EXP(twopi * ci * scalfac) 3442 IF (.NOT. noncolin) THEN 3443 IF (reduce_unk) THEN 3444 wann_func(qxx, qyy, qzz, 1) = wann_func(qxx, qyy, qzz, 1) + & 3445 psic_small(ipoint, 1) * catmp 3446 ELSE 3447 wann_func(qxx, qyy, qzz, 1) = wann_func(qxx, qyy, qzz, 1) + & 3448 psic(ipoint) * catmp 3449 ENDIF 3450 ELSE 3451 IF (reduce_unk) THEN 3452 wann_func(qxx, qyy, qzz, 1) = wann_func(qxx, qyy, qzz, 1) + & 3453 psic_small(ipoint, 1) * catmp 3454 wann_func(qxx, qyy, qzz, 2) = wann_func(qxx, qyy, qzz, 2) + & 3455 psic_small(ipoint, 2) * catmp 3456 ELSE 3457 wann_func(qxx, qyy, qzz, 1) = wann_func(qxx, qyy, qzz, 1) + & 3458 psic_nc(ipoint, 1) * catmp 3459 wann_func(qxx, qyy, qzz, 2) = wann_func(qxx, qyy, qzz, 2) + & 3460 psic_nc(ipoint, 2) * catmp 3461 ENDIF 3462 ENDIF 3463 ENDDO ! nxx 3464 ENDDO ! nyy 3465 ENDDO ! nzz 3466 ENDDO ! ik 3467 ! 3468#if defined(__MPI) 3469 IF (meta_ionode) THEN 3470 CALL MPI_REDUCE( MPI_IN_PLACE, wann_func, 2 * npol * ngridwf_max, MPI_DOUBLE_PRECISION, & 3471 MPI_SUM, meta_ionode_id, world_comm, ierr ) 3472 ELSE 3473 CALL MPI_REDUCE( wann_func, wann_func, 2 * npol * ngridwf_max, MPI_DOUBLE_PRECISION, & 3474 MPI_SUM, meta_ionode_id, world_comm, ierr ) 3475 ENDIF 3476 IF (ierr /= 0) CALL errore('write_plot', 'mpi_reduce', ierr) 3477#endif 3478 ! 3479 IF (meta_ionode) THEN 3480 ! 3481 !! [HL] For spinor Wannier functions, the step below is not necessary. 3482 ! 3483 IF (.NOT. noncolin) THEN 3484 ! fix the global phase by setting the wannier to 3485 ! be real at the point where it has max. modulus 3486 tmaxx = 0.0d0 3487 wmod = CMPLX(1.0d0, 0.0d0, KIND = DP) 3488 DO nzz = 1, ilength(3) 3489 DO nyy = 1, ilength(2) 3490 DO nxx = 1, ilength(1) 3491 wann_func(nxx, nyy, nzz, 1) = wann_func(nxx, nyy, nzz, 1)/REAL(iknum, KIND = DP) 3492 tmax = REAL(wann_func(nxx, nyy, nzz, 1) * & 3493 CONJG(wann_func(nxx, nyy, nzz, 1)), KIND = DP) 3494 IF (tmax > tmaxx) THEN 3495 tmaxx = tmax 3496 wmod = wann_func(nxx, nyy, nzz, 1) 3497 ENDIF 3498 ENDDO 3499 ENDDO 3500 ENDDO 3501 wmod = wmod / DSQRT(DBLE(wmod)**2 + AIMAG(wmod)**2) 3502 wann_func(:, :, :, :) = wann_func(:, :, :, :) / wmod 3503 ! 3504 ! Check the 'reality' of the WF 3505 ! 3506 ratmax = 0.0d0 3507 DO nzz = 1, ilength(3) 3508 DO nyy = 1, ilength(2) 3509 DO nxx = 1, ilength(1) 3510 IF (ABS(REAL(wann_func(nxx, nyy, nzz, 1), KIND = DP)) >= 0.01d0) THEN 3511 ratio = ABS(AIMAG(wann_func(nxx, nyy, nzz, 1))) / & 3512 ABS(REAL(wann_func(nxx, nyy, nzz, 1), KIND = DP)) 3513 ratmax = MAX(ratmax, ratio) 3514 ENDIF 3515 ENDDO 3516 ENDDO 3517 ENDDO 3518 WRITE (stdout, '(6x,a,i4,7x,a,f11.6)') 'Wannier Function Num: ', & 3519 wanplotlist(loop_w), 'Maximum Im/Re Ratio = ', ratmax 3520 ELSE 3521 DO nzz = 1, ilength(3) 3522 DO nyy = 1, ilength(2) 3523 DO nxx = 1, ilength(1) 3524 wann_func(nxx, nyy, nzz, 1) = CMPLX(DSQRT( & 3525 REAL(wann_func(nxx, nyy, nzz, 1) * CONJG(wann_func(nxx, nyy, nzz, 1)), KIND = DP) & 3526 + REAL(wann_func(nxx, nyy, nzz, 2) * CONJG(wann_func(nxx, nyy, nzz, 2)), KIND = DP)) & 3527 / DBLE(iknum), 0.0d0, KIND = DP) 3528 ENDDO 3529 ENDDO 3530 ENDDO 3531 ENDIF 3532 CALL internal_cube_format(loop_w) 3533 ENDIF 3534 CALL mp_barrier(world_comm) 3535 ENDDO ! loop_w 3536 ! 3537 DEALLOCATE(wann_func, STAT = ierr) 3538 IF (ierr /= 0) CALL errore('write_plot', 'Error deallocating wann_func', 1) 3539 ! 3540 WRITE(stdout, '(/)') 3541 WRITE(stdout, *) ' cube files written' 3542 ! 3543 IF (reduce_unk) THEN 3544 DEALLOCATE(psic_small, STAT = ierr) 3545 IF (ierr /= 0) CALL errore('write_plot', 'Error deallocating psic_small', 1) 3546 ENDIF 3547 DEALLOCATE(u_kc, STAT = ierr) 3548 IF (ierr /= 0) CALL errore('write_plot', 'Error deallocating u_kc', 1) 3549 DEALLOCATE(rotwf, STAT = ierr) 3550 IF (ierr /= 0) CALL errore('write_plot', 'Error deallocating rotwf', 1) 3551 DEALLOCATE(wanplotlist, STAT = ierr) 3552 IF (ierr /= 0) CALL errore('write_plot', 'Error deallocating wanplotlist', 1) 3553 ! 3554 CALL stop_clock('write_plot') 3555 ! 3556 RETURN 3557 ! 3558 CONTAINS 3559 ! 3560 !----------------------------------------------------------------------- 3561 SUBROUTINE internal_cube_format(loop_w) 3562 !----------------------------------------------------------------------- 3563 !! 3564 !! Write WFs in Gaussian cube format. 3565 !! 3566 !! HL 03/2020: 3567 !! Imported and adapted from the same name of subroutine in plot.F90 3568 !! in the directory of src in Wannier90 3569 !! 3570 USE ions_base, ONLY : nat, tau, ityp, atm, nsp 3571 USE io_var, ONLY : iun_plot 3572 USE io_files, ONLY : prefix 3573 ! 3574 IMPLICIT NONE 3575 ! 3576 INTEGER, INTENT(in) :: loop_w 3577 !! Index of Wannier functions to plot 3578 CHARACTER(LEN = 60) :: wancube 3579 !! Name of Wannier function cube files 3580 CHARACTER(LEN = 2), DIMENSION(109) :: periodic_table = (/ & 3581 & 'H ', 'He', & 3582 & 'Li', 'Be', 'B ', 'C ', 'N ', 'O ', 'F ', 'Ne', & 3583 & 'Na', 'Mg', 'Al', 'Si', 'P ', 'S ', 'Cl', 'Ar', & 3584 & 'K ', 'Ca', 'Sc', 'Ti', 'V ', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', & 3585 & 'Rb', 'Sr', 'Y ', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I ', 'Xe', & 3586 & 'Cs', 'Ba', & 3587 & 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', & 3588 & 'Hf', 'Ta', 'W ', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', & 3589 & 'Fr', 'Ra', & 3590 & 'Ac', 'Th', 'Pa', 'U ', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', & 3591 & 'Rf', 'Db', 'Sg', 'Bh', 'Hs', 'Mt'/) 3592 CHARACTER(LEN = 9) :: cdate 3593 !! Current date 3594 CHARACTER(LEN = 9) :: ctime 3595 !! Current time 3596 INTEGER :: iname 3597 !! Counter on elements 3598 INTEGER :: max_elements 3599 !! Max. of atomic number 3600 INTEGER :: isp 3601 !! Counter on species 3602 INTEGER :: iat 3603 !! Counter on atoms 3604 INTEGER :: icount 3605 !! Counter on atoms within a given radius of Wannier centre 3606 INTEGER, ALLOCATABLE :: atomic_Z(:) 3607 !! Atomic number 3608 REAL(KIND = DP) :: val_Q 3609 !! Dummy value for cube file 3610 REAL(KIND = DP) :: dist 3611 !! Distance from Wannier centre to atoms 3612 REAL(KIND = DP) :: wcf(3) 3613 !! Wannier centre in fractional coords. 3614 REAL(KIND = DP) :: diff(3) 3615 !! Vector (in fractional coords.) from Wannier centre to "centre of mass" 3616 REAL(KIND = DP) :: difc(3) 3617 !! Temporary difference vector 3618 REAL(KIND = DP) :: xau(3, nat) 3619 !! Atomic positions (in fractional coords.) 3620 ! 3621 ALLOCATE(atomic_Z(nsp), STAT = ierr) 3622 IF (ierr /= 0) CALL errore('internal_cube_format', 'Error allocating atomic_Z', 1) 3623 ! 3624 val_Q = 1.0d0 ! dummy value for cube file 3625 ! 3626 ! Assign atomic numbers to species 3627 ! 3628 max_elements = SIZE(periodic_table) 3629 DO isp = 1, nsp 3630 DO iname = 1, max_elements 3631 IF (atm(isp) == periodic_table(iname)) THEN 3632 atomic_Z(isp) = iname 3633 EXIT 3634 ENDIF 3635 ENDDO 3636 ENDDO 3637 ! 3638202 FORMAT(a, '_', i5.5, '.cube') 3639 ! 3640 ! Compute the coordinates of each atom in the basis of the 3641 ! direct lattice vectors 3642 ! 3643 DO iat = 1, nat 3644 DO ipol = 1, 3 3645 xau(ipol,iat) = bg(1,ipol)*tau(1,iat) + bg(2,ipol)*tau(2,iat) + & 3646 bg(3,ipol)*tau(3,iat) 3647 ENDDO 3648 ENDDO 3649 ! 3650 wann_index = wanplotlist(loop_w) 3651 WRITE(wancube, 202) TRIM(prefix), wann_index 3652 IF (iverbosity == 1) THEN 3653 WRITE(stdout, '(a,i12)') 'loop_w =', loop_w 3654 WRITE(stdout, '(a,3i12)') 'ngi =', ngx, ngy, ngz 3655 WRITE(stdout, '(a,3f12.6)') 'dgrid =', (dgrid(i), i=1, 3) 3656 WRITE(stdout, '(a,3f12.6)') 'rstart =', (rstart(i), i=1, 3) 3657 WRITE(stdout, '(a,3f12.6)') 'rend =', (rend(i), i=1, 3) 3658 WRITE(stdout, '(a,3f12.6)') 'rlength =', (rlength(i), i=1, 3) 3659 WRITE(stdout, '(a,3i12)') 'istart =', (istart(i), i=1, 3) 3660 WRITE(stdout, '(a,3i12)') 'iend =', (iend(i), i=1, 3) 3661 WRITE(stdout, '(a,3i12)') 'ilength =', (ilength(i), i=1, 3) 3662 WRITE(stdout, '(a,3f12.6)') 'orig =', (orig(i), i=1, 3) 3663 WRITE(stdout, '(a,3f12.6)') 'wann_cen=', (wann_centers(i, wann_index), i=1, 3) 3664 ENDIF 3665 ! 3666 ! WF centre in fractional coordinates 3667 ! 3668 wcf(:) = wann_centers(:, wann_index) / bohr / alat 3669 CALL cryst_to_cart(1, wcf(:), bg, -1) 3670 ! 3671 IF (iverbosity == 1) THEN 3672 WRITE(stdout, '(a,3f12.6)') 'wcf =', (wcf(i), i=1, 3) 3673 ENDIF 3674 ! 3675 icount = 0 3676 DO iat = 1, nat 3677 DO nzz = -ngs(3)/2, (ngs(3) + 1)/2 3678 DO nyy = -ngs(2)/2, (ngs(2) + 1)/2 3679 DO nxx = -ngs(1)/2, (ngs(1) + 1)/2 3680 diff(:) = xau(:, iat) - wcf(:) & 3681 + (/REAL(nxx, KIND = DP), REAL(nyy, KIND = DP), REAL(nzz, KIND = DP)/) 3682 difc(:) = diff(:) 3683 CALL cryst_to_cart(1, difc(:), at, 1) 3684 dist = DSQRT(difc(1)*difc(1) + difc(2)*difc(2) + difc(3)*difc(3)) * alat 3685 IF (dist < (wannier_plot_scale * wannier_plot_radius / bohr)) THEN 3686 icount = icount + 1 3687 ENDIF 3688 ENDDO 3689 ENDDO 3690 ENDDO 3691 ENDDO ! iat 3692 IF (iverbosity == 1) WRITE(stdout, '(a,i12)') 'icount =', icount 3693 ! 3694 ! Write cube file (everything in Bohr) 3695 ! 3696 OPEN(UNIT = iun_plot, FILE=TRIM(wancube), FORM='formatted', STATUS='unknown') 3697 CALL date_and_tim(cdate, ctime ) 3698 ! First two lines are comments 3699 WRITE(iun_plot, *) ' Generated by EPW code' 3700 WRITE(iun_plot, *) ' On ', cdate, ' at ', ctime 3701 ! Number of atoms, origin of cube (Cartesians) wrt simulation (home) cell 3702 WRITE(iun_plot, '(i4,3f13.5)') icount, orig(1), orig(2), orig(3) 3703 ! 3704 ! Number of grid points in each direction, lattice vector 3705 ! 3706 WRITE(iun_plot, '(i4,3f13.5)') ilength(1), & 3707 at(1, 1) * alat / (REAL(ngx, KIND = DP)), & 3708 at(2, 1) * alat / (REAL(ngx, KIND = DP)), & 3709 at(3, 1) * alat / (REAL(ngx, KIND = DP)) 3710 WRITE(iun_plot, '(i4,3f13.5)') ilength(2), & 3711 at(1, 2) * alat / (REAL(ngy, KIND = DP)), & 3712 at(2, 2) * alat / (REAL(ngy, KIND = DP)), & 3713 at(3, 2) * alat / (REAL(ngy, KIND = DP)) 3714 WRITE(iun_plot, '(i4,3f13.5)') ilength(3), & 3715 at(1, 3) * alat / (REAL(ngz, KIND = DP)), & 3716 at(2, 3) * alat / (REAL(ngz, KIND = DP)), & 3717 at(3, 3) * alat / (REAL(ngz, KIND = DP)) 3718 ! 3719 DO iat = 1, nat 3720 DO nzz = -ngs(3)/2, (ngs(3) + 1)/2 3721 DO nyy = -ngs(2)/2, (ngs(2) + 1)/2 3722 DO nxx = -ngs(1)/2, (ngs(1) + 1)/2 3723 diff(:) = xau(:, iat) - wcf(:) & 3724 + (/REAL(nxx, KIND = DP), REAL(nyy, KIND = DP), REAL(nzz, KIND = DP)/) 3725 difc(:) = diff(:) 3726 CALL cryst_to_cart(1, difc(:), at, 1) 3727 dist = DSQRT(difc(1) * difc(1) + difc(2) * difc(2) + difc(3) * difc(3)) * alat 3728 IF (dist < (wannier_plot_scale * wannier_plot_radius / bohr)) THEN 3729 diff(:) = xau(:, iat) & 3730 + (/REAL(nxx, KIND = DP), REAL(nyy, KIND = DP), REAL(nzz, KIND = DP)/) 3731 difc(:) = diff(:) 3732 CALL cryst_to_cart(1, difc(:), at, 1) 3733 difc(:) = difc(:) * alat 3734 WRITE(iun_plot, '(i4,4f13.5)') atomic_Z(ityp(iat)), val_Q, (difc(i), i=1, 3) 3735 ENDIF 3736 ENDDO 3737 ENDDO 3738 ENDDO 3739 ENDDO ! iat 3740 ! 3741 ! Volumetric data in batches of 6 values per line, 'z'-direction first. 3742 ! 3743 DO nxx = 1, ilength(1) 3744 DO nyy = 1, ilength(2) 3745 DO nzz = 1, ilength(3) 3746 WRITE(iun_plot, '(E13.5)', ADVANCE = 'no') REAL(wann_func(nxx,nyy,nzz,1), KIND = DP) 3747 IF ((MOD(nzz, 6) == 0) .OR. (nzz == ilength(3))) WRITE(iun_plot, '(a)') '' 3748 ENDDO 3749 ENDDO 3750 ENDDO 3751 ! 3752 CLOSE(iun_plot) 3753 ! 3754 DEALLOCATE(atomic_Z, STAT = ierr) 3755 IF (ierr /= 0) CALL errore('internal_cube_format', 'Error deallocating atomic_Z', 1) 3756 ! 3757 RETURN 3758 ! 3759 END SUBROUTINE internal_cube_format 3760 ! 3761 !------------------------------------------------------------------------ 3762 END SUBROUTINE write_plot 3763 !----------------------------------------------------------------------- 3764 ! 3765 !----------------------------------------------------------------------- 3766 SUBROUTINE ylm_expansion() 3767 !----------------------------------------------------------------------- 3768 !! 3769 !! Expansion 3770 !! 3771 USE kinds, ONLY : DP 3772 USE random_numbers,ONLY : randy 3773 USE wannierEPW, ONLY : n_proj, xaxis, zaxis, csph, l_w, mr_w 3774 USE matrix_inversion, ONLY: invmat 3775 USE constants_epw, ONLY : zero 3776 ! 3777 IMPLICIT NONE 3778 ! 3779 ! local variables 3780 INTEGER, PARAMETER :: lmax2 = 16 3781 !! 3782 INTEGER :: i 3783 !! Couter on polarizations 3784 INTEGER :: ir 3785 !! Counter on spherical harmonics 3786 INTEGER :: iw 3787 !! Counter on wannier projections 3788 INTEGER :: ierr 3789 !! Error status 3790 ! 3791 REAL(KIND = DP) :: u(3, 3) 3792 !! 3793 REAL(KIND = DP), ALLOCATABLE :: r(:, :) 3794 !! 3795 REAL(KIND = DP), ALLOCATABLE :: rr(:) 3796 !! 3797 REAL(KIND = DP), ALLOCATABLE :: rp(:, :) 3798 !! 3799 REAL(KIND = DP), ALLOCATABLE :: ylm_w(:) 3800 !! 3801 REAL(KIND = DP), ALLOCATABLE :: ylm(:, :) 3802 !! 3803 REAL(KIND = DP), ALLOCATABLE :: mly(:, :) 3804 !! 3805 ! 3806 ALLOCATE(r(3, lmax2), STAT = ierr) 3807 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error allocating r', 1) 3808 ALLOCATE(rp(3, lmax2), STAT = ierr) 3809 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error allocating rp', 1) 3810 ALLOCATE(rr(lmax2), STAT = ierr) 3811 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error allocating rr', 1) 3812 ALLOCATE(ylm_w(lmax2), STAT = ierr) 3813 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error allocating ylm_w', 1) 3814 ALLOCATE(ylm(lmax2, lmax2), STAT = ierr) 3815 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error allocating ylm', 1) 3816 ALLOCATE(mly(lmax2, lmax2), STAT = ierr) 3817 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error allocating mly', 1) 3818 r(:, :) = zero 3819 rr(:) = zero 3820 rp(:, :) = zero 3821 ylm_w(:) = zero 3822 ylm(:, :) = zero 3823 mly(:, :) = zero 3824 ! 3825 ! generate a set of nr=lmax2 random vectors 3826 DO ir = 1, lmax2 3827 DO i = 1, 3 3828 r(i, ir) = randy() - 0.5d0 3829 ENDDO 3830 ENDDO 3831 rr(:) = r(1, :) * r(1, :) + r(2, :) * r(2, :) + r(3, :) * r(3, :) 3832 !- compute ylm(ir,lm) 3833 CALL ylmr2(lmax2, lmax2, r, rr, ylm) 3834 !- store the inverse of ylm(ir,lm) in mly(lm,ir) 3835 CALL invmat(lmax2, ylm, mly) 3836 !- check that r points are independent 3837 CALL check_inverse(lmax2, ylm, mly) 3838 ! 3839 DO iw = 1, n_proj 3840 ! 3841 !- define the u matrix that rotate the reference frame 3842 CALL set_u_matrix(xaxis(:, iw), zaxis(:, iw), u) 3843 !- find rotated r-vectors 3844 rp(:, :) = MATMUL(u(:, :), r(:, :)) 3845 !- set ylm funtion according to wannier90 (l,mr) indexing in the rotaterd points 3846 CALL ylm_wannier(ylm_w, l_w(iw), mr_w(iw), rp, lmax2) 3847 ! 3848 csph(:, iw) = MATMUL(mly(:, :), ylm_w(:)) 3849 ! 3850 ! 3851 ENDDO 3852 ! 3853 DEALLOCATE(r, STAT = ierr) 3854 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error deallocating r', 1) 3855 DEALLOCATE(rp, STAT = ierr) 3856 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error deallocating rp', 1) 3857 DEALLOCATE(rr, STAT = ierr) 3858 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error deallocating rr', 1) 3859 DEALLOCATE(ylm_w, STAT = ierr) 3860 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error deallocating ylm_w', 1) 3861 DEALLOCATE(ylm, STAT = ierr) 3862 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error deallocating ylm', 1) 3863 DEALLOCATE(mly, STAT = ierr) 3864 IF (ierr /= 0) CALL errore('ylm_expansion', 'Error deallocating mly', 1) 3865 ! 3866 RETURN 3867 ! 3868 !----------------------------------------------------------------------- 3869 END SUBROUTINE ylm_expansion 3870 !----------------------------------------------------------------------- 3871 ! 3872 !----------------------------------------------------------------------- 3873 SUBROUTINE check_inverse(lmax2, ylm, mly) 3874 !----------------------------------------------------------------------- 3875 !! 3876 !! Check the inverse 3877 !! 3878 USE kinds, ONLY : DP 3879 USE constants_epw, ONLY : zero, eps8 3880 ! 3881 IMPLICIT NONE 3882 ! 3883 ! I/O variables 3884 INTEGER, INTENT(in) :: lmax2 3885 !! 3886 REAL(KIND = DP), INTENT(in) :: ylm(lmax2, lmax2) 3887 !! 3888 REAL(KIND = DP), INTENT(in) :: mly(lmax2, lmax2) 3889 !! 3890 ! 3891 ! local variables 3892 INTEGER :: lm 3893 !! Counter on spherical harmonics 3894 INTEGER :: ierr 3895 !! Error status 3896 REAL(KIND = DP) :: capel 3897 !! 3898 REAL(KIND = DP), ALLOCATABLE :: uno(:, :) 3899 !! 3900 ! 3901 ALLOCATE(uno(lmax2, lmax2), STAT = ierr) 3902 IF (ierr /= 0) CALL errore('check_inverse', 'Error allocating uno', 1) 3903 uno(:, :) = zero 3904 ! 3905 uno = MATMUL(mly, ylm) 3906 capel = zero 3907 DO lm = 1, lmax2 3908 uno(lm, lm) = uno(lm, lm) - 1.d0 3909 ENDDO 3910 capel = capel + SUM(ABS(uno(1:lmax2, 1:lmax2))) 3911 IF (capel > eps8) CALL errore('ylm_expansion', & 3912 'Inversion failed: r(*, 1:nr) are not all independent!!', 1) 3913 DEALLOCATE(uno, STAT = ierr) 3914 IF (ierr /= 0) CALL errore('check_inverse', 'Error deallocating uno', 1) 3915 ! 3916 RETURN 3917 ! 3918 !------------------------------------------------------------------------ 3919 END SUBROUTINE check_inverse 3920 !----------------------------------------------------------------------- 3921 ! 3922 !----------------------------------------------------------------------- 3923 SUBROUTINE set_u_matrix(x, z, u) 3924 !----------------------------------------------------------------------- 3925 !! 3926 !! Set the U matrix 3927 !! 3928 USE kinds, ONLY : DP 3929 USE constants_epw, ONLY : eps8 3930 ! 3931 IMPLICIT NONE 3932 ! 3933 ! I/O variables 3934 REAL(KIND = DP), INTENT(in) :: x(3) 3935 !! 3936 REAL(KIND = DP), INTENT(in) :: z(3) 3937 !! 3938 REAL(KIND = DP), INTENT(out) :: u(3, 3) 3939 !! 3940 ! local variables 3941 REAL(KIND = DP) :: coseno 3942 !! cosine between x(3) and z(3) 3943 REAL(KIND = DP) :: xx 3944 !! Norm of x(3) 3945 REAL(KIND = DP) :: zz 3946 !! Norm of z(3) 3947 REAL(KIND = DP) :: y(3) 3948 !! 3949 ! 3950 xx = DSQRT(SUM(x(:) * x(:))) 3951 IF (xx < eps8) CALL errore('set_u_matrix', '|xaxis| < eps8', 1) 3952 ! 3953 zz = DSQRT(SUM(z(:) * z(:))) 3954 IF (zz < eps8) CALL errore ('set_u_matrix', '|zaxis| < eps8', 1) 3955 ! 3956 coseno = SUM(x(:) * z(:)) / xx / zz 3957 IF (ABS(coseno) > eps8) CALL errore('set_u_matrix', 'xaxis and zaxis are not orthogonal!', 1) 3958 ! 3959 y(1) = (z(2) * x(3) - x(2) * z(3)) / xx / zz 3960 y(2) = (z(3) * x(1) - x(3) * z(1)) / xx / zz 3961 y(3) = (z(1) * x(2) - x(1) * z(2)) / xx / zz 3962 ! 3963 u(1, :) = x(:) / xx 3964 u(2, :) = y(:) 3965 u(3, :) = z(:) / zz 3966 ! 3967 ! 3968 RETURN 3969 ! 3970 !------------------------------------------------------------------------ 3971 END SUBROUTINE set_u_matrix 3972 !----------------------------------------------------------------------- 3973 ! 3974 !----------------------------------------------------------------------- 3975 SUBROUTINE ylm_wannier(ylm, l, mr, r, nr) 3976 !----------------------------------------------------------------------- 3977 !! 3978 !! This routine returns in ylm(r) the values at the nr points r(1:3,1:nr) 3979 !! of the spherical harmonic identified by indices (l,mr) 3980 !! in table 3.1 of the wannierf90 specification. 3981 !! 3982 !! No reference to the particular ylm ordering internal to quantum-espresso 3983 !! is assumed. 3984 !! 3985 !! If ordering in wannier90 code is changed or extended this should be the 3986 !! only place to be modified accordingly 3987 !! 3988 USE kinds, ONLY : DP 3989 USE constants_epw, ONLY : pibytwo, eps8 3990 USE constants, ONLY : pi 3991 USE low_lvl, ONLY : fy3x2my2, fxx2m3y2, fxx2m3y2, fxyz, fzx2my2, & 3992 fyz2, fxz2, fz3, dxy, dx2my2, dyz, dxz, dz2, & 3993 p_z, py, px, s 3994 ! 3995 IMPLICIT NONE 3996 ! 3997 ! I/O variables 3998 INTEGER, INTENT(in) :: l, mr 3999 !! quantum numbers 4000 INTEGER, INTENT(in) :: nr 4001 !! 4002 ! 4003 REAL(KIND = DP), INTENT(in) :: r(3, nr) 4004 !! 4005 REAL(KIND = DP), INTENT(out) :: ylm(nr) 4006 !! spherical harmonics 4007 ! 4008 ! local variables 4009 INTEGER :: ir 4010 !! 4011 REAL(KIND = DP) :: rr 4012 !! 4013 REAL(KIND = DP) :: cost 4014 !! 4015 REAL(KIND = DP) :: phi 4016 !! 4017 REAL(KIND = DP) :: bs2, bs3, bs6, bs12 4018 !! Inverse DSQRT of 2, 3, 6, and 12 4019 ! 4020 bs2 = 1.d0 / DSQRT(2.d0) 4021 bs3 = 1.d0 / DSQRT(3.d0) 4022 bs6 = 1.d0 / DSQRT(6.d0) 4023 bs12 = 1.d0 / DSQRT(12.d0) 4024 ! 4025 IF (l > 3 .OR. l < -5) CALL errore('ylm_wannier', 'l out of range ', 1) 4026 IF (l >= 0) THEN 4027 IF (mr < 1 .OR. mr > 2 * l + 1) CALL errore('ylm_wannier', 'mr out of range' ,1) 4028 ELSE 4029 IF (mr < 1 .OR. mr > ABS(l) + 1 ) CALL errore('ylm_wannier', 'mr out of range', 1) 4030 ENDIF 4031 ! 4032 DO ir = 1, nr 4033 rr = DSQRT(SUM(r(:, ir) * r(:, ir))) 4034 IF (rr < eps8) CALL errore('ylm_wannier', 'rr too small', 1) 4035 ! 4036 cost = r(3, ir) / rr 4037 ! 4038 ! beware the arc tan, it is defined modulo pi 4039 ! 4040 IF (r(1, ir) > eps8) THEN 4041 phi = ATAN(r(2, ir) / r(1, ir)) 4042 ELSEIF (r(1, ir) < -eps8) THEN 4043 phi = ATAN(r(2, ir) / r(1, ir)) + pi 4044 ELSE 4045 phi = SIGN(pibytwo, r(2, ir)) 4046 ENDIF 4047 ! 4048 IF (l == 0) THEN ! s orbital 4049 ylm(ir) = s() 4050 ENDIF 4051 IF (l == 1) THEN ! p orbitals 4052 IF (mr == 1) ylm(ir) = p_z(cost) 4053 IF (mr == 2) ylm(ir) = px(cost, phi) 4054 IF (mr == 3) ylm(ir) = py(cost, phi) 4055 ENDIF 4056 IF (l == 2) THEN ! d orbitals 4057 IF (mr == 1) ylm(ir) = dz2(cost) 4058 IF (mr == 2) ylm(ir) = dxz(cost, phi) 4059 IF (mr == 3) ylm(ir) = dyz(cost, phi) 4060 IF (mr == 4) ylm(ir) = dx2my2(cost, phi) 4061 IF (mr == 5) ylm(ir) = dxy(cost, phi) 4062 ENDIF 4063 IF (l == 3) THEN ! f orbitals 4064 IF (mr == 1) ylm(ir) = fz3(cost) 4065 IF (mr == 2) ylm(ir) = fxz2(cost, phi) 4066 IF (mr == 3) ylm(ir) = fyz2(cost, phi) 4067 IF (mr == 4) ylm(ir) = fzx2my2(cost, phi) 4068 IF (mr == 5) ylm(ir) = fxyz(cost, phi) 4069 IF (mr == 6) ylm(ir) = fxx2m3y2(cost, phi) 4070 IF (mr == 7) ylm(ir) = fy3x2my2(cost, phi) 4071 ENDIF 4072 IF (l == -1) THEN ! sp hybrids 4073 IF (mr == 1) ylm(ir) = bs2 * (s() + px(cost, phi)) 4074 IF (mr == 2) ylm(ir) = bs2 * (s() - px(cost, phi)) 4075 ENDIF 4076 IF (l == -2) THEN ! sp2 hybrids 4077 IF (mr == 1) ylm(ir) = bs3 * s() - bs6 * px(cost, phi) + bs2 * py(cost, phi) 4078 IF (mr == 2) ylm(ir) = bs3 * s() - bs6 * px(cost, phi) - bs2 * py(cost, phi) 4079 IF (mr == 3) ylm(ir) = bs3 * s() + 2.d0 *bs6 * px(cost, phi) 4080 ENDIF 4081 IF (l == -3) THEN ! sp3 hybrids 4082 IF (mr == 1) ylm(ir) = 0.5d0 *(s() + px(cost, phi) + py(cost, phi) + p_z(cost)) 4083 IF (mr == 2) ylm(ir) = 0.5d0 *(s() + px(cost, phi) - py(cost, phi) - p_z(cost)) 4084 IF (mr == 3) ylm(ir) = 0.5d0 *(s() - px(cost, phi) + py(cost, phi) - p_z(cost)) 4085 IF (mr == 4) ylm(ir) = 0.5d0 *(s() - px(cost, phi) - py(cost, phi) + p_z(cost)) 4086 ENDIF 4087 IF (l == -4) THEN ! sp3d hybrids 4088 IF (mr == 1) ylm(ir) = bs3 * s() -bs6 * px(cost, phi) + bs2 * py(cost, phi) 4089 IF (mr == 2) ylm(ir) = bs3 * s() -bs6 * px(cost, phi) - bs2 * py(cost, phi) 4090 IF (mr == 3) ylm(ir) = bs3 * s() +2.d0 * bs6 * px(cost, phi) 4091 IF (mr == 4) ylm(ir) = bs2 * p_z(cost) + bs2 * dz2(cost) 4092 IF (mr == 5) ylm(ir) =-bs2 * p_z(cost) + bs2 * dz2(cost) 4093 ENDIF 4094 IF (l == -5) THEN ! sp3d2 hybrids 4095 IF (mr == 1) ylm(ir) = bs6 *s() - bs2 * px(cost, phi) - bs12 * dz2(cost) + .5 * dx2my2(cost, phi) 4096 IF (mr == 2) ylm(ir) = bs6 *s() + bs2 * px(cost, phi) - bs12 * dz2(cost) + .5 * dx2my2(cost, phi) 4097 IF (mr == 3) ylm(ir) = bs6 *s() - bs2 * py(cost, phi) - bs12 * dz2(cost) - .5 * dx2my2(cost, phi) 4098 IF (mr == 4) ylm(ir) = bs6 *s() + bs2 * py(cost, phi) - bs12 * dz2(cost) - .5 * dx2my2(cost, phi) 4099 IF (mr == 5) ylm(ir) = bs6 *s() - bs2 * p_z(cost) + bs3 * dz2(cost) 4100 IF (mr == 6) ylm(ir) = bs6 *s() + bs2 * p_z(cost) + bs3 * dz2(cost) 4101 ENDIF 4102 ! 4103 ENDDO 4104 ! 4105 RETURN 4106 ! 4107 !----------------------------------------------------------------------- 4108 END SUBROUTINE ylm_wannier 4109 !----------------------------------------------------------------------- 4110 ! 4111 !----------------------------------------------------------------------- 4112 SUBROUTINE radialpart(ng, q, alfa, rvalue, lmax, radial) 4113 !----------------------------------------------------------------------- 4114 !! 4115 !! This routine computes a table with the radial Fourier transform 4116 !! of the radial functions. 4117 !! 4118 USE kinds, ONLY : DP 4119 USE constants, ONLY : fpi 4120 USE cell_base, ONLY : omega 4121 USE constants_epw, ONLY : zero 4122 ! 4123 IMPLICIT NONE 4124 ! 4125 ! I/O 4126 INTEGER, INTENT(in) :: ng 4127 !! Number of plane waves 4128 INTEGER, INTENT(in) :: rvalue 4129 !! 4130 INTEGER, INTENT(in) :: lmax 4131 !! Maxium angular momentum 4132 ! 4133 REAL(KIND = DP), INTENT(in) :: alfa 4134 !! 4135 REAL(KIND = DP), INTENT(in) :: q(ng) 4136 !! 4137 REAL(KIND = DP), INTENT(out) :: radial(ng, 0:lmax) 4138 !! 4139 ! 4140 ! local variables 4141 INTEGER :: l 4142 !! Counter on angular momentum 4143 INTEGER :: ig 4144 !! Counter on plane waves 4145 INTEGER :: ir 4146 !! Counter on mesh_r 4147 INTEGER :: mesh_r 4148 !! Size of the radial FFT mesh 4149 INTEGER :: ierr 4150 !! Error status 4151 ! 4152 REAL(KIND = DP), PARAMETER :: xmin = -6.d0 4153 !! 4154 REAL(KIND = DP), PARAMETER :: dx = 0.025d0 4155 !! 4156 REAL(KIND = DP), PARAMETER :: rmax = 10.d0 4157 !! 4158 REAL(KIND = DP) :: rad_int 4159 !! 4160 REAL(KIND = DP) :: pref 4161 !! 4162 REAL(KIND = DP) :: x 4163 !! 4164 REAL(KIND = DP), ALLOCATABLE :: bes(:) 4165 !! 4166 REAL(KIND = DP), ALLOCATABLE :: func_r(:) 4167 !! 4168 REAL(KIND = DP), ALLOCATABLE :: r(:) 4169 !! 4170 REAL(KIND = DP), ALLOCATABLE :: rij(:) 4171 !! 4172 REAL(KIND = DP), ALLOCATABLE :: aux(:) 4173 !! 4174 ! 4175 mesh_r = NINT((log(rmax) - xmin) / dx + 1) 4176 ! 4177 ALLOCATE(bes(mesh_r), STAT = ierr) 4178 IF (ierr /= 0) CALL errore('radialpart', 'Error allocating bes', 1) 4179 ALLOCATE(func_r(mesh_r), STAT = ierr) 4180 IF (ierr /= 0) CALL errore('radialpart', 'Error allocating func_r', 1) 4181 ALLOCATE(r(mesh_r), STAT = ierr) 4182 IF (ierr /= 0) CALL errore('radialpart', 'Error allocating r', 1) 4183 ALLOCATE(rij(mesh_r), STAT = ierr) 4184 IF (ierr /= 0) CALL errore('radialpart', 'Error allocating rij', 1) 4185 ALLOCATE(aux(mesh_r), STAT = ierr) 4186 IF (ierr /= 0) CALL errore('radialpart', 'Error allocating aux', 1) 4187 bes(:) = zero 4188 func_r(:) = zero 4189 r(:) = zero 4190 rij(:) = zero 4191 aux(:) = zero 4192 ! 4193 ! compute the radial mesh 4194 ! 4195 DO ir = 1, mesh_r 4196 x = xmin + DBLE(ir - 1) * dx 4197 r(ir) = EXP(x) / alfa 4198 rij(ir) = dx * r(ir) 4199 ENDDO 4200 ! 4201 IF (rvalue == 1) func_r(:) = 2.d0 * alfa**(3.d0 / 2.d0) * EXP(-alfa * r(:)) 4202 IF (rvalue == 2) func_r(:) = 1.d0 / DSQRT(8.d0) * alfa**(3.d0 / 2.d0) * & 4203 (2.0d0 - alfa * r(:)) * EXP(-alfa * r(:) * 0.5d0) 4204 IF (rvalue == 3) func_r(:) = DSQRT(4.d0 / 27.d0) * alfa**(2.0d0 / 3.0d0) * & 4205 (1.d0 - 1.5d0 * alfa * r(:) + & 4206 2.d0 * (alfa * r(:))**2 / 27.d0) * & 4207 EXP(-alfa * r(:) / 3.0d0) 4208 pref = fpi / DSQRT(omega) 4209 ! 4210 DO l = 0, lmax 4211 DO ig = 1, ng 4212 CALL sph_bes(mesh_r, r(1), q(ig), l, bes) 4213 aux(:) = bes(:) * func_r(:) * r(:) * r(:) 4214 CALL simpson(mesh_r, aux, rij, rad_int) 4215 radial(ig, l) = rad_int * pref 4216 ENDDO 4217 ENDDO 4218 ! 4219 DEALLOCATE(bes, STAT = ierr) 4220 IF (ierr /= 0) CALL errore('radialpart', 'Error deallocating bes', 1) 4221 DEALLOCATE(func_r, STAT = ierr) 4222 IF (ierr /= 0) CALL errore('radialpart', 'Error deallocating func_r', 1) 4223 DEALLOCATE(r, STAT = ierr) 4224 IF (ierr /= 0) CALL errore('radialpart', 'Error deallocating r', 1) 4225 DEALLOCATE(rij, STAT = ierr) 4226 IF (ierr /= 0) CALL errore('radialpart', 'Error deallocating rij', 1) 4227 DEALLOCATE(aux, STAT = ierr) 4228 IF (ierr /= 0) CALL errore('radialpart', 'Error deallocating aux', 1) 4229 ! 4230 RETURN 4231 ! 4232 !----------------------------------------------------------------------- 4233 END SUBROUTINE radialpart 4234 !----------------------------------------------------------------------- 4235 !------------------------------------------------------------------------------ 4236 END MODULE pw2wan2epw 4237 !------------------------------------------------------------------------------ 4238