1!!! this file contains post-espresso-5.1 commits: r11037, r11038 !!! 2 3!------------------------------------------------------------------------------- 4! 5! pw2bgw.f90 6! $MPIRUN pw2bgw.x -in pw2bgw.in > pw2bgw.out 7! converts Quantum ESPRESSO output files to BerkeleyGW input files 8! based on espresso-4.3.2/PP/pw_export.f90 9! written by G. Samsonidze (April 2008) 10! 11! The format of the input file for pw2bgw is described in the sample file pw2bgw.inp. 12! 13! NOTE: you CANNOT use USPP, PAW, or spinors for BerkeleyGW! 14! 15! real wavefunctions are constructed by applying the Gram-Schmidt process 16! based on paratecSGL/src/para/gwreal.f90 17! 18! make sure to obtain the whole sphere of G-vectors from espresso 19! do not use "K_POINTS gamma" as it produces only half the sphere 20! use "K_POINTS { tpiba | automatic | crystal }" even for the 21! Gamma-point calculation 22! 23!------------------------------------------------------------------------------- 24! 25! Compilation: use install.sh, or manually follow these steps: 26! 27! copy pw2bgw.f90 and make.bgw to espresso-4.3.2/PP 28! add this line to the end of espresso-4.3.2/PP/Makefile: "include make.bgw" 29! make as follows : cd espresso-4.3.2 ; make pw ph ; cd PP ; make pw2bgw.x 30! 31!------------------------------------------------------------------------------- 32! 33! subroutines 34! 35! write_wfng - generates complex wavefunctions in G-space (normalized to 1) 36! real_wfng - constructs real wavefunctions by applying the Gram-Schmidt 37! process (called from write_wfng) 38! write_rhog - generates complex charge density in G-space 39! (units of the number of electronic states per unit cell) 40! write_vxcg - generates complex exchange-correlation potential in G-space 41! (units of Rydberg) [only local part of Vxc] 42! write_vxc0 - prints complex exchange-correlation potential at G=0 43! (units of eV) [only local part of Vxc] 44! write_vxc - calculates matrix elements of exchange-correlation potential 45! in R-space [only local part of Vxc] or G-space [supports 46! non-local Vxc]. In units of eV. 47! write_vnlg - generates Kleinman-Bylander projectors in G-space 48! (units of Rydberg) 49! check_inversion - checks whether real/complex version is appropriate 50! write_header - writes header for VXC and RHO files 51! k_pools - determines start and end k-points for a pool 52! set_spin - sets parameters associated with LSDA and spinor cases 53! 54!------------------------------------------------------------------------------- 55! 56! Quantum ESPRESSO stores the wavefunctions in is-ik-ib-ig order 57! BerkeleyGW stores the wavefunctions in ik-ib-is-ig order 58! the outer loop is over is(QE)/ik(BGW) and the inner loop is over ig 59! ik = k-point index, is = spin index, ib = band index, ig = G-vector index 60! 61! write_wfng reverts the order of is and ik using smap and kmap arrays, 62! distributes wavefunctions over processors by ig (either real case or 63! spin-polarized case), calls real_wfng that applies the Gram-Schmidt 64! process (real case), reverts the order of is and ib (spin-polarized 65! case), and writes wavefunctions to disk 66! 67!------------------------------------------------------------------------------- 68 69#include "f_defs.h" 70 71module pw2bgw_m 72 73 implicit none 74 75 private 76 77 public :: & 78 write_wfng, & 79 real_wfng, & 80 write_rhog, & 81 write_vxcg, & 82 write_vxc0, & 83 write_vxc, & 84 write_vnlg, & 85 check_inversion, & 86 check_dfpt_modes, & 87 initialize_phonon, & 88 initialize_electric, & 89 dfpt_dwf 90 91contains 92 93!------------------------------------------------------------------------------- 94 95SUBROUTINE write_wfng ( output_file_name, real_or_complex, & 96 wfng_kgrid, wfng_nk1, wfng_nk2, wfng_nk3, wfng_dk1, wfng_dk2, & 97 wfng_dk3, wfng_occupation, wfng_nvmin, wfng_nvmax ) 98 99 USE cell_base, ONLY : omega, alat, tpiba, tpiba2, at, bg, & 100 ibrav, symm_type 101 USE constants, ONLY : pi, tpi, eps6 102 USE grid_dimensions, ONLY : nr1, nr2, nr3 103 USE gvect, ONLY : ngm, ngm_g, ig_l2g, g, mill, ecutrho 104 USE io_files, ONLY : iunwfc, nwordwfc 105 USE io_global, ONLY : ionode, ionode_id 106 USE ions_base, ONLY : nat, atm, ityp, tau 107 USE kinds, ONLY : DP 108 USE klist, ONLY : xk, wk, ngk, nks, nkstot 109 USE lsda_mod, ONLY : nspin, isk 110 USE mp, ONLY : mp_sum, mp_max, mp_get, mp_bcast, mp_barrier 111 USE mp_global, ONLY : mpime, nproc, world_comm, me_pool, & 112 root_pool, npool, nproc_pool, intra_pool_comm 113 USE mp_wave, ONLY : mergewf 114 USE start_k, ONLY : nk1, nk2, nk3, k1, k2, k3 115 USE symm_base, ONLY : s, ftau, nsym 116 USE wavefunctions_module, ONLY : evc 117 USE wvfct, ONLY : npwx, nbnd, npw, et, wg, g2kin, ecutwfc 118#ifdef __PARA 119 USE parallel_include, ONLY : MPI_DOUBLE_COMPLEX 120#endif 121 122 IMPLICIT NONE 123 124 character ( len = 256 ), intent (in) :: output_file_name 125 integer, intent (in) :: real_or_complex 126 logical, intent (in) :: wfng_kgrid 127 integer, intent (in) :: wfng_nk1 128 integer, intent (in) :: wfng_nk2 129 integer, intent (in) :: wfng_nk3 130 real (DP), intent (in) :: wfng_dk1 131 real (DP), intent (in) :: wfng_dk2 132 real (DP), intent (in) :: wfng_dk3 133 logical, intent (in) :: wfng_occupation 134 integer, intent (in) :: wfng_nvmin 135 integer, intent (in) :: wfng_nvmax 136 137 character :: cdate*9, ctime*9, sdate*32, stime*32, stitle*32 138 logical :: proc_wf, bad_kgrid 139 integer :: unit, i, j, k, cell_symmetry, nrecord 140 integer :: id, ib, ik, iks, ike, is, ig, ierr 141 integer :: nd, ntran, nb, nk_l, nk_g, ns, nst, nsf, ng_l, ng_g 142 integer :: ngg, npw_g, npwx_g 143 integer :: local_pw, ipsour, igwx, ngkdist_g, ngkdist_l, npw_g2 144 real (DP) :: alat2, recvol, dr1, t1 ( 3 ), t2 ( 3 ) 145 real (DP) :: r1 ( 3, 3 ), r2 ( 3, 3 ), adot ( 3, 3 ) 146 real (DP) :: bdot ( 3, 3 ), translation ( 3, 48 ) 147 integer, allocatable :: kmap ( : ) 148 integer, allocatable :: smap ( : ) 149 integer, allocatable :: ifmin ( : ) 150 integer, allocatable :: ifmax ( : ) 151 integer, allocatable :: itmp ( : ) 152 integer, allocatable :: ngk_g ( : ) 153 integer, allocatable :: ipmask ( : ) 154 integer, allocatable :: igwk ( : ) 155 integer, allocatable :: igwf_l2g ( : ) 156 integer, allocatable :: g_g ( :, : ) 157 integer, allocatable :: igk_l2g ( :, : ) 158 real (DP), allocatable :: et_g ( :, : ) 159 real (DP), allocatable :: wg_g ( :, : ) 160 real (DP), allocatable :: energy ( :, : ) 161 complex (DP), allocatable :: wfng ( : ) 162 complex (DP), allocatable :: wfng_buf ( :, : ) 163 complex (DP), allocatable :: wfng_dist ( :, :, : ) 164 165 INTEGER, EXTERNAL :: atomic_number 166 167 CALL check_inversion ( real_or_complex, nsym, s, nspin, .true., .true. ) 168 169 IF ( real_or_complex .EQ. 1 .OR. nspin .GT. 1 ) THEN 170 proc_wf = .TRUE. 171 ELSE 172 proc_wf = .FALSE. 173 ENDIF 174 175 bad_kgrid = .FALSE. 176 IF ( wfng_kgrid ) THEN 177 IF ( wfng_nk1 .LE. 0 .OR. wfng_nk2 .LE. 0 .OR. wfng_nk3 .LE. 0 ) & 178 bad_kgrid = .TRUE. 179 ELSE 180 IF ( nk1 .LE. 0 .OR. nk2 .LE. 0 .OR. nk3 .LE. 0 ) & 181 bad_kgrid = .TRUE. 182 ENDIF 183 IF ( bad_kgrid .AND. ionode ) THEN 184 WRITE ( 0, 101 ) 185 ENDIF 186 187 CALL date_and_tim ( cdate, ctime ) 188 WRITE ( sdate, '(A2,"-",A3,"-",A4,21X)' ) cdate(1:2), cdate(3:5), cdate(6:9) 189 WRITE ( stime, '(A8,24X)' ) ctime(1:8) 190 IF ( real_or_complex .EQ. 1 ) THEN 191 WRITE ( stitle, '("WFN-Real",24X)' ) 192 ELSE 193 WRITE ( stitle, '("WFN-Complex",21X)' ) 194 ENDIF 195 196 unit = 4 197 nrecord = 1 198 nd = 3 199 200 nb = nbnd 201 nk_l = nks 202 nk_g = nkstot 203 call set_spin(ns,nst,nsf,nspin) 204 205 ng_l = ngm 206 ng_g = ngm_g 207 208 CALL k_pools(iks, ike) 209 210 ALLOCATE ( kmap ( nk_g ) ) 211 ALLOCATE ( smap ( nk_g ) ) 212 213 DO i = 1, nk_g 214 j = ( i - 1 ) / ns 215 k = i - 1 - j * ns 216 kmap ( i ) = j + k * ( nk_g / ns ) + 1 217 smap ( i ) = k + 1 218 ENDDO 219 ierr = 0 220 DO i = 1, nk_g 221 ik = kmap ( i ) 222 is = smap ( i ) 223 IF ( ik .GE. iks .AND. ik .LE. ike .AND. is .NE. isk ( ik ) ) & 224 ierr = ierr + 1 225 ENDDO 226 CALL mp_max ( ierr ) 227 IF ( ierr .GT. 0 ) & 228 CALL errore ( 'write_wfng', 'smap', ierr ) 229 230 alat2 = alat ** 2 231 recvol = 8.0D0 * pi**3 / omega 232 233 DO i = 1, nd 234 DO j = 1, nd 235 adot ( j, i ) = 0.0D0 236 ENDDO 237 ENDDO 238 DO i = 1, nd 239 DO j = 1, nd 240 DO k = 1, nd 241 adot ( j, i ) = adot ( j, i ) + & 242 at ( k, j ) * at ( k, i ) * alat2 243 ENDDO 244 ENDDO 245 ENDDO 246 247 DO i = 1, nd 248 DO j = 1, nd 249 bdot ( j, i ) = 0.0D0 250 ENDDO 251 ENDDO 252 DO i = 1, nd 253 DO j = 1, nd 254 DO k = 1, nd 255 bdot ( j, i ) = bdot ( j, i ) + & 256 bg ( k, j ) * bg ( k, i ) * tpiba2 257 ENDDO 258 ENDDO 259 ENDDO 260 261 ierr = 0 262 IF ( ibrav .EQ. 0 ) THEN 263 IF ( TRIM ( symm_type ) .EQ. 'cubic' ) THEN 264 cell_symmetry = 0 265 ELSEIF ( TRIM ( symm_type ) .EQ. 'hexagonal' ) THEN 266 cell_symmetry = 1 267 ELSE 268 ierr = 1 269 ENDIF 270 ELSEIF ( ibrav .GE. 1 .AND. ibrav .LE. 3 ) THEN 271 cell_symmetry = 0 272 ELSEIF ( ibrav .GE. 4 .AND. ibrav .LE. 5 ) THEN 273 cell_symmetry = 1 274 ELSEIF ( ibrav .GE. 6 .AND. ibrav .LE. 14 ) THEN 275 cell_symmetry = 0 276 ELSE 277 ierr = 1 278 ENDIF 279 IF ( ierr .GT. 0 ) & 280 CALL errore ( 'write_wfng', 'cell_symmetry', ierr ) 281 282 ntran = nsym 283 DO i = 1, ntran 284 DO j = 1, nd 285 DO k = 1, nd 286 r1 ( k, j ) = dble ( s ( k, j, i ) ) 287 ENDDO 288 ENDDO 289 CALL invmat ( 3, r1, r2, dr1 ) 290 t1 ( 1 ) = dble ( ftau ( 1, i ) ) / dble ( nr1 ) 291 t1 ( 2 ) = dble ( ftau ( 2, i ) ) / dble ( nr2 ) 292 t1 ( 3 ) = dble ( ftau ( 3, i ) ) / dble ( nr3 ) 293 DO j = 1, nd 294 t2 ( j ) = 0.0D0 295 DO k = 1, nd 296 t2 ( j ) = t2 ( j ) + r2 ( k, j ) * t1 ( k ) 297 ENDDO 298 IF ( t2 ( j ) .GE. eps6 + 0.5D0 ) & 299 t2 ( j ) = t2 ( j ) - dble ( int ( t2 ( j ) + 0.5D0 ) ) 300 IF ( t2 ( j ) .LT. eps6 - 0.5D0 ) & 301 t2 ( j ) = t2 ( j ) - dble ( int ( t2 ( j ) - 0.5D0 ) ) 302 ENDDO 303 DO j = 1, nd 304 translation ( j, i ) = t2 ( j ) * tpi 305 ENDDO 306 ENDDO 307 308 ALLOCATE ( et_g ( nb, nk_g ) ) 309 310 DO ik = 1, nk_l 311 DO ib = 1, nb 312 et_g ( ib, ik ) = et ( ib, ik ) 313 ENDDO 314 ENDDO 315#ifdef __PARA 316 CALL poolrecover ( et_g, nb, nk_g, nk_l ) 317 CALL mp_bcast ( et_g, ionode_id ) 318#endif 319 320 ALLOCATE ( wg_g ( nb, nk_g ) ) 321 ALLOCATE ( ifmin ( nk_g ) ) 322 ALLOCATE ( ifmax ( nk_g ) ) 323 324 IF ( wfng_occupation ) THEN 325 326 DO ik = 1, nk_g 327 DO ib = 1, nb 328 IF ( ib .GE. wfng_nvmin .AND. ib .LE. wfng_nvmax ) THEN 329 wg_g ( ib, ik ) = 1.0D0 330 ELSE 331 wg_g ( ib, ik ) = 0.0D0 332 ENDIF 333 ENDDO 334 ENDDO 335 DO ik = 1, nk_g 336 ifmin ( ik ) = wfng_nvmin 337 ENDDO 338 DO ik = 1, nk_g 339 ifmax ( ik ) = wfng_nvmax 340 ENDDO 341 342 ELSE 343 344 DO ik = 1, nk_l 345 DO ib = 1, nb 346 wg_g ( ib, ik ) = wg ( ib, ik ) 347 IF ( abs ( wk ( ik ) ) .GT. eps6 ) THEN 348 wg_g ( ib, ik ) = wg_g ( ib, ik ) / wk ( ik ) 349 ENDIF 350 ENDDO 351 ENDDO 352#ifdef __PARA 353 CALL poolrecover ( wg_g, nb, nk_g, nk_l ) 354#endif 355 DO ik = 1, nk_g 356 ifmin ( ik ) = 0 357 ENDDO 358 DO ik = 1, nk_g 359 ifmax ( ik ) = 0 360 ENDDO 361 DO ik = 1, nk_g 362 DO ib = 1, nb 363 IF ( wg_g( ib, ik ) .GT. 0.5D0 ) THEN 364 IF ( ifmin ( ik ) .EQ. 0 ) ifmin ( ik ) = ib 365 ifmax ( ik ) = ib 366 ENDIF 367 ENDDO 368 ENDDO 369 370 ENDIF 371 372 ALLOCATE ( g_g ( nd, ng_g ) ) 373 374 DO ig = 1, ng_g 375 DO id = 1, nd 376 g_g ( id, ig ) = 0 377 ENDDO 378 ENDDO 379 DO ig = 1, ng_l 380 g_g ( 1, ig_l2g ( ig ) ) = mill ( 1, ig ) 381 g_g ( 2, ig_l2g ( ig ) ) = mill ( 2, ig ) 382 g_g ( 3, ig_l2g ( ig ) ) = mill ( 3, ig ) 383 ENDDO 384 CALL mp_sum ( g_g, intra_pool_comm ) 385 386 ALLOCATE ( igk_l2g ( npwx, nk_l ) ) 387 388 ALLOCATE ( itmp ( npwx ) ) 389 DO ik = 1, nk_l 390 DO i = 1, npwx 391 itmp ( i ) = 0 392 ENDDO 393 npw = npwx 394 CALL gk_sort ( xk ( 1, ik + iks - 1 ), ng_l, g, ecutwfc / tpiba2, & 395 npw, itmp ( 1 ), g2kin ) 396 DO ig = 1, npw 397 igk_l2g ( ig, ik ) = ig_l2g ( itmp ( ig ) ) 398 ENDDO 399 DO ig = npw + 1, npwx 400 igk_l2g ( ig, ik ) = 0 401 ENDDO 402 ngk ( ik ) = npw 403 ENDDO 404 DEALLOCATE ( itmp ) 405 406 ALLOCATE ( ngk_g ( nk_g ) ) 407 408 DO ik = 1, nk_g 409 ngk_g ( ik ) = 0 410 ENDDO 411 DO ik = 1, nk_l 412 ngk_g ( ik + iks - 1 ) = ngk ( ik ) 413 ENDDO 414 CALL mp_sum ( ngk_g ) 415 416 npw_g = MAXVAL ( igk_l2g ( :, : ) ) 417 CALL mp_max ( npw_g ) 418 419 npwx_g = MAXVAL ( ngk_g ( : ) ) 420 421 CALL cryst_to_cart ( nk_g / ns, xk, at, - 1 ) 422 423 IF ( ionode ) THEN 424 OPEN ( unit = unit, file = TRIM ( output_file_name ), & 425 form = 'unformatted', status = 'replace' ) 426 WRITE ( unit ) stitle, sdate, stime 427 WRITE ( unit ) nspin, ng_g, ntran, cell_symmetry, nat, ecutrho, & 428 nk_g / ns, nb, npwx_g, ecutwfc 429 IF ( wfng_kgrid ) THEN 430 WRITE ( unit ) nr1, nr2, nr3, wfng_nk1, wfng_nk2, wfng_nk3, & 431 wfng_dk1, wfng_dk2, wfng_dk3 432 ELSE 433 WRITE ( unit ) nr1, nr2, nr3, nk1, nk2, nk3, & 434 0.5D0 * dble ( k1 ), 0.5D0 * dble ( k2 ), 0.5D0 * dble ( k3 ) 435 ENDIF 436 WRITE ( unit ) omega, alat, ( ( at ( j, i ), j = 1, nd ), i = 1, nd ), & 437 ( ( adot ( j, i ), j = 1, nd ), i = 1, nd ) 438 WRITE ( unit ) recvol, tpiba, ( ( bg ( j, i ), j = 1, nd ), i = 1, nd ), & 439 ( ( bdot ( j, i ), j = 1, nd ), i = 1, nd ) 440 WRITE ( unit ) ( ( ( s ( k, j, i ), k = 1, nd ), j = 1, nd ), i = 1, ntran ) 441 WRITE ( unit ) ( ( translation ( j, i ), j = 1, nd ), i = 1, ntran ) 442 WRITE ( unit ) ( ( tau ( j, i ), j = 1, nd ), atomic_number ( atm ( ityp ( i ) ) ), i = 1, nat ) 443 WRITE ( unit ) ( ngk_g ( ik ), ik = 1, nk_g / ns ) 444 WRITE ( unit ) ( wk ( ik ) * dble ( nst ) / 2.0D0, ik = 1, nk_g / ns ) 445 WRITE ( unit ) ( ( xk ( id, ik ), id = 1, nd ), ik = 1, nk_g / ns ) 446 WRITE ( unit ) ( ifmin ( ik ), ik = 1, nk_g ) 447 WRITE ( unit ) ( ifmax ( ik ), ik = 1, nk_g ) 448 WRITE ( unit ) ( ( et_g ( ib, ik ), ib = 1, nb ), ik = 1, nk_g ) 449 WRITE ( unit ) ( ( wg_g ( ib, ik ), ib = 1, nb ), ik = 1, nk_g ) 450 WRITE ( unit ) nrecord 451 WRITE ( unit ) ng_g 452 WRITE ( unit ) ( ( g_g ( id, ig ), id = 1, nd ), ig = 1, ng_g ) 453 ENDIF 454 455 CALL cryst_to_cart ( nk_g / ns, xk, bg, 1 ) 456 457 DEALLOCATE ( wg_g ) 458 DEALLOCATE ( ifmax ) 459 DEALLOCATE ( ifmin ) 460 461 ALLOCATE ( igwk ( npwx_g ) ) 462 463 IF ( proc_wf ) THEN 464 IF ( nspin == 4) THEN 465 ELSE 466 IF ( MOD ( npwx_g, nproc ) .EQ. 0 ) THEN 467 ngkdist_l = npwx_g / nproc 468 ELSE 469 ngkdist_l = npwx_g / nproc + 1 470 ENDIF 471 ENDIF 472 ngkdist_g = ngkdist_l * nproc 473 IF ( real_or_complex .EQ. 1 ) & 474 ALLOCATE ( energy ( nb, ns ) ) 475 ALLOCATE ( wfng_buf ( ngkdist_g, ns ) ) 476 ALLOCATE ( wfng_dist ( ngkdist_l, nb, ns ) ) 477 ENDIF 478 479 DO i = 1, nk_g 480 481 ik = kmap ( i ) 482 is = smap ( i ) 483 484 IF ( real_or_complex .EQ. 1 ) THEN 485 DO ib = 1, nb 486 energy ( ib, is ) = et_g ( ib, i ) 487 ENDDO 488 ENDIF 489 490 DO j = 1, npwx_g 491 igwk ( j ) = 0 492 ENDDO 493 ALLOCATE ( itmp ( npw_g ) ) 494 DO j = 1, npw_g 495 itmp ( j ) = 0 496 ENDDO 497 IF ( ik .GE. iks .AND. ik .LE. ike ) THEN 498 DO ig = 1, ngk ( ik - iks + 1 ) 499 itmp ( igk_l2g ( ig, ik - iks + 1 ) ) = igk_l2g ( ig, ik - iks + 1 ) 500 ENDDO 501 ENDIF 502 CALL mp_sum ( itmp ) 503 ngg = 0 504 DO ig = 1, npw_g 505 IF ( itmp ( ig ) .EQ. ig ) THEN 506 ngg = ngg + 1 507 igwk ( ngg ) = ig 508 ENDIF 509 ENDDO 510 DEALLOCATE ( itmp ) 511 512 IF ( ionode ) THEN 513 IF ( is .EQ. 1 ) THEN 514 WRITE ( unit ) nrecord 515 WRITE ( unit ) ngk_g ( ik ) 516 WRITE ( unit ) ( ( g_g ( id, igwk ( ig ) ), id = 1, nd ), & 517 ig = 1, ngk_g ( ik ) ) 518 ENDIF 519 ENDIF 520 521 local_pw = 0 522 IF ( ik .GE. iks .AND. ik .LE. ike ) THEN 523 CALL davcio ( evc, nwordwfc, iunwfc, ik - iks + 1, - 1 ) 524 local_pw = ngk ( ik - iks + 1 ) 525 ENDIF 526 527 ALLOCATE ( igwf_l2g ( local_pw ) ) 528 529 DO ig = 1, local_pw 530 igwf_l2g ( ig ) = 0 531 ENDDO 532 DO ig = 1, local_pw 533 ngg = igk_l2g ( ig, ik - iks + 1 ) 534 DO j = 1, ngk_g ( ik ) 535 IF ( ngg .EQ. igwk ( j ) ) THEN 536 igwf_l2g ( ig ) = j 537 EXIT 538 ENDIF 539 ENDDO 540 ENDDO 541 542 ALLOCATE ( ipmask ( nproc ) ) 543 DO j = 1, nproc 544 ipmask ( j ) = 0 545 ENDDO 546 ipsour = ionode_id 547 IF ( npool .GT. 1 ) THEN 548 IF ( ( ik .GE. iks ) .AND. ( ik .LE. ike ) ) THEN 549 IF ( me_pool .EQ. root_pool ) ipmask ( mpime + 1 ) = 1 550 ENDIF 551 CALL mp_sum ( ipmask ) 552 DO j = 1, nproc 553 IF ( ipmask ( j ) .EQ. 1 ) ipsour = j - 1 554 ENDDO 555 ENDIF 556 DEALLOCATE ( ipmask ) 557 558 igwx = 0 559 ierr = 0 560 IF ( ik .GE. iks .AND. ik .LE. ike ) & 561 igwx = MAXVAL ( igwf_l2g ( 1 : local_pw ) ) 562 CALL mp_max ( igwx, intra_pool_comm ) 563 IF ( ipsour .NE. ionode_id ) & 564 CALL mp_get ( igwx, igwx, mpime, ionode_id, ipsour, 1 ) 565 ierr = 0 566 IF ( ik .GE. iks .AND. ik .LE. ike .AND. igwx .NE. ngk_g ( ik ) ) & 567 ierr = 1 568 CALL mp_max ( ierr ) 569 IF ( ierr .GT. 0 ) & 570 CALL errore ( 'write_wfng', 'igwx ngk_g', ierr ) 571 572 ALLOCATE ( wfng ( MAX ( 1, igwx ) ) ) 573 574 575 DO ib = 1, nb 576 577 DO j = 1, igwx 578 wfng ( j ) = ( 0.0D0, 0.0D0 ) 579 ENDDO 580 IF ( npool .GT. 1 ) THEN 581 IF ( ( ik .GE. iks ) .AND. ( ik .LE. ike ) ) THEN 582 CALL mergewf ( evc ( :, ib ), wfng, local_pw, igwf_l2g, & 583 me_pool, nproc_pool, root_pool, intra_pool_comm ) 584 ENDIF 585 IF ( ipsour .NE. ionode_id ) THEN 586 CALL mp_get ( wfng, wfng, mpime, ionode_id, ipsour, ib, & 587 world_comm ) 588 ENDIF 589 ELSE 590 CALL mergewf ( evc ( :, ib ), wfng, local_pw, igwf_l2g, & 591 mpime, nproc, ionode_id, world_comm ) 592 593 ENDIF 594 595 IF ( proc_wf ) THEN 596 DO ig = 1, igwx 597 wfng_buf ( ig, is ) = wfng ( ig ) 598 ENDDO 599 DO ig = igwx + 1, ngkdist_g 600 wfng_buf ( ig, is ) = ( 0.0D0, 0.0D0 ) 601 ENDDO 602#ifdef __PARA 603 CALL mp_barrier ( ) 604 CALL MPI_Scatter ( wfng_buf ( :, is ), ngkdist_l, MPI_DOUBLE_COMPLEX, & 605 wfng_dist ( :, ib, is ), ngkdist_l, MPI_DOUBLE_COMPLEX, & 606 ionode_id, world_comm, ierr ) 607 IF ( ierr .GT. 0 ) & 608 CALL errore ( 'write_wfng', 'mpi_scatter', ierr ) 609#else 610 DO ig = 1, ngkdist_g 611 wfng_dist ( ig, ib, is ) = wfng_buf ( ig, is ) 612 ENDDO 613#endif 614 ELSE 615 IF ( ionode ) THEN 616 WRITE ( unit ) nrecord 617 WRITE ( unit ) ngk_g ( ik ) 618 WRITE ( unit ) ( wfng ( ig ), ig = 1, igwx ) 619 ENDIF 620 ENDIF 621 622 ENDDO 623 624 DEALLOCATE ( wfng ) 625 DEALLOCATE ( igwf_l2g ) 626 627 IF ( proc_wf .AND. is .EQ. ns ) THEN 628 IF ( real_or_complex .EQ. 1 ) THEN 629 CALL start_clock ( 'real_wfng' ) 630 CALL real_wfng ( ik, ngkdist_l, nb, ns, energy, wfng_dist ) 631 CALL stop_clock ( 'real_wfng' ) 632 ENDIF 633 DO ib = 1, nb 634 DO is = 1, ns 635#ifdef __PARA 636 CALL mp_barrier ( ) 637 CALL MPI_Gather ( wfng_dist ( :, ib, is ), ngkdist_l, & 638 MPI_DOUBLE_COMPLEX, wfng_buf ( :, is ), ngkdist_l, & 639 MPI_DOUBLE_COMPLEX, ionode_id, world_comm, ierr ) 640 IF ( ierr .GT. 0 ) & 641 CALL errore ( 'write_wfng', 'mpi_gather', ierr ) 642#else 643 DO ig = 1, ngkdist_g 644 wfng_buf ( ig, is ) = wfng_dist ( ig, ib, is ) 645 ENDDO 646#endif 647 ENDDO 648 IF ( ionode ) THEN 649 WRITE ( unit ) nrecord 650 WRITE ( unit ) ngk_g ( ik ) 651 IF ( real_or_complex .EQ. 1 ) THEN 652 WRITE ( unit ) ( ( dble ( wfng_buf ( ig, is ) ), & 653 ig = 1, igwx ), is = 1, ns ) 654 ELSE 655 WRITE ( unit ) ( ( wfng_buf ( ig, is ), & 656 ig = 1, igwx ), is = 1, ns ) 657 ENDIF 658 ENDIF 659 ENDDO 660 ENDIF 661 662 ENDDO 663 664 DEALLOCATE ( igwk ) 665 DEALLOCATE ( ngk_g ) 666 DEALLOCATE ( igk_l2g ) 667 DEALLOCATE ( et_g ) 668 669 IF ( proc_wf ) THEN 670 IF ( real_or_complex .EQ. 1 ) & 671 DEALLOCATE ( energy ) 672 DEALLOCATE ( wfng_buf ) 673 DEALLOCATE ( wfng_dist ) 674 ENDIF 675 676 IF ( ionode ) THEN 677 CLOSE ( unit = unit, status = 'keep' ) 678 ENDIF 679 680 DEALLOCATE ( g_g ) 681 DEALLOCATE ( smap ) 682 DEALLOCATE ( kmap ) 683 684 CALL mp_barrier ( ) 685 686 RETURN 687 688101 FORMAT ( /, 5X, "WARNING: kgrid is set to zero in the wavefunction file.", & 689 /, 14X, "The resulting file will only be usable as the fine grid in inteqp.", / ) 690 691END SUBROUTINE write_wfng 692 693!------------------------------------------------------------------------------- 694 695SUBROUTINE real_wfng ( ik, ngkdist_l, nb, ns, energy, wfng_dist ) 696 697 USE constants, ONLY : eps6 698 USE kinds, ONLY : DP 699 USE io_global, ONLY : ionode 700 USE mp, ONLY : mp_sum 701 702 IMPLICIT NONE 703 704 integer, intent (in) :: ik, ngkdist_l, nb, ns 705 real (DP), intent (in) :: energy ( nb, ns ) 706 complex (DP), intent (inout) :: wfng_dist ( ngkdist_l, nb, ns ) 707 708 real (DP), PARAMETER :: eps2 = 1.0D-2 709 real (DP), PARAMETER :: eps5 = 1.0D-5 710 711 character :: tmpstr*80 712 integer :: i, j, k, is, ib, jb, ig, inum, deg, mdeg, inc 713 integer :: dimension_span, reduced_span, ierr 714 real (DP) :: x 715 integer, allocatable :: imap ( :, : ) 716 integer, allocatable :: inums ( : ) 717 integer, allocatable :: inull ( : ) 718 integer, allocatable :: null_map ( :, : ) 719 real (DP), allocatable :: psi ( :, : ) 720 real (DP), allocatable :: phi ( :, : ) 721 real (DP), allocatable :: vec ( : ) 722 complex (DP), allocatable :: wfc ( : ) 723 724 mdeg = 1 725 DO is = 1, ns 726 DO ib = 1, nb - 1 727 deg = 1 728 DO jb = ib + 1, nb 729 IF ( abs ( energy ( ib, is ) - energy ( jb, is ) ) & 730 .LT. eps5 * dble ( jb - ib + 1 ) ) deg = deg + 1 731 ENDDO 732 IF ( deg .GT. mdeg ) mdeg = deg 733 ENDDO 734 ENDDO 735 mdeg = mdeg * 2 736 737 ALLOCATE ( imap ( nb, ns ) ) 738 ALLOCATE ( inums ( ns ) ) 739 ALLOCATE ( inull ( nb ) ) 740 ALLOCATE ( null_map ( mdeg, nb ) ) 741 742 DO is = 1, ns 743 inum = 1 744 DO ib = 1, nb 745 IF ( ib .EQ. nb ) THEN 746 imap ( inum, is ) = ib 747 inum = inum + 1 748 ELSEIF ( abs ( energy ( ib, is ) - & 749 energy ( ib + 1, is ) ) .GT. eps5 ) THEN 750 imap ( inum, is ) = ib 751 inum = inum + 1 752 ENDIF 753 ENDDO 754 inum = inum - 1 755 inums ( is ) = inum 756 ENDDO 757 758 ALLOCATE ( wfc ( ngkdist_l ) ) 759 ALLOCATE ( psi ( ngkdist_l, mdeg ) ) 760 ALLOCATE ( phi ( ngkdist_l, mdeg ) ) 761 ALLOCATE ( vec ( ngkdist_l ) ) 762 763 DO is = 1, ns 764 inc = 1 765 inum = inums ( is ) 766 DO i = 1, inum 767 inull ( i ) = 1 768 DO ib = inc, imap ( i, is ) 769 DO ig = 1, ngkdist_l 770 wfc ( ig ) = wfng_dist ( ig, ib, is ) 771 ENDDO 772 x = 0.0D0 773 DO ig = 1, ngkdist_l 774 x = x + dble ( wfc ( ig ) ) **2 775 ENDDO 776 CALL mp_sum ( x ) 777 IF ( x .LT. eps2 ) null_map ( inull ( i ), i ) = 0 778 IF ( x .GT. eps2 ) null_map ( inull ( i ), i ) = 1 779 inull ( i ) = inull ( i ) + 1 780 x = 0.0D0 781 DO ig = 1, ngkdist_l 782 x = x + aimag ( wfc ( ig ) ) **2 783 ENDDO 784 CALL mp_sum ( x ) 785 IF ( x .LT. eps2 ) null_map ( inull ( i ), i ) = 0 786 IF ( x .GT. eps2 ) null_map ( inull ( i ), i ) = 1 787 inull ( i ) = inull ( i ) + 1 788 ENDDO 789 inull ( i ) = inull ( i ) - 1 790 inc = imap ( i, is ) + 1 791 ENDDO 792 inc = 1 793 ib = 1 794 DO i = 1, inum 795 k = 1 796 DO j = 1, 2 * ( imap ( i, is ) - inc ) + 1, 2 797 IF ( null_map ( j, i ) .EQ. 1 .OR. & 798 null_map ( j + 1, i ) .EQ. 1 ) THEN 799 DO ig = 1, ngkdist_l 800 wfc ( ig ) = wfng_dist ( ig, ib, is ) 801 ENDDO 802 IF ( null_map ( j, i ) .EQ. 1 ) THEN 803 DO ig = 1, ngkdist_l 804 phi ( ig, k ) = dble ( wfc ( ig ) ) 805 ENDDO 806 k = k + 1 807 ENDIF 808 IF ( null_map ( j + 1, i ) .EQ. 1 ) THEN 809 DO ig = 1, ngkdist_l 810 phi ( ig, k ) = aimag ( wfc ( ig ) ) 811 ENDDO 812 k = k + 1 813 ENDIF 814 ib = ib + 1 815 ENDIF 816 ENDDO 817 dimension_span = k - 1 818 IF ( dimension_span .EQ. 0 ) THEN 819 ierr = 201 820 WRITE ( tmpstr, 201 ) ik, is, inc 821 CALL errore ( 'real_wfng', tmpstr, ierr ) 822 ENDIF 823 DO j = 1, dimension_span 824 x = 0.0D0 825 DO ig = 1, ngkdist_l 826 x = x + phi ( ig, j ) **2 827 ENDDO 828 CALL mp_sum ( x ) 829 x = sqrt ( x ) 830 DO ig = 1, ngkdist_l 831 phi ( ig, j ) = phi ( ig, j ) / x 832 ENDDO 833 ENDDO 834! 835! the Gram-Schmidt process begins 836! 837 reduced_span = 1 838 DO ig = 1, ngkdist_l 839 psi ( ig, 1 ) = phi ( ig, 1 ) 840 ENDDO 841 DO j = 1, dimension_span - 1 842 DO ig = 1, ngkdist_l 843 vec ( ig ) = phi ( ig, j + 1 ) 844 ENDDO 845 DO k = 1, reduced_span 846 x = 0.0D0 847 DO ig = 1, ngkdist_l 848 x = x + phi ( ig, j + 1 ) * psi ( ig, k ) 849 ENDDO 850 CALL mp_sum ( x ) 851 DO ig = 1, ngkdist_l 852 vec ( ig ) = vec ( ig ) - psi ( ig, k ) * x 853 ENDDO 854 ENDDO 855 x = 0.0D0 856 DO ig = 1, ngkdist_l 857 x = x + vec ( ig ) **2 858 ENDDO 859 CALL mp_sum ( x ) 860 x = sqrt ( x ) 861 IF ( x .GT. eps6 ) THEN 862 reduced_span = reduced_span + 1 863 DO ig = 1, ngkdist_l 864 psi ( ig, reduced_span ) = vec ( ig ) / x 865 ENDDO 866 ENDIF 867 ENDDO 868! 869! the Gram-Schmidt process ends 870! 871 IF ( reduced_span .LT. imap ( i, is ) - inc + 1 ) THEN 872 ierr = 202 873 WRITE ( tmpstr, 202 ) ik, is, inc 874 CALL errore ( 'real_wfng', tmpstr, ierr ) 875 ENDIF 876 DO ib = inc, imap ( i, is ) 877 DO ig = 1, ngkdist_l 878 wfng_dist ( ig, ib, is ) = & 879 CMPLX ( psi ( ig, ib - inc + 1 ), 0.0D0 ) 880 ENDDO 881 ENDDO 882 inc = imap ( i, is ) + 1 883 ENDDO 884 ENDDO 885 886 DEALLOCATE ( vec ) 887 DEALLOCATE ( phi ) 888 DEALLOCATE ( psi ) 889 DEALLOCATE ( wfc ) 890 DEALLOCATE ( null_map ) 891 DEALLOCATE ( inull ) 892 DEALLOCATE ( inums ) 893 DEALLOCATE ( imap ) 894 895 RETURN 896 897201 FORMAT("failed Gram-Schmidt dimension span for kpoint =",i6," spin =",i2," band =",i6) 898202 FORMAT("failed Gram-Schmidt reduced span for kpoint =",i6," spin =",i2," band =",i6) 899 900END SUBROUTINE real_wfng 901 902!------------------------------------------------------------------------------- 903 904SUBROUTINE write_rhog ( output_file_name, real_or_complex ) 905 906 USE cell_base, ONLY : omega 907 USE gvect, ONLY : ngm, ngm_g, ig_l2g, mill, ecutrho 908 USE io_global, ONLY : ionode 909 USE ions_base, ONLY : nat, atm, ityp, tau 910 USE kinds, ONLY : DP 911 USE lsda_mod, ONLY : nspin 912 USE mp, ONLY : mp_sum 913 USE mp_global, ONLY : intra_pool_comm 914 USE scf, ONLY : rho 915 USE symm_base, ONLY : s, ftau, nsym 916 917 IMPLICIT NONE 918 919 character ( len = 256 ), intent (in) :: output_file_name 920 integer, intent (in) :: real_or_complex 921 922 character :: stitle*32 923 integer :: unit, is, ig 924 integer :: ns, nst, nsf, ng_l, ng_g 925 integer :: nrecord 926 complex (DP), allocatable :: rhog_g ( :, : ) 927 928 INTEGER, EXTERNAL :: atomic_number 929 930 CALL check_inversion ( real_or_complex, nsym, s, nspin, .true., .true. ) 931 932 IF ( real_or_complex .EQ. 1 ) THEN 933 WRITE ( stitle, '("RHO-Real",24X)' ) 934 ELSE 935 WRITE ( stitle, '("RHO-Complex",21X)' ) 936 ENDIF 937 938 unit = 4 939 nrecord = 1 940 941 ng_l = ngm 942 ng_g = ngm_g 943 call set_spin(ns,nst,nsf,nspin) 944 945 ALLOCATE ( rhog_g ( ng_g, ns ) ) 946 947 DO is = 1, ns 948 DO ig = 1, ng_g 949 rhog_g ( ig, is ) = ( 0.0D0, 0.0D0 ) 950 ENDDO 951 ENDDO 952 953 DO is = 1, ns 954 DO ig = 1, ng_l 955 rhog_g ( ig_l2g ( ig ), is ) = rho%of_g ( ig, is ) 956 ENDDO 957 ENDDO 958 959 CALL mp_sum ( rhog_g, intra_pool_comm ) 960 961 DO is = 1, ns 962 DO ig = 1, ng_g 963 rhog_g ( ig, is ) = rhog_g ( ig, is ) * CMPLX ( omega, 0.0D0 ) 964 ENDDO 965 ENDDO 966 967 IF ( ionode ) OPEN ( unit = unit, file = TRIM ( output_file_name ), & 968 form = 'unformatted', status = 'replace' ) 969 970 CALL write_header(unit, stitle) 971 972 IF ( ionode ) THEN 973 WRITE ( unit ) nrecord 974 WRITE ( unit ) ng_g 975 IF ( real_or_complex .EQ. 1 ) THEN 976 WRITE ( unit ) ( ( dble ( rhog_g ( ig, is ) ), & 977 ig = 1, ng_g ), is = 1, ns ) 978 ELSE 979 WRITE ( unit ) ( ( rhog_g ( ig, is ), & 980 ig = 1, ng_g ), is = 1, ns ) 981 ENDIF 982 CLOSE ( unit = unit, status = 'keep' ) 983 ENDIF 984 985 DEALLOCATE ( rhog_g ) 986 987 RETURN 988 989END SUBROUTINE write_rhog 990 991!------------------------------------------------------------------------------- 992 993SUBROUTINE write_vxcg ( output_file_name, real_or_complex, & 994 vxc_zero_rho_core, input_dft, exx_flag ) 995 996 USE cell_base, ONLY : omega 997 USE ener, ONLY : etxc, vtxc 998 USE fft_base, ONLY : dfftp 999 USE fft_interfaces, ONLY : fwfft 1000 USE funct, ONLY : enforce_input_dft 1001 USE grid_dimensions, ONLY : nrxx 1002 USE gvect, ONLY : ngm, ngm_g, ig_l2g, nl, mill, ecutrho 1003 USE io_global, ONLY : ionode 1004 USE ions_base, ONLY : nat, atm, ityp, tau 1005 USE kinds, ONLY : DP 1006 USE lsda_mod, ONLY : nspin 1007 USE mp, ONLY : mp_sum 1008 USE mp_global, ONLY : intra_pool_comm 1009 USE scf, ONLY : rho, rho_core, rhog_core 1010 USE symm_base, ONLY : s, ftau, nsym 1011 USE wavefunctions_module, ONLY : psic 1012#ifdef EXX 1013 USE funct, ONLY : start_exx, stop_exx 1014#endif 1015 1016 IMPLICIT NONE 1017 1018 character ( len = 256 ), intent (in) :: output_file_name 1019 integer, intent (in) :: real_or_complex 1020 logical, intent (in) :: vxc_zero_rho_core 1021 character ( len = 20 ), intent (in) :: input_dft 1022 logical, intent (in) :: exx_flag 1023 1024 character :: stitle*32 1025 integer :: unit, is, ir, ig 1026 integer :: ns, nst, nsf, nr, ng_l, ng_g 1027 integer :: nrecord 1028 real (DP), allocatable :: vxcr_g ( :, : ) 1029 complex (DP), allocatable :: vxcg_g ( :, : ) 1030 1031 INTEGER, EXTERNAL :: atomic_number 1032 1033 CALL check_inversion ( real_or_complex, nsym, s, nspin, .true., .true. ) 1034 1035 IF ( real_or_complex .EQ. 1 ) THEN 1036 WRITE ( stitle, '("VXC-Real",24X)' ) 1037 ELSE 1038 WRITE ( stitle, '("VXC-Complex",21X)' ) 1039 ENDIF 1040 1041 unit = 4 1042 nrecord = 1 1043 1044 call set_spin(ns,nst,nsf,nspin) 1045 nr = nrxx 1046 ng_l = ngm 1047 ng_g = ngm_g 1048 1049 ALLOCATE ( vxcr_g ( nr, nsf ) ) 1050 ALLOCATE ( vxcg_g ( ng_g, nsf ) ) 1051 1052 DO is = 1, nsf 1053 DO ig = 1, ng_g 1054 vxcg_g ( ig, is ) = ( 0.0D0, 0.0D0 ) 1055 ENDDO 1056 ENDDO 1057 1058 CALL enforce_input_dft ( input_dft ) 1059#ifdef EXX 1060 IF ( exx_flag ) CALL start_exx ( ) 1061#endif 1062 vxcr_g ( :, : ) = 0.0D0 1063 IF ( vxc_zero_rho_core ) THEN 1064 rho_core ( : ) = 0.0D0 1065 rhog_core ( : ) = ( 0.0D0, 0.0D0 ) 1066 ENDIF 1067 CALL v_xc ( rho, rho_core, rhog_core, etxc, vtxc, vxcr_g ) 1068#ifdef EXX 1069 IF ( exx_flag ) CALL stop_exx ( ) 1070#endif 1071 DO is = 1, ns 1072 DO ir = 1, nr 1073 psic ( ir ) = CMPLX ( vxcr_g ( ir, is ), 0.0D0 ) 1074 ENDDO 1075 CALL fwfft ( 'Dense', psic, dfftp ) 1076 DO ig = 1, ng_l 1077 vxcg_g ( ig_l2g ( ig ), is ) = psic ( nl ( ig ) ) 1078 ENDDO 1079 ENDDO 1080 1081 CALL mp_sum ( vxcg_g, intra_pool_comm ) 1082 1083 IF ( ionode ) OPEN ( unit = unit, file = TRIM ( output_file_name ), & 1084 form = 'unformatted', status = 'replace' ) 1085 1086 CALL write_header(unit, stitle) 1087 1088 IF ( ionode ) THEN 1089 WRITE ( unit ) nrecord 1090 WRITE ( unit ) ng_g 1091 IF ( real_or_complex .EQ. 1 ) THEN 1092 WRITE ( unit ) ( ( dble ( vxcg_g ( ig, is ) ), & 1093 ig = 1, ng_g ), is = 1, ns ) 1094 ELSE 1095 WRITE ( unit ) ( ( vxcg_g ( ig, is ), ig = 1, ng_g ), is = 1, nsf ) 1096 ENDIF 1097 CLOSE ( unit = unit, status = 'keep' ) 1098 ENDIF 1099 1100 DEALLOCATE ( vxcg_g ) 1101 DEALLOCATE ( vxcr_g ) 1102 1103 RETURN 1104 1105END SUBROUTINE write_vxcg 1106 1107!------------------------------------------------------------------------------- 1108 1109SUBROUTINE write_vxc0 ( output_file_name, vxc_zero_rho_core, input_dft, & 1110 exx_flag ) 1111 1112 USE constants, ONLY : RYTOEV 1113 USE ener, ONLY : etxc, vtxc 1114 USE fft_base, ONLY : dfftp 1115 USE fft_interfaces, ONLY : fwfft 1116 USE funct, ONLY : enforce_input_dft 1117 USE grid_dimensions, ONLY : nrxx 1118 USE gvect, ONLY : ngm, nl, mill 1119 USE io_global, ONLY : ionode 1120 USE kinds, ONLY : DP 1121 USE lsda_mod, ONLY : nspin 1122 USE mp, ONLY : mp_sum 1123 USE mp_global, ONLY : intra_pool_comm 1124 USE scf, ONLY : rho, rho_core, rhog_core 1125 USE wavefunctions_module, ONLY : psic 1126#ifdef EXX 1127 USE funct, ONLY : start_exx, stop_exx 1128#endif 1129 1130 IMPLICIT NONE 1131 1132 character ( len = 256 ), intent (in) :: output_file_name 1133 logical, intent (in) :: vxc_zero_rho_core 1134 character ( len = 20 ), intent (in) :: input_dft 1135 logical, intent (in) :: exx_flag 1136 1137 integer :: unit 1138 integer :: is, ir, ig 1139 integer :: ns, nst, nsf, nr, ng_l 1140 real (DP), allocatable :: vxcr_g ( :, : ) 1141 complex (DP), allocatable :: vxc0_g ( : ) 1142 1143 unit = 4 1144 1145 nr = nrxx 1146 ng_l = ngm 1147 call set_spin(ns,nst,nsf,nspin) 1148 1149 ALLOCATE ( vxcr_g ( nr, nsf ) ) 1150 ALLOCATE ( vxc0_g ( nsf ) ) 1151 1152 DO is = 1, nsf 1153 vxc0_g ( is ) = ( 0.0D0, 0.0D0 ) 1154 ENDDO 1155 1156 CALL enforce_input_dft ( input_dft ) 1157#ifdef EXX 1158 IF ( exx_flag ) CALL start_exx ( ) 1159#endif 1160 vxcr_g ( :, : ) = 0.0D0 1161 IF ( vxc_zero_rho_core ) THEN 1162 rho_core ( : ) = 0.0D0 1163 rhog_core ( : ) = ( 0.0D0, 0.0D0 ) 1164 ENDIF 1165 CALL v_xc ( rho, rho_core, rhog_core, etxc, vtxc, vxcr_g ) 1166#ifdef EXX 1167 IF ( exx_flag ) CALL stop_exx ( ) 1168#endif 1169 DO is = 1, nsf 1170 DO ir = 1, nr 1171 psic ( ir ) = CMPLX ( vxcr_g ( ir, is ), 0.0D0 ) 1172 ENDDO 1173 CALL fwfft ( 'Dense', psic, dfftp ) 1174 DO ig = 1, ng_l 1175 IF ( mill ( 1, ig ) .EQ. 0 .AND. mill ( 2, ig ) .EQ. 0 .AND. & 1176 mill ( 3, ig ) .EQ. 0 ) vxc0_g ( is ) = psic ( nl ( ig ) ) 1177 ENDDO 1178 ENDDO 1179 1180 CALL mp_sum ( vxc0_g, intra_pool_comm ) 1181 1182 DO is = 1, nsf 1183 vxc0_g ( is ) = vxc0_g ( is ) * CMPLX ( RYTOEV, 0.0D0 ) 1184 ENDDO 1185 1186 IF ( ionode ) THEN 1187 OPEN (unit = unit, file = TRIM (output_file_name), & 1188 form = 'formatted', status = 'replace') 1189 WRITE ( unit, 101 ) 1190 DO is = 1, nsf 1191 WRITE ( unit, 102 ) is, vxc0_g ( is ) 1192 ENDDO 1193 WRITE ( unit, 103 ) 1194 CLOSE (unit = unit, status = 'keep') 1195 ENDIF 1196 1197 DEALLOCATE ( vxcr_g ) 1198 DEALLOCATE ( vxc0_g ) 1199 1200 RETURN 1201 1202101 FORMAT ( /, 5X, "--------------------------------------------", & 1203 /, 5X, "spin Re Vxc(G=0) (eV) Im Vxc(G=0) (eV)", & 1204 /, 5X, "--------------------------------------------" ) 1205102 FORMAT ( 5X, I1, 3X, 2F20.15 ) 1206103 FORMAT ( 5X, "--------------------------------------------", / ) 1207 1208END SUBROUTINE write_vxc0 1209 1210!------------------------------------------------------------------------------- 1211 1212SUBROUTINE write_vxc (vxc_integral, output_file_name, diag_nmin, diag_nmax, & 1213 offdiag_nmin, offdiag_nmax, vxc_zero_rho_core, input_dft, exx_flag, & 1214 dfpt_type, dfpt_mode) 1215 1216 USE becmod, ONLY : bec_type, allocate_bec_type, deallocate_bec_type 1217 USE constants, ONLY : rytoev, eps6 1218 USE control_ph, ONLY : nbnd_occ ! initialize! 1219 USE cell_base, ONLY : tpiba2, at, bg 1220 USE ener, ONLY : etxc, vtxc 1221 USE eqv, ONLY : dpsi, dvpsi 1222 USE fft_base, ONLY : dfftp 1223 USE fft_interfaces, ONLY : fwfft, invfft 1224 USE funct, ONLY : enforce_input_dft 1225 USE grid_dimensions, ONLY : nr1x, nr2x, nr3x, nrxx 1226 USE gvect, ONLY : ngm, g, nl 1227 USE io_files, ONLY : nwordwfc, iunwfc 1228 USE io_global, ONLY : ionode 1229 USE kinds, ONLY : DP 1230 USE klist, ONLY : xk, nkstot 1231 USE lsda_mod, ONLY : nspin, isk 1232 USE modes, ONLY : u 1233 USE mp, ONLY : mp_sum 1234 USE mp_global, ONLY : intra_pool_comm, inter_pool_comm 1235 USE phus, ONLY : becp1 1236 USE units_ph, ONLY : iudvscf, lrdrho 1237 USE scf, ONLY : rho, rho_core, rhog_core 1238 USE uspp, ONLY : nkb 1239 USE wavefunctions_module, ONLY : evc, psic 1240 USE wvfct, ONLY : npwx, npw, nbnd, igk, g2kin, ecutwfc 1241#ifdef EXX 1242 USE funct, ONLY : start_exx, stop_exx, exx_is_active 1243 USE exx, ONLY : vexx 1244#endif 1245 1246 IMPLICIT NONE 1247 1248 character (len = 1), intent (in) :: vxc_integral ! 'g' = G-space, 'r' = R-space 1249 character (len = 256), intent (in) :: output_file_name 1250 integer, intent (inout) :: diag_nmin 1251 integer, intent (inout) :: diag_nmax 1252 integer, intent (inout) :: offdiag_nmin 1253 integer, intent (inout) :: offdiag_nmax 1254 logical, intent (in) :: vxc_zero_rho_core 1255 character (len = 20), intent (in) :: input_dft 1256 logical, intent (in) :: exx_flag 1257 integer, intent (in) :: dfpt_type ! 0 = none, 1 = ionic, 2 = electric 1258 integer, intent (in) :: dfpt_mode 1259 1260 integer :: ik, is, ib, ig, ir, unit, iks, ike, ndiag, noffdiag, ib2 1261 integer :: ns, nst, nsf 1262 real(DP) :: dummyr, maxvaldvscf 1263 complex (DP) :: dummy 1264 complex (DP), allocatable :: mtxeld (:, :) 1265 complex (DP), allocatable :: mtxelo (:, :, :) 1266 real (DP), allocatable :: vxcr (:, :) 1267 complex (DP), allocatable :: dvscf (:, :) 1268 complex (DP), allocatable :: psic2 (:) 1269 complex (DP), allocatable :: hpsi (:) 1270 type(bec_type) :: becp2 ! we do not care, but need to pass this 1271 1272 IF ( vxc_integral /= 'r' .AND. vxc_integral /= 'g' ) & 1273 CALL errore ( 'write_vxc', 'vxc_integral', 1 ) 1274 1275 IF ( nspin .EQ. 4) CALL errore ( 'write_vxc', 'nspin', 4 ) 1276 1277 if(diag_nmin > diag_nmax) then 1278 call errore ( 'write_vxc', 'diag_nmin > diag_nmax', diag_nmin ) 1279 endif 1280 IF (diag_nmin .LT. 1) diag_nmin = 1 1281 IF (diag_nmax .GT. nbnd) then 1282 if(ionode) write(0,'(a,i6)') 'WARNING: resetting diag_nmax to max number of bands', nbnd 1283 diag_nmax = nbnd 1284 ENDIF 1285 1286 if(offdiag_nmin > offdiag_nmax) then 1287 call errore ( 'write_vxc', 'offdiag_nmin > offdiag_nmax', offdiag_nmin ) 1288 endif 1289 ndiag = MAX (diag_nmax - diag_nmin + 1, 0) 1290 IF (offdiag_nmin .LT. 1) offdiag_nmin = 1 1291 IF (offdiag_nmax .GT. nbnd) then 1292 if(ionode) write(0,'(a,i6)') 'WARNING: resetting offdiag_nmax to max number of bands', nbnd 1293 offdiag_nmax = nbnd 1294 ENDIF 1295 1296 noffdiag = MAX (offdiag_nmax - offdiag_nmin + 1, 0) 1297 IF (ndiag .EQ. 0 .AND. noffdiag .EQ. 0) RETURN 1298 1299 unit = 4 1300 1301 call set_spin(ns, nst, nsf, nspin) 1302 CALL k_pools(iks, ike) 1303 1304 IF (ndiag .GT. 0) THEN 1305 ALLOCATE (mtxeld (ndiag, nkstot)) 1306 mtxeld (:, :) = (0.0D0, 0.0D0) 1307 ENDIF 1308 IF (noffdiag .GT. 0) THEN 1309 ALLOCATE (mtxelo (noffdiag, noffdiag, nkstot)) 1310 mtxelo (:, :, :) = (0.0D0, 0.0D0) 1311 ENDIF 1312 1313 ALLOCATE (vxcr (nrxx, nsf)) 1314 IF (noffdiag .GT. 0) ALLOCATE (psic2 (nrxx)) 1315 if(vxc_integral == 'g') ALLOCATE (hpsi (nrxx)) 1316 1317 if(dfpt_type == 0) then 1318 1319 CALL enforce_input_dft (input_dft) 1320#ifdef EXX 1321 IF (exx_flag) CALL start_exx () 1322#endif 1323 1324 vxcr (:, :) = 0.0D0 1325 IF ( vxc_zero_rho_core ) THEN 1326 rho_core ( : ) = 0.0D0 1327 rhog_core ( : ) = ( 0.0D0, 0.0D0 ) 1328 ENDIF 1329 CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxcr) 1330 1331 else 1332 1333 IF(dfpt_type == 1) WRITE(6,'(a,i6)') 'Calculating ionic matrix elements for mode ', dfpt_mode 1334 IF(dfpt_type == 2) WRITE(6,'(a,i6)') 'Calculating electric matrix elements for polarization ', dfpt_mode 1335 IF(vxc_integral == 'r') WRITE(6,'(a)') 'Only the Hartree and XC parts of the perturbation are being included.' 1336 1337 ! iudvscf, lrdrho set in check_dfpt_modes 1338 ALLOCATE (dvscf (nrxx, nsf)) 1339 CALL davcio_drho(dvscf, lrdrho, iudvscf, dfpt_mode, - 1 ) 1340 write(6,'(a)') 'Read dvscf from disk.' 1341 1342 ! dvscf is real at q=0 1343 if(any(aimag(dvscf(:,:)) > eps6)) then 1344 write(6,'(a)') 'This mode has imaginary dvscf (max = ', maxval(abs(aimag(dvscf(:,:)))), ') so not calculating.' 1345 RETURN 1346 endif 1347 vxcr(:,:) = dble(dvscf(:,:)) 1348 1349 maxvaldvscf = maxval(abs(dvscf(:,:))) 1350 DEALLOCATE(dvscf) 1351 1352 if(maxvaldvscf < 1e-9) then 1353 WRITE(6,'(a)') 'This mode has dvscf = 0 (max abs = ', maxvaldvscf, ') so not calculating.' 1354 RETURN 1355 endif 1356 1357 if(dfpt_type == 2) call allocate_bec_type (nkb, nbnd, becp2) 1358 endif 1359 1360 DO ik = iks, ike 1361 CALL gk_sort (xk (1, ik - iks + 1), ngm, g, ecutwfc & 1362 / tpiba2, npw, igk, g2kin) 1363 CALL davcio (evc, nwordwfc, iunwfc, ik - iks + 1, -1) 1364 1365 if(vxc_integral == 'g') then 1366 if(dfpt_type == 1) then 1367 ! initialize dvpsi: bare ionic pert applied to all the wfns 'evc' 1368 ! gradient of local and non-local part of pseudopotential 1369 dvpsi(:,:) = (0d0, 0d0) 1370 call dvqpsi_us(ik, u (1, dfpt_mode), .false.) 1371 else if(dfpt_type == 2) then 1372 ! for electric field, from PH/dvpsi_e. check initialization of this stuff! 1373 ! calculate the commutator [H,x_ipol] psi > and store it in dvpsi 1374 ! dpsi used as workspace 1375 dpsi(:,:) = (0d0, 0d0) 1376 dvpsi(:,:) = (0d0, 0d0) 1377 call commutator_Hx_psi (ik, nbnd_occ(ik), becp1(ik), becp2, dfpt_mode, dvpsi, dpsi ) 1378 endif 1379 endif 1380 1381 IF (ndiag .GT. 0) THEN 1382 DO ib = diag_nmin, diag_nmax 1383 psic (:) = (0.0D0, 0.0D0) 1384 DO ig = 1, npw 1385 psic (nl (igk (ig))) = evc (ig, ib) 1386 ENDDO 1387 CALL invfft ('Dense', psic, dfftp) 1388 1389 if(vxc_integral == 'g') then 1390 DO ir = 1, nrxx 1391 psic (ir) = psic (ir) * vxcr (ir, isk (ik)) 1392 ENDDO 1393 CALL fwfft ('Dense', psic, dfftp) 1394 1395 if(dfpt_type == 1) then 1396 ! add bare part of ionic perturbation 1397 DO ig = 1, npw 1398 psic (nl (igk(ig))) = psic (nl (igk(ig))) + dvpsi(ig, ib) 1399 ENDDO 1400 endif 1401 1402 hpsi (:) = (0.0D0, 0.0D0) 1403 DO ig = 1, npw 1404 hpsi (ig) = psic (nl (igk (ig))) 1405 ENDDO 1406 psic (:) = (0.0D0, 0.0D0) 1407 DO ig = 1, npw 1408 psic (ig) = evc (ig, ib) 1409 ENDDO 1410#ifdef EXX 1411 IF (exx_is_active () .and. dfpt_type == 0) & 1412 CALL vexx (npwx, npw, 1, psic, hpsi) 1413#endif 1414 dummy = (0.0D0, 0.0D0) 1415 DO ig = 1, npw 1416 dummy = dummy + conjg (psic (ig)) * hpsi (ig) 1417 ENDDO 1418 dummy = dummy * rytoev 1419 CALL mp_sum (dummy, intra_pool_comm) 1420 mtxeld (ib - diag_nmin + 1, ik) = dummy 1421 else 1422 dummyr = 0.0D0 1423 DO ir = 1, nrxx 1424 dummyr = dummyr + vxcr (ir, isk (ik)) & 1425 * (dble (psic (ir)) **2 + aimag (psic (ir)) **2) 1426 ENDDO 1427 dummyr = dummyr * rytoev / dble (nr1x * nr2x * nr3x) 1428 CALL mp_sum (dummyr, intra_pool_comm) 1429 mtxeld (ib - diag_nmin + 1, ik) = dummyr 1430 endif 1431 ENDDO 1432 ENDIF 1433 IF (noffdiag .GT. 0) THEN 1434 DO ib = offdiag_nmin, offdiag_nmax 1435 psic (:) = (0.0D0, 0.0D0) 1436 DO ig = 1, npw 1437 psic (nl (igk (ig))) = evc (ig, ib) 1438 ENDDO 1439 CALL invfft ('Dense', psic, dfftp) 1440 1441 if(vxc_integral == 'g') then 1442 DO ir = 1, nrxx 1443 psic (ir) = psic (ir) * vxcr (ir, isk (ik)) 1444 ENDDO 1445 CALL fwfft ('Dense', psic, dfftp) 1446 1447 if(dfpt_type == 1) then 1448 ! add bare part of ionic perturbation 1449 DO ig = 1, npw 1450 psic (nl (igk(ig))) = psic (nl (igk(ig))) + dvpsi(ig, ib) 1451 ENDDO 1452 endif 1453 1454 hpsi (:) = (0.0D0, 0.0D0) 1455 DO ig = 1, npw 1456 hpsi (ig) = psic (nl (igk (ig))) 1457 ENDDO 1458 psic (:) = (0.0D0, 0.0D0) 1459 DO ig = 1, npw 1460 psic (ig) = evc (ig, ib) 1461 ENDDO 1462#ifdef EXX 1463 IF (exx_is_active () .and. dfpt_type == 0) & 1464 CALL vexx (npwx, npw, 1, psic, hpsi) 1465#endif 1466 endif 1467 DO ib2 = offdiag_nmin, offdiag_nmax 1468 psic2 (:) = (0.0D0, 0.0D0) 1469 dummy = (0.0D0, 0.0D0) 1470 if(vxc_integral == 'g') then 1471 DO ig = 1, npw 1472 psic2 (ig) = evc (ig, ib2) 1473 ENDDO 1474 DO ig = 1, npw 1475 dummy = dummy + conjg (psic2 (ig)) * hpsi (ig) 1476 ENDDO 1477 else 1478 DO ig = 1, npw 1479 psic2 (nl (igk (ig))) = evc (ig, ib2) 1480 ENDDO 1481 CALL invfft ('Dense', psic2, dfftp) 1482 DO ir = 1, nrxx 1483 dummy = dummy + CMPLX (vxcr (ir, isk (ik)), 0.0D0) & 1484 * conjg (psic2 (ir)) * psic (ir) 1485 ENDDO 1486 dummy = dummy / dble (nr1x * nr2x * nr3x) 1487 endif 1488 dummy = dummy * rytoev 1489 CALL mp_sum (dummy, intra_pool_comm) 1490 mtxelo (ib2 - offdiag_nmin + 1, ib - offdiag_nmin + 1, ik) = dummy 1491 ENDDO 1492 ENDDO 1493 ENDIF 1494 ENDDO 1495 1496 ! from dvpsi_e 1497 if(dfpt_type == 2 .and. nkb > 0) call deallocate_bec_type (becp2) 1498 1499#ifdef EXX 1500 IF (exx_flag) CALL stop_exx () 1501#endif 1502 1503 DEALLOCATE (vxcr) 1504 IF (noffdiag .GT. 0) DEALLOCATE (psic2) 1505 if(vxc_integral == 'g') DEALLOCATE (hpsi) 1506 1507 IF (ndiag .GT. 0) CALL mp_sum (mtxeld, inter_pool_comm) 1508 IF (noffdiag .GT. 0) CALL mp_sum (mtxelo, inter_pool_comm) 1509 1510 CALL cryst_to_cart (nkstot, xk, at, -1) 1511 1512 IF (ionode) THEN 1513 OPEN (unit = unit, file = TRIM (output_file_name), & 1514 form = 'formatted', status = 'replace') 1515 DO ik = 1, nkstot / ns 1516 WRITE (unit, 101) xk(:, ik), ns * ndiag, & 1517 ns * noffdiag **2 1518 DO is = 1, ns 1519 IF (ndiag .GT. 0) THEN 1520 DO ib = diag_nmin, diag_nmax 1521 WRITE (unit, 102) is, ib, mtxeld & 1522 (ib - diag_nmin + 1, ik + (is - 1) * nkstot / ns) 1523 ENDDO 1524 ENDIF 1525 IF (noffdiag .GT. 0) THEN 1526 DO ib = offdiag_nmin, offdiag_nmax 1527 DO ib2 = offdiag_nmin, offdiag_nmax 1528 WRITE (unit, 103) is, ib2, ib, mtxelo & 1529 (ib2 - offdiag_nmin + 1, ib - offdiag_nmin + 1, & 1530 ik + (is - 1) * nkstot / ns) 1531 ENDDO 1532 ENDDO 1533 ENDIF 1534 ENDDO 1535 ENDDO 1536 CLOSE (unit = unit, status = 'keep') 1537 ENDIF 1538 1539 CALL cryst_to_cart (nkstot, xk, bg, 1) 1540 1541 IF (ndiag .GT. 0) DEALLOCATE (mtxeld) 1542 IF (noffdiag .GT. 0) DEALLOCATE (mtxelo) 1543 1544 RETURN 1545 1546 101 FORMAT (3F13.9, 2I8) 1547 102 FORMAT (2I8, 2F15.9) 1548 103 FORMAT (3I8, 2F15.9) 1549 1550END SUBROUTINE write_vxc 1551 1552!------------------------------------------------------------------------------- 1553 1554SUBROUTINE write_vnlg (output_file_name, wfng_kgrid, wfng_nk1, wfng_nk2, & 1555 wfng_nk3, wfng_dk1, wfng_dk2, wfng_dk3) 1556 1557 USE cell_base, ONLY : omega, alat, tpiba, tpiba2, at, bg, & 1558 ibrav, symm_type 1559 USE constants, ONLY : pi, tpi, eps6 1560 USE grid_dimensions, ONLY : nr1, nr2, nr3 1561 USE gvect, ONLY : ngm, ngm_g, ig_l2g, g, mill, ecutrho 1562 USE io_global, ONLY : ionode, ionode_id 1563 USE ions_base, ONLY : nat, atm, ityp, tau, nsp 1564 USE kinds, ONLY : DP 1565 USE klist, ONLY : xk, wk, ngk, nks, nkstot 1566 USE lsda_mod, ONLY : nspin, isk 1567 USE mp, ONLY : mp_sum, mp_max, mp_get, mp_barrier 1568 USE mp_global, ONLY : mpime, nproc, world_comm, me_pool, & 1569 root_pool, npool, nproc_pool, intra_pool_comm 1570 USE mp_wave, ONLY : mergewf 1571 USE start_k, ONLY : nk1, nk2, nk3, k1, k2, k3 1572 USE symm_base, ONLY : s, ftau, nsym 1573 USE uspp, ONLY : nkb, vkb, deeq 1574 USE uspp_param, ONLY : nhm, nh 1575 USE wvfct, ONLY : npwx, npw, g2kin, ecutwfc 1576 1577 IMPLICIT NONE 1578 1579 character (len = 256), intent (in) :: output_file_name 1580 logical, intent (in) :: wfng_kgrid 1581 integer, intent (in) :: wfng_nk1 1582 integer, intent (in) :: wfng_nk2 1583 integer, intent (in) :: wfng_nk3 1584 real (DP), intent (in) :: wfng_dk1 1585 real (DP), intent (in) :: wfng_dk2 1586 real (DP), intent (in) :: wfng_dk3 1587 1588 character :: cdate*9, ctime*9, sdate*32, stime*32, stitle*32 1589 integer :: i, j, k, ierr, ik, is, ig, ikb, iat, isp, ih, jh, & 1590 unit, iks, ike, npw_g, npwx_g, ngg, ipsour, & 1591 igwx, local_pw, id, nd, ntran, cell_symmetry, nrecord 1592 real (DP) :: alat2, recvol, dr1, t1 ( 3 ), t2 ( 3 ) 1593 real (DP) :: r1 ( 3, 3 ), r2 ( 3, 3 ), adot ( 3, 3 ) 1594 real (DP) :: bdot ( 3, 3 ), translation ( 3, 48 ) 1595 integer, allocatable :: kmap ( : ) 1596 integer, allocatable :: smap ( : ) 1597 integer, allocatable :: gvec ( :, : ) 1598 integer, allocatable :: ngk_g ( : ) 1599 integer, allocatable :: igk_l2g ( :, : ) 1600 integer, allocatable :: itmp ( : ) 1601 integer, allocatable :: igwk ( : ) 1602 integer, allocatable :: igwf_l2g ( : ) 1603 integer, allocatable :: ipmask ( : ) 1604 complex (DP), allocatable :: vkb_g ( : ) 1605 1606 INTEGER, EXTERNAL :: atomic_number 1607 1608 IF ( nkb == 0 ) RETURN 1609 1610 CALL date_and_tim ( cdate, ctime ) 1611 WRITE ( sdate, '(A2,"-",A3,"-",A4,21X)' ) cdate(1:2), cdate(3:5), cdate(6:9) 1612 WRITE ( stime, '(A8,24X)' ) ctime(1:8) 1613 WRITE ( stitle, '("VKB-Complex",21X)' ) 1614 1615 IF ( nspin .EQ. 4) CALL errore ( 'write_vnlg', 'nspin', 4) 1616 1617 unit = 4 1618 nrecord = 1 1619 nd = 3 1620 1621 CALL k_pools(iks, ike) 1622 1623 ierr = 0 1624 IF ( ibrav .EQ. 0 ) THEN 1625 IF ( TRIM ( symm_type ) .EQ. 'cubic' ) THEN 1626 cell_symmetry = 0 1627 ELSEIF ( TRIM ( symm_type ) .EQ. 'hexagonal' ) THEN 1628 cell_symmetry = 1 1629 ELSE 1630 ierr = 1 1631 ENDIF 1632 ELSEIF ( ibrav .GE. 1 .AND. ibrav .LE. 3 ) THEN 1633 cell_symmetry = 0 1634 ELSEIF ( ibrav .GE. 4 .AND. ibrav .LE. 5 ) THEN 1635 cell_symmetry = 1 1636 ELSEIF ( ibrav .GE. 6 .AND. ibrav .LE. 14 ) THEN 1637 cell_symmetry = 0 1638 ELSE 1639 ierr = 1 1640 ENDIF 1641 IF ( ierr .GT. 0 ) & 1642 CALL errore ( 'write_vnlg', 'cell_symmetry', ierr ) 1643 1644 ntran = nsym 1645 DO i = 1, ntran 1646 DO j = 1, nd 1647 DO k = 1, nd 1648 r1 ( k, j ) = dble ( s ( k, j, i ) ) 1649 ENDDO 1650 ENDDO 1651 CALL invmat ( 3, r1, r2, dr1 ) 1652 t1 ( 1 ) = dble ( ftau ( 1, i ) ) / dble ( nr1 ) 1653 t1 ( 2 ) = dble ( ftau ( 2, i ) ) / dble ( nr2 ) 1654 t1 ( 3 ) = dble ( ftau ( 3, i ) ) / dble ( nr3 ) 1655 DO j = 1, nd 1656 t2 ( j ) = 0.0D0 1657 DO k = 1, nd 1658 t2 ( j ) = t2 ( j ) + r2 ( k, j ) * t1 ( k ) 1659 ENDDO 1660 IF ( t2 ( j ) .GE. eps6 + 0.5D0 ) & 1661 t2 ( j ) = t2 ( j ) - dble ( int ( t2 ( j ) + 0.5D0 ) ) 1662 IF ( t2 ( j ) .LT. eps6 - 0.5D0 ) & 1663 t2 ( j ) = t2 ( j ) - dble ( int ( t2 ( j ) - 0.5D0 ) ) 1664 ENDDO 1665 DO j = 1, nd 1666 translation ( j, i ) = t2 ( j ) * tpi 1667 ENDDO 1668 ENDDO 1669 1670 alat2 = alat ** 2 1671 recvol = 8.0D0 * pi**3 / omega 1672 1673 DO i = 1, nd 1674 DO j = 1, nd 1675 adot ( j, i ) = 0.0D0 1676 ENDDO 1677 ENDDO 1678 DO i = 1, nd 1679 DO j = 1, nd 1680 DO k = 1, nd 1681 adot ( j, i ) = adot ( j, i ) + & 1682 at ( k, j ) * at ( k, i ) * alat2 1683 ENDDO 1684 ENDDO 1685 ENDDO 1686 1687 DO i = 1, nd 1688 DO j = 1, nd 1689 bdot ( j, i ) = 0.0D0 1690 ENDDO 1691 ENDDO 1692 DO i = 1, nd 1693 DO j = 1, nd 1694 DO k = 1, nd 1695 bdot ( j, i ) = bdot ( j, i ) + & 1696 bg ( k, j ) * bg ( k, i ) * tpiba2 1697 ENDDO 1698 ENDDO 1699 ENDDO 1700 1701 ALLOCATE ( kmap ( nkstot ) ) 1702 ALLOCATE ( smap ( nkstot ) ) 1703 1704 DO i = 1, nkstot 1705 j = ( i - 1 ) / nspin 1706 k = i - 1 - j * nspin 1707 kmap ( i ) = j + k * ( nkstot / nspin ) + 1 1708 smap ( i ) = k + 1 1709 ENDDO 1710 ierr = 0 1711 DO i = 1, nkstot 1712 ik = kmap ( i ) 1713 is = smap ( i ) 1714 IF ( ik .GE. iks .AND. ik .LE. ike .AND. is .NE. isk ( ik ) ) & 1715 ierr = ierr + 1 1716 ENDDO 1717 CALL mp_max ( ierr ) 1718 IF ( ierr .GT. 0 ) & 1719 CALL errore ( 'write_vnlg', 'smap', ierr ) 1720 1721 ALLOCATE ( gvec ( 3, ngm_g ) ) 1722 gvec = 0 1723 DO ig = 1, ngm 1724 gvec ( 1, ig_l2g ( ig ) ) = mill ( 1, ig ) 1725 gvec ( 2, ig_l2g ( ig ) ) = mill ( 2, ig ) 1726 gvec ( 3, ig_l2g ( ig ) ) = mill ( 3, ig ) 1727 ENDDO 1728 CALL mp_sum ( gvec, intra_pool_comm ) 1729 1730 ALLOCATE ( ngk_g ( nkstot ) ) 1731 ALLOCATE ( igk_l2g ( npwx, nks ) ) 1732 ngk_g = 0 1733 igk_l2g = 0 1734 ALLOCATE ( itmp ( npwx ) ) 1735 DO ik = 1, nks 1736 itmp = 0 1737 npw = npwx 1738 CALL gk_sort ( xk ( 1, ik + iks - 1 ), ngm, g, ecutwfc / tpiba2, & 1739 npw, itmp ( 1 ), g2kin ) 1740 DO ig = 1, npw 1741 igk_l2g ( ig, ik ) = ig_l2g ( itmp ( ig ) ) 1742 ENDDO 1743 ngk ( ik ) = npw 1744 ENDDO 1745 DEALLOCATE ( itmp ) 1746 DO ik = 1, nks 1747 ngk_g ( ik + iks - 1 ) = ngk ( ik ) 1748 ENDDO 1749 CALL mp_sum ( ngk_g ) 1750 npw_g = MAXVAL ( igk_l2g ( :, : ) ) 1751 CALL mp_max ( npw_g ) 1752 npwx_g = MAXVAL ( ngk_g ( : ) ) 1753 1754 CALL cryst_to_cart (nkstot, xk, at, -1) 1755 1756 IF (ionode) THEN 1757 OPEN (unit = unit, file = TRIM (output_file_name), & 1758 form = 'unformatted', status = 'replace') 1759 WRITE ( unit ) stitle, sdate, stime 1760 WRITE ( unit ) nspin, ngm_g, ntran, cell_symmetry, nat, ecutrho, & 1761 nkstot / nspin, nsp, nkb, nhm, npwx_g, ecutwfc 1762 IF ( wfng_kgrid ) THEN 1763 WRITE ( unit ) nr1, nr2, nr3, wfng_nk1, wfng_nk2, wfng_nk3, & 1764 wfng_dk1, wfng_dk2, wfng_dk3 1765 ELSE 1766 WRITE ( unit ) nr1, nr2, nr3, nk1, nk2, nk3, & 1767 0.5D0 * dble ( k1 ), 0.5D0 * dble ( k2 ), 0.5D0 * dble ( k3 ) 1768 ENDIF 1769 WRITE ( unit ) omega, alat, ( ( at ( j, i ), j = 1, nd ), i = 1, nd ), & 1770 ( ( adot ( j, i ), j = 1, nd ), i = 1, nd ) 1771 WRITE ( unit ) recvol, tpiba, ( ( bg ( j, i ), j = 1, nd ), i = 1, nd ), & 1772 ( ( bdot ( j, i ), j = 1, nd ), i = 1, nd ) 1773 WRITE ( unit ) ( ( ( s ( k, j, i ), k = 1, nd ), j = 1, nd ), i = 1, ntran ) 1774 WRITE ( unit ) ( ( translation ( j, i ), j = 1, nd ), i = 1, ntran ) 1775 WRITE ( unit ) ( ( tau ( j, i ), j = 1, nd ), atomic_number ( atm ( ityp ( i ) ) ), i = 1, nat ) 1776 WRITE ( unit ) ( ngk_g ( ik ), ik = 1, nkstot / nspin ) 1777 WRITE ( unit ) ( wk ( ik ) * dble ( nspin ) / 2.0D0, ik = 1, nkstot / nspin ) 1778 WRITE ( unit ) ( ( xk ( id, ik ), id = 1, nd ), ik = 1, nkstot / nspin ) 1779 WRITE ( unit ) ( ityp ( iat ), iat = 1, nat ) 1780 WRITE ( unit ) ( nh ( isp ), isp = 1, nsp ) 1781 WRITE ( unit ) ( ( ( ( deeq ( jh, ih, iat, is ), & 1782 jh = 1, nhm ), ih = 1, nhm ), iat = 1, nat ), is = 1, nspin ) 1783 WRITE ( unit ) nrecord 1784 WRITE ( unit ) ngm_g 1785 WRITE ( unit ) ( ( gvec ( id, ig ), id = 1, nd ), ig = 1, ngm_g ) 1786 ENDIF 1787 1788 CALL cryst_to_cart (nkstot, xk, bg, 1) 1789 1790 ALLOCATE ( igwk ( npwx_g ) ) 1791 1792 DO i = 1, nkstot 1793 1794 ik = kmap ( i ) 1795 is = smap ( i ) 1796 1797 igwk = 0 1798 1799 ALLOCATE ( itmp ( npw_g ) ) 1800 itmp = 0 1801 IF ( ik .GE. iks .AND. ik .LE. ike ) THEN 1802 DO ig = 1, ngk ( ik - iks + 1 ) 1803 itmp ( igk_l2g ( ig, ik - iks + 1 ) ) = igk_l2g ( ig, ik - iks + 1 ) 1804 ENDDO 1805 ENDIF 1806 CALL mp_sum ( itmp ) 1807 ngg = 0 1808 DO ig = 1, npw_g 1809 IF ( itmp ( ig ) .EQ. ig ) THEN 1810 ngg = ngg + 1 1811 igwk ( ngg ) = ig 1812 ENDIF 1813 ENDDO 1814 DEALLOCATE ( itmp ) 1815 1816 IF ( ionode ) THEN 1817 IF ( is .EQ. 1 ) THEN 1818 WRITE ( unit ) nrecord 1819 WRITE ( unit ) ngk_g ( ik ) 1820 WRITE ( unit ) ( ( gvec ( j, igwk ( ig ) ), j = 1, 3 ), & 1821 ig = 1, ngk_g ( ik ) ) 1822 ENDIF 1823 ENDIF 1824 1825 local_pw = 0 1826 IF ( ik .GE. iks .AND. ik .LE. ike ) THEN 1827 ALLOCATE ( itmp ( npwx ) ) 1828 npw = npwx 1829 CALL gk_sort ( xk ( 1, ik ), ngm, g, ecutwfc / tpiba2, & 1830 npw, itmp ( 1 ), g2kin ) 1831 CALL init_us_2 ( npw, itmp, xk ( 1, ik ), vkb ) 1832 local_pw = ngk ( ik - iks + 1 ) 1833 DEALLOCATE ( itmp ) 1834 ENDIF 1835 1836 ALLOCATE ( igwf_l2g ( local_pw ) ) 1837 igwf_l2g = 0 1838 DO ig = 1, local_pw 1839 ngg = igk_l2g ( ig, ik - iks + 1 ) 1840 DO j = 1, ngk_g ( ik ) 1841 IF ( ngg .EQ. igwk ( j ) ) THEN 1842 igwf_l2g ( ig ) = j 1843 EXIT 1844 ENDIF 1845 ENDDO 1846 ENDDO 1847 1848 ALLOCATE ( ipmask ( nproc ) ) 1849 ipmask = 0 1850 ipsour = ionode_id 1851 IF ( npool .GT. 1 ) THEN 1852 IF ( ( ik .GE. iks ) .AND. ( ik .LE. ike ) ) THEN 1853 IF ( me_pool .EQ. root_pool ) ipmask ( mpime + 1 ) = 1 1854 ENDIF 1855 CALL mp_sum ( ipmask ) 1856 DO j = 1, nproc 1857 IF ( ipmask ( j ) .EQ. 1 ) ipsour = j - 1 1858 ENDDO 1859 ENDIF 1860 DEALLOCATE ( ipmask ) 1861 1862 igwx = 0 1863 ierr = 0 1864 IF ( ik .GE. iks .AND. ik .LE. ike ) & 1865 igwx = MAXVAL ( igwf_l2g ( 1 : local_pw ) ) 1866 CALL mp_max ( igwx, intra_pool_comm ) 1867 IF ( ipsour .NE. ionode_id ) & 1868 CALL mp_get ( igwx, igwx, mpime, ionode_id, ipsour, 1 ) 1869 ierr = 0 1870 IF ( ik .GE. iks .AND. ik .LE. ike .AND. igwx .NE. ngk_g ( ik ) ) & 1871 ierr = 1 1872 CALL mp_max ( ierr ) 1873 IF ( ierr .GT. 0 ) & 1874 CALL errore ( 'write_vnlg', 'igwx ngk_g', ierr ) 1875 1876 ALLOCATE ( vkb_g ( MAX ( 1, igwx ) ) ) 1877 1878 DO ikb = 1, nkb 1879 1880 vkb_g = ( 0.0D0, 0.0D0 ) 1881 IF ( npool .GT. 1 ) THEN 1882 IF ( ( ik .GE. iks ) .AND. ( ik .LE. ike ) ) THEN 1883 CALL mergewf ( vkb ( :, ikb ), vkb_g, local_pw, igwf_l2g, & 1884 me_pool, nproc_pool, root_pool, intra_pool_comm ) 1885 ENDIF 1886 IF ( ipsour .NE. ionode_id ) THEN 1887 CALL mp_get ( vkb_g, vkb_g, mpime, ionode_id, ipsour, ikb, & 1888 world_comm ) 1889 ENDIF 1890 ELSE 1891 CALL mergewf ( vkb ( :, ikb ), vkb_g, local_pw, igwf_l2g, & 1892 mpime, nproc, ionode_id, world_comm ) 1893 ENDIF 1894 1895 IF ( ionode ) THEN 1896 WRITE ( unit ) nrecord 1897 WRITE ( unit ) igwx 1898 WRITE ( unit ) ( vkb_g ( ig ), ig = 1, igwx ) 1899 ENDIF 1900 1901 ENDDO 1902 1903 DEALLOCATE ( vkb_g ) 1904 DEALLOCATE ( igwf_l2g ) 1905 1906 ENDDO 1907 1908 IF ( ionode ) THEN 1909 CLOSE ( unit = unit, status = 'keep' ) 1910 ENDIF 1911 1912 DEALLOCATE ( igwk ) 1913 DEALLOCATE ( igk_l2g ) 1914 DEALLOCATE ( ngk_g ) 1915 DEALLOCATE ( gvec ) 1916 DEALLOCATE ( smap ) 1917 DEALLOCATE ( kmap ) 1918 1919 RETURN 1920 1921END SUBROUTINE write_vnlg 1922 1923!------------------------------------------------------------------------------- 1924 1925SUBROUTINE check_dfpt_modes (imode, dfpt_type, nmode) 1926 1927 USE control_ph, ONLY : tmp_dir_ph 1928 USE grid_dimensions, ONLY : nr1x, nr2x, nr3x 1929 USE io_files, ONLY : diropn 1930 USE io_global, ONLY : ionode, ionode_id 1931 USE ions_base, ONLY : nat 1932 USE lsda_mod, ONLY : nspin 1933 USE mp, ONLY : mp_bcast 1934 USE units_ph, ONLY : iudvscf, lrdrho 1935 1936 IMPLICIT NONE 1937 1938 integer, intent(in) :: imode 1939 integer, intent(in) :: dfpt_type 1940 integer, intent(out) :: nmode 1941 1942 integer :: statb(13) 1943 logical :: exst 1944 integer*8 :: unf_recl 1945 character*10 :: filename 1946 1947 iudvscf = 27 1948 lrdrho = 2 * nr1x * nr2x * nr3x * nspin ! from PH/openfilq.f90 1949 1950 IF(ionode) THEN 1951 1952 IF(dfpt_type == 1) THEN 1953 filename = "dvscf" 1954 ELSE 1955 filename = "dvscf.E" 1956 ENDIF 1957 1958 CALL diropn (iudvscf, filename, lrdrho, exst, tmp_dir_ = tmp_dir_ph ) 1959 IF(.NOT. exst) CALL errore ('check_dfpt_modes', 'could not open file', 7) 1960 1961 ! from io_files.f90, diropn 1962#if defined(__SX6) 1963# define DIRECT_IO_FACTOR 1 1964#else 1965# define DIRECT_IO_FACTOR 8 1966#endif 1967 1968 ! NOTE: ifort needs -assume byterecl or else nrec will be 4 times too large 1969 unf_recl = DIRECT_IO_FACTOR * int(lrdrho, kind=kind(unf_recl)) 1970 CALL fstat ( iudvscf, statb) ! see EPW, readdvscf.f90 1971 nmode = statb(8) / unf_recl 1972 IF(nmode * unf_recl /= statb(8) .and. ionode) & 1973 WRITE(0,*) 'WARNING: dvscf file does not contain an integral number of records.', & 1974 statb(8) * 1d0 / unf_recl 1975 IF(ionode) WRITE(6,*) 'Number of modes in dvscf file: ', nmode 1976 IF(ionode) WRITE(6,*) 'Starting with mode number: ', imode 1977 IF (imode > nmode) CALL errore('pw2bgw', 'starting mode is not present in dvscf file', iudvscf) 1978 1979 IF(dfpt_type == 1) then 1980 IF(nmode < 3 * nat) WRITE(0,'(a)') 'WARNING: dvscf file does not contain all modes.' 1981 IF(nmode > 3 * nat) WRITE(0,'(a)') 'WARNING: dvscf file contains too many modes.' 1982 ENDIF 1983 IF(dfpt_type == 2) then 1984 IF(nmode < 3) WRITE(0,'(a)') 'WARNING: dvscf_e file does not contain all polarizations.' 1985 IF(nmode > 3) WRITE(0,'(a)') 'WARNING: dvscf_e file contains too many polarizations.' 1986 ENDIF 1987 ENDIF 1988 1989 CALL mp_bcast ( nmode, ionode_id ) 1990 1991 RETURN 1992 1993END SUBROUTINE check_dfpt_modes 1994 1995!------------------------------------------------------------------------------- 1996 1997subroutine check_inversion(real_or_complex, ntran, mtrx, nspin, warn, real_need_inv) 1998 1999! check_inversion Originally By D. Strubbe Last Modified 10/14/2010 2000! Check whether our choice of real/complex version is appropriate given the 2001! presence or absence of inversion symmetry. 2002 2003 USE io_global, ONLY : ionode 2004 2005 implicit none 2006 2007 integer, intent(in) :: real_or_complex 2008 integer, intent(in) :: ntran 2009 integer, intent(in) :: mtrx(3, 3, 48) 2010 integer, intent(in) :: nspin 2011 logical, intent(in) :: warn ! set to false to suppress warnings, for converters 2012 logical, intent(in) :: real_need_inv ! use for generating routines to block real without inversion 2013 2014 integer :: invflag, isym, ii, jj, itest 2015 2016 invflag = 0 2017 do isym = 1, ntran 2018 itest = 0 2019 do ii = 1, 3 2020 do jj = 1, 3 2021 if(ii .eq. jj) then 2022 itest = itest + (mtrx(ii, jj, isym) + 1)**2 2023 else 2024 itest = itest + mtrx(ii, jj, isym)**2 2025 endif 2026 enddo 2027 enddo 2028 if(itest .eq. 0) invflag = invflag + 1 2029 if(invflag .gt. 1) call errore('check_inversion', 'More than one inversion symmetry operation is present.', invflag) 2030 enddo 2031 2032 if(real_or_complex .eq. 2) then 2033 if(invflag .ne. 0 .and. warn .and. nspin == 1) then 2034 if(ionode) write(0, '(a)') 'WARNING: Inversion symmetry is present. The real version would be faster.' 2035 endif 2036 else 2037 if(invflag .eq. 0) then 2038 if(real_need_inv) then 2039 call errore('check_inversion', 'The real version cannot be used without inversion symmetry.', -1) 2040 endif 2041 if(ionode) then 2042 write(0, '(a)') 'WARNING: Inversion symmetry is absent in symmetries used to reduce k-grid.' 2043 write(0, '(a)') 'Be sure inversion is still a spatial symmetry, or you must use complex version instead.' 2044 endif 2045 endif 2046 if(nspin > 1) then 2047 call errore('check_inversion', & 2048 'Real version may only be used for spin-unpolarized calculations.', nspin) 2049 endif 2050 endif 2051 2052 return 2053 2054end subroutine check_inversion 2055 2056!------------------------------------------------------------------------------- 2057! Set up some things for DFPT electron-phonon matrix elements. 2058subroutine initialize_phonon() 2059 2060 USE atom, ONLY : msh, rgrid 2061 USE becmod, ONLY : calbec, allocate_bec_type 2062 USE cell_base, ONLY : omega, tpiba, tpiba2 2063 USE control_ph, ONLY : current_iq, trans, tmp_dir_ph 2064 USE eqv, ONLY : dpsi, dvpsi, vlocq 2065 USE gvect, ONLY : g, ngm 2066 USE io_files, ONLY : tmp_dir 2067 USE io_global, ONLY : ionode 2068 USE ions_base, ONLY : nat, ntyp => nsp 2069 USE kinds, ONLY : DP 2070 USE klist, ONLY : nks, xk, nkstot 2071 USE modes, ONLY : u, npert, name_rap_mode 2072 USE noncollin_module, ONLY : npol 2073 USE ph_restart, ONLY : ph_readfile 2074 USE phus, ONLY : alphap, becp1 2075 USE qpoint, ONLY : npwq, igkq, xq, eigqts, ikks, ikqs, nksq ! phcom 2076 USE uspp, ONLY : nkb, vkb 2077 USE uspp_param, ONLY : upf 2078 USE wavefunctions_module, ONLY : evc 2079 USE wvfct, ONLY : igk, nbnd, npw, npwx, g2kin, ecutwfc 2080 2081 IMPLICIT NONE 2082 2083 integer :: nt, ik, ipol, ibnd, ig, ierr, ii, jj, iks, ike 2084 COMPLEX(DP), ALLOCATABLE :: aux1(:,:) 2085 2086 if (ionode) then 2087 WRITE ( 6, '(a)') 'Performing initialization for ionic perturbation.' 2088 endif 2089 2090 tmp_dir_ph = trim(tmp_dir) // '_ph0/' 2091 2092 if ( ANY (upf(1:ntyp)%nlcc) ) then 2093 CALL errore ( 'initialize_phonon', 'NLCC not implemented.', ierr ) 2094 endif 2095 2096 ! if gamma_gamma tricks are used, the dvscf file and wfns do not mean what we think they do... 2097 if(nkstot == 1 .and. ionode) write(0,'(a)') & 2098 "WARNING: Be sure you are not using a phonon run with gamma_gamma tricks. Must be nogg = .true." 2099 2100! allocate_phq 2101 allocate (u ( 3 * nat, 3 * nat)) 2102 allocate (npert ( 3 * nat)) 2103 allocate (name_rap_mode( 3 * nat)) 2104 2105 current_iq = 1 2106 trans = .true. 2107! lgamma = .true. ! using this causes crash from iotk in ph_readfile?? 2108 CALL ph_readfile ('data_u', ierr) ! reads from XML, for uact 2109 IF ( ierr /= 0 ) CALL errore ( 'initialize_phonon', 'failed to open phonon XML', ierr ) 2110 if(ionode) then 2111 OPEN (unit = 16, file = 'displacements.dat', form = 'formatted', status = 'replace') 2112 write(6,'(a)') 'displacements read from XML file' 2113 do ii = 1, 3 * nat 2114 do jj = 1, 3 * nat 2115 ! we must use (Re, Im) format for this to recognized as complex by a Fortran read 2116 write(16,'(a,f12.7,a,f12.7,a)',advance='no') ' (', dble(u(jj, ii)), ',', aimag(u(jj, ii)), ')' 2117 enddo 2118 write(16,*) 2119 enddo 2120 CLOSE(16) 2121 endif 2122 2123! allocate_phq 2124 allocate (dvpsi ( npwx*npol , nbnd)) 2125 2126! if (lgamma) then 2127 igkq => igk 2128 2129! phq_init 2130! IF ( lgamma ) THEN 2131 npwq = npw 2132 2133 allocate(eigqts(nat)) 2134 eigqts(:) = CMPLX(1d0, 0d0) 2135 2136 xq(1:3) = 0d0 2137 2138! allocate_phq.f90 2139 2140 if(ionode) write(6,'(a,i6)') 'Number of species: ', ntyp 2141 allocate (vlocq ( ngm , ntyp)) 2142 2143! phq_init.f90: 2144 2145 ! the fourier components of the local potential at q+G 2146 vlocq(:,:) = 0.D0 2147 ! 2148 DO nt = 1, ntyp 2149 ! 2150 IF (upf(nt)%tcoulombp) then 2151 CALL setlocq_coul ( xq, upf(nt)%zp, tpiba2, ngm, g, omega, vlocq(1,nt) ) 2152 ELSE 2153 CALL setlocq( xq, rgrid(nt)%mesh, msh(nt), rgrid(nt)%rab, rgrid(nt)%r,& 2154 upf(nt)%vloc(1), upf(nt)%zp, tpiba2, ngm, g, omega, & 2155 vlocq(1,nt) ) 2156 END IF 2157 ! 2158! write(6,*) 'init setlocq nt = ', nt 2159! write(6,*) vlocq(:, nt) 2160 END DO 2161 2162! deallocate_phq.f90:53: if(allocated(vlocq)) deallocate (vlocq) 2163 2164! initialize_ph 2165! IF ( lgamma ) THEN 2166 nksq = nks 2167 ALLOCATE(ikks(nksq), ikqs(nksq)) 2168 DO ik=1,nksq 2169 ikks(ik) = ik 2170 ikqs(ik) = ik 2171 ENDDO 2172 2173 ALLOCATE (becp1(nksq)) 2174 ALLOCATE (alphap(3,nksq)) 2175 2176 DO ik=1,nksq 2177 call allocate_bec_type ( nkb, nbnd, becp1(ik) ) 2178 DO ipol=1,3 2179 call allocate_bec_type ( nkb, nbnd, alphap(ipol,ik) ) 2180 ENDDO 2181 END DO 2182! CALL allocate_bec_type ( nkb, nbnd, becp ) 2183 2184! phq_init 2185 ! ... e) we compute the becp terms which are used in the rest of 2186 ! ... the code 2187 ! 2188 ALLOCATE( aux1( npwx*npol, nbnd ) ) 2189 2190 CALL k_pools(iks, ike) 2191 2192 do ik = 1, nksq 2193 2194 CALL gk_sort (xk (1, ik - iks + 1), ngm, g, ecutwfc & 2195 / tpiba2, npw, igk, g2kin) 2196 CALL init_us_2 ( npw, igk, xk ( 1, ik - iks + 1), vkb ) 2197 2198! npw = npwx 2199 2200! CALL davcio (evc, nwordwfc, iunwfc, ik - iks + 1, -1) ? 2201 CALL calbec (npw, vkb, evc, becp1(ik) ) 2202 ! 2203 ! ... e') we compute the derivative of the becp term with respect to an 2204 ! atomic displacement 2205 ! 2206 DO ipol = 1, 3 2207 aux1=(0.d0,0.d0) 2208 DO ibnd = 1, nbnd 2209 DO ig = 1, npw 2210 aux1(ig,ibnd) = evc(ig,ibnd) * tpiba * ( 0.D0, 1.D0 ) * & 2211 ( xk(ipol,ik) + g(ipol,igk(ig)) ) 2212 END DO 2213! IF (noncolin) THEN 2214! DO ig = 1, npw 2215! aux1(ig+npwx,ibnd)=evc(ig+npwx,ibnd)*tpiba*(0.D0,1.D0)*& 2216! ( xk(ipol,ikk) + g(ipol,igk(ig)) ) 2217! END DO 2218! END IF 2219 END DO 2220 CALL calbec (npw, vkb, aux1, alphap(ipol,ik) ) 2221 END DO 2222 enddo 2223 2224 deallocate(aux1) 2225 2226 return 2227end subroutine initialize_phonon 2228 2229!------------------------------------------------------------------------------- 2230! Set up some things for DFPT electron-electric matrix elements. 2231subroutine initialize_electric() 2232 2233 USE atom, ONLY : msh, rgrid 2234 USE becmod, ONLY : calbec, allocate_bec_type 2235 USE cell_base, ONLY : omega, tpiba, tpiba2 2236 USE control_ph, ONLY : current_iq, trans, tmp_dir_ph 2237 USE eqv, ONLY : dpsi, dvpsi, vlocq 2238 USE gvect, ONLY : g, ngm 2239 USE io_files, ONLY : tmp_dir 2240 USE io_global, ONLY : ionode 2241 USE ions_base, ONLY : nat, ntyp => nsp 2242 USE kinds, ONLY : DP 2243 USE klist, ONLY : nks, xk, nkstot 2244 USE modes, ONLY : u, npert, name_rap_mode 2245 USE noncollin_module, ONLY : npol 2246 USE ph_restart, ONLY : ph_readfile 2247 USE phus, ONLY : alphap, becp1 2248 USE qpoint, ONLY : npwq, igkq, xq, eigqts, ikks, ikqs, nksq ! phcom 2249 USE uspp, ONLY : nkb, vkb 2250 USE uspp_param, ONLY : upf 2251 USE wavefunctions_module, ONLY : evc 2252 USE wvfct, ONLY : igk, nbnd, npw, npwx, g2kin, ecutwfc 2253 2254 IMPLICIT NONE 2255 2256 if (ionode) then 2257 WRITE ( 6, '(a)') 'Performing initialization for electric perturbation.' 2258 endif 2259 2260 tmp_dir_ph = trim(tmp_dir) // '_ph0/' 2261 2262 ! if gamma_gamma tricks are used, the dvscf file and wfns do not mean what we think they do... 2263 if(nkstot == 1 .and. ionode) write(0,'(a)') & 2264 "WARNING: Be sure you are not using a phonon run with gamma_gamma tricks. Must be nogg = .true." 2265 2266 allocate ( dpsi ( npwx*npol , nbnd)) 2267 2268 ! more is needed! 2269 2270 return 2271 2272end subroutine initialize_electric 2273 2274!------------------------------------------------------------------------------- 2275 2276subroutine write_header(unit, stitle) 2277 2278 USE cell_base, ONLY : omega, alat, tpiba, tpiba2, at, bg, & 2279 ibrav, symm_type 2280 USE constants, ONLY : pi, tpi, eps6 2281 USE grid_dimensions, ONLY : nr1, nr2, nr3 2282 USE gvect, ONLY : ngm, ngm_g, ig_l2g, mill, ecutrho 2283 USE io_global, ONLY : ionode 2284 USE ions_base, ONLY : nat, atm, ityp, tau 2285 USE kinds, ONLY : DP 2286 USE lsda_mod, ONLY : nspin 2287 USE mp, ONLY : mp_sum 2288 USE mp_global, ONLY : intra_pool_comm 2289 USE scf, ONLY : rho 2290 USE symm_base, ONLY : s, ftau, nsym 2291 2292 IMPLICIT NONE 2293 2294 integer, intent(in) :: unit 2295 character(len=32), intent(in) :: stitle 2296 2297 character :: cdate*9, ctime*9, sdate*32, stime*32 2298 integer :: id, ig, i, j, k, ierr 2299 integer :: nd, ns, nst, nsf, ng_l, ng_g 2300 integer :: ntran, cell_symmetry, nrecord 2301 real (DP) :: alat2, recvol, dr1, t1 ( 3 ), t2 ( 3 ) 2302 real (DP) :: r1 ( 3, 3 ), r2 ( 3, 3 ), adot ( 3, 3 ) 2303 real (DP) :: bdot ( 3, 3 ), translation ( 3, 48 ) 2304 integer, allocatable :: g_g ( :, : ) 2305 2306 INTEGER, EXTERNAL :: atomic_number 2307 2308 CALL date_and_tim ( cdate, ctime ) 2309 WRITE ( sdate, '(A2,"-",A3,"-",A4,21X)' ) cdate(1:2), cdate(3:5), cdate(6:9) 2310 WRITE ( stime, '(A8,24X)' ) ctime(1:8) 2311 2312 nrecord = 1 2313 nd = 3 2314 2315 ng_l = ngm 2316 ng_g = ngm_g 2317 call set_spin(ns,nst,nsf,nspin) 2318 2319 ierr = 0 2320 IF ( ibrav .EQ. 0 ) THEN 2321 IF ( TRIM ( symm_type ) .EQ. 'cubic' ) THEN 2322 cell_symmetry = 0 2323 ELSEIF ( TRIM ( symm_type ) .EQ. 'hexagonal' ) THEN 2324 cell_symmetry = 1 2325 ELSE 2326 ierr = 1 2327 ENDIF 2328 ELSEIF ( ibrav .GE. 1 .AND. ibrav .LE. 3 ) THEN 2329 cell_symmetry = 0 2330 ELSEIF ( ibrav .GE. 4 .AND. ibrav .LE. 5 ) THEN 2331 cell_symmetry = 1 2332 ELSEIF ( ibrav .GE. 6 .AND. ibrav .LE. 14 ) THEN 2333 cell_symmetry = 0 2334 ELSE 2335 ierr = 1 2336 ENDIF 2337 IF ( ierr .GT. 0 ) & 2338 CALL errore ( 'write_header', 'cell_symmetry', ierr ) 2339 2340 ntran = nsym 2341 DO i = 1, ntran 2342 DO j = 1, nd 2343 DO k = 1, nd 2344 r1 ( k, j ) = dble ( s ( k, j, i ) ) 2345 ENDDO 2346 ENDDO 2347 CALL invmat ( 3, r1, r2, dr1 ) 2348 t1 ( 1 ) = dble ( ftau ( 1, i ) ) / dble ( nr1 ) 2349 t1 ( 2 ) = dble ( ftau ( 2, i ) ) / dble ( nr2 ) 2350 t1 ( 3 ) = dble ( ftau ( 3, i ) ) / dble ( nr3 ) 2351 DO j = 1, nd 2352 t2 ( j ) = 0.0D0 2353 DO k = 1, nd 2354 t2 ( j ) = t2 ( j ) + r2 ( k, j ) * t1 ( k ) 2355 ENDDO 2356 IF ( t2 ( j ) .GE. eps6 + 0.5D0 ) & 2357 t2 ( j ) = t2 ( j ) - dble ( int ( t2 ( j ) + 0.5D0 ) ) 2358 IF ( t2 ( j ) .LT. eps6 - 0.5D0 ) & 2359 t2 ( j ) = t2 ( j ) - dble ( int ( t2 ( j ) - 0.5D0 ) ) 2360 ENDDO 2361 DO j = 1, nd 2362 translation ( j, i ) = t2 ( j ) * tpi 2363 ENDDO 2364 ENDDO 2365 2366 alat2 = alat ** 2 2367 recvol = 8.0D0 * pi**3 / omega 2368 2369 DO i = 1, nd 2370 DO j = 1, nd 2371 adot ( j, i ) = 0.0D0 2372 ENDDO 2373 ENDDO 2374 DO i = 1, nd 2375 DO j = 1, nd 2376 DO k = 1, nd 2377 adot ( j, i ) = adot ( j, i ) + & 2378 at ( k, j ) * at ( k, i ) * alat2 2379 ENDDO 2380 ENDDO 2381 ENDDO 2382 2383 DO i = 1, nd 2384 DO j = 1, nd 2385 bdot ( j, i ) = 0.0D0 2386 ENDDO 2387 ENDDO 2388 DO i = 1, nd 2389 DO j = 1, nd 2390 DO k = 1, nd 2391 bdot ( j, i ) = bdot ( j, i ) + & 2392 bg ( k, j ) * bg ( k, i ) * tpiba2 2393 ENDDO 2394 ENDDO 2395 ENDDO 2396 2397 ALLOCATE ( g_g ( nd, ng_g ) ) 2398 2399 DO ig = 1, ng_g 2400 DO id = 1, nd 2401 g_g ( id, ig ) = 0 2402 ENDDO 2403 ENDDO 2404 2405 DO ig = 1, ng_l 2406 g_g ( 1, ig_l2g ( ig ) ) = mill ( 1, ig ) 2407 g_g ( 2, ig_l2g ( ig ) ) = mill ( 2, ig ) 2408 g_g ( 3, ig_l2g ( ig ) ) = mill ( 3, ig ) 2409 ENDDO 2410 2411 CALL mp_sum ( g_g, intra_pool_comm ) 2412 2413 IF(ionode) then 2414 WRITE ( unit ) stitle, sdate, stime 2415 WRITE ( unit ) nspin, ng_g, ntran, cell_symmetry, nat, ecutrho 2416 WRITE ( unit ) nr1, nr2, nr3 2417 WRITE ( unit ) omega, alat, ( ( at ( j, i ), j = 1, nd ), i = 1, nd ), & 2418 ( ( adot ( j, i ), j = 1, nd ), i = 1, nd ) 2419 WRITE ( unit ) recvol, tpiba, ( ( bg ( j, i ), j = 1, nd ), i = 1, nd ), & 2420 ( ( bdot ( j, i ), j = 1, nd ), i = 1, nd ) 2421 WRITE ( unit ) ( ( ( s ( k, j, i ), k = 1, nd ), j = 1, nd ), i = 1, ntran ) 2422 WRITE ( unit ) ( ( translation ( j, i ), j = 1, nd ), i = 1, ntran ) 2423 WRITE ( unit ) ( ( tau ( j, i ), j = 1, nd ), atomic_number ( atm ( ityp ( i ) ) ), i = 1, nat ) 2424 WRITE ( unit ) nrecord 2425 WRITE ( unit ) ng_g 2426 WRITE ( unit ) ( ( g_g ( id, ig ), id = 1, nd ), ig = 1, ng_g ) 2427 ENDIF 2428 2429 DEALLOCATE ( g_g ) 2430 RETURN 2431 2432end subroutine write_header 2433 2434!------------------------------------------------------------------------------- 2435 2436! Alternate method of calculating DFPT matrix elements, for testing. 2437! Project variation of occupied wavefunctions (dwf files) onto wavefunctions. 2438! Must hack ph code to preserve dwf's. Must be done one mode at a time. 2439! Only gives occ-unocc mtxels, which are actually not the ones needed in forces. 2440subroutine dfpt_dwf() 2441 2442 USE cell_base, ONLY : tpiba2 2443 USE constants, ONLY : rytoev 2444 USE control_ph, ONLY : tmp_dir_ph 2445 USE eqv, ONLY : dpsi 2446 USE gvect, ONLY : g, ngm 2447 USE io_files, ONLY : diropn, prefix, nwordwfc, iunwfc, nd_nmbr 2448 USE kinds, ONLY : DP 2449 USE klist, ONLY : xk, nelec 2450 USE lsda_mod, ONLY : nspin 2451 USE mp, ONLY : mp_sum 2452 USE mp_global, ONLY : intra_pool_comm 2453 USE noncollin_module, ONLY : npol 2454 USE qpoint, ONLY : nksq 2455 USE units_ph, ONLY : iudwf, lrdwf 2456 USE wavefunctions_module, ONLY : evc 2457 USE wvfct, ONLY : nbnd, npw, npwx, igk, g2kin, ecutwfc, et 2458 2459 IMPLICIT NONE 2460 2461 integer :: nrec, ipert, ik, ib, ib2, nbnd_occ, iks, ike 2462 logical :: exst 2463 complex(DP), allocatable :: mtxel(:,:) 2464 2465 nbnd_occ = nelec 2466 if(nspin == 1) nbnd_occ = nbnd_occ / 2d0 2467 allocate (dpsi (npwx*npol, nbnd_occ / nspin)) 2468 allocate (mtxel (nbnd_occ, nbnd)) 2469 2470 ! from openfilq.f90 2471 iudwf = 22 2472 ! we do not actually know the number of bands in scf/ph. we guess... 2473 lrdwf = 2 * nbnd_occ * npwx * npol 2474 CALL diropn (iudwf, 'dwf', lrdwf, exst, tmp_dir_ = tmp_dir_ph) 2475 IF (.NOT.exst) CALL errore ('dfpt_dwf','file '//trim(prefix)//'.dwf not found', 1) 2476 2477 ! this would be the index of the perturbation within the irrep. 2478 ! we do not have this info though, so we will just deal with the first one 2479 ipert = 1 2480 2481 CALL k_pools(iks, ike) 2482 do ik = iks, ike 2483 nrec = (ipert - 1) * nksq + ik 2484 CALL gk_sort (xk (1, ik - iks + 1), ngm, g, ecutwfc & 2485 / tpiba2, npw, igk, g2kin) 2486 CALL davcio (evc, nwordwfc, iunwfc, ik - iks + 1, -1) 2487 call davcio (dpsi, lrdwf, iudwf, nrec, -1) 2488! write(37,*) evc 2489! write(38,*) dpsi 2490 do ib = 1, nbnd_occ 2491 do ib2 = 1, nbnd 2492 ! zgemm? 2493! mtxel(ib, ib2) = sum(conjg(evc(1:npw, ib)) * evc(1:npw, ib2)) 2494 mtxel(ib, ib2) = sum(conjg(dpsi(1:npw, ib)) * evc(1:npw, ib2)) 2495 CALL mp_sum (mtxel(ib, ib2), intra_pool_comm) 2496 write(6,'(2i6,2f12.6)') ib, ib2, mtxel(ib, ib2) * (et(ib2, ik) - et(ib, ik)) * rytoev 2497 enddo 2498 enddo 2499 enddo 2500 2501 deallocate(dpsi) 2502 deallocate(mtxel) 2503 2504 RETURN 2505 2506end subroutine dfpt_dwf 2507 2508!------------------------------------------------------------------------------- 2509 2510! determine start and end k-points for this pool 2511! taken from PW/pw_restart.f90 2512subroutine k_pools(iks, ike) 2513 2514 USE klist, ONLY : nkstot 2515 USE mp_global, ONLY : kunit, my_pool_id, npool 2516 2517 IMPLICIT NONE 2518 2519 integer, intent(out) :: iks, ike 2520 2521 integer :: nkbl, nkl, nkr 2522 2523! ... find out the number of pools 2524! npool = nproc / nproc_pool 2525! apparently this is already initialized. --DAS 2526 2527! ... find out number of k-points blocks 2528 nkbl = nkstot / kunit 2529! ... k-points per pool 2530 nkl = kunit * ( nkbl / npool ) 2531! ... find out the remainder 2532 nkr = ( nkstot - nkl * npool ) / kunit 2533! ... Assign the remainder to the first nkr pools 2534 IF ( my_pool_id .LT. nkr ) nkl = nkl + kunit 2535! ... find out the index of the first k-point in this pool 2536 iks = nkl * my_pool_id + 1 2537 IF ( my_pool_id .GE. nkr ) iks = iks + nkr * kunit 2538! ... find out the index of the last k-point in this pool 2539 ike = iks + nkl - 1 2540 2541 RETURN 2542 2543end subroutine k_pools 2544 2545!------------------------------------------------------------------------------- 2546 2547! Set spin-related parameters ns, nst, nsf 2548! 2549subroutine set_spin(ns, nst, nsf, nspin) 2550 2551 IMPLICIT NONE 2552 2553 integer, intent(out) :: ns, nst, nsf 2554 integer, intent(in) :: nspin 2555 2556 IF ( nspin == 4 ) THEN 2557 ELSE 2558 ns = nspin 2559 nst = nspin 2560 nsf = nspin 2561 ENDIF 2562 2563 RETURN 2564 2565end subroutine set_spin 2566 2567!------------------------------------------------------------------------------- 2568 2569end module pw2bgw_m 2570 2571PROGRAM pw2bgw 2572 2573 USE constants, ONLY : eps12 2574 USE control_flags, ONLY : gamma_only 2575 USE environment, ONLY : environment_start, environment_end 2576 USE io_files, ONLY : prefix, tmp_dir, outdir 2577 USE io_global, ONLY : ionode, ionode_id 2578 USE kinds, ONLY : DP 2579 USE klist, ONLY : nkstot 2580 USE mp, ONLY : mp_bcast 2581 USE mp_global, ONLY : mp_startup 2582 USE paw_variables, ONLY : okpaw 2583 use scf, ONLY: rho_core, rhog_core 2584 USE uspp, ONLY : okvan 2585 USE wvfct, ONLY : nbnd 2586 2587 USE pw2bgw_m 2588 2589 IMPLICIT NONE 2590 2591 character(len=6) :: codename = 'PW2BGW' 2592 2593 integer :: real_or_complex 2594 logical :: wfng_flag 2595 character ( len = 256 ) :: wfng_file 2596 logical :: wfng_kgrid 2597 integer :: wfng_nk1 2598 integer :: wfng_nk2 2599 integer :: wfng_nk3 2600 real (DP) :: wfng_dk1 2601 real (DP) :: wfng_dk2 2602 real (DP) :: wfng_dk3 2603 logical :: wfng_occupation 2604 integer :: wfng_nvmin 2605 integer :: wfng_nvmax 2606 logical :: rhog_flag 2607 character ( len = 256 ) :: rhog_file 2608 logical :: vxcg_flag 2609 character ( len = 256 ) :: vxcg_file 2610 logical :: vxc0_flag 2611 character ( len = 256 ) :: vxc0_file 2612 logical :: vxc_flag 2613 character ( len = 256 ) :: vxc_file 2614 integer :: dfpt_type ! 0 = none, 1 = ionic, 2 = electric 2615 logical :: dfpt_from_dwf 2616 character ( len = 256 ) :: dfpt_file 2617 character :: vxc_integral 2618 integer :: vxc_diag_nmin 2619 integer :: vxc_diag_nmax 2620 integer :: vxc_offdiag_nmin 2621 integer :: vxc_offdiag_nmax 2622 integer :: dfpt_mode_start 2623 integer :: dfpt_mode_end 2624 logical :: vxc_zero_rho_core 2625 character ( len = 20 ) :: input_dft 2626 logical :: exx_flag 2627 logical :: vnlg_flag 2628 character ( len = 256 ) :: vnlg_file 2629 2630 NAMELIST / input_pw2bgw / prefix, outdir, & 2631 real_or_complex, wfng_flag, wfng_file, wfng_kgrid, & 2632 wfng_nk1, wfng_nk2, wfng_nk3, wfng_dk1, wfng_dk2, wfng_dk3, & 2633 wfng_occupation, wfng_nvmin, wfng_nvmax, rhog_flag, rhog_file, & 2634 vxcg_flag, vxcg_file, vxc0_flag, vxc0_file, vxc_flag, & 2635 vxc_file, vxc_integral, vxc_diag_nmin, vxc_diag_nmax, & 2636 vxc_offdiag_nmin, vxc_offdiag_nmax, vxc_zero_rho_core, & 2637 input_dft, exx_flag, vnlg_flag, vnlg_file 2638 2639 integer :: ii, ios, imode, nmode 2640 character ( len = 256 ) :: output_file_name, tmpstr 2641 logical :: die_spinors 2642 2643 character (len=256), external :: trimcheck 2644 character (len=1), external :: lowercase 2645 2646#ifdef __PARA 2647 CALL mp_startup ( ) 2648#endif 2649 CALL environment_start ( codename ) 2650 2651 ! assign defaults 2652 prefix = 'prefix' 2653 CALL get_env ( 'ESPRESSO_TMPDIR', outdir ) 2654 IF ( TRIM ( outdir ) == ' ' ) outdir = './' 2655 real_or_complex = 2 2656 wfng_flag = .FALSE. 2657 wfng_file = 'WFN' 2658 wfng_kgrid = .FALSE. 2659 wfng_nk1 = 0 2660 wfng_nk2 = 0 2661 wfng_nk3 = 0 2662 wfng_dk1 = 0.0D0 2663 wfng_dk2 = 0.0D0 2664 wfng_dk3 = 0.0D0 2665 wfng_occupation = .FALSE. 2666 wfng_nvmin = 0 2667 wfng_nvmax = 0 2668 rhog_flag = .FALSE. 2669 rhog_file = 'RHO' 2670 vxcg_flag = .FALSE. 2671 vxcg_file = 'VXC' 2672 vxc0_flag = .FALSE. 2673 vxc0_file = 'vxc0.dat' 2674 vxc_flag = .FALSE. 2675 vxc_file = 'vxc.dat' 2676 vxc_integral = 'g' 2677 vxc_diag_nmin = 0 2678 vxc_diag_nmax = 0 2679 vxc_offdiag_nmin = 0 2680 vxc_offdiag_nmax = 0 2681 vxc_zero_rho_core = .TRUE. 2682 input_dft = 'sla+pz' 2683 exx_flag = .FALSE. 2684 vnlg_flag = .FALSE. 2685 vnlg_file = 'VNL' 2686 dfpt_type = 0 2687 dfpt_from_dwf = .FALSE. 2688 dfpt_file = 'dfpt.dat' 2689 dfpt_mode_start = 1 2690 dfpt_mode_end = 0 ! do all 2691 2692 IF ( ionode ) THEN 2693 CALL input_from_file ( ) 2694 READ ( 5, input_pw2bgw, iostat = ios ) 2695 IF ( ios /= 0 ) CALL errore ( codename, 'input_pw2bgw', abs ( ios ) ) 2696 2697 DO ii = 1, LEN_TRIM (vxc_integral) 2698 vxc_integral(ii:ii) = lowercase (vxc_integral(ii:ii)) 2699 END DO 2700 2701 IF ( real_or_complex /= 1 .AND. real_or_complex /= 2 ) & 2702 CALL errore ( codename, 'real_or_complex', 1 ) 2703 IF ( vxc_integral /= 'r' .AND. vxc_integral /= 'g' ) & 2704 CALL errore ( codename, 'vxc_integral', 1 ) 2705 ENDIF 2706 2707 tmp_dir = trimcheck ( outdir ) 2708 CALL mp_bcast ( outdir, ionode_id ) 2709 CALL mp_bcast ( tmp_dir, ionode_id ) 2710 CALL mp_bcast ( prefix, ionode_id ) 2711 CALL mp_bcast ( real_or_complex, ionode_id ) 2712 CALL mp_bcast ( wfng_flag, ionode_id ) 2713 CALL mp_bcast ( wfng_file, ionode_id ) 2714 CALL mp_bcast ( wfng_kgrid, ionode_id ) 2715 CALL mp_bcast ( wfng_nk1, ionode_id ) 2716 CALL mp_bcast ( wfng_nk2, ionode_id ) 2717 CALL mp_bcast ( wfng_nk3, ionode_id ) 2718 CALL mp_bcast ( wfng_dk1, ionode_id ) 2719 CALL mp_bcast ( wfng_dk2, ionode_id ) 2720 CALL mp_bcast ( wfng_dk3, ionode_id ) 2721 CALL mp_bcast ( wfng_occupation, ionode_id ) 2722 CALL mp_bcast ( wfng_nvmin, ionode_id ) 2723 CALL mp_bcast ( wfng_nvmax, ionode_id ) 2724 CALL mp_bcast ( rhog_flag, ionode_id ) 2725 CALL mp_bcast ( rhog_file, ionode_id ) 2726 CALL mp_bcast ( vxcg_flag, ionode_id ) 2727 CALL mp_bcast ( vxcg_file, ionode_id ) 2728 CALL mp_bcast ( vxc0_flag, ionode_id ) 2729 CALL mp_bcast ( vxc0_file, ionode_id ) 2730 CALL mp_bcast ( vxc_flag, ionode_id ) 2731 CALL mp_bcast ( vxc_integral, ionode_id ) 2732 CALL mp_bcast ( vxc_file, ionode_id ) 2733 CALL mp_bcast ( vxc_diag_nmin, ionode_id ) 2734 CALL mp_bcast ( vxc_diag_nmax, ionode_id ) 2735 CALL mp_bcast ( vxc_offdiag_nmin, ionode_id ) 2736 CALL mp_bcast ( vxc_offdiag_nmax, ionode_id ) 2737 CALL mp_bcast ( vxc_zero_rho_core, ionode_id ) 2738 CALL mp_bcast ( input_dft, ionode_id ) 2739 CALL mp_bcast ( exx_flag, ionode_id ) 2740 CALL mp_bcast ( vnlg_flag, ionode_id ) 2741 CALL mp_bcast ( vnlg_file, ionode_id ) 2742 CALL mp_bcast ( dfpt_type, ionode_id ) 2743 CALL mp_bcast ( dfpt_from_dwf, ionode_id ) 2744 CALL mp_bcast ( dfpt_file, ionode_id ) 2745 CALL mp_bcast ( dfpt_mode_start, ionode_id ) 2746 CALL mp_bcast ( dfpt_mode_end, ionode_id ) 2747 2748 if(dfpt_type < 0 .or. dfpt_type > 2) & 2749 call errore ( 'pw2bgw', 'unknown dfpt_type', dfpt_type) 2750 if(dfpt_type /= 0 .and. dfpt_mode_start < 1) & 2751 call errore ( 'pw2bgw', 'dfpt_mode_start < 1', dfpt_mode_start) 2752 2753 if (ionode) then 2754 WRITE ( 6, '(a)') 'Reading basic information from files.' 2755 endif 2756 CALL read_file ( ) 2757 2758 if (ionode) then 2759 if (MAX (MAXVAL (ABS (rho_core (:) ) ), MAXVAL (ABS (rhog_core (:) ) ) ) & 2760 .LT. eps12) then 2761 WRITE ( 6, '(/,5x,"NLCC is absent")' ) 2762 else 2763 WRITE ( 6, '(/,5x,"NLCC is present")' ) 2764 endif 2765 endif 2766 if (okvan) call errore ( 'pw2bgw', 'BGW cannot use USPP.', 3 ) 2767 if (okpaw) call errore ( 'pw2bgw', 'BGW cannot use PAW.', 4 ) 2768 if (gamma_only) call errore ( 'pw2bgw', 'BGW cannot use gamma-only run.', 5 ) 2769 die_spinors = (nspin == 4) 2770 if (die_spinors) call errore ( 'pw2bgw', 'BGW cannot use spinors.', 6 ) 2771 if (real_or_complex == 1 .AND. vxc_flag .AND. vxc_offdiag_nmax > 0) & 2772 call errore ( 'pw2bgw', 'Off-diagonal matrix elements of Vxc ' // & 2773 'with real wavefunctions are not implemented, compute them in ' // & 2774 'Sigma using VXC.', 7) 2775 2776 CALL openfil_pp ( ) 2777 2778 if ( ionode ) WRITE ( 6, '("")' ) 2779 2780 if ( ionode ) write ( 6, '(a, i8)' ) "Total number of bands available: ", nbnd 2781 if ( ionode ) write ( 6, '(a, i8)' ) "Total number of k-points available: ", nkstot 2782 2783 IF ( wfng_flag ) THEN 2784 output_file_name = TRIM ( outdir ) // '/' // TRIM ( wfng_file ) 2785 IF ( ionode ) WRITE ( 6, '(5x,"call write_wfng")' ) 2786 CALL start_clock ( 'write_wfng' ) 2787 CALL write_wfng ( output_file_name, real_or_complex, & 2788 wfng_kgrid, wfng_nk1, wfng_nk2, wfng_nk3, wfng_dk1, wfng_dk2, & 2789 wfng_dk3, wfng_occupation, wfng_nvmin, wfng_nvmax ) 2790 CALL stop_clock ( 'write_wfng' ) 2791 IF ( ionode ) WRITE ( 6, '(5x,"done write_wfng",/)' ) 2792 ENDIF 2793 2794 IF ( rhog_flag ) THEN 2795 output_file_name = TRIM ( outdir ) // '/' // TRIM ( rhog_file ) 2796 IF ( ionode ) WRITE ( 6, '(5x,"call write_rhog")' ) 2797 CALL start_clock ( 'write_rhog' ) 2798 CALL write_rhog ( output_file_name, real_or_complex ) 2799 CALL stop_clock ( 'write_rhog' ) 2800 IF ( ionode ) WRITE ( 6, '(5x,"done write_rhog",/)' ) 2801 ENDIF 2802 2803 IF ( vxcg_flag ) THEN 2804 output_file_name = TRIM ( outdir ) // '/' // TRIM ( vxcg_file ) 2805 IF ( ionode ) WRITE ( 6, '(5x,"call write_vxcg")' ) 2806 CALL start_clock ( 'write_vxcg' ) 2807 CALL write_vxcg ( output_file_name, real_or_complex, & 2808 vxc_zero_rho_core, input_dft, exx_flag ) 2809 CALL stop_clock ( 'write_vxcg' ) 2810 IF ( ionode ) WRITE ( 6, '(5x,"done write_vxcg",/)' ) 2811 ENDIF 2812 2813 IF ( vxc0_flag ) THEN 2814 output_file_name = TRIM ( outdir ) // '/' // TRIM ( vxc0_file ) 2815 IF ( ionode ) WRITE ( 6, '(5x,"call write_vxc0")' ) 2816 CALL start_clock ( 'write_vxc0' ) 2817 CALL write_vxc0 ( output_file_name, vxc_zero_rho_core, input_dft, & 2818 exx_flag ) 2819 CALL stop_clock ( 'write_vxc0' ) 2820 IF ( ionode ) WRITE ( 6, '(5x,"done write_vxc0",/)' ) 2821 ENDIF 2822 2823 IF ( vxc_flag ) THEN 2824 output_file_name = TRIM ( outdir ) // '/' // TRIM ( vxc_file ) 2825 IF ( ionode ) WRITE ( 6, '(5x,"call write_vxc")' ) 2826 CALL start_clock ( 'write_vxc' ) 2827 CALL write_vxc ( vxc_integral, output_file_name, & 2828 vxc_diag_nmin, vxc_diag_nmax, & 2829 vxc_offdiag_nmin, vxc_offdiag_nmax, & 2830 vxc_zero_rho_core, input_dft, exx_flag, 0, 0 ) 2831 CALL stop_clock ( 'write_vxc' ) 2832 IF ( ionode ) WRITE ( 6, '(5x,"done write_vxc",/)' ) 2833 ENDIF 2834 2835 IF ( vnlg_flag ) THEN 2836 output_file_name = TRIM ( outdir ) // '/' // TRIM ( vnlg_file ) 2837 IF ( ionode ) WRITE ( 6, '(5x,"call write_vnlg")' ) 2838 CALL start_clock ( 'write_vnlg' ) 2839 CALL write_vnlg ( output_file_name, wfng_kgrid, wfng_nk1, & 2840 wfng_nk2, wfng_nk3, wfng_dk1, wfng_dk2, wfng_dk3 ) 2841 CALL stop_clock ( 'write_vnlg' ) 2842 IF ( ionode ) WRITE ( 6, '(5x,"done write_vnlg",/)' ) 2843 ENDIF 2844 2845 IF ( dfpt_type /= 0 ) THEN 2846 IF(dfpt_type == 1) THEN 2847 call initialize_phonon() 2848 ELSE 2849 call initialize_electric() 2850 ENDIF 2851 CALL check_dfpt_modes ( dfpt_mode_start, dfpt_type, nmode ) 2852 IF(dfpt_mode_end == 0) dfpt_mode_end = nmode 2853 2854 IF(dfpt_from_dwf) CALL dfpt_dwf() 2855 2856 DO imode = dfpt_mode_start, dfpt_mode_end 2857 WRITE(tmpstr,'(i6.6)') imode 2858 IF ( ionode ) WRITE ( 6, '(5x,"call write_dfpt")' ) 2859 CALL start_clock ( 'write_dfpt' ) 2860 2861! output_file_name = TRIM ( outdir ) // '/' // TRIM ( dfpt_file ) // "_dvscf_mode" // TRIM ( tmpstr ) 2862! CALL write_vxc ( 'r', output_file_name, & 2863! vxc_diag_nmin, vxc_diag_nmax, & 2864! vxc_diag_nmin, vxc_diag_nmax, & ! we need the whole matrix here 2865! input_dft, .false., dfpt_type, imode) 2866 2867 output_file_name = TRIM ( outdir ) // '/' // TRIM ( dfpt_file ) // "_mode" // TRIM ( tmpstr ) 2868 CALL write_vxc ( 'g', output_file_name, & 2869 vxc_diag_nmin, vxc_diag_nmax, & 2870 vxc_diag_nmin, vxc_diag_nmax, & ! we need the whole matrix here 2871 vxc_zero_rho_core, input_dft, .false., dfpt_type, imode) 2872 CALL stop_clock ( 'write_dfpt' ) 2873 IF ( ionode ) WRITE ( 6, '(5x,"done write_dfpt",/)' ) 2874 ENDDO 2875 ENDIF 2876 2877 IF ( ionode ) WRITE ( 6, * ) 2878 IF ( wfng_flag ) CALL print_clock ( 'write_wfng' ) 2879 IF ( rhog_flag ) CALL print_clock ( 'write_rhog' ) 2880 IF ( vxcg_flag ) CALL print_clock ( 'write_vxcg' ) 2881 IF ( vxc0_flag ) CALL print_clock ( 'write_vxc0' ) 2882 IF ( vxc_flag ) CALL print_clock ( 'write_vxc' ) 2883 IF ( dfpt_type /= 0 ) CALL print_clock ( 'write_dfpt' ) 2884 IF ( vnlg_flag ) CALL print_clock ( 'write_vnlg' ) 2885 IF ( wfng_flag .AND. real_or_complex .EQ. 1 ) THEN 2886 IF ( ionode ) WRITE ( 6, '(/,5x,"Called by write_wfng:")' ) 2887 CALL print_clock ( 'real_wfng' ) 2888 ENDIF 2889 2890 CALL environment_end ( codename ) 2891 CALL stop_pp ( ) 2892 2893 STOP 2894 2895END PROGRAM pw2bgw 2896