1 ! 2 ! Copyright (C) 2010-2016 Samuel Ponce', Roxana Margine, Carla Verdi, Feliciano Giustino 3 ! 4 ! This file is distributed under the terms of the GNU General Public 5 ! License. See the file `LICENSE' in the root directory of the 6 ! present distribution, or http://www.gnu.org/copyleft.gpl.txt . 7 ! 8 ! 9 !---------------------------------------------------------------------- 10 MODULE io_epw 11 !---------------------------------------------------------------------- 12 !! 13 !! This module contains various writing or reading routines to files. 14 !! Most of them are for restart purposes. 15 !! 16 IMPLICIT NONE 17 ! 18 CONTAINS 19 ! 20 !----------------------------------------------------------------- 21 SUBROUTINE rwepmatw(epmatw, nbnd, np, nmodes, nrec, iun, iop) 22 !----------------------------------------------------------------- 23 !! 24 !! A simple wrapper to the davcio routine to read/write arrays 25 !! instead of vectors 26 !! 27 !----------------------------------------------------------------- 28 USE kinds, ONLY : DP 29 USE mp, ONLY : mp_barrier 30 ! 31 IMPLICIT NONE 32 ! 33 INTEGER, INTENT(in) :: nbnd 34 !! Total number of bands 35 INTEGER, INTENT(in) :: np 36 !! np is either nrr_k or nq (epmatwe and epmatwp have the same structure) 37 INTEGER, INTENT(in) :: nmodes 38 !! Number of modes 39 INTEGER, INTENT(in) :: nrec 40 !! Place where to start reading/writing 41 INTEGER, INTENT(in) :: iun 42 !! Record number 43 INTEGER, INTENT(in) :: iop 44 !! If -1, read and if +1 write the matrix 45 COMPLEX(KIND = DP), INTENT(inout) :: epmatw(nbnd, nbnd, np, nmodes) 46 !! El-ph matrix to read or write 47 ! 48 ! Local variables 49 INTEGER :: lrec 50 !! Record length 51 INTEGER :: i 52 !! Index number 53 INTEGER :: ibnd 54 !! Band index 55 INTEGER :: jbnd 56 !! Band index 57 INTEGER :: imode 58 !! Mode index 59 INTEGER :: ip 60 !! REal space index (either nrr_k or nq) 61 COMPLEX(KIND = DP):: aux(nbnd * nbnd * np * nmodes) 62 !! 1-D vector to store the matrix elements. 63 ! 64 lrec = 2 * nbnd * nbnd * np * nmodes 65 ! 66 IF (iop == -1) then 67 ! 68 ! read matrix 69 ! 70 CALL davcio(aux, lrec, iun, nrec, -1) 71 ! 72 i = 0 73 DO imode = 1, nmodes 74 DO ip = 1, np 75 DO jbnd = 1, nbnd 76 DO ibnd = 1, nbnd 77 i = i + 1 78 epmatw(ibnd, jbnd, ip, imode) = aux(i) 79 ! 80 ENDDO 81 ENDDO 82 ENDDO 83 ENDDO 84 ! 85 ELSEIF (iop == 1) THEN 86 ! 87 ! write matrix 88 ! 89 i = 0 90 DO imode = 1, nmodes 91 DO ip = 1, np 92 DO jbnd = 1, nbnd 93 DO ibnd = 1, nbnd 94 i = i + 1 95 aux(i) = epmatw(ibnd, jbnd, ip, imode) 96 ENDDO 97 ENDDO 98 ENDDO 99 ENDDO 100 ! 101 CALL davcio(aux, lrec, iun, nrec, +1) 102 ! 103 ELSE 104 ! 105 CALL errore('rwepmatw', 'iop not permitted', 1) 106 ! 107 ENDIF 108 ! 109 !---------------------------------------------------------------------- 110 END SUBROUTINE rwepmatw 111 !---------------------------------------------------------------------- 112 ! 113 !---------------------------------------------------------------------- 114 SUBROUTINE epw_write(nrr_k, nrr_q, nrr_g, w_centers) 115 !---------------------------------------------------------------------- 116 !! 117 !! Routine to write files on real-space grid for fine grid interpolation 118 !! 119 USE kinds, ONLY : DP 120 USE epwcom, ONLY : nbndsub, vme, eig_read, etf_mem 121 USE pwcom, ONLY : ef, nelec 122 USE elph2, ONLY : chw, rdw, cdmew, cvmew, chw_ks, & 123 zstar, epsi, epmatwp 124 USE ions_base, ONLY : amass, ityp, nat, tau 125 USE cell_base, ONLY : at, bg, omega, alat 126 USE modes, ONLY : nmodes 127 USE io_var, ONLY : epwdata, iundmedata, iunvmedata, iunksdata, iunepmatwp, & 128 crystal 129 USE noncollin_module, ONLY : noncolin 130 USE io_files, ONLY : prefix, diropn, tmp_dir 131 USE mp, ONLY : mp_barrier 132 USE mp_world, ONLY : mpime 133 USE io_global, ONLY : ionode_id, stdout 134 ! 135 IMPLICIT NONE 136 ! 137 INTEGER, INTENT(in) :: nrr_k 138 !! Number of WS vectors for the electrons 139 INTEGER, INTENT(in) :: nrr_q 140 !! Number of WS vectors for the phonons 141 INTEGER, INTENT(in) :: nrr_g 142 !! Number of WS vectors for the electron-phonons 143 REAL(KIND = DP), INTENT(in) :: w_centers(3, nbndsub) 144 !! Wannier center 145 ! 146 ! Local variables 147 CHARACTER(LEN = 256) :: filint 148 !! Name of the file 149 LOGICAL :: exst 150 !! The file exists 151 INTEGER :: ibnd, jbnd 152 !! Band index 153 INTEGER :: jmode, imode 154 !! Mode index 155 INTEGER :: irk, irq, irg 156 !! WS vector looping index on electron, phonons and el-ph 157 INTEGER :: ipol 158 !! Cartesian direction (polarison direction) 159 INTEGER :: lrepmatw 160 !! Record length 161 INTEGER*8 :: unf_recl 162 !! Record length 163 INTEGER :: direct_io_factor 164 !! Type of IOlength 165 INTEGER :: ierr 166 !! Error index 167 REAL(KIND = DP) :: dummy 168 !! Dummy variable 169 ! 170 WRITE(stdout,'(/5x,"Writing Hamiltonian, Dynamical matrix and EP vertex in Wann rep to file"/)') 171 ! 172 IF (mpime == ionode_id) THEN 173 ! 174 OPEN(UNIT = epwdata, FILE = 'epwdata.fmt') 175 OPEN(UNIT = crystal, FILE = 'crystal.fmt') 176 IF (vme) THEN 177 OPEN(UNIT = iunvmedata, FILE = 'vmedata.fmt') 178 ELSE 179 OPEN(UNIT = iundmedata, FILE = 'dmedata.fmt') 180 ENDIF 181 IF (eig_read) OPEN(UNIT = iunksdata, FILE = 'ksdata.fmt') 182 WRITE(crystal,*) nat 183 WRITE(crystal,*) nmodes 184 WRITE(crystal,*) nelec 185 WRITE(crystal,*) at 186 WRITE(crystal,*) bg 187 WRITE(crystal,*) omega 188 WRITE(crystal,*) alat 189 WRITE(crystal,*) tau 190 WRITE(crystal,*) amass 191 WRITE(crystal,*) ityp 192 WRITE(crystal,*) noncolin 193 WRITE(crystal,*) w_centers 194 ! 195 WRITE(epwdata,*) ef 196 WRITE(epwdata,*) nbndsub, nrr_k, nmodes, nrr_q, nrr_g 197 WRITE(epwdata,*) zstar, epsi 198 ! 199 DO ibnd = 1, nbndsub 200 DO jbnd = 1, nbndsub 201 DO irk = 1, nrr_k 202 WRITE (epwdata,*) chw(ibnd, jbnd, irk) 203 IF (eig_read) WRITE (iunksdata,*) chw_ks(ibnd, jbnd, irk) 204 DO ipol = 1, 3 205 IF (vme) THEN 206 WRITE(iunvmedata,*) cvmew(ipol, ibnd, jbnd, irk) 207 ELSE 208 WRITE(iundmedata,*) cdmew(ipol, ibnd, jbnd, irk) 209 ENDIF 210 ENDDO 211 ENDDO 212 ENDDO 213 ENDDO 214 ! 215 DO imode = 1, nmodes 216 DO jmode = 1, nmodes 217 DO irq = 1, nrr_q 218 WRITE(epwdata,*) rdw(imode, jmode, irq) 219 ENDDO 220 ENDDO 221 ENDDO 222 ! 223 IF (etf_mem == 0) THEN 224 ! SP: The call to epmatwp is now inside the loop 225 ! This is important as otherwise the lrepmatw INTEGER 226 ! could become too large for integer(kind=4). 227 ! Note that in Fortran the record length has to be a integer 228 ! of kind 4. 229 lrepmatw = 2 * nbndsub * nbndsub * nrr_k * nmodes 230 filint = TRIM(tmp_dir) // TRIM(prefix)//'.epmatwp' 231 INQUIRE(IOLENGTH = direct_io_factor) dummy 232 unf_recl = direct_io_factor * INT(lrepmatw, KIND = KIND(unf_recl)) 233 IF (unf_recl <= 0) CALL errore('epw_write', 'wrong record length', 3) 234 OPEN(iunepmatwp, FILE = TRIM(ADJUSTL(filint)), IOSTAT = ierr, form='unformatted', & 235 STATUS = 'unknown', ACCESS = 'direct', RECL = unf_recl) 236 IF (ierr /= 0) CALL errore('epw_write', 'error opening ' // TRIM(filint), 1) 237 ! 238 !CALL diropn(iunepmatwp, 'epmatwp', lrepmatw, exst) 239 DO irg = 1, nrr_g 240 CALL davcio(epmatwp(:, :, :, :, irg), lrepmatw, iunepmatwp, irg, +1) 241 ENDDO 242 ! 243 CLOSE(iunepmatwp) 244 ENDIF 245 ! 246 CLOSE(epwdata) 247 CLOSE(crystal) 248 IF (vme) THEN 249 CLOSE(iunvmedata) 250 ELSE 251 CLOSE(iundmedata) 252 ENDIF 253 IF (eig_read) CLOSE(iunksdata) 254 ! 255 ENDIF 256 !-------------------------------------------------------------------------------- 257 END SUBROUTINE epw_write 258 !-------------------------------------------------------------------------------- 259 ! 260 !-------------------------------------------------------------------------------- 261 SUBROUTINE epw_read(nrr_k, nrr_q, nrr_g) 262 !-------------------------------------------------------------------------------- 263 !! 264 !! Routine to read the real space quantities for fine grid interpolation 265 !! 266 USE kinds, ONLY : DP 267 USE epwcom, ONLY : nbndsub, vme, eig_read, etf_mem, lifc 268 USE pwcom, ONLY : ef 269 USE elph2, ONLY : chw, rdw, epmatwp, cdmew, cvmew, chw_ks, zstar, epsi 270 USE ions_base, ONLY : nat 271 USE modes, ONLY : nmodes 272 USE io_global, ONLY : stdout 273 USE io_files, ONLY : prefix, diropn, tmp_dir 274 USE io_var, ONLY : epwdata, iundmedata, iunvmedata, iunksdata, iunepmatwp 275 USE constants_epw, ONLY : czero, zero 276#if defined(__NAG) 277 USE f90_unix_io,ONLY : flush 278#endif 279 USE io_global, ONLY : ionode_id 280 USE mp, ONLY : mp_barrier, mp_bcast 281 USE mp_global, ONLY : world_comm 282 USE mp_world, ONLY : mpime 283 ! 284 IMPLICIT NONE 285 ! 286 INTEGER, INTENT(out) :: nrr_k 287 !! Number of WS vectors for the electrons 288 INTEGER, INTENT(out) :: nrr_q 289 !! Number of WS vectors for the phonons 290 INTEGER, INTENT(out) :: nrr_g 291 !! Number of WS vectors for the electron-phonons 292 ! 293 ! Local variables 294 ! 295 CHARACTER(LEN = 256) :: filint 296 !! Name of the file 297 LOGICAL :: exst 298 !! The file exists 299 INTEGER :: ibnd, jbnd 300 !! Band index 301 INTEGER :: jmode, imode 302 !! Mode index 303 INTEGER :: irk, irq, irg 304 !! WS vector looping index on electron, phonons and el-ph 305 INTEGER :: ipol 306 !! Cartesian direction (polarison direction) 307 INTEGER :: lrepmatw 308 !! Record length 309 INTEGER :: ios 310 !! Status of files 311 INTEGER :: ierr 312 !! Error status 313 INTEGER*8 :: unf_recl 314 !! Record length 315 INTEGER :: direct_io_factor 316 !! Type of IOlength 317 REAL(KIND = DP) :: dummy 318 !! Dummy variable 319 320 ! 321 WRITE(stdout,'(/5x,"Reading Hamiltonian, Dynamical matrix and EP vertex in Wann rep from file"/)') 322 FLUSH(stdout) 323 ! 324 ! This is important in restart mode as zstar etc has not been allocated 325 ALLOCATE(zstar(3, 3, nat), STAT = ierr) 326 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating zstar', 1) 327 ALLOCATE(epsi(3, 3), STAT = ierr) 328 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating epsi', 1) 329 ! 330 IF (mpime == ionode_id) THEN 331 ! 332 OPEN(UNIT = epwdata, FILE = 'epwdata.fmt', STATUS = 'old', IOSTAT = ios) 333 IF (ios /= 0) CALL errore ('epw_read', 'error opening epwdata.fmt', epwdata) 334 IF (eig_read) OPEN(UNIT = iunksdata, FILE = 'ksdata.fmt', STATUS = 'old', IOSTAT = ios) 335 IF (eig_read .AND. ios /= 0) CALL errore ('epw_read', 'error opening ksdata.fmt', iunksdata) 336 IF (vme) THEN 337 OPEN(UNIT = iunvmedata, FILE = 'vmedata.fmt', STATUS = 'old', IOSTAT = ios) 338 IF (ios /= 0) CALL errore ('epw_read', 'error opening vmedata.fmt', iunvmedata) 339 ELSE 340 OPEN(UNIT = iundmedata, FILE = 'dmedata.fmt', STATUS = 'old', IOSTAT = ios) 341 IF (ios /= 0) CALL errore ('epw_read', 'error opening dmedata.fmt', iundmedata) 342 ENDIF 343 READ(epwdata,*) ef 344 READ(epwdata,*) nbndsub, nrr_k, nmodes, nrr_q, nrr_g 345 READ(epwdata,*) zstar, epsi 346 ! 347 ENDIF 348 CALL mp_bcast(ef, ionode_id, world_comm) 349 CALL mp_bcast(nbndsub, ionode_id, world_comm) 350 CALL mp_bcast(nrr_k, ionode_id, world_comm) 351 CALL mp_bcast(nmodes, ionode_id, world_comm) 352 CALL mp_bcast(nrr_q, ionode_id, world_comm) 353 CALL mp_bcast(nrr_g, ionode_id, world_comm) 354 CALL mp_bcast(zstar, ionode_id, world_comm) 355 CALL mp_bcast(epsi, ionode_id, world_comm) 356 ! 357 ALLOCATE(chw(nbndsub, nbndsub, nrr_k), STAT = ierr) 358 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating chw', 1) 359 ALLOCATE(chw_ks(nbndsub, nbndsub, nrr_k), STAT = ierr) 360 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating chw_ks', 1) 361 ALLOCATE(rdw(nmodes, nmodes, nrr_q), STAT = ierr) 362 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating rdw', 1) 363 IF (vme) THEN 364 ALLOCATE(cvmew(3, nbndsub, nbndsub, nrr_k), STAT = ierr) 365 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating cvmew', 1) 366 ELSE 367 ALLOCATE(cdmew(3, nbndsub, nbndsub, nrr_k), STAT = ierr) 368 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating cdmew', 1) 369 ENDIF 370 ! 371 IF (mpime == ionode_id) THEN 372 ! 373 DO ibnd = 1, nbndsub 374 DO jbnd = 1, nbndsub 375 DO irk = 1, nrr_k 376 READ(epwdata,*) chw(ibnd, jbnd, irk) 377 IF (eig_read) READ(iunksdata,*) chw_ks(ibnd, jbnd, irk) 378 DO ipol = 1,3 379 IF (vme) THEN 380 READ(iunvmedata,*) cvmew(ipol, ibnd, jbnd, irk) 381 ELSE 382 READ(iundmedata,*) cdmew(ipol, ibnd, jbnd, irk) 383 ENDIF 384 ENDDO 385 ENDDO 386 ENDDO 387 ENDDO 388 ! 389 IF (.NOT. lifc) THEN 390 DO imode = 1, nmodes 391 DO jmode = 1, nmodes 392 DO irq = 1, nrr_q 393 READ(epwdata,*) rdw(imode, jmode, irq) 394 ENDDO 395 ENDDO 396 ENDDO 397 ENDIF 398 ! 399 ENDIF 400 ! 401 CALL mp_bcast(chw, ionode_id, world_comm) 402 ! 403 IF (eig_read) CALL mp_bcast(chw_ks, ionode_id, world_comm) 404 IF (.NOT. lifc) CALL mp_bcast(rdw, ionode_id, world_comm) 405 ! 406 IF (vme) THEN 407 CALL mp_bcast(cvmew, ionode_id, world_comm) 408 ELSE 409 CALL mp_bcast(cdmew, ionode_id, world_comm) 410 ENDIF 411 ! 412 IF (lifc) THEN 413 CALL read_ifc_epw 414 ENDIF 415 ! 416 IF (etf_mem == 0) THEN 417 ALLOCATE(epmatwp(nbndsub, nbndsub, nrr_k, nmodes, nrr_g), STAT = ierr) 418 IF (ierr /= 0) CALL errore('epw_read', 'Error allocating epmatwp', 1) 419 epmatwp = czero 420 IF (mpime == ionode_id) THEN 421 ! SP: The call to epmatwp is now inside the loop 422 ! This is important as otherwise the lrepmatw INTEGER 423 ! could become too large for integer(kind=4). 424 ! Note that in Fortran the record length has to be a integer 425 ! of kind 4. 426 lrepmatw = 2 * nbndsub * nbndsub * nrr_k * nmodes 427 filint = TRIM(tmp_dir) // TRIM(prefix)//'.epmatwp' 428 ! 429 INQUIRE(IOLENGTH = direct_io_factor) dummy 430 unf_recl = direct_io_factor * INT(lrepmatw, KIND = KIND(unf_recl)) 431 IF (unf_recl <= 0) CALL errore('epw_read', 'wrong record length', 3) 432 OPEN(iunepmatwp, FILE = TRIM(ADJUSTL(filint)), IOSTAT = ierr, FORM = 'unformatted', & 433 STATUS = 'unknown', ACCESS = 'direct', RECL = unf_recl) 434 IF (ierr /= 0) CALL errore('epw_read', 'error opening ' // TRIM(filint), 1) 435 ! 436 DO irg = 1, nrr_g 437 CALL davcio(epmatwp(:, :, :, :, irg), lrepmatw, iunepmatwp, irg, -1) 438 ENDDO 439 ! 440 CLOSE(iunepmatwp) 441 ENDIF 442 ! 443 CALL mp_bcast(epmatwp, ionode_id, world_comm) 444 ! 445 ENDIF 446 ! 447 !CALL mp_barrier(inter_pool_comm) 448 IF (mpime == ionode_id) THEN 449 CLOSE(epwdata) 450 IF (vme) THEN 451 CLOSE(iunvmedata) 452 ELSE 453 CLOSE(iundmedata) 454 ENDIF 455 ENDIF 456 ! 457 WRITE(stdout, '(/5x,"Finished reading Wann rep data from file"/)') 458 ! 459 !------------------------------------------------------------------------------ 460 END SUBROUTINE epw_read 461 !------------------------------------------------------------------------------ 462 ! 463 !--------------------------------------------------------------------------------- 464 SUBROUTINE read_ifc_epw() 465 !--------------------------------------------------------------------------------- 466 !! 467 !! Read IFC in real space from the file generated by q2r. 468 !! Adapted from PH/matdyn.x by C. Verdi and S. Ponce 469 !! 470 ! 471 USE kinds, ONLY : DP 472 USE elph2, ONLY : ifc, zstar, epsi 473 USE epwcom, ONLY : asr_typ, dvscf_dir, nqc1, nqc2, nqc3 474 USE ions_base, ONLY : nat 475 USE cell_base, ONLY : ibrav, omega, at, bg, celldm, alat 476 USE io_global, ONLY : stdout 477 USE io_var, ONLY : iunifc 478 USE noncollin_module, ONLY : nspin_mag 479 USE io_global, ONLY : ionode_id, ionode 480 USE mp, ONLY : mp_barrier, mp_bcast 481 USE mp_global, ONLY : intra_pool_comm, inter_pool_comm, root_pool 482 USE mp_world, ONLY : mpime, world_comm 483 USE io_dyn_mat,ONLY : read_dyn_mat_param, read_dyn_mat_header, & 484 read_ifc_param, read_ifc 485#if defined(__NAG) 486 USE f90_unix_io, ONLY : flush 487#endif 488 ! 489 IMPLICIT NONE 490 ! 491 ! Local variables 492 LOGICAL :: lpolar_ 493 !! Polar flag 494 LOGICAL :: has_zstar 495 !! Does it has Born effective charges 496 LOGICAL :: is_xml_file 497 !! Is the file XML 498 CHARACTER(LEN = 80) :: line 499 !! 500 CHARACTER(LEN = 256) :: tempfile 501 !! 502 CHARACTER(LEN = 3), ALLOCATABLE :: atm(:) 503 !! 504 INTEGER :: ios 505 !! 506 INTEGER :: i, j 507 !! 508 INTEGER :: m1, m2, m3 509 !! 510 INTEGER :: na, nb 511 !! 512 INTEGER :: idum 513 !! 514 INTEGER :: ibid, jbid 515 !! 516 INTEGER :: nabid 517 !! 518 INTEGER :: nbbid 519 !! 520 INTEGER :: m1bid, m2bid, m3bid 521 !! 522 INTEGER :: ntyp_ 523 !! 524 INTEGER :: nat_ 525 !! 526 INTEGER :: ibrav_ 527 !! 528 INTEGER :: ityp_(nat) 529 !! 530 INTEGER :: nqs 531 !! 532 INTEGER :: ierr 533 !! Error status 534 INTEGER, PARAMETER :: ntypx = 10 535 !! 536 REAL(KIND = DP):: tau_(3, nat) 537 !! 538 REAL(KIND = DP):: amass2(ntypx) 539 !! 540 REAL(KIND = DP), ALLOCATABLE :: m_loc(:, :) 541 !! 542 ! 543 WRITE(stdout, '(/5x,"Reading interatomic force constants"/)') 544 FLUSH(stdout) 545 ! 546 ! Generic name for the ifc.q2r file. If it is xml, the file will be named ifc.q2r.xml instead 547 tempfile = TRIM(dvscf_dir) // 'ifc.q2r' 548 ! The following function will check if the file exists in xml format 549 CALL check_is_xml_file(tempfile, is_xml_file) 550 ! 551 IF (is_xml_file) THEN 552 !IF (mpime == 0) THEN 553 ! ionode = .TRUE. 554 !ELSE 555 ! ionode = .FALSE. 556 !ENDIF 557 ! 558 ! pass the 'tempfile' as the '.xml' extension is added in the next routine 559 CALL read_dyn_mat_param(tempfile, ntyp_, nat_) 560 ALLOCATE(m_loc(3, nat_), STAT = ierr) 561 IF (ierr /= 0) CALL errore('read_ifc_epw', 'Error allocating m_loc', 1) 562 ALLOCATE(atm(ntyp_), STAT = ierr) 563 IF (ierr /= 0) CALL errore('read_ifc_epw', 'Error allocating atm', 1) 564 CALL read_dyn_mat_header(ntyp_, nat_, ibrav, nspin_mag, & 565 celldm, at, bg, omega, atm, amass2, & 566 tau_, ityp_, m_loc, nqs, has_zstar, epsi, zstar) 567 CALL volume(alat, at(1, 1), at(1, 2), at(1, 3), omega) 568 CALL read_ifc_param(nqc1, nqc2, nqc3) 569 CALL read_ifc(nqc1, nqc2, nqc3, nat_, ifc) 570 DEALLOCATE(m_loc, STAT = ierr) 571 IF (ierr /= 0) CALL errore('read_ifc_epw', 'Error deallocating m_loc', 1) 572 DEALLOCATE(atm, STAT = ierr) 573 IF (ierr /= 0) CALL errore('read_ifc_epw', 'Error deallocating atm', 1) 574 ! 575 ELSE ! is_xml_file 576 IF (mpime == ionode_id) THEN 577 ! 578 OPEN(UNIT = iunifc, FILE = tempfile, STATUS = 'old', IOSTAT = ios) 579 IF (ios /= 0) call errore ('read_ifc_epw', 'error opening ifc.q2r', iunifc) 580 ! 581 ! read real-space interatomic force constants 582 ! 583 READ(iunifc,'(3i4)') ntyp_ , nat_ , ibrav_ 584 IF (ibrav_ == 0) THEN 585 DO i = 1, 3 586 READ(iunifc, *) line 587 ENDDO 588 ENDIF 589 DO i = 1, ntyp_ 590 READ(iunifc, '(a)') line 591 ENDDO 592 DO na = 1, nat 593 READ(iunifc, *) idum, idum, (tau_(j, na), j = 1, 3) 594 ENDDO 595 READ(iunifc, *) lpolar_ 596 ! 597 IF (lpolar_) THEN 598 READ (iunifc,*) ((epsi(i, j), j = 1, 3), i = 1, 3) 599 DO na = 1, nat 600 READ(iunifc, *) idum 601 READ(iunifc, *) ((zstar(i, j, na), j = 1, 3), i = 1, 3) 602 ENDDO 603 WRITE(stdout, '(5x,a)') "Read Z* and epsilon" 604 ENDIF 605 ! 606 READ (iunifc,*) idum 607 ! 608 ifc = 0.d0 609 DO i = 1, 3 610 DO j = 1, 3 611 DO na = 1, nat 612 DO nb = 1, nat 613 READ(iunifc, *) ibid, jbid, nabid, nbbid 614 IF (i /= ibid .OR. j /= jbid .OR. na /= nabid .OR. nb /= nbbid) & 615 CALL errore('read_epw', 'error in reading ifc', 1) 616 READ(iunifc, *) (((m1bid, m2bid, m3bid, ifc(m1, m2, m3, i, j, na, nb), & 617 m1 = 1, nqc1), m2 = 1, nqc2), m3 = 1, nqc3) 618 ENDDO 619 ENDDO 620 ENDDO 621 ENDDO 622 CLOSE(iunifc) 623 ENDIF ! mpime 624 ! It has to be casted like this because mpi cannot cast 7 indices 625 DO i = 1, 3 626 DO j = 1, 3 627 DO na = 1, nat 628 DO nb = 1, nat 629 CALL mp_bcast(ifc(:, :, :, i, j, na, nb), ionode_id, inter_pool_comm) 630 CALL mp_bcast(ifc(:, :, :, i, j, na, nb), root_pool, intra_pool_comm) 631 ENDDO 632 ENDDO 633 ENDDO 634 ENDDO 635 CALL mp_bcast(zstar, ionode_id, world_comm) 636 CALL mp_bcast(epsi, ionode_id, world_comm) 637 CALL mp_bcast(tau_, ionode_id, world_comm) 638 CALL mp_bcast(ibrav_, ionode_id, world_comm) 639 ENDIF ! has_xml 640 ! 641 WRITE(stdout,'(5x,"IFC last ", 1f12.7)') ifc(nqc1, nqc2, nqc3, 3, 3, nat, nat) 642 ! 643 CALL set_asr2(asr_typ, nqc1, nqc2, nqc3, ifc, zstar, nat, ibrav_, tau_) 644 ! 645 WRITE(stdout, '(/5x,"Finished reading ifcs"/)') 646 ! 647 !------------------------------------------------------------------------------- 648 END SUBROUTINE read_ifc_epw 649 !------------------------------------------------------------------------------- 650 ! 651 !---------------------------------------------------------------------- 652 SUBROUTINE set_asr2(asr, nr1, nr2, nr3, frc, zeu, nat, ibrav, tau) 653 !----------------------------------------------------------------------- 654 !! 655 !! Set the acoustic sum rule. 656 !! Taken directly from PHonon/PH/q2trans.f90 657 !! It would be better to take the set_asr for /Modules/. 658 !! However they are different (frc) and to be consitent with q2r.x we take this one. 659 ! 660 USE kinds, ONLY : DP 661 USE io_global, ONLY : stdout 662 ! 663 IMPLICIT NONE 664 ! 665 CHARACTER(LEN = 10), INTENT(in) :: asr 666 !! Acoustic sum rule 667 INTEGER, INTENT(in) :: nr1, nr2, nr3 668 !! 669 INTEGER, INTENT(in) :: nat 670 !! Number of atoms 671 INTEGER, INTENT(in) :: ibrav 672 !! Bravais lattice 673 REAL(KIND = DP), INTENT(in) :: tau(3, nat) 674 !! Atomic position 675 REAL(KIND = DP), INTENT(inout) :: frc(nr1, nr2, nr3, 3, 3, nat, nat) 676 !! Real-space IFC 677 REAL(KIND = DP), INTENT(inout) :: zeu(3, 3, nat) 678 !! Z* 679 ! 680 ! Local variables 681 TYPE vector 682 REAL(KIND = DP), pointer :: vec(:, :, :, :, :, :, :) 683 END TYPE vector 684 ! 685 TYPE (vector) u(6 * 3 * nat) 686 ! These are the "vectors" associated with the sum rules on force-constants 687 INTEGER :: axis 688 !! 689 INTEGER :: n 690 !! 691 INTEGER :: i, j 692 !! 693 INTEGER :: na, nb 694 !! 695 INTEGER :: n1, n2, n3 696 !! 697 INTEGER :: m, p, k, l, q, r 698 !! 699 INTEGER :: i1, j1 700 !! 701 INTEGER :: na1 702 !! 703 INTEGER :: ierr 704 !! Error status 705 INTEGER :: u_less(6 * 3 * nat) 706 !! indices of the vectors u that are not independent to the preceding ones 707 INTEGER :: n_less 708 !! Number of index that are no independent 709 INTEGER :: i_less 710 !! temporary parameter 711 INTEGER :: zeu_less(6 * 3) 712 !! indices of the vectors zeu_u that are not independent to the preceding ones, 713 INTEGER :: nzeu_less 714 !! number of such vectors 715 INTEGER :: izeu_less 716 !! temporary parameter 717 INTEGER, ALLOCATABLE :: ind_v(:, :, :) 718 !! 719 REAL(KIND = DP) :: zeu_new(3, 3, nat) 720 !! 721 REAL(KIND = DP) :: scal 722 !! 723 REAL(KIND = DP) :: norm2 724 !! 725 REAL(KIND = DP) :: sum 726 !! 727 REAL(KIND = DP) :: zeu_u(6 * 3, 3, 3, nat) 728 !! These are the "vectors" associated with the sum rules on effective charges 729 REAL(KIND = DP) :: zeu_w(3, 3, nat) 730 !! 731 REAL(KIND = DP) :: zeu_x(3, 3, nat) 732 !! 733 REAL(KIND = DP), ALLOCATABLE :: w(:, :, :, :, :, :, :) 734 !! temporary vectors and parameters 735 REAL(KIND = DP), ALLOCATABLE :: x(:, :, :, :, :, :, :) 736 !! temporary vectors and parameters 737 REAL(KIND = DP), ALLOCATABLE :: frc_new(:,:,:,:,:,:,:) 738 ! 739 REAL(KIND = DP), ALLOCATABLE :: v(:, :) 740 !! These are the "vectors" associated with symmetry conditions, coded by 741 !! indicating the positions (i.e. the seven indices) of the non-zero elements (there 742 !! should be only 2 of them) and the value of that element. We do so in order 743 !! to limit the amount of memory used. 744 ! 745 ! Initialization. n is the number of sum rules to be considered (if 746 ! asr/='simple') 747 ! and 'axis' is the rotation axis in the case of a 1D system 748 ! (i.e. the rotation axis is (Ox) if axis='1', (Oy) if axis='2' and (Oz) if 749 ! axis='3') 750 ! 751 IF ((asr /= 'simple') .AND. (asr /= 'crystal') .AND. (asr /= 'one-dim') .AND. (asr /= 'zero-dim')) THEN 752 CALL errore('set_asr','invalid Acoustic Sum Rule:' // asr, 1) 753 ENDIF 754 ! 755 IF (asr == 'simple') THEN 756 ! 757 ! Simple Acoustic Sum Rule on effective charges 758 ! 759 DO i = 1, 3 760 DO j = 1, 3 761 sum = 0.0d0 762 DO na = 1, nat 763 sum = sum + zeu(i, j, na) 764 ENDDO 765 DO na = 1, nat 766 zeu(i, j, na) = zeu(i, j, na) - sum / nat 767 ENDDO 768 ENDDO 769 ENDDO 770 ! 771 ! Simple Acoustic Sum Rule on force constants in real space 772 ! 773 DO i = 1, 3 774 DO j = 1, 3 775 DO na = 1, nat 776 sum = 0.0d0 777 DO nb = 1, nat 778 DO n1 = 1, nr1 779 DO n2 = 1, nr2 780 DO n3 = 1, nr3 781 sum = sum + frc(n1, n2, n3, i, j, na, nb) 782 ENDDO 783 ENDDO 784 ENDDO 785 ENDDO 786 frc(1, 1, 1, i, j, na, na) = frc(1, 1, 1, i, j, na, na) - sum 787 ENDDO 788 ENDDO 789 ENDDO 790 WRITE (stdout,'(5x,a)') " Imposed simple ASR" 791 RETURN 792 ! 793 ENDIF ! simple 794 ! 795 IF (asr == 'crystal') n = 3 796 IF (asr == 'one-dim') THEN 797 ! the direction of periodicity is the rotation axis 798 ! It will work only if the crystal axis considered is one of 799 ! the cartesian axis (typically, ibrav = 1, 6 or 8, or 4 along the z-direction) 800 IF (nr1 * nr2 * nr3 == 1) axis = 3 801 IF ((nr1 /= 1) .AND. (nr2 * nr3 == 1)) axis = 1 802 IF ((nr2 /= 1) .AND. (nr1 * nr3 == 1)) axis = 2 803 IF ((nr3 /= 1) .AND. (nr1 * nr2 == 1)) axis = 3 804 IF (((nr1 /= 1) .AND. (nr2 /=1 )) .OR. ((nr2 /= 1) .AND. (nr3 /= 1)) .OR. ((nr1 /= 1) .AND. (nr3 /=1 ))) THEN 805 CALL errore('set_asr', 'too many directions of periodicity in 1D system', axis) 806 ENDIF 807 IF ((ibrav /= 1) .AND. (ibrav /= 6) .AND. (ibrav /= 8) .AND. ((ibrav /= 4) .OR. (axis /= 3))) THEN 808 WRITE(stdout, *) 'asr: rotational axis may be wrong' 809 ENDIF 810 WRITE(stdout, '("asr rotation axis in 1D system= ", I4)') axis 811 n = 4 812 ENDIF 813 IF(asr == 'zero-dim') n = 6 814 ! 815 ! Acoustic Sum Rule on effective charges 816 ! 817 ! generating the vectors of the orthogonal of the subspace to project the effective charges matrix on 818 ! 819 zeu_u(:, :, :, :) = 0.0d0 820 DO i = 1, 3 821 DO j = 1, 3 822 DO na = 1, nat 823 zeu_new(i, j, na) = zeu(i, j, na) 824 ENDDO 825 ENDDO 826 ENDDO 827 ! 828 p = 0 829 DO i = 1, 3 830 DO j = 1, 3 831 ! These are the 3*3 vectors associated with the 832 ! translational acoustic sum rules 833 p = p + 1 834 zeu_u(p, i, j, :) = 1.0d0 835 ! 836 ENDDO 837 ENDDO 838 ! 839 IF (n == 4) THEN 840 DO i = 1, 3 841 ! These are the 3 vectors associated with the 842 ! single rotational sum rule (1D system) 843 p = p + 1 844 DO na = 1, nat 845 zeu_u(p, i, MOD(axis, 3) + 1, na) = -tau(MOD(axis + 1, 3) + 1,na) 846 zeu_u(p, i, MOD(axis + 1, 3) + 1, na) = tau(MOD(axis, 3) + 1, na) 847 ENDDO 848 ENDDO 849 ENDIF 850 ! 851 IF (n == 6) THEN 852 DO i = 1, 3 853 DO j = 1, 3 854 ! These are the 3*3 vectors associated with the 855 ! three rotational sum rules (0D system - typ. molecule) 856 p = p + 1 857 DO na = 1, nat 858 zeu_u(p, i, MOD(j, 3) + 1, na) = -tau(MOD(j + 1, 3) + 1, na) 859 zeu_u(p, i, MOD(j + 1, 3) + 1, na) = tau(MOD(j, 3) + 1, na) 860 ENDDO 861 ENDDO 862 ENDDO 863 ENDIF 864 ! 865 ! Gram-Schmidt orthonormalization of the set of vectors created. 866 ! 867 nzeu_less = 0 868 DO k = 1, p 869 zeu_w(:, :, :) = zeu_u(k, :, :, :) 870 zeu_x(:, :, :) = zeu_u(k, :, :, :) 871 DO q = 1, k - 1 872 r = 1 873 DO izeu_less = 1, nzeu_less 874 IF (zeu_less(izeu_less) == q) r = 0 875 ENDDO 876 IF (r /= 0) THEN 877 CALL sp_zeu(zeu_x, zeu_u(q, :, :, :), nat,scal) 878 zeu_w(:, :, :) = zeu_w(:, :, :) - scal * zeu_u(q, :, :, :) 879 ENDIF 880 ENDDO 881 CALL sp_zeu(zeu_w, zeu_w, nat, norm2) 882 IF (norm2 > 1.0d-16) THEN 883 zeu_u(k,:,:,:) = zeu_w(:, :, :) / DSQRT(norm2) 884 ELSE 885 nzeu_less = nzeu_less + 1 886 zeu_less(nzeu_less) = k 887 ENDIF 888 ENDDO 889 ! 890 ! Projection of the effective charge "vector" on the orthogonal of the 891 ! subspace of the vectors verifying the sum rules 892 ! 893 zeu_w(:, :, :) = 0.0d0 894 DO k = 1, p 895 r = 1 896 DO izeu_less = 1,nzeu_less 897 IF (zeu_less(izeu_less) == k) r = 0 898 ENDDO 899 IF (r /= 0) THEN 900 zeu_x(:, :, :) = zeu_u(k, :, :, :) 901 CALL sp_zeu(zeu_x, zeu_new, nat, scal) 902 zeu_w(:, :, :) = zeu_w(:, :, :) + scal * zeu_u(k, :, :, :) 903 ENDIF 904 ENDDO 905 ! 906 ! Final substraction of the former projection to the initial zeu, to get 907 ! the new "projected" zeu 908 ! 909 zeu_new(:, :, :) = zeu_new(:, :, :) - zeu_w(:, :, :) 910 CALL sp_zeu(zeu_w, zeu_w, nat, norm2) 911 WRITE(stdout, '(5x,"Norm of the difference between old and new effective charges: ", 1f12.7)') DSQRT(norm2) 912 ! 913 DO i = 1, 3 914 DO j = 1, 3 915 DO na = 1, nat 916 zeu(i, j, na) = zeu_new(i, j, na) 917 ENDDO 918 ENDDO 919 ENDDO 920 ! 921 ! Acoustic Sum Rule on force constants 922 ! 923 ! generating the vectors of the orthogonal of the subspace to project 924 ! the force-constants matrix on 925 ! 926 DO k = 1, 18 * nat 927 ALLOCATE(u(k) % vec(nr1, nr2, nr3, 3, 3, nat, nat), STAT = ierr) 928 IF (ierr /= 0) CALL errore('set_asr2', 'Error allocating u(k) % vec', 1) 929 u(k) % vec (:, :, :, :, :, :, :) = 0.0d0 930 ENDDO 931 ALLOCATE(frc_new(nr1, nr2, nr3, 3, 3, nat, nat), STAT = ierr) 932 IF (ierr /= 0) CALL errore('set_asr2', 'Error allocating frc_new', 1) 933 DO i = 1, 3 934 DO j = 1, 3 935 DO na = 1, nat 936 DO nb = 1, nat 937 DO n1 = 1, nr1 938 DO n2 = 1, nr2 939 DO n3 = 1, nr3 940 frc_new(n1, n2, n3, i, j, na, nb) = frc(n1, n2, n3, i, j, na, nb) 941 ENDDO 942 ENDDO 943 ENDDO 944 ENDDO 945 ENDDO 946 ENDDO 947 ENDDO 948 ! 949 p = 0 950 DO i = 1, 3 951 DO j = 1, 3 952 DO na = 1, nat 953 ! These are the 3*3*nat vectors associated with the 954 ! translational acoustic sum rules 955 p = p + 1 956 u(p) % vec(:, :, :, i, j, na, :) = 1.0d0 957 ! 958 ENDDO 959 ENDDO 960 ENDDO 961 ! 962 IF (n == 4) THEN 963 DO i = 1, 3 964 DO na = 1, nat 965 ! These are the 3*nat vectors associated with the 966 ! single rotational sum rule (1D system) 967 p = p + 1 968 DO nb = 1, nat 969 u(p) % vec(: ,:, :, i, MOD(axis, 3) + 1, na, nb) = -tau(MOD(axis + 1, 3) + 1, nb) 970 u(p) % vec(:, :, :, i, MOD(axis + 1, 3) + 1, na, nb) = tau(MOD(axis, 3) + 1, nb) 971 ENDDO 972 ENDDO 973 ENDDO 974 ENDIF 975 ! 976 IF (n == 6) THEN 977 DO i = 1, 3 978 DO j = 1, 3 979 DO na = 1, nat 980 ! These are the 3*3*nat vectors associated with the 981 ! three rotational sum rules (0D system - typ. molecule) 982 p = p + 1 983 DO nb = 1, nat 984 u(p) % vec(:, :, :, i, MOD(j, 3) + 1, na, nb) = -tau(MOD(j + 1, 3) + 1, nb) 985 u(p) % vec(:, :, :, i, MOD(j + 1, 3) + 1, na, nb) = tau(MOD(j, 3) + 1, nb) 986 ENDDO 987 ENDDO 988 ENDDO 989 ENDDO 990 ENDIF 991 ! 992 ALLOCATE(ind_v(9 * nat * nat * nr1 * nr2 * nr3, 2, 7), STAT = ierr) 993 IF (ierr /= 0) CALL errore('set_asr2', 'Error allocating ind_v', 1) 994 ALLOCATE(v(9 * nat * nat * nr1 * nr2 * nr3, 2), STAT = ierr) 995 IF (ierr /= 0) CALL errore('set_asr2', 'Error allocating v', 1) 996 m = 0 997 DO i = 1, 3 998 DO j = 1, 3 999 DO na = 1, nat 1000 DO nb = 1, nat 1001 DO n1 = 1, nr1 1002 DO n2 = 1, nr2 1003 DO n3 = 1, nr3 1004 ! These are the vectors associated with the symmetry 1005 ! constraints 1006 q = 1 1007 l = 1 1008 DO WHILE((l <= m) .AND. (q /= 0)) 1009 IF ((ind_v(l, 1, 1) == n1) .AND. (ind_v(l, 1, 2) == n2) .AND. & 1010 (ind_v(l, 1, 3) == n3) .AND. (ind_v(l, 1, 4) == i) .AND. & 1011 (ind_v(l, 1, 5) == j) .AND. (ind_v(l, 1, 6) == na) .AND. & 1012 (ind_v(l, 1, 7) == nb)) q = 0 1013 IF ((ind_v(l, 2, 1) == n1) .AND. (ind_v(l, 2, 2) == n2) .AND. & 1014 (ind_v(l, 2, 3) == n3) .AND. (ind_v(l, 2, 4) == i) .AND. & 1015 (ind_v(l, 2, 5) == j) .AND. (ind_v(l, 2, 6) == na) .AND. & 1016 (ind_v(l, 2, 7) == nb)) q = 0 1017 l = l + 1 1018 ENDDO 1019 IF ((n1 == MOD(nr1 + 1 - n1, nr1) + 1) .AND. (n2 == MOD(nr2 + 1 - n2, nr2) + 1) & 1020 .AND. (n3 == MOD(nr3 + 1 - n3, nr3) + 1) .AND. (i == j) .AND. (na == nb)) q = 0 1021 IF (q /= 0) THEN 1022 m = m + 1 1023 ind_v(m, 1, 1) = n1 1024 ind_v(m, 1, 2) = n2 1025 ind_v(m, 1, 3) = n3 1026 ind_v(m, 1, 4) = i 1027 ind_v(m, 1, 5) = j 1028 ind_v(m, 1, 6) = na 1029 ind_v(m, 1, 7) = nb 1030 v(m, 1) = 1.0d0 / DSQRT(2.0d0) 1031 ind_v(m, 2, 1) = MOD(nr1 + 1 - n1, nr1) + 1 1032 ind_v(m, 2, 2) = MOD(nr2 + 1 - n2, nr2) + 1 1033 ind_v(m, 2, 3) = MOD(nr3 + 1 - n3, nr3) + 1 1034 ind_v(m, 2, 4) = j 1035 ind_v(m, 2, 5) = i 1036 ind_v(m, 2, 6) = nb 1037 ind_v(m, 2, 7) = na 1038 v(m, 2) = -1.0d0 / DSQRT(2.0d0) 1039 ENDIF 1040 ENDDO 1041 ENDDO 1042 ENDDO 1043 ENDDO 1044 ENDDO 1045 ENDDO 1046 ENDDO 1047 ! 1048 ! Gram-Schmidt orthonormalization of the set of vectors created. 1049 ! Note that the vectors corresponding to symmetry constraints are already 1050 ! orthonormalized by construction. 1051 ! 1052 n_less = 0 1053 ALLOCATE(w(nr1, nr2, nr3, 3, 3, nat, nat), STAT = ierr) 1054 IF (ierr /= 0) CALL errore('set_asr2', 'Error allocating w', 1) 1055 ALLOCATE(x(nr1, nr2, nr3, 3, 3, nat, nat), STAT = ierr) 1056 IF (ierr /= 0) CALL errore('set_asr2', 'Error allocating x', 1) 1057 DO k = 1, p 1058 w(:, :, :, :, :, :, :) = u(k) % vec(:, :, :, :, :, :, :) 1059 x(:, :, :, :, :, :, :) = u(k) % vec(:, :, :, :, :, :, :) 1060 DO l = 1, m 1061 ! 1062 CALL sp2(x, v(l, :), ind_v(l, :, :), nr1, nr2, nr3, nat, scal) 1063 DO r = 1, 2 1064 n1 = ind_v(l, r, 1) 1065 n2 = ind_v(l, r, 2) 1066 n3 = ind_v(l, r, 3) 1067 i = ind_v(l, r, 4) 1068 j = ind_v(l, r, 5) 1069 na = ind_v(l, r, 6) 1070 nb = ind_v(l, r, 7) 1071 w(n1, n2, n3, i, j, na, nb) = w(n1, n2, n3, i, j, na, nb) - scal * v(l, r) 1072 ENDDO 1073 ENDDO 1074 IF (k <= (9 * nat)) THEN 1075 na1 = MOD(k, nat) 1076 IF (na1 == 0) na1 = nat 1077 j1 = MOD((k - na1) / nat, 3) + 1 1078 i1 = MOD((((k - na1) / nat) - j1 + 1) / 3, 3) + 1 1079 ELSE 1080 q = k - 9 * nat 1081 IF (n == 4) THEN 1082 na1 = MOD(q, nat) 1083 IF (na1 == 0) na1 = nat 1084 i1 = MOD((q - na1) / nat, 3) + 1 1085 ELSE 1086 na1 = MOD(q, nat) 1087 IF (na1 == 0) na1 = nat 1088 j1 = MOD((q - na1) / nat, 3) + 1 1089 i1 = MOD((((q - na1) / nat) - j1 + 1) / 3, 3) + 1 1090 ENDIF 1091 ENDIF 1092 DO q = 1, k - 1 1093 r = 1 1094 DO i_less = 1, n_less 1095 IF (u_less(i_less) == q) r = 0 1096 ENDDO 1097 IF (r /= 0) THEN 1098 CALL sp3(x, u(q) % vec(:, :, :, :, :, :, :), i1, na1, nr1, nr2, nr3, nat, scal) 1099 w(:, :, :, :, :, :, :) = w(:, :, :, :, :, :, :) - scal * u(q) % vec(:, :, :, :, :, :, :) 1100 ENDIF 1101 ENDDO 1102 CALL sp1(w, w, nr1, nr2, nr3, nat, norm2) 1103 IF (norm2 > 1.0d-16) THEN 1104 u(k) % vec(:, :, :, :, :, :, :) = w(:, :, :, :, :, :, :) / DSQRT(norm2) 1105 ELSE 1106 n_less = n_less + 1 1107 u_less(n_less) = k 1108 ENDIF 1109 ENDDO 1110 ! 1111 ! Projection of the force-constants "vector" on the orthogonal of the 1112 ! subspace of the vectors verifying the sum rules and symmetry contraints 1113 ! 1114 w(:, :, :, :, :, :, :) = 0.0d0 1115 DO l = 1, m 1116 CALL sp2(frc_new, v(l, :), ind_v(l, :, :), nr1, nr2, nr3, nat, scal) 1117 DO r = 1, 2 1118 n1 = ind_v(l, r, 1) 1119 n2 = ind_v(l, r, 2) 1120 n3 = ind_v(l, r, 3) 1121 i = ind_v(l, r, 4) 1122 j = ind_v(l, r, 5) 1123 na = ind_v(l, r, 6) 1124 nb = ind_v(l, r, 7) 1125 w(n1, n2, n3, i, j, na, nb) = w(n1, n2, n3, i, j, na, nb) + scal * v(l, r) 1126 ENDDO 1127 ENDDO 1128 DO k = 1, p 1129 r = 1 1130 DO i_less = 1,n_less 1131 IF (u_less(i_less) == k) r = 0 1132 ENDDO 1133 IF (r /= 0) THEN 1134 x(:, :, :, :, :, :, :) = u(k) % vec(:, :, :, :, :, :, :) 1135 call sp1(x, frc_new, nr1, nr2, nr3, nat, scal) 1136 w(:, :, :, :, :, :, :) = w(:, :, :, :, :, :, :) + scal * u(k)%vec(:, :, :, :, :, :, :) 1137 ENDIF 1138 DEALLOCATE(u(k)%vec, STAT = ierr) 1139 IF (ierr /= 0) CALL errore('set_asr2', 'Error deallocating u(k)%vec', 1) 1140 ENDDO 1141 ! 1142 ! Final substraction of the former projection to the initial frc, to get 1143 ! the new "projected" frc 1144 ! 1145 frc_new(:, :, :, :, :, :, :) = frc_new(:, :, :, :, :, :, :) - w(:, :, :, :, :, :, :) 1146 CALL sp1(w, w, nr1, nr2, nr3, nat, norm2) 1147 ! 1148 WRITE(stdout, '(5x,"Norm of the difference between old and new force-constants: ", 1f12.7)') DSQRT(norm2) 1149 ! 1150 DO i = 1, 3 1151 DO j = 1, 3 1152 DO na = 1, nat 1153 DO nb = 1, nat 1154 DO n1 = 1, nr1 1155 DO n2 = 1, nr2 1156 DO n3 = 1, nr3 1157 frc(n1, n2, n3, i, j, na, nb) = frc_new(n1, n2, n3, i, j, na, nb) 1158 ENDDO 1159 ENDDO 1160 ENDDO 1161 ENDDO 1162 ENDDO 1163 ENDDO 1164 ENDDO 1165 DEALLOCATE(x, STAT = ierr) 1166 IF (ierr /= 0) CALL errore('set_asr2', 'Error deallocating x', 1) 1167 DEALLOCATE(w, STAT = ierr) 1168 IF (ierr /= 0) CALL errore('set_asr2', 'Error deallocating w', 1) 1169 DEALLOCATE(v, STAT = ierr) 1170 IF (ierr /= 0) CALL errore('set_asr2', 'Error deallocating v', 1) 1171 DEALLOCATE(ind_v, STAT = ierr) 1172 IF (ierr /= 0) CALL errore('set_asr2', 'Error deallocating ind_v', 1) 1173 DEALLOCATE(frc_new, STAT = ierr) 1174 IF (ierr /= 0) CALL errore('set_asr2', 'Error deallocating frc_new', 1) 1175 WRITE(stdout, '(5x,a)') "Imposed crystal ASR" 1176 ! 1177 RETURN 1178 !---------------------------------------------------------------------- 1179 END SUBROUTINE set_asr2 1180 !---------------------------------------------------------------------- 1181 ! 1182 !---------------------------------------------------------------------- 1183 SUBROUTINE sp_zeu(zeu_u, zeu_v, nat, scal) 1184 !----------------------------------------------------------------------- 1185 !! 1186 !! Does the scalar product of two effective charges matrices zeu_u and zeu_v 1187 !! (considered as vectors in the R^(3*3*nat) space, and coded in the usual way) 1188 !! 1189 USE kinds, ONLY: DP 1190 ! 1191 IMPLICIT NONE 1192 ! 1193 INTEGER, INTENT(in) :: nat 1194 !! 1195 REAL(KIND = DP), INTENT(in) :: zeu_u(3, 3, nat) 1196 !! 1197 REAL(KIND = DP), INTENT(in) :: zeu_v(3, 3, nat) 1198 !! 1199 REAL(KIND = DP), INTENT(inout) :: scal 1200 !! 1201 ! Local variables 1202 INTEGER :: i 1203 !! 1204 INTEGER :: j 1205 !! 1206 INTEGER :: na 1207 !! 1208 ! 1209 scal = 0.0d0 1210 DO i = 1, 3 1211 DO j = 1, 3 1212 DO na = 1, nat 1213 scal = scal + zeu_u(i, j, na) * zeu_v(i, j, na) 1214 ENDDO 1215 ENDDO 1216 ENDDO 1217 ! 1218 RETURN 1219 ! 1220 !----------------------------------------------------------------------- 1221 END SUBROUTINE sp_zeu 1222 !----------------------------------------------------------------------- 1223 ! 1224 !---------------------------------------------------------------------- 1225 SUBROUTINE sp1(u, v, nr1, nr2, nr3, nat, scal) 1226 !----------------------------------------------------------------------- 1227 !! 1228 !! Does the scalar product of two force-constants matrices u and v (considered as 1229 !! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space, and coded in the usual way) 1230 !! 1231 USE kinds, ONLY: DP 1232 ! 1233 IMPLICIT NONE 1234 ! 1235 INTEGER, INTENT(in) :: nr1, nr2, nr3 1236 !! Supercell dims 1237 INTEGER, INTENT(in) :: nat 1238 !! Number of atoms 1239 REAL(KIND = DP), INTENT(in) :: u(nr1, nr2, nr3, 3, 3, nat, nat) 1240 !! First force-constent matrix 1241 REAL(KIND = DP), INTENT(in) :: v(nr1, nr2, nr3, 3, 3, nat, nat) 1242 !! Second force-constent matrix 1243 REAL(KIND = DP), INTENT(out) :: scal 1244 !! Scalar product 1245 ! 1246 ! Local variables 1247 INTEGER :: i, j 1248 !! Cartesian direction 1249 INTEGER :: na, nb 1250 !! Atoms index 1251 INTEGER :: n1, n2, n3 1252 !! Supercell index 1253 ! 1254 scal = 0.0d0 1255 DO i = 1, 3 1256 DO j = 1, 3 1257 DO na = 1, nat 1258 DO nb = 1, nat 1259 DO n1 = 1, nr1 1260 DO n2 = 1, nr2 1261 DO n3 = 1, nr3 1262 scal = scal + u(n1, n2, n3, i, j, na, nb) * v(n1, n2, n3, i, j, na, nb) 1263 ENDDO 1264 ENDDO 1265 ENDDO 1266 ENDDO 1267 ENDDO 1268 ENDDO 1269 ENDDO 1270 ! 1271 RETURN 1272 ! 1273 !---------------------------------------------------------------------- 1274 END SUBROUTINE sp1 1275 !---------------------------------------------------------------------- 1276 ! 1277 !---------------------------------------------------------------------- 1278 SUBROUTINE sp2(u, v, ind_v, nr1, nr2, nr3, nat, scal) 1279 !----------------------------------------------------------------------- 1280 !! 1281 !! Does the scalar product of two force-constants matrices u and v (considered as 1282 !! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space). u is coded in the usual way 1283 !! but v is coded as explained when defining the vectors corresponding to the symmetry constraints 1284 !! 1285 USE kinds, ONLY: DP 1286 ! 1287 IMPLICIT NONE 1288 ! 1289 INTEGER, INTENT(in) :: nr1, nr2, nr3 1290 !! Supercell dims 1291 INTEGER, INTENT(in) :: nat 1292 !! Number of atoms 1293 INTEGER, INTENT(in) :: ind_v(2, 7) 1294 !! Index vector 1295 REAL(KIND = DP), INTENT(in) :: u(nr1, nr2, nr3, 3, 3, nat, nat) 1296 !! Input vector 1297 REAL(KIND = DP), INTENT(in) :: v(2) 1298 !! input vector 1299 REAL(KIND = DP), INTENT(out) :: scal 1300 !! Output vector 1301 ! 1302 ! Local variables 1303 INTEGER :: i 1304 !! Index 1305 ! 1306 scal = 0.0d0 1307 DO i = 1, 2 1308 scal = scal + u(ind_v(i, 1), ind_v(i, 2), ind_v(i, 3), ind_v(i, 4), ind_v(i, 5), ind_v(i, 6), & 1309 ind_v(i, 7)) * v(i) 1310 ENDDO 1311 ! 1312 RETURN 1313 ! 1314 !----------------------------------------------------------------------- 1315 END SUBROUTINE sp2 1316 !----------------------------------------------------------------------- 1317 ! 1318 !---------------------------------------------------------------------- 1319 SUBROUTINE sp3(u, v, i, na, nr1, nr2, nr3, nat, scal) 1320 !----------------------------------------------------------------------- 1321 ! 1322 ! like sp1, but in the particular case when u is one of the u(k)%vec 1323 ! defined in set_asr (before orthonormalization). In this case most of the 1324 ! terms are zero (the ones that are not are characterized by i and na), so 1325 ! that a lot of computer time can be saved (during Gram-Schmidt). 1326 ! 1327 USE kinds, ONLY: DP 1328 ! 1329 IMPLICIT NONE 1330 ! 1331 INTEGER, INTENT(in) :: nr1, nr2, nr3 1332 !! Supercell dims 1333 INTEGER, INTENT(in) :: nat 1334 !! Number of atoms 1335 REAL(KIND = DP), INTENT(in) :: u(nr1, nr2, nr3, 3, 3, nat, nat) 1336 !! First force-constent matrix 1337 REAL(KIND = DP), INTENT(in) :: v(nr1, nr2, nr3, 3, 3, nat, nat) 1338 !! Second force-constent matrix 1339 REAL(KIND = DP), INTENT(out) :: scal 1340 !! Scalar product 1341 ! 1342 ! Local variables 1343 INTEGER :: i, j 1344 !! Cartesian direction 1345 INTEGER :: na, nb 1346 !! Atoms index 1347 INTEGER :: n1, n2, n3 1348 !! Supercell index 1349 ! 1350 scal = 0.0d0 1351 DO j = 1,3 1352 DO nb = 1, nat 1353 DO n1 = 1, nr1 1354 DO n2 = 1, nr2 1355 DO n3 = 1, nr3 1356 scal = scal + u(n1, n2, n3, i, j, na, nb) * v(n1, n2, n3, i, j, na, nb) 1357 ENDDO 1358 ENDDO 1359 ENDDO 1360 ENDDO 1361 ENDDO 1362 ! 1363 RETURN 1364 ! 1365 !------------------------------------------------------------------------------- 1366 END SUBROUTINE sp3 1367 !------------------------------------------------------------------------------- 1368 ! 1369 !------------------------------------------------------------------------------- 1370 SUBROUTINE check_is_xml_file(filename, is_xml_file) 1371 !------------------------------------------------------------------------------- 1372 !! 1373 !! This routine checks if a file is formatted in XML. It does so by 1374 !! checking if the file exists and if the file + '.xml' in its name exists. 1375 !! If both of them or none of them exists, an error is raised. If only one of 1376 !! them exists, it sets the 'is_xml_file' to .TRUE. of .FALSE. depending of 1377 !! the file found. 1378 !! 1379 IMPLICIT NONE 1380 ! 1381 CHARACTER(LEN = 256), INTENT(in) :: filename 1382 !! The name of the file to check if formatted in XML format 1383 !! This string is assumed to be trimmed 1384 LOGICAL, INTENT(out) :: is_xml_file 1385 !! Is .TRUE. if the file is in xml format. .FALSE. otherwise. 1386 ! 1387 ! Local variables 1388 CHARACTER(LEN = 256) :: filename_xml 1389 !! File name 1390 CHARACTER(LEN = 1024) :: errmsg 1391 !! Error message 1392 LOGICAL :: is_plain_text_file 1393 !! Plain tex 1394 ! 1395 filename_xml = TRIM(filename) // '.xml' 1396 filename_xml = TRIM(filename_xml) 1397 INQUIRE(FILE = filename, EXIST = is_plain_text_file) 1398 INQUIRE(FILE = filename_xml, EXIST = is_xml_file) 1399 ! Tell user if any inconsistencies 1400 IF (is_xml_file .AND. is_plain_text_file) THEN 1401 ! 2 different type of files exist => warn user 1402 errmsg = "Detected both: '" // filename // "' and '" // filename_xml // & 1403 &"' which one to choose?" 1404 CALL errore('check_is_xml_file', errmsg, 1) 1405 ELSEIF (.NOT. is_xml_file .AND. .NOT. is_plain_text_file) THEN 1406 errmsg = "Expected a file named either '" // filename //"' or '"& 1407 &// filename_xml // "' but none was found." 1408 CALL errore('check_is_xml_file', errmsg, 1) 1409 ENDIF 1410 ! else one of the file in exists 1411 !------------------------------------------------------------------------------ 1412 END SUBROUTINE check_is_xml_file 1413 !------------------------------------------------------------------------------ 1414 ! 1415 !------------------------------------------------------------- 1416 SUBROUTINE readdvscf(dvscf, recn, iq, nqc) 1417 !------------------------------------------------------------- 1418 !! 1419 !! Open dvscf files as direct access, read, and close again 1420 !! 1421 !! SP - Nov 2017 1422 !! Replaced fstat by Fortran instric function inquire. 1423 !! 1424 !! RM - Nov/Dec 2014 1425 !! Imported the noncolinear case implemented by xlzhang 1426 !! 1427 !------------------------------------------------------------- 1428 USE kinds, ONLY : DP 1429 USE io_files, ONLY : prefix 1430 USE units_ph, ONLY : lrdrho 1431 USE fft_base, ONLY : dfftp 1432 USE epwcom, ONLY : dvscf_dir 1433 USE io_var, ONLY : iudvscf 1434 USE low_lvl, ONLY : set_ndnmbr 1435 USE noncollin_module, ONLY : nspin_mag 1436 ! 1437 IMPLICIT NONE 1438 ! 1439 INTEGER, INTENT(in) :: recn 1440 !! perturbation number 1441 INTEGER, INTENT(in) :: iq 1442 !! the current q-point 1443 INTEGER, INTENT(in) :: nqc 1444 !! the total number of q-points in the list 1445 COMPLEX(KIND = DP), INTENT(inout) :: dvscf(dfftp%nnr, nspin_mag) 1446 !! dVscf potential is read from file 1447 ! 1448 ! Local variables 1449 ! 1450 CHARACTER(LEN = 256) :: tempfile 1451 !! Temp file 1452 CHARACTER(LEN = 4) :: filelab 1453 !! File number 1454 INTEGER :: unf_recl 1455 !! Rcl unit 1456 INTEGER :: ios 1457 !! Error number 1458 INTEGER(KIND = 8) :: mult_unit 1459 !! Record length 1460 INTEGER(KIND = 8) :: file_size 1461 !! File size 1462 REAL(KIND = DP) :: dummy 1463 !! Dummy variable 1464 ! 1465 ! the call to set_ndnmbr is just a trick to get quickly 1466 ! a file label by exploiting an existing subroutine 1467 ! (if you look at the sub you will find that the original 1468 ! purpose was for pools and nodes) 1469 ! 1470 CALL set_ndnmbr(0, iq, 1, nqc, filelab) 1471 tempfile = TRIM(dvscf_dir) // TRIM(prefix) // '.dvscf_q' // filelab 1472 INQUIRE(IOLENGTH = unf_recl) dummy 1473 unf_recl = unf_recl * lrdrho 1474 mult_unit = unf_recl 1475 mult_unit = recn * mult_unit 1476 ! 1477 ! open the dvscf file, read and close 1478 ! 1479 OPEN(iudvscf, FILE = tempfile, FORM = 'unformatted', & 1480 ACCESS = 'direct', IOSTAT = ios, RECL = unf_recl, STATUS = 'old') 1481 IF (ios /= 0) CALL errore('readdvscf', 'error opening ' // tempfile, iudvscf) 1482 ! 1483 ! check that the binary file is long enough 1484 INQUIRE(FILE = tempfile, SIZE = file_size) 1485 IF (mult_unit > file_size) CALL errore('readdvscf', & 1486 TRIM(tempfile) //' too short, check ecut', iudvscf) 1487 ! 1488 READ(iudvscf, REC = recn) dvscf 1489 CLOSE(iudvscf, STATUS = 'keep') 1490 ! 1491 RETURN 1492 ! 1493 !------------------------------------------------------------- 1494 END SUBROUTINE readdvscf 1495 !------------------------------------------------------------- 1496 ! 1497 !------------------------------------------------------------- 1498 SUBROUTINE readint3paw(int3paw, recn, iq, nqc) 1499 !------------------------------------------------------------- 1500 !! 1501 !! Open int3paw files as direct access, read, and close again 1502 !! 1503 !! HL - Mar 2020 based on the subroutine of readdvscf 1504 !! 1505 !------------------------------------------------------------- 1506 USE kinds, ONLY : DP 1507 USE io_files, ONLY : prefix 1508 USE units_ph, ONLY : lint3paw 1509 USE epwcom, ONLY : dvscf_dir 1510 USE io_var, ONLY : iuint3paw 1511 USE low_lvl, ONLY : set_ndnmbr 1512 USE uspp_param, ONLY : nhm 1513 USE ions_base, ONLY : nat 1514 USE noncollin_module, ONLY : nspin_mag 1515 ! 1516 IMPLICIT NONE 1517 ! 1518 INTEGER, INTENT(in) :: recn 1519 !! perturbation number 1520 INTEGER, INTENT(in) :: iq 1521 !! the current q-point 1522 INTEGER, INTENT(in) :: nqc 1523 !! the total number of q-points in the list 1524 COMPLEX(KIND = DP), INTENT(inout) :: int3paw(nhm, nhm, nat, nspin_mag) 1525 !! int3paw is read from file 1526 ! 1527 ! Local variables 1528 ! 1529 CHARACTER(LEN = 256) :: tempfile 1530 !! Temp file 1531 CHARACTER(LEN = 4) :: filelab 1532 !! File number 1533 INTEGER :: unf_recl 1534 !! Rcl unit 1535 INTEGER :: ios 1536 !! Error number 1537 INTEGER(KIND = 8) :: mult_unit 1538 !! Record length 1539 INTEGER(KIND = 8) :: file_size 1540 !! File size 1541 REAL(KIND = DP) :: dummy 1542 !! Dummy variable 1543 ! 1544 ! the call to set_ndnmbr is just a trick to get quickly 1545 ! a file label by exploiting an existing subroutine 1546 ! (if you look at the sub you will find that the original 1547 ! purpose was for pools and nodes) 1548 ! 1549 CALL set_ndnmbr(0, iq, 1, nqc, filelab) 1550 tempfile = TRIM(dvscf_dir) // TRIM(prefix) // '.dvscf_paw_q' // filelab 1551 INQUIRE(IOLENGTH = unf_recl) dummy 1552 unf_recl = unf_recl * lint3paw 1553 mult_unit = unf_recl 1554 mult_unit = recn * mult_unit 1555 ! 1556 ! open the dvscf_paw file, read and close 1557 ! 1558 OPEN(iuint3paw, FILE = tempfile, FORM = 'unformatted', & 1559 ACCESS = 'direct', IOSTAT = ios, RECL = unf_recl, STATUS = 'old') 1560 IF (ios /= 0) CALL errore('readint3paw', 'error opening ' // tempfile, iuint3paw) 1561 ! 1562 ! check that the binary file is long enough 1563 INQUIRE(FILE = tempfile, SIZE = file_size) 1564 IF (mult_unit > file_size) CALL errore('readint3paw', & 1565 TRIM(tempfile) //' too short', iuint3paw) 1566 ! 1567 READ(iuint3paw, REC = recn) int3paw 1568 CLOSE(iuint3paw, STATUS = 'keep') 1569 ! 1570 RETURN 1571 ! 1572 !------------------------------------------------------------- 1573 END SUBROUTINE readint3paw 1574 !------------------------------------------------------------- 1575 ! 1576 !------------------------------------------------------------ 1577 SUBROUTINE readwfc(ipool, recn, evc0) 1578 !------------------------------------------------------------ 1579 !! 1580 !! Open wfc files as direct access, read, and close again 1581 !! 1582 !! RM - Nov/Dec 2014 1583 !! Imported the noncolinear case implemented by xlzhang 1584 !! 1585 ! 1586 USE kinds, ONLY : DP 1587 USE io_files, ONLY : prefix, tmp_dir 1588 USE units_lr, ONLY : lrwfc, iuwfc 1589 USE wvfct, ONLY : npwx 1590 USE pwcom, ONLY : nbnd 1591 USE low_lvl, ONLY : set_ndnmbr 1592 USE noncollin_module, ONLY : npol 1593 USE mp_global, ONLY : nproc_pool, me_pool, npool 1594 ! 1595 IMPLICIT NONE 1596 ! 1597 INTEGER, INTENT(in) :: recn 1598 !! kpoint number 1599 INTEGER, INTENT(in) :: ipool 1600 !! poolfile number to be read (not used in serial case) 1601 COMPLEX(KIND = DP), INTENT(out) :: evc0(npwx * npol, nbnd) 1602 !! wavefunction is read from file 1603 ! 1604 ! Local variables 1605 CHARACTER(LEN = 256) :: tempfile 1606 !! Temp file 1607 CHARACTER(LEN = 4) :: nd_nmbr0 1608 !! File number 1609 INTEGER :: unf_recl 1610 !! Rcl unit 1611 INTEGER :: ios 1612 !! Error number 1613 REAL(KIND = DP) :: dummy 1614 !! Dummy variable 1615 ! 1616 ! Open the wfc file, read and close 1617 CALL set_ndnmbr(ipool, me_pool, nproc_pool, npool, nd_nmbr0) 1618 !print*,'nd_nmbr0 ',nd_nmbr0 1619 ! 1620#if defined(__MPI) 1621 tempfile = TRIM(tmp_dir) // TRIM(prefix) // '.wfc' // nd_nmbr0 1622# else 1623 tempfile = TRIM(tmp_dir) // TRIM(prefix) // '.wfc' 1624#endif 1625 INQUIRE(IOLENGTH = unf_recl) dummy 1626 unf_recl = unf_recl * lrwfc 1627 ! 1628 OPEN(iuwfc, FILE = tempfile, FORM = 'unformatted', ACCESS = 'direct', IOSTAT = ios, RECL = unf_recl) 1629 IF (ios /= 0) CALL errore('readwfc', 'error opening wfc file', iuwfc) 1630 READ(iuwfc, REC = recn) evc0 1631 CLOSE(iuwfc, STATUS = 'keep') 1632 ! 1633 RETURN 1634 ! 1635 !------------------------------------------------------------ 1636 END SUBROUTINE readwfc 1637 !------------------------------------------------------------ 1638 ! 1639 !-------------------------------------------------------------- 1640 SUBROUTINE readgmap(nkstot) 1641 !-------------------------------------------------------------- 1642 !! 1643 !! read the map of G vectors G -> G-G_0 for a given q point 1644 !! (this is used for the folding of k+q into the first BZ) 1645 !! 1646 ! 1647 USE kinds, ONLY : DP 1648 USE mp_global,ONLY : world_comm 1649 USE mp, ONLY : mp_bcast, mp_max 1650 USE io_global,ONLY : meta_ionode, meta_ionode_id 1651 USE io_var, ONLY : iukgmap 1652 USE elph2, ONLY : ngxxf, ngxx, ng0vec, shift, gmap, g0vec_all_r 1653 USE io_files, ONLY : prefix 1654 ! 1655 IMPLICIT NONE 1656 ! 1657 INTEGER, INTENT(in) :: nkstot 1658 !! Total number of k-points 1659 ! 1660 ! Lork variables 1661 INTEGER :: ik 1662 !! Counter on k-points 1663 INTEGER :: ik1 1664 !! Temporary indices when reading kgmap files 1665 INTEGER :: ig0 1666 !! Counter on G_0 vectors 1667 INTEGER :: ishift 1668 !! Counter on G_0 vectors 1669 INTEGER :: ig 1670 !! Counter on G vectors 1671 INTEGER :: ios 1672 !! Integer variable for I/O control 1673 INTEGER :: ierr 1674 !! Error status 1675 ! 1676 IF (meta_ionode) THEN 1677 ! 1678 OPEN(iukgmap, FILE = TRIM(prefix)//'.kgmap', FORM = 'formatted', STATUS = 'old', IOSTAT = ios) 1679 IF (ios /=0) CALL errore('readgmap', 'error opening kgmap file', iukgmap) 1680 ! 1681 READ(iukgmap, *) ngxxf 1682 ! 1683 !! HL: The part below should be removed later since it is useless. 1684 ! 1685 DO ik = 1, nkstot 1686 READ(iukgmap, *) ik1, shift(ik1) 1687 ENDDO 1688 ! 1689 READ(iukgmap, *) ng0vec 1690 ! 1691 ENDIF 1692 ! 1693 ! first node broadcasts ng0vec to all nodes for allocation of gmap 1694 ! 1695 CALL mp_bcast(ngxxf, meta_ionode_id, world_comm) 1696 CALL mp_bcast(ng0vec, meta_ionode_id, world_comm) 1697 ! 1698 ALLOCATE(gmap(ngxx * ng0vec), STAT = ierr) 1699 IF (ierr /= 0) CALL errore('readgmap', 'Error allocating gmap', 1) 1700 gmap(:) = 0 1701 ! 1702 IF (meta_ionode) THEN 1703 ! 1704 DO ig0 = 1, ng0vec 1705 READ(iukgmap,*) g0vec_all_r(:,ig0) 1706 ENDDO 1707 DO ig = 1, ngxx 1708 ! 1709 ! at variance with the nscf calculation, here gmap is read as a vector, 1710 ! 1711 READ(iukgmap,*) (gmap(ng0vec * (ig - 1) + ishift), ishift = 1, ng0vec) 1712 ENDDO 1713 ! 1714 CLOSE(iukgmap) 1715 ! 1716 ENDIF 1717 ! 1718 ! first node broadcasts arrays to all nodes 1719 ! 1720 CALL mp_bcast(g0vec_all_r, meta_ionode_id, world_comm) 1721 CALL mp_bcast(gmap, meta_ionode_id, world_comm) 1722 ! 1723 !----------------------------------------------------------------------- 1724 END SUBROUTINE readgmap 1725 !----------------------------------------------------------------------- 1726 ! 1727 !-------------------------------------------------------------- 1728 SUBROUTINE readkmap(nkstot) 1729 !-------------------------------------------------------------- 1730 !! 1731 !! Read the index of G_0 such that k+q+G_0 belongs to the 1st BZ 1732 !! 1733 ! 1734 USE kinds, ONLY : DP 1735 USE mp_global,ONLY : world_comm 1736 USE mp, ONLY : mp_bcast 1737 USE io_global,ONLY : meta_ionode, meta_ionode_id 1738 USE io_var, ONLY : iukmap 1739 USE io_files, ONLY : prefix 1740 USE elph2, ONLY : shift 1741 ! 1742 IMPLICIT NONE 1743 ! 1744 INTEGER, INTENT(in) :: nkstot 1745 !! Total number of k-points 1746 ! 1747 ! Local variables 1748 INTEGER :: ik 1749 !! Counter on k-points 1750 INTEGER :: ik1 1751 !! Temporary indices when reading kmap files 1752 INTEGER :: itmp 1753 !! Temporary indices when reading kmap files 1754 INTEGER :: ios 1755 !! Integer variable for I/O control 1756 ! 1757 IF (meta_ionode) THEN 1758 ! 1759 OPEN(iukmap, FILE = TRIM(prefix)//'.kmap', FORM = 'formatted', STATUS = 'old', IOSTAT = ios) 1760 IF (ios /= 0) CALL errore('readkmap', 'error opening kmap file', iukmap) 1761 DO ik = 1, nkstot 1762 READ(iukmap,*) ik1, itmp, shift(ik1) 1763 ENDDO 1764 CLOSE(iukmap) 1765 ! 1766 ENDIF 1767 ! 1768 ! first node broadcasts shift to all nodes 1769 ! 1770 CALL mp_bcast(shift, meta_ionode_id, world_comm) 1771 ! 1772 !----------------------------------------------------------------------- 1773 END SUBROUTINE readkmap 1774 !----------------------------------------------------------------------- 1775 ! 1776 !----------------------------------------------------------------------- 1777 SUBROUTINE openfilepw() 1778 !----------------------------------------------------------------------- 1779 !! 1780 !! Adapted from the code PH/openfilq - Quantum-ESPRESSO group 1781 !! This routine opens the WF files necessary for the EPW 1782 !! calculation. 1783 !! 1784 !! RM - Nov/Dec 2014 1785 !! Imported the noncolinear case implemented by xlzhang 1786 !! 1787 !----------------------------------------------------------------------- 1788 USE io_files, ONLY : prefix, diropn, seqopn 1789 USE units_lr, ONLY : iuwfc, lrwfc 1790 USE wvfct, ONLY : nbnd, npwx 1791 USE noncollin_module, ONLY : npol, nspin_mag 1792 USE units_ph, ONLY : lrdrho, lint3paw 1793 USE fft_base, ONLY : dfftp 1794 USE uspp_param, ONLY : nhm 1795 USE ions_base, ONLY : nat 1796 ! 1797 IMPLICIT NONE 1798 ! 1799 ! Local variables 1800 LOGICAL :: exst 1801 !! logical variable to check file existe 1802 ! 1803 IF (len_TRIM(prefix) == 0) CALL errore('openfilepw', 'wrong prefix', 1) 1804 ! 1805 ! The file with the wavefunctions 1806 ! 1807 iuwfc = 20 1808 lrwfc = 2 * nbnd * npwx * npol 1809 CALL diropn(iuwfc, 'wfc', lrwfc, exst) 1810 IF (.NOT. exst) CALL errore ('openfilepw', 'file ' // TRIM(prefix) // '.wfc' // ' not found', 1) 1811 ! 1812 ! file for setting unitary gauges of eigenstates 1813 ! 1814 lrdrho = 2 * dfftp%nr1x * dfftp%nr2x * dfftp%nr3x * nspin_mag 1815 lint3paw = 2 * nhm * nhm * nat * nspin_mag 1816 ! 1817 RETURN 1818 ! 1819 !---------------------------------------------------------------------------- 1820 END SUBROUTINE openfilepw 1821 !---------------------------------------------------------------------------- 1822 ! 1823 !---------------------------------------------------------------------------- 1824 SUBROUTINE param_get_range_vector(field, length, lcount, i_value) 1825 !---------------------------------------------------------------------------- 1826 !! 1827 !! Read a range vector eg. 1,2,3,4-10 or 1 3 400:100 1828 !! if(lcount) we return the number of states in length 1829 !! 1830 !! HL - April 2020 1831 !! Imported and adapted from the same name of subroutine in parameters.F90 1832 !! in the directory of src in Wannier90 1833 !! 1834 !---------------------------------------------------------------------------- 1835 ! 1836 IMPLICIT NONE 1837 ! 1838 CHARACTER(*), INTENT(in) :: field 1839 !! Field read for parsing 1840 INTEGER, INTENT(inout) :: length 1841 !! Number of states 1842 LOGICAL, INTENT(in) :: lcount 1843 !! If T only count states 1844 INTEGER, OPTIONAL, INTENT(out) :: i_value(length) 1845 !! States specified in range vector 1846 ! 1847 INTEGER :: loop 1848 !! Loop index 1849 INTEGER :: num1 1850 !! Integer number read 1851 INTEGER :: num2 1852 !! Integer number read 1853 INTEGER :: i_punc 1854 !! Position returned after scanning punctuation marks 1855 INTEGER :: counter 1856 !! Counter index 1857 INTEGER :: i_digit 1858 !! Position returned after scanning numbers 1859 INTEGER :: loop_r 1860 !! Loop index 1861 INTEGER :: range_size 1862 !! Size of range 1863 CHARACTER(LEN = 255) :: dummy 1864 !! Copy of field read for parsing 1865 CHARACTER(LEN = 10), PARAMETER :: c_digit = "0123456789" 1866 CHARACTER(LEN = 2), PARAMETER :: c_range = "-:" 1867 CHARACTER(LEN = 3), PARAMETER :: c_sep = " ,;" 1868 CHARACTER(LEN = 5), PARAMETER :: c_punc = " ,;-:" 1869 CHARACTER(LEN = 5) :: c_num1 1870 !! Number read 1871 CHARACTER(LEN = 5) :: c_num2 1872 !! Number read 1873 ! 1874 IF (lcount .AND. PRESENT(i_value)) & 1875 CALL errore('param_get_range_vector', 'incorrect call', 1) 1876 ! 1877 dummy = field 1878 dummy = ADJUSTL(dummy) 1879 ! 1880 counter = 0 1881 IF (LEN_TRIM(dummy) == 0) THEN 1882 length = counter 1883 RETURN 1884 ENDIF 1885 ! 1886 DO 1887 i_punc = SCAN(dummy, c_punc) 1888 IF (i_punc == 0) & 1889 CALL errore('param_get_range_vector', 'Error parsing field', 1) 1890 c_num1 = dummy(1:i_punc - 1) 1891 READ(c_num1, *, ERR = 101, END = 101) num1 1892 dummy = ADJUSTL(dummy(i_punc:)) 1893 !look for range 1894 IF (SCAN(dummy, c_range) == 1) THEN 1895 i_digit = SCAN(dummy, c_digit) 1896 dummy = ADJUSTL(dummy(i_digit:)) 1897 i_punc = SCAN(dummy, c_punc) 1898 c_num2 = dummy(1:i_punc - 1) 1899 READ(c_num2, *, ERR = 101, END = 101) num2 1900 dummy = ADJUSTL(dummy(i_punc:)) 1901 range_size = ABS(num2 - num1) + 1 1902 DO loop_r = 1, range_size 1903 counter = counter + 1 1904 IF (.NOT. lcount) i_value(counter) = MIN(num1, num2) + loop_r - 1 1905 ENDDO 1906 ELSE 1907 counter = counter + 1 1908 IF (.NOT. lcount) i_value(counter) = num1 1909 ENDIF 1910 IF (SCAN(dummy, c_sep) == 1) dummy = ADJUSTL(dummy(2:)) 1911 IF (SCAN(dummy, c_range) == 1) & 1912 CALL errore('param_get_range_vector', 'Error parsing field: incorrect range', 1) 1913 IF (INDEX(dummy, ' ') == 1) EXIT 1914 ENDDO 1915 ! 1916 IF (lcount) length = counter 1917 IF (.NOT. lcount) THEN 1918 DO loop = 1, counter - 1 1919 DO loop_r = loop + 1, counter 1920 IF (i_value(loop) == i_value(loop_r)) & 1921 CALL errore('param_get_range_vector', 'Error parsing field: duplicate values', 1) 1922 ENDDO 1923 ENDDO 1924 ENDIF 1925 ! 1926 RETURN 1927 ! 1928101 CALL errore('param_get_range_vector', 'Error parsing field', 1) 1929 ! 1930 !---------------------------------------------------------------------------- 1931 END SUBROUTINE param_get_range_vector 1932 !---------------------------------------------------------------------------- 1933 !------------------------------------------------------------------------------ 1934 END MODULE io_epw 1935 !------------------------------------------------------------------------------ 1936