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