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