1 2! Copyright (C) 2005-2009 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 13Aprile2005 8! GENERATES INPUT for GW code 9!tested on: Silicon bulk, Germanium Bulk, Na4, InP bulk 10! Please note just symmorphic symm. op. have to be used 11! Use input option of pw.x: force_symmorphic=.TRUE. 12 13! Update 20 November 2017 (Davide Grassano, Adriano Mosca Conte) 14! Added spin-orbit case for nspin == 4 15! Added flag to input file: 16! Emin: Starting energy of eps spectra 17! Emax: Max energy of eps spectra 18! DeltaE: Energy steps for eps sectra 19! qplda: disable/enable priting of qplda file (default = false) 20! vkb: disable/enable priting of fort.15 file (default = false) 21! vxcdiag: disable/enable priting of vxcdiag.dat file (default = false) 22! Using proper units(tpiba) for k-vectors in k+G sum before calculating matrixelements 23 24!----------------------------------------------------------------------- 25PROGRAM pw2gw 26 !----------------------------------------------------------------------- 27 28 ! This subroutine writes files containing plane wave coefficients 29 ! and other stuff needed by GW codes 30 31 USE io_files, ONLY : prefix, tmp_dir 32 USE io_global, ONLY : ionode, ionode_id 33 USE mp, ONLY : mp_bcast 34 USE mp_global, ONLY : mp_startup 35 USE mp_images, ONLY : intra_image_comm 36 USE mp_pools, ONLY : kunit 37 USE environment,ONLY : environment_start, environment_end 38 USE us, ONLY : spline_ps 39 USE kinds, ONLY : DP 40 ! 41 IMPLICIT NONE 42 ! 43 CHARACTER(LEN=256), EXTERNAL :: trimcheck 44 CHARACTER(LEN=256) :: outdir 45 ! 46 INTEGER :: ios 47 INTEGER :: kunittmp 48 LOGICAL :: use_gmaps, qplda, vkb, vxcdiag 49 CHARACTER(len=20) :: what 50 REAL(kind=DP) :: Emin, Emax, DeltaE 51 52 NAMELIST / inputpp / prefix, outdir, what, use_gmaps, Emin, Emax, DeltaE, & 53 qplda, vkb, vxcdiag 54 ! 55 ! initialise environment 56 ! 57#if defined(__MPI) 58 CALL mp_startup ( ) 59#endif 60 CALL environment_start ( 'PW2GW' ) 61 ! 62 ! set default values for variables in namelist 63 ! 64 prefix = 'pwscf' 65 CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) 66 IF ( trim( outdir ) == ' ' ) outdir = './' 67 what = 'gw' 68 qplda = .false. 69 vkb = .false. 70 vxcdiag = .false. 71 use_gmaps = .false. 72 Emin = 0.0 73 Emax = 30.0 74 DeltaE = 0.05 75 76 ios = 0 77 IF ( ionode ) THEN 78 ! 79 READ (5, inputpp, iostat=ios) 80 tmp_dir = trimcheck (outdir) 81 ! 82 ENDIF 83 ! 84 CALL mp_bcast( ios, ionode_id, intra_image_comm ) 85 IF (ios /= 0) CALL errore('pw2gw', 'reading inputpp namelist', abs(ios)) 86 ! 87 ! ... Broadcast variables 88 ! 89 CALL mp_bcast( prefix, ionode_id, intra_image_comm ) 90 CALL mp_bcast(tmp_dir, ionode_id, intra_image_comm ) 91 CALL mp_bcast( what, ionode_id, intra_image_comm ) 92 CALL mp_bcast( use_gmaps, ionode_id, intra_image_comm ) 93 CALL mp_bcast( qplda, ionode_id, intra_image_comm ) 94 CALL mp_bcast( vkb, ionode_id, intra_image_comm ) 95 CALL mp_bcast( vxcdiag, ionode_id, intra_image_comm ) 96 CALL mp_bcast( Emin, ionode_id, intra_image_comm ) 97 CALL mp_bcast( Emax, ionode_id, intra_image_comm ) 98 CALL mp_bcast( DeltaE, ionode_id, intra_image_comm ) 99 100 ! 101 102 spline_ps = .false. 103 104 CALL read_file 105 CALL openfil_pp 106 ! 107 CALL mp_bcast(spline_ps, ionode_id, intra_image_comm) 108#if defined __MPI 109 kunittmp = kunit 110#else 111 kunittmp = 1 112#endif 113 ! 114 !WRITE(*,*) what, qplda, vxcdiag, Emin, Emax, DeltaE 115 IF( trim( what ) == 'gw' ) THEN 116 CALL compute_gw ( Emin, Emax, DeltaE, use_gmaps, qplda, vkb, vxcdiag ) 117 ELSEIF( trim( what ) == 'gmaps' ) THEN 118 CALL write_gmaps ( kunittmp ) 119 ENDIF 120 ! 121 CALL environment_end ( 'PW2GW' ) 122 ! 123 CALL stop_pp 124 125END PROGRAM pw2gw 126 127 128SUBROUTINE compute_gw( omegamin, omegamax, d_omega, use_gmaps, qplda, vkb, vxcdiag_f ) 129 130 ! This routine creates the QPLDA and the matrixelements 131 ! tform = .false. UNFORMATTED QPLDA 132 ! tform = .true. FORMATTED QPLDA 133 ! tsingle must be always true 134 135 USE kinds, ONLY : DP, sgl 136 USE constants, ONLY : eps8, pi, AUTOEV, rytoev 137 USE cell_base, ONLY : alat, tpiba2, at, bg, omega 138 USE symm_base, ONLY : s, nsym 139 USE wvfct, ONLY : npwx, nbnd, wg, et 140 USE gvecw, ONLY : gcutw 141 USE control_flags, ONLY : gamma_only 142 USE gvect, ONLY : ngm, g, gg, ig_l2g 143 USE fft_base, ONLY: dfftp 144 USE fft_interfaces, ONLY : fwfft, invfft 145 USE klist , ONLY : nks, xk, wk, ngk, igk_k 146 USE lsda_mod, ONLY : nspin 147 USE io_files, ONLY : nwordwfc, iunwfc 148 USE wavefunctions, ONLY : evc, psic 149 USE io_global, ONLY : ionode, ionode_id 150 USE mp, ONLY : mp_sum , mp_max 151 USE mp_pools, ONLY : npool 152 USE mp_images, ONLY : intra_image_comm 153 USE mp_world, ONLY : mpime, nproc 154 USE mp_wave, ONLY : mergewf 155 USE parallel_include 156 USE scf, ONLY : rho, rho_core, rhog_core 157 USE ener, ONLY : etxc, vtxc 158 159 USE uspp_param, ONLY : upf, nh 160 USE uspp, ONLY : nhtol 161 USE us, ONLY : tab, tab_d2y, spline_ps 162 USE ions_base, ONLY : ntyp => nsp 163 USE klist, ONLY : ngk 164 165 IMPLICIT NONE 166 167 LOGICAL, INTENT(in) :: use_gmaps, qplda, vkb, vxcdiag_f 168 REAL(kind=DP), INTENT( in) :: omegamax, omegamin, d_omega 169 170 INTEGER :: ii(16), ngw, nkpt, ig, ik, ir, n, i,j,k, io = 98, iband1, iband2 171 INTEGER :: npw, omax, o, iproc 172 INTEGER, ALLOCATABLE :: in1(:), in2(:), in3(:) 173 INTEGER, ALLOCATABLE :: in1_tmp(:), in2_tmp(:), in3_tmp(:) 174 INTEGER, ALLOCATABLE :: inx_rcv(:), ig_l2g_rcv(:) 175 LOGICAL :: t_form = .false., t_single = .true. 176 REAL(kind=sgl) :: a1_s(3), a2_s(3), a3_s(3) 177 REAL(kind=sgl), ALLOCATABLE :: xk_s(:,:), eig_s(:,:), focc_s(:,:) 178 REAL(kind=DP):: g2max, a1(3), a2(3), a3(3),norm, xkgk(3), rrhotwx(3), delta 179 REAL(kind=DP):: alpha, egap, halfalpha, Df, const, dummy 180 !REAL(kind=DP), PARAMETER :: omegamax = 30.0 181 REAL(kind=DP), ALLOCATABLE:: gsort(:), eig(:,:), focc(:,:), kpg(:,:), omegatt(:), omeg(:) 182 REAL(kind=DP), ALLOCATABLE:: pp1(:,:), pp2(:,:), pp3(:,:) 183 REAL(kind=DP), ALLOCATABLE:: epsx(:,:), epsy(:,:), epsz(:,:) 184 REAL(kind=DP), ALLOCATABLE:: epstx(:), epsty(:), epstz(:) 185 REAL(kind=DP) :: epsxx, epsyy, epszz 186 REAL(kind=DP) :: vxcdiag 187 REAL(kind=DP), ALLOCATABLE :: vxc(:,:) 188 COMPLEX(kind=DP):: rhotwx(3), ctemp, dasomma(3) 189 COMPLEX(kind=DP), ALLOCATABLE:: c0(:), c0_m(:,:), c0_tmp_dp(:) !, c0_tmp(:) !, c0_gamma(:) 190 COMPLEX(kind=sgl), ALLOCATABLE:: c0_s(:), c0_tmp(:) !, c0_gamma_s(:) 191 CHARACTER(len=80) :: titleo(2) 192 INTEGER :: igwx, igwxx, comm, ierr, ig_max, igwx_r 193 INTEGER :: igwx_p(nproc) 194 INTEGER, ALLOCATABLE :: igk_l2g(:) 195 INTEGER :: npol 196 ! 197 REAL(kind=DP), ALLOCATABLE :: vkb0(:), djl(:), vec_tab(:), vec_tab_d2y(:) 198 INTEGER :: nb, nt, size_tab, size_tab_d2y, ipw, l 199 ! 200 ! REAL(kind=DP) :: norma ! Variable needed only for DEBUG 201 ! 202#if defined __MPI 203 INTEGER :: istatus( MPI_STATUS_SIZE ) 204#endif 205 ! 206 IF( nspin == 4 ) WRITE(*,*) "nspin = 4" 207 IF( nspin == 2 ) CALL errore('pw2gw','Spin polarization not implemented',1) 208 IF( npool > 1 ) CALL errore('pw2gw','parallel run with pools not allowed yet',1) 209 ! 210 ! 211 IF( mpime == 0 ) THEN 212 IF (t_form) THEN 213 WRITE (6,'(//" writing LDA info on unit 98 FORMATTED")') 214 IF( qplda) OPEN (io, FILE='QPLDA',STATUS='unknown',FORM='FORMATTED') 215 ELSE 216 WRITE (6,'(//" writing LDA info on unit io UNFORMATTED")') 217 IF( qplda) OPEN (io, FILE='QPLDA',STATUS='unknown',FORM='UNFORMATTED') 218 ENDIF 219 WRITE (6,'(//" writing matrixelements on unit 98 FORMATTED")') 220 OPEN (90, FILE='matrixelements',STATUS='unknown',FORM='FORMATTED') 221 ENDIF 222 ! 223 ! file's title [2 lines] 224 ! 225 titleo(1)='pw2gw' 226 titleo(2)='test version' 227 IF( mpime == 0 ) THEN 228 IF (t_form) THEN 229 IF( qplda) WRITE (io,'(A80/A80)') titleo(1), titleo(2) 230 ELSE 231 IF( qplda) WRITE (io) titleo(1) 232 IF( qplda) WRITE (io) titleo(2) 233 ENDIF 234 ! 235 WRITE(6,*) 'qplda title' 236 WRITE(6,*) titleo(1) 237 WRITE(6,*) titleo(2) 238 ENDIF 239 ! 240 ! Read 16 integers (reserved for future flags) 241 ! Flags used so far: 242 ! I1 = 0 if QPLDA file is formatted, 1 if unformatted 243 ! I2 = 0 if RWG format, 1 if BF format 244 ! I3 = 1 if non-symmorphic operations (+vectors) included, otherwise 0 245 ! 246 ii(:) = 0 247 IF (t_form) THEN 248 ii(1)=0 249 IF( mpime == 0 .AND. qplda ) WRITE (io,'(16I5)') ii 250 ELSE 251 ii(1)=1 252 IF( mpime == 0 .AND. qplda ) WRITE (io) ii 253 ENDIF 254 ! 255 WRITE(6,'(16I5)') ii 256 ! 257 ! write real-space lattice vectors (Cartesian, in au) [3 lines] 258 ! 259 a1(:)=at(:,1)*alat 260 a2(:)=at(:,2)*alat 261 a3(:)=at(:,3)*alat 262 a1_s(:) = a1(:) 263 a2_s(:) = a2(:) 264 a3_s(:) = a3(:) 265 266 IF( mpime == 0 ) THEN 267 ! 268 IF (t_form) THEN 269 IF( qplda) WRITE (io,'(3E26.18)') a1, a2, a3 270 ELSE 271 IF (t_single) THEN 272 IF( qplda) WRITE (io) a1_s, a2_s, a3_s 273 ELSE 274 IF( qplda) WRITE (io) a1, a2, a3 275 ENDIF 276 ENDIF 277 ! 278 WRITE(6,*) 'Vettori di reticolo diretto' 279 WRITE(6,'(a,3E26.18)') 'a1', a1_s 280 WRITE(6,'(a,3E26.18)') 'a2', a2_s 281 WRITE(6,'(a,3E26.18)') 'a3', a3_s 282 ! 283 ENDIF 284 ! 285 ! Write symmetry operations. 286 ! The matrix s is the transpose of the symmetry matrix in direct space, 287 ! in units of a_i. But the transpose of the symmetry matrix in real space 288 ! is the symmetry matrix in reciprocal space so "s" is already the symmetry 289 ! matrix in reciprocal space in units of b_i 290 ! The gw code will read row by row a matrix and will treat it as symmetry 291 ! matrix in reciprocal space in units of b_i 292 ! In other words, the gw code expects as input the direct space symmetry 293 ! matrix, in units of a_i, written columnwise 294 ! 295 IF( mpime == 0 ) THEN 296 WRITE(6,*)'nrot=',nsym 297 WRITE(6,'(3E26.18)') (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) 298 IF (t_form) THEN 299 IF( qplda) WRITE (io,'(I2)') nsym 300 IF( qplda) WRITE (io,'(3E26.18)') (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) 301 IF (ii(3) == 1) THEN 302 ! READ (10,1020) ((VOFFSET(I,J),I=1,3),J=1,NOP) 303 ! WRITE (6,'(//" Run program CNVNSY to convert QPLDA file first.")') 304 CALL errore('pw2gw','non-symmorphic translation vectors',ii(3)) 305 ENDIF 306 ELSE 307 IF( qplda) WRITE (io) nsym 308 IF (t_single) THEN 309 IF( qplda) WRITE (io) (((float(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) 310 ELSE 311 WRITE (io) (((dble(s(i,j,k)),j=1,3),i=1,3),k=1,nsym) 312 ENDIF 313 IF (ii(3) == 1) THEN 314 ! READ (10,1020) ((VOFFSET(I,J),I=1,3),J=1,NOP) 315 CALL errore('pw2gw','non-symmorphic translation vectors',ii(3)) 316 ENDIF 317 ENDIF 318 ENDIF 319 ! 320 ! write reciprocal lattice vectors (in reciprocal lattice units; 321 ! ie in the basis of the reciprocal lattice basis vectors) 322 ! 323 ! PWscf stores psi(k+G), using |k+G| to order the components; 324 ! GW codes require on input psi_k(G), using the same set of G 325 ! 326 g2max = 0.0d0 327 !DEBUG 328 IF (ionode) WRITE(6,*) ' nks ', nks 329 IF (ionode) WRITE(6,*) ' k points in cartesian coordinates' 330 IF (ionode) WRITE(6,'(1x,3f10.6)') ( (xk(i,ik),i=1,3), ik=1,nks) 331 !DEBUG 332 igwx = 0 ! maximum G vector index 333 g2max = gcutw ! RAGGIO DELLA SFERA |G+k|<cut 334 DO ik = 1, nks 335 npw = ngk(ik) 336 ! WRITE( 6, * ) 'DEBUG g2max ', g2max 337 ! g2max, g2kin = RAGGIO DELLA SFERA |G+k|<cut, non MASSIMO |G| nella sfera 338 ! g2max <= gcutw PER COSTRUZIONE 339 igwx = max( igwx, maxval( igk_k(1:npw,ik) ) ) 340 ENDDO 341 !IF (ionode) write(*,*) "igwx = ", igwx 342 ! 343 ! ngw = number of G-vectors (complete shells) such that G2 <= G2max 344 ! ovvero <= RAGGIO della SFERA, in pratica trova i G2 relativi a GAMMA 345 ! 346 DO ngw = 1, ngm 347 IF ( gg(ngw) > g2max + eps8) GOTO 100 348 ENDDO 349 CALL errore ( 'pw2gw','max G in PW not found?!?',ngw) 350100 ngw = ngw - 1 351 352 ! Pongo NGW pari al massimo indice tra i vettori G che fanno parte delle 353 ! sfere |G+k|<cut per qualsiasi k 354 ! 355 !IF (ionode) write( 6, * ) ' igwx= ', igwx 356 357 ngw = igwx 358 359 ! 360 ! PWscf stores G in order of increasing module 361 ! 362 ALLOCATE (in1(ngw), in2(ngw), in3(ngw)) 363 DO ig=1,ngw 364 in1(ig) = nint ( at(1,1)*g(1,ig) + at(2,1)*g(2,ig) + at(3,1)*g(3,ig) ) 365 in2(ig) = nint ( at(1,2)*g(1,ig) + at(2,2)*g(2,ig) + at(3,2)*g(3,ig) ) 366 in3(ig) = nint ( at(1,3)*g(1,ig) + at(2,3)*g(2,ig) + at(3,3)*g(3,ig) ) 367 ENDDO 368 369 igwxx = maxval( ig_l2g( 1:ngw ) ) 370 CALL mp_max( igwxx, intra_image_comm ) 371 IF (ionode) WRITE(*,*) "NDIMCP = ", igwxx 372 373 igwx_p = 0 374 igwx_p( mpime + 1 ) = igwx 375 CALL mp_sum( igwx_p, intra_image_comm ) 376 377 IF( mpime == 0 ) THEN 378 ! 379 ! allocate arrays 380 ! 381 ALLOCATE ( in1_tmp(igwxx), in2_tmp(igwxx), in3_tmp(igwxx) ) 382 ! copy local data of the root proc into global vector 383 DO ig = 1, ngw 384 in1_tmp( ig_l2g(ig) ) = in1(ig) 385 in2_tmp( ig_l2g(ig) ) = in2(ig) 386 in3_tmp( ig_l2g(ig) ) = in3(ig) 387 ENDDO 388 ! 389#if defined __MPI 390 ALLOCATE( ig_l2g_rcv( igwxx ) ) 391 ALLOCATE( inx_rcv( igwxx ) ) 392 ! 393 DO iproc = 2, nproc 394 CALL MPI_RECV( ig_l2g_rcv, igwx_p( iproc ), MPI_INTEGER, (iproc-1), iproc, intra_image_comm, istatus, IERR ) 395 ! 396 CALL MPI_RECV( inx_rcv, igwx_p( iproc ), MPI_INTEGER, (iproc-1), iproc+NPROC, intra_image_comm, istatus, IERR ) 397 DO ig = 1, igwx_p( iproc ) 398 in1_tmp( ig_l2g_rcv( ig ) ) = inx_rcv( ig ) 399 ENDDO 400 CALL MPI_RECV( inx_rcv, igwx_p( iproc ), MPI_INTEGER, (iproc-1), iproc+2*NPROC, intra_image_comm, istatus, IERR ) 401 DO ig = 1, igwx_p( iproc ) 402 in2_tmp( ig_l2g_rcv( ig ) ) = inx_rcv( ig ) 403 ENDDO 404 CALL MPI_RECV( inx_rcv, igwx_p( iproc ), MPI_INTEGER, (iproc-1), iproc+3*NPROC, intra_image_comm, istatus, IERR ) 405 DO ig = 1, igwx_p( iproc ) 406 in3_tmp( ig_l2g_rcv( ig ) ) = inx_rcv( ig ) 407 ENDDO 408 ENDDO 409 ! 410 DEALLOCATE( ig_l2g_rcv ) 411 DEALLOCATE( inx_rcv ) 412#endif 413 ! 414 ELSE 415 ! 416#if defined __MPI 417 CALL MPI_SEND( ig_l2g, igwx, MPI_INTEGER, 0, mpime+1, intra_image_comm, IERR ) 418 CALL MPI_SEND( in1(1), igwx, MPI_INTEGER, 0, mpime+1+NPROC, intra_image_comm, IERR ) 419 CALL MPI_SEND( in2(1), igwx, MPI_INTEGER, 0, mpime+1+2*NPROC, intra_image_comm, IERR ) 420 CALL MPI_SEND( in3(1), igwx, MPI_INTEGER, 0, mpime+1+3*NPROC, intra_image_comm, IERR ) 421#endif 422 ! 423 ENDIF 424 425 IF (mpime == 0) WRITE(*,*) "fine debug sui punti g" 426 427 428 429 IF (t_form) THEN 430 IF( mpime == 0 ) THEN 431 IF( qplda) WRITE (io,'(I12)') igwxx 432 IF( qplda) WRITE (io,'(3I5)') (in1_tmp(ig),in2_tmp(ig),in3_tmp(ig),ig=1,igwxx) 433 ENDIF 434 ELSE 435 IF( mpime == 0 ) THEN 436 IF( qplda) WRITE (io) igwxx 437 IF( qplda) WRITE (io) (in1_tmp(ig),in2_tmp(ig),in3_tmp(ig),ig=1,igwxx) 438 ENDIF 439 ENDIF 440 ! 441 DEALLOCATE ( in1, in2, in3 ) 442 IF( mpime == 0 ) THEN 443 DEALLOCATE ( in1_tmp, in2_tmp, in3_tmp ) 444 ENDIF 445 ! 446 ! WRITE k-points (in RL units) 447 ! 448 ! transformation in relative coordinates with respect to b1,b2,b3 449 ! 450 CALL cryst_to_cart (nks, xk, at, -1) 451 ! xk(3,nkpt) in input deve essere in coordinate cartesiane! 452 nkpt = nks 453 ALLOCATE (xk_s(3,nkpt)) 454 xk_s(:,:) = xk(:,1:nkpt) 455 ! 456 IF( mpime == 0 ) THEN 457 OPEN(65,file='k.dat') 458 WRITE(65,'(1x,3f10.6,1x,f10.6)') ( xk_s(:,ik), wk(ik)*0.5, ik=1,nks ) 459 CLOSE(unit=65) 460 ENDIF 461 462 IF( mpime == 0 ) THEN 463 IF (t_form) THEN 464 IF( qplda) WRITE (io,'(I12)') nkpt 465 IF( qplda) WRITE (io,'(3E26.18)') ((xk_s(i,ik),i=1,3),ik=1,nkpt) 466 ELSE 467 IF( qplda) WRITE (io) nkpt 468 IF( qplda) WRITE(6,*) 'nkpt',nkpt 469 IF(t_single) THEN 470 IF( qplda) WRITE (io) ((xk_s(i,ik),i=1,3),ik=1,nkpt) 471 ELSE 472 IF( qplda) WRITE (io) ((xk(i,ik),i=1,3),ik=1,nkpt) 473 ENDIF 474 ENDIF 475 WRITE(6,'(1x,3f10.6)') ( (xk_s(i,ik),i=1,3), ik=1,nkpt) 476 ENDIF 477 478 IF( vkb) THEN 479! -------------------------- 480! vkb0 481! -------------------------- 482 DO ik=1,nkpt 483 npw = ngk(ik) 484 WRITE(15,*) "npw", npw 485 ALLOCATE(vkb0(1:npw)) 486 487 size_tab=size(tab,1) 488 size_tab_d2y=size(tab_d2y,1) 489 490 ALLOCATE(vec_tab(1:size_tab)) 491 if(.not.allocated(vec_tab_d2y)) ALLOCATE(vec_tab_d2y(1:size_tab_d2y)) 492 493 DO nt = 1, ntyp 494 DO nb = 1, upf(nt)%nbeta 495 vkb0(:) = 0.0_dp 496 vec_tab(:) = 0.0_dp 497 vec_tab_d2y(:) = 0.0_dp 498 vec_tab(:) = tab(:,nb,nt) 499 IF(spline_ps) vec_tab_d2y(:) = tab_d2y(:,nb,nt) 500 CALL gen_us_vkb0(ik,npw,vkb0,size_tab,vec_tab,spline_ps,vec_tab_d2y) 501 WRITE(15,*) "---------------DEBUG-VKB0----------------------" 502 WRITE(15,*) "ik= ", ik 503 WRITE(15,*) "nt= ", nt 504 WRITE(15,*) "nb= ", nb 505 WRITE(15,*) "l= ", upf(nt)%lll(nb) 506 WRITE (15,'(8f15.9)') vkb0 507 WRITE(15,*) "--------------END-DEBUG------------------------" 508 ! WRITE(io) vkb0 509 ENDDO 510 ENDDO 511 512 DEALLOCATE(vkb0) 513 DEALLOCATE(vec_tab) 514 IF(allocated(vec_tab_d2y)) DEALLOCATE(vec_tab_d2y) 515 516 ENDDO 517!--------------------------- 518! djl 519!--------------------------- 520 DO ik=1,nkpt 521 npw = ngk(ik) 522 523 ALLOCATE(djl(1:npw)) 524 525 size_tab=size(tab,1) 526 size_tab_d2y=size(tab_d2y,1) 527 528 ALLOCATE(vec_tab(1:size_tab)) 529 IF(.not. allocated(vec_tab_d2y)) ALLOCATE(vec_tab_d2y(1:size_tab_d2y)) 530 DO nt = 1, ntyp 531 DO nb = 1, upf(nt)%nbeta 532 djl(:) = 0.0_dp 533 vec_tab(:) = 0.0_dp 534 vec_tab_d2y(:) = 0.0_dp 535 vec_tab(:) = tab(:,nb,nt) 536 IF(spline_ps) vec_tab_d2y(:) = tab_d2y(:,nb,nt) 537 CALL gen_us_djl(ik,npw,djl,size_tab,vec_tab,spline_ps,vec_tab_d2y) 538 ! WRITE(0,*) "---------------DEBUG-----------------------" 539 ! WRITE(0,*) "spline: ", spline_ps 540 ! WRITE(0,*) "ik= ", ik 541 ! WRITE(0,*) "nt= ", nt 542 ! WRITE(0,*) "nb= ", nb 543 ! WRITE(0,*) "l= ", upf(nt)%lll(nb) 544 ! WRITE (0,'(8f15.9)') djl 545 ! WRITE(0,*) "--------------END-DEBUG------------------------" 546 ! WRITE(io) djl 547 ENDDO 548 ENDDO 549 550 DEALLOCATE(djl) 551 DEALLOCATE(vec_tab) 552 IF(allocated(vec_tab_d2y)) DEALLOCATE(vec_tab_d2y) 553 554 ENDDO 555 556 ENDIF !vkb 557!----------------------- 558!----------------------- 559 560 ! 561 ! WRITE energies (Hartrees) (in ascending order, hopefully they are ordered) 562 ! 563 n = nbnd 564 ALLOCATE (eig(n,nkpt), eig_s(n,nkpt)) 565 eig(:,:) = et(:,1:nkpt)*0.5d0 566 567 IF (t_form) THEN 568 IF( mpime == 0 .AND. qplda) WRITE (io,'(I12)') n 569 IF( mpime == 0 .AND. qplda) WRITE (io,'(3E26.18)') ((eig(i,ik),ik=1,nkpt),i=1,n) 570 ELSE 571 IF( mpime == 0 .AND. qplda) WRITE (io) n 572 IF(t_single) THEN 573 DO ik=1,nkpt 574 DO i=1,n 575 eig_s(i,ik)=eig(i,ik) 576 ENDDO 577 ENDDO 578 WRITE(6,*) 'nbndsi=',n 579 IF( mpime == 0 .AND. qplda) WRITE (io) ((eig_s(i,ik),ik=1,nkpt),i=1,n) 580 ELSE 581 WRITE(6,*) 'nbndsi=',n 582 IF( mpime == 0 .AND. qplda) WRITE (io) ((eig(i,ik),ik=1,nkpt),i=1,n) 583 ENDIF 584 ENDIF 585 586! write(6,*) 'autovalori energia per 10bande e tutti kpt' 587! WRITE(6,'(10F10.7)') ( ( eig(i,ik)*27.21 ,ik=1,nkpt), i=1,10 ) 588! DEALLOCATE (eig_s, eig) 589 ! 590 ! occupation numbers 591 ! 592 ALLOCATE (focc(n,nkpt), focc_s(n,nkpt)) 593! focc(:,:) = wg(:,1:nkpt) 594! focc_s(:,:) = wg(:,1:nkpt) 595 DO j=1,n 596 DO ik=1,nkpt 597 focc(j,ik)=wg(j,ik)*2.0d0/wk(ik) 598 ENDDO 599 ENDDO 600 601 focc_s(:,:) = focc(:,:) 602 603 IF( mpime == 0 ) THEN 604 IF (t_form) THEN 605 IF( qplda) WRITE (io,'(3E26.18)') ((focc(i,ik), ik=1,nkpt), i=1,n) 606 ELSE 607 IF(t_single) THEN 608 IF( qplda) WRITE (io) ((focc_s(i,ik),ik=1,nkpt),i=1,n) 609 ELSE 610 IF( qplda) WRITE (io) ((focc(i,ik),ik=1,nkpt),i=1,n) 611 ENDIF 612 ENDIF 613 ENDIF 614 615 WRITE (6,*) nkpt 616 WRITE (6,*) 'weights:' 617 WRITE (6,'(10f10.7)') (wk(ik), ik=1,nkpt) 618 619 DO ik = 1, nkpt 620 WRITE (6,*) 'ik=', ik 621 WRITE (6,'(10f10.7)') (focc_s(j,ik), j=1,n) 622 ENDDO 623 624! DEALLOCATE (focc_s, focc) 625 ! 626 627 ! omax = nbnd*6 628 omax = floor((omegamax - omegamin)/d_omega) 629 WRITE(6,*) 'io sono omax = ', omax 630 alpha = d_omega 631 WRITE(6,*) 'alpha = ', alpha 632 halfalpha= alpha*.5 633 WRITE(6,*) 'halfalpha = ', halfalpha 634 const = 4.d0*pi**2*AUTOEV**3/omega 635 WRITE(6,*) 'const = ', const 636 637 WRITE(*,*) "sono qui 6" 638 ALLOCATE( omeg(omax+1)) 639 ALLOCATE( epsx(nkpt,omax+1), epsy(nkpt,omax+1), epsz(nkpt,omax+1) ) 640 ALLOCATE( epstx(omax+1), epsty(omax+1), epstz(omax+1) ) 641 ALLOCATE( pp1(nkpt,omax+1), pp2(nkpt,omax+1), pp3(nkpt,omax+1), omegatt(omax+1) ) 642 643 DO o = 1, omax + 1 644 omeg(o) = omegamin + (o-1)*alpha 645 ENDDO 646 647 epstx(:)=0.0 648 epsty(:)=0.0 649 epstz(:)=0.0 650 epsx(:,:)=0.0 651 epsy(:,:)=0.0 652 epsz(:,:)=0.0 653 pp1(:,:)=0.0 654 pp2(:,:)=0.0 655 pp3(:,:)=0.0 656 657 WRITE(6,*) pp1(1,1), epstx(1) 658 659 npol=1 660 IF (nspin == 4) npol=2 661 write(*,*) "igwx, npw", igwx, npw 662 ALLOCATE ( c0(igwx*npol), c0_s(igwx*npol), kpg(3,igwx), c0_m(igwx*npol,n), & 663 c0_tmp(igwxx*npol) ) 664 c0 = 0.0 665 c0_s = 0.0 666 kpg = 0.0 667 c0_m = 0.0 668 c0_tmp = 0.0 669 !IF (gamma_only) ALLOCATE ( c0_gamma(2*igwx-1), c0_gamma_s(2*igwx-1) ) 670 CALL cryst_to_cart (nks, xk, bg, +1) 671 IF (ionode) WRITE(6,*) 'Costruisco le psi ed il matrixelements' 672 IF (ionode) WRITE(6,*) 'Controllo: I punti k ora devo essere in coordinate cartesiane!' 673 IF (ionode) WRITE(6,'(1x,3f10.6)') ( (xk(i,ik),i=1,3), ik=1,nkpt) 674 675 676 DO ik = 1, nkpt 677 ! 678 npw = ngk(ik) 679 ALLOCATE( igk_l2g( npw ) ) 680 ! 681 DO ig = 1, npw 682 ! 683 igk_l2g(ig) = ig_l2g(igk_k(ig,ik)) 684 ! 685 ENDDO 686 ! 687 IF( use_gmaps ) THEN 688 ! 689 c0_m = 0.0d0 690 ! 691 CALL read_and_collect( c0_m, size( c0_m, 1 ), n, ik ) 692 ! 693 ELSE 694 ! 695 CALL davcio( evc, 2*nwordwfc, iunwfc, ik, -1 ) 696 ! 697 ! copy coefficient from array evc(:,n) (ordered as |k+G|) 698 ! into array c0 with |G| ordering 699 ! 700 DO ig = 1, npw 701 IF( igk_k(ig,ik) < 1 .or. igk_k(ig,ik) > size( c0 ) ) & 702 CALL errore(' pw2gw ', ' c0 too small ', 1 ) 703 ENDDO 704 705 ! read wavefunctions and write the matrixelemnts 706 707 DO i = 1, n 708 709 ALLOCATE( c0_tmp_dp( igwxx ) ) 710 711 CALL mergewf( evc(:,i), c0_tmp_dp, npw, igk_l2g(:), mpime, nproc, 0, intra_image_comm ) 712 ! 713 ! important: missing components must be set to zero 714 c0 (:) = 0.d0 715 DO ig=1,npw 716 c0(igk_k(ig,ik)) = evc(ig,i) 717 IF (nspin == 4) c0(igk_k(ig,ik)+igwx) = evc(ig+npwx,i) 718 ENDDO 719 c0_m(:,i)=c0(:) 720 721 c0_tmp = c0_tmp_dp 722 IF( mpime == 0 .AND. qplda) WRITE(io) c0_tmp ! c0_s 723 724 DEALLOCATE( c0_tmp_dp ) 725 726 ENDDO 727 ENDIF 728 729 DEALLOCATE( igk_l2g ) 730 731 ! k + g thet must be in 2piba units 732 kpg(:,:) = 0.d0 733 DO ig=1,npw 734 kpg(:,igk_k(ig,ik))= xk(:,ik)+g(:,igk_k(ig,ik)) 735 ENDDO 736 737 DO iband1 = 1,n 738 IF ( focc(iband1,ik)>=1e-4) THEN 739 DO iband2 = 1,n 740 delta=2.0d0-focc(iband2,ik) 741 IF (delta>1e-4) THEN 742 743 rhotwx = 0.0 744 DO ig=1,igwx 745 xkgk(1)= kpg(1,ig) 746 xkgk(2)= kpg(2,ig) 747 xkgk(3)= kpg(3,ig) 748 ctemp= conjg(c0_m(ig,iband1))*c0_m(ig,iband2) 749 rhotwx(1) = rhotwx(1) + xkgk(1) * ctemp 750 rhotwx(2) = rhotwx(2) + xkgk(2) * ctemp 751 rhotwx(3) = rhotwx(3) + xkgk(3) * ctemp 752 IF (nspin == 4) THEN 753 ctemp= conjg(c0_m(ig+igwx,iband1))*c0_m(ig+igwx,iband2) 754 rhotwx(1) = rhotwx(1) + xkgk(1) * ctemp 755 rhotwx(2) = rhotwx(2) + xkgk(2) * ctemp 756 rhotwx(3) = rhotwx(3) + xkgk(3) * ctemp 757 ENDIF 758 ENDDO 759 760 CALL mp_sum( rhotwx, intra_image_comm ) 761 762 IF (mpime == 0) THEN 763 rrhotwx(1)=tpiba2* real(rhotwx(1)*conjg(rhotwx(1))) 764 rrhotwx(2)=tpiba2* real(rhotwx(2)*conjg(rhotwx(2))) 765 rrhotwx(3)=tpiba2* real(rhotwx(3)*conjg(rhotwx(3))) 766 WRITE (90,'(1x,3i5,3e16.8,2f8.4)') ik,iband1,iband2,rrhotwx(1),rrhotwx(2), & 767 rrhotwx(3),(eig(iband2,ik)-eig(iband1,ik))*AUTOEV, (focc(iband1,ik)-focc(iband2,ik)) 768 egap = (eig(iband2,ik)-eig(iband1,ik))*AUTOEV 769 Df = focc(iband1,ik)-focc(iband2,ik) 770 IF (egap>1e-3.and.Df>1e-4) THEN 771 DO o=1, omax+1 772 dummy = abs(egap - omeg(o)) 773 IF (dummy<halfalpha) THEN 774 pp1(ik,o)=pp1(ik,o)+rrhotwx(1)*Df/egap**2 775 pp2(ik,o)=pp2(ik,o)+rrhotwx(2)*Df/egap**2 776 pp3(ik,o)=pp3(ik,o)+rrhotwx(3)*Df/egap**2 777 ENDIF 778 ENDDO 779 ENDIF 780 ENDIF 781 ENDIF 782 ENDDO 783 ENDIF 784 ENDDO 785 786 ENDDO 787 788! CALCULATE POTENTIAL MATRIX ELEMNTS 789 790 IF( vxcdiag_f) THEN 791 OPEN (UNIT=313,FILE="vxcdiag.dat",STATUS="UNKNOWN") 792 WRITE(313,*) "# K BND <Vxc>" 793 ALLOCATE ( vxc(dfftp%nnr,nspin) ) 794 ! 795 CALL v_xc (rho, rho_core, rhog_core, etxc, vtxc, vxc) 796 ! 797 DO ik=1,nkpt 798 npw = ngk(ik) 799 CALL davcio( evc, 2*nwordwfc, iunwfc, ik, -1 ) 800 DO iband1 = 1, nbnd 801 psic(:) = (0.d0, 0.d0) 802 DO ig = 1, npw 803 psic(dfftp%nl(igk_k(ig,ik))) = evc(ig,iband1) 804 ENDDO 805 806 CALL invfft ('Rho', psic, dfftp) 807 vxcdiag = 0.0d0 808 !norma = 0.0d0 809 DO ir = 1, dfftp%nnr 810 vxcdiag = vxcdiag + vxc(ir,nspin) * & 811 ( dble(psic (ir) ) **2 + aimag(psic (ir) ) **2) 812 ! norma = norma + ( DBLE(psic (ir) ) **2 + AIMAG(psic (ir) ) **2) / (nr1*nr2*nr3) 813 ENDDO 814 ! PG: this is the correct integral - 27/8/2010 815 vxcdiag = vxcdiag * rytoev / (dfftp%nr1*dfftp%nr2*dfftp%nr3) 816 CALL mp_sum( vxcdiag, intra_image_comm ) !, intra_pool_comm ) 817 ! ONLY FOR DEBUG! 818 !IF (norma /= 1.0) THEN 819 ! WRITE(*,*) "norma =", norma 820 ! WRITE(*,*) "nrxx =", nrxx 821 ! STOP 822 !ENDIF 823 WRITE(313,"(i1,2x,i1,4x,f18.14)") ik, iband1, vxcdiag 824 ENDDO 825 ENDDO 826 DEALLOCATE ( vxc ) 827 CLOSE (313) 828 ENDIF 829 830! 831 ! 832 IF ( mpime == 0 ) THEN 833 834 WRITE(6, * ) ' out from k-points loop' 835 WRITE(6, * ) ' Starting writing epsx,y,z,tot' 836 837 WRITE(6,*) pp1(1,100) 838 WRITE(6,*) pp1(1,350) 839 840 OPEN (91, FILE='epsX.dat',STATUS='unknown',FORM='FORMATTED') 841 OPEN (92, FILE='epsY.dat',STATUS='unknown',FORM='FORMATTED') 842 OPEN (93, FILE='epsZ.dat',STATUS='unknown',FORM='FORMATTED') 843 OPEN (94, FILE='epsTOT.dat',STATUS='unknown',FORM='FORMATTED') 844 845 DO ik = 1, nkpt 846 DO o =2, omax+1 847 epsx(ik,o) = const * pp1(ik,o)*wk(ik)*0.5/ alpha 848 epsy(ik,o) = const * pp2(ik,o)*wk(ik)*0.5/ alpha 849 epsz(ik,o) = const * pp3(ik,o)*wk(ik)*0.5/ alpha 850 ENDDO 851 ENDDO 852 853 WRITE(6, * ) epsx(1,150),epsx(1,300) 854 855 DO o = 2, omax + 1 856 omegatt(o) = (omeg(o-1)+omeg(o))*0.5 857 DO ik = 1, nkpt 858 epsxx= (epsx(ik,o-1)+epsx(ik,o))*0.5 859 epsyy= (epsy(ik,o-1)+epsy(ik,o))*0.5 860 epszz= (epsz(ik,o-1)+epsz(ik,o))*0.5 861 epstx(o)=epstx(o)+epsxx 862 epsty(o)=epsty(o)+epsyy 863 epstz(o)=epstz(o)+epszz 864 ENDDO 865 WRITE(91,"(f15.6,1x,f15.6)") omegatt(o), epstx(o) 866 WRITE(92,"(f15.6,1x,f15.6)") omegatt(o), epsty(o) 867 WRITE(93,"(f15.6,1x,f15.6)") omegatt(o), epstz(o) 868 WRITE(94,"(f15.6,1x,f15.6)") omegatt(o), (epstx(o)+ epsty(o)+ epstz(o))/3.0 869 ENDDO 870 871 WRITE(6, * ) ' Hey bello sto a fini' 872 873 CLOSE(91) 874 CLOSE(92) 875 CLOSE(93) 876 CLOSE(94) 877 878 ENDIF 879 880 DEALLOCATE (xk_s) 881 DEALLOCATE (eig_s, eig) 882 DEALLOCATE (focc_s, focc) 883 DEALLOCATE (c0_s, c0, kpg, c0_m) 884 DEALLOCATE (omeg, pp1,pp2, pp3, omegatt) 885 DEALLOCATE ( epsx, epsy, epsz ) 886 DEALLOCATE ( epstx, epsty, epstz ) 887 ! 888 IF( mpime == 0 .AND. qplda) CLOSE(io) 889 IF( mpime == 0 ) CLOSE(90) 890 ! 891END SUBROUTINE compute_gw 892 893 894!----------------------------------------------------------------------- 895SUBROUTINE write_gmaps ( kunit) 896 !----------------------------------------------------------------------- 897 ! 898 USE io_global, ONLY : stdout 899 USE cell_base, ONLY : at, bg, alat 900 USE ions_base, ONLY : atm, nat 901 USE gvect, ONLY : ngm, ngm_g, ig_l2g, g 902 USE lsda_mod, ONLY : nspin, isk 903 USE ions_base, ONLY : ntyp => nsp, tau, ityp 904 USE wvfct, ONLY : nbnd, npwx, et 905 USE gvecw, ONLY : gcutw 906 USE klist, ONLY : nkstot, ngk, nks, xk 907 USE wavefunctions, ONLY : evc 908 USE io_files, ONLY : nd_nmbr, tmp_dir, prefix, iunwfc, nwordwfc 909 USE io_global, ONLY : ionode 910 USE mp_images, ONLY : intra_image_comm, my_image_id 911 USE mp_pools, ONLY : nproc_pool, my_pool_id, intra_pool_comm 912 USE mp, ONLY : mp_sum, mp_max 913 USE mp_world, ONLY : nproc, mpime 914 915 916 IMPLICIT NONE 917 INTEGER :: kunit 918 919 INTEGER :: npw, i, j, k, ig, ik, ibnd, na, ngg, ikw 920 INTEGER, ALLOCATABLE :: kisort(:) 921 INTEGER :: ike, iks, npw_g, npwx_g, ispin 922 INTEGER, EXTERNAL :: global_kpoint_index 923 INTEGER, ALLOCATABLE :: ngk_g( : ) 924 INTEGER, ALLOCATABLE :: ngk_gw( : ) 925 INTEGER, ALLOCATABLE :: itmp( :, : ) 926 INTEGER, ALLOCATABLE :: igwk( : ) 927 INTEGER, ALLOCATABLE :: igk_l2g( :, : ) 928 REAL(kind=8), ALLOCATABLE :: gk(:) 929 930 real(kind=8) :: wfc_scal 931 LOGICAL :: twf0, twfm, twrite_wfc 932 933 ! 934 ! 935 IF( ionode ) WRITE( stdout, fmt="(//,'WRITING G-MAPS for each processor' )" ) 936 937 IF( nkstot > 0 ) THEN 938 939 IF( ( kunit < 1 ) .or. ( mod( nkstot, kunit ) /= 0 ) ) & 940 CALL errore( ' write_wannier ',' wrong kunit ', 1 ) 941 942 IF( ( nproc_pool > nproc ) .or. ( mod( nproc, nproc_pool ) /= 0 ) ) & 943 CALL errore( ' write_wannier ',' nproc_pool ', 1 ) 944 945 iks = global_kpoint_index (nkstot, 1) 946 ike = iks + nks - 1 947 948 ENDIF 949 950 ! find out the global number of G vectors: ngm_g 951 ngm_g = ngm 952 CALL mp_sum( ngm_g, intra_pool_comm ) 953 954 955 ! build the G+k array indexes 956 ALLOCATE ( igk_l2g( npwx, ik ) ) 957 ALLOCATE ( kisort( npwx ) ) 958 ALLOCATE ( gk( npwx ) ) 959 DO ik = 1, nks 960 kisort = 0 961 CALL gk_sort (xk (1, ik+iks-1), ngm, g, gcutw, npw, kisort(1), gk) 962 DO ig = 1, npw 963 igk_l2g(ig,ik) = ig_l2g(kisort(ig)) 964 ENDDO 965 ngk (ik) = npw 966 ENDDO 967 DEALLOCATE (gk, kisort) 968 969 ! compute the global number of G+k vectors for each k point 970 ALLOCATE( ngk_g( nkstot ) ) 971 ALLOCATE( ngk_gw( nkstot/nspin ) ) 972 ngk_g = 0 973 ngk_g( iks:ike ) = ngk( 1:nks ) 974 CALL mp_sum( ngk_g, intra_image_comm ) 975 976 ! compute the Maximum G vector index among all G+k an processors 977 npw_g = maxval( ig_l2g(:) ) ! ( igk_l2g(:,:) ) 978 CALL mp_max( npw_g, intra_image_comm ) 979 980 ! compute the Maximum number of G vector among all k points 981 npwx_g = maxval( ngk_g( 1:nkstot ) ) 982 983 984 ALLOCATE( igwk( npwx_g ) ) 985 986 DO ik = 1, nkstot 987 igwk = 0 988 ALLOCATE( itmp( npw_g, 1 ) ) 989 itmp = 0 990 IF( ik >= iks .and. ik <= ike ) THEN 991 DO ig = 1, ngk( ik-iks+1 ) 992 itmp( ig_l2g( ig ), 1 ) = ig_l2g( ig ) 993 ENDDO 994 ENDIF 995 CALL mp_sum( itmp, intra_image_comm ) 996 ngg = 0 997 DO ig = 1, npw_g 998 IF( itmp( ig, 1 ) == ig ) THEN 999 ngg = ngg + 1 1000 igwk( ngg ) = ig 1001 ENDIF 1002 ENDDO 1003 IF( ngg /= ngk_g( ik ) ) THEN 1004 WRITE( stdout,*) ' ik, ngg, ngk_g = ', ik, ngg, ngk_g( ik ) 1005 ENDIF 1006 DEALLOCATE( itmp ) 1007 IF( ionode ) THEN 1008 ! write (40)( igwk(ig), ig = 1, npwx_g ) 1009 ENDIF 1010 ENDDO 1011 1012 DEALLOCATE( igwk ) 1013 1014 DO ik = 1, nkstot 1015 IF( (ik >= iks) .and. (ik <= ike) ) THEN 1016 ispin = isk( ik ) 1017 WRITE( 100 + mpime ) ik, iks, ike, nkstot, kunit, nproc, ispin, nspin, npw_g, & 1018 nbnd, ngk(ik-iks+1), 2*nwordwfc, npwx, iunwfc, nd_nmbr 1019 WRITE( 100 + mpime ) ( igk_l2g( i, ik-iks+1 ), i = 1, ngk(ik-iks+1) ) 1020 ENDIF 1021 ENDDO 1022 1023 DEALLOCATE ( ngk_g ) 1024 DEALLOCATE ( ngk_gw ) 1025 DEALLOCATE (igk_l2g) 1026 1027END SUBROUTINE write_gmaps 1028 1029 1030SUBROUTINE read_and_collect( c, ldc, n, ik ) 1031 USE io_global, ONLY : stdout 1032 USE io_files, ONLY : prefix 1033 USE kinds, ONLY : DP, sgl 1034 1035 IMPLICIT NONE 1036 1037 INTEGER :: ldc, n, ik 1038 COMPLEX(DP) :: c( ldc, n ) 1039 INTEGER :: ik_ , iks, ike, nkstot, kunit, nproc_ , ispin, nspin, npw_g , nbnd , ngk 1040 INTEGER :: nwordwfc, npwx, iunwfc 1041 INTEGER :: nfile, ip, i, j 1042 COMPLEX(DP), ALLOCATABLE :: evc( :, : ) 1043 INTEGER, ALLOCATABLE :: igk_l2g( : ) 1044 LOGICAL :: exst 1045 CHARACTER(len=3) :: nd_nmbr 1046 1047 READ( 100 ) ik_ , iks, ike, nkstot, kunit, nproc_ , ispin, nspin, npw_g , & 1048 nbnd , ngk, nwordwfc, npwx, iunwfc, nd_nmbr 1049 1050 REWIND( 100 ) 1051 1052 nfile = nproc_ 1053 1054 CLOSE( iunwfc ) 1055 1056 DO ip = 0, nfile - 1 1057 READ( 100 + ip ) ik_ , iks, ike, nkstot, kunit, nproc_ , ispin, nspin, npw_g , & 1058 nbnd , ngk, nwordwfc, npwx, iunwfc, nd_nmbr 1059 WRITE( stdout, * ) 'DEBUG nd_nmbr ', nd_nmbr 1060 IF( ( ik_ == ik ) .and. ( ik_ >= iks ) .and. ( ik_ <= ike ) ) THEN 1061 ALLOCATE( evc( npwx, nbnd ) ) 1062 ALLOCATE( igk_l2g( ngk ) ) 1063 READ( 100 + ip ) ( igk_l2g( i ), i = 1, ngk ) 1064 CALL diropn_gw ( 99, trim( prefix )//'.wfc', 2*nwordwfc, exst, ip, nd_nmbr ) 1065 CALL davcio ( evc, 2*nwordwfc, 99, (ik-iks+1), - 1 ) 1066 CLOSE( 99 ) 1067 DO j = 1, n 1068 DO i = 1, ngk 1069 c( igk_l2g( i ), j ) = evc( i, j ) 1070 ENDDO 1071 ENDDO 1072 DEALLOCATE( evc ) 1073 DEALLOCATE( igk_l2g ) 1074 ENDIF 1075 REWIND( 100 + ip ) 1076 ENDDO 1077 1078 RETURN 1079END SUBROUTINE 1080 1081! 1082! Copyright (C) 2001-2003 PWSCF group 1083! This file is distributed under the terms of the 1084! GNU General Public License. See the file `License' 1085! in the root directory of the present distribution, 1086! or http://www.gnu.org/copyleft/gpl.txt . 1087! 1088! 1089!----------------------------------------------------------------------- 1090SUBROUTINE diropn_gw (unit, filename, recl, exst, mpime, nd_nmbr_ ) 1091 !----------------------------------------------------------------------- 1092 ! 1093 ! this routine opens a file in tmp_dir for direct I/O access 1094 ! If appropriate, the node number is added to the file name 1095 ! 1096 USE kinds 1097 USE io_files 1098 IMPLICIT NONE 1099 1100 ! 1101 ! first the input variables 1102 ! 1103 CHARACTER(len=*) :: filename 1104 ! input: name of the file to open 1105 INTEGER :: unit, recl 1106 ! input: unit of the file to open 1107 ! input: length of the records 1108 LOGICAL :: exst 1109 ! output: if true the file exists 1110 INTEGER :: mpime 1111 ! input: processor index 1112 CHARACTER(len=3) :: nd_nmbr_ 1113 ! 1114 ! local variables 1115 ! 1116 CHARACTER(len=256) :: tempfile 1117 ! complete file name 1118 CHARACTER(len=80) :: assstr 1119 INTEGER :: ios, unf_recl, ierr, direct_io_factor 1120 ! used to check I/O operations 1121 ! length of the record 1122 ! error code 1123 LOGICAL :: opnd 1124 ! if true the file is already opened 1125 REAL(dp):: dummy 1126 1127 IF (unit < 0) CALL errore ('diropn', 'wrong unit', 1) 1128 ! 1129 ! we first check that the file is not already openend 1130 ! 1131 ios = 0 1132 INQUIRE (unit = unit, opened = opnd) 1133 IF (opnd) CALL errore ('diropn', "can't open a connected unit", abs(unit)) 1134 ! 1135 ! then we check the filename 1136 ! 1137 1138 IF (filename == ' ') CALL errore ('diropn', 'filename not given', 2) 1139 tempfile = trim(tmp_dir) // trim(filename) // trim( nd_nmbr_ ) 1140 1141 INQUIRE (file = tempfile, exist = exst) 1142 ! 1143 ! the unit for record length is unfortunately machine-dependent 1144 ! 1145 INQUIRE (IOLENGTH=direct_io_factor) dummy 1146 unf_recl = direct_io_factor * recl 1147 IF (unf_recl <= 0) CALL errore ('diropn', 'wrong record length', 3) 1148 ! 1149 OPEN ( unit, file = trim(tempfile), iostat = ios, form = 'unformatted', & 1150 status = 'unknown', access = 'direct', recl = unf_recl ) 1151 1152 IF (ios /= 0) CALL errore ('diropn', 'error opening '//filename, unit) 1153 RETURN 1154END SUBROUTINE diropn_gw 1155 1156!---------------------------------------------------------------------- 1157subroutine gen_us_djl (ik,npw,djl,size_tab,vec_tab, spline_ps, vec_tab_d2y) 1158 !---------------------------------------------------------------------- 1159 ! 1160 ! Calculates the kleinman-bylander pseudopotentials with the 1161 ! derivative of the spherical harmonics projected on vector u 1162 ! 1163 USE kinds, ONLY : DP 1164 USE io_global, ONLY : stdout 1165 USE constants, ONLY : tpi 1166 USE cell_base, ONLY : tpiba 1167 USE klist, ONLY : xk, igk_k 1168 USE gvect, ONLY : g 1169 USE us, ONLY : nqx, dq 1170 USE splinelib, ONLY : splint_deriv 1171 USE uspp_param, ONLY : upf 1172 ! 1173 implicit none 1174 ! 1175 integer, intent(in) :: ik, npw 1176 real(DP), intent(inout) ::djl(1:npw) 1177 integer, intent(in) :: size_tab 1178 real(DP), intent(in) :: vec_tab(1:size_tab) 1179 real(DP), intent(in) :: vec_tab_d2y(1:size_tab) 1180 logical :: spline_ps 1181 ! 1182 integer :: i0, i1, i2, & 1183 i3, ig 1184 real(DP), allocatable :: gk(:,:), q (:) 1185 real(DP) :: px, ux, vx, wx 1186 complex(DP), allocatable :: sk (:) 1187 1188 integer :: iq 1189 real(DP), allocatable :: xdata(:) 1190 real(DP) :: qt 1191 1192 1193 allocate ( gk(3,npw) ) 1194 allocate ( q(npw) ) 1195 1196 do ig = 1, npw 1197 gk (1, ig) = xk (1, ik) + g (1, igk_k (ig,ik) ) 1198 gk (2, ig) = xk (2, ik) + g (2, igk_k (ig,ik) ) 1199 gk (3, ig) = xk (3, ik) + g (3, igk_k (ig,ik) ) 1200 q (ig) = gk(1, ig)**2 + gk(2, ig)**2 + gk(3, ig)**2 1201 enddo 1202 1203 do ig = 1, npw 1204 q (ig) = sqrt ( q(ig) ) * tpiba 1205 end do 1206 1207 if (spline_ps) then 1208 allocate(xdata(nqx)) 1209 do iq = 1, nqx 1210 xdata(iq) = (iq - 1) * dq 1211 enddo 1212 endif 1213 1214 ! calculate beta in G-space using an interpolation table 1215 do ig = 1, npw 1216 qt = sqrt(q(ig)) * tpiba 1217 if (spline_ps) then 1218 djl(ig) = splint_deriv(xdata, vec_tab(:), & 1219 vec_tab_d2y(:), qt) 1220 else 1221 px = qt / dq - int (qt / dq) 1222 ux = 1.d0 - px 1223 vx = 2.d0 - px 1224 wx = 3.d0 - px 1225 i0 = qt / dq + 1 1226 i1 = i0 + 1 1227 i2 = i0 + 2 1228 i3 = i0 + 3 1229 djl (ig) = vec_tab (i0) * (-vx*wx-ux*wx-ux*vx) / 6.d0 + & 1230 vec_tab (i1) * (+vx*wx-px*wx-px*vx) / 2.d0 - & 1231 vec_tab (i2) * (+ux*wx-px*wx-px*ux) / 2.d0 + & 1232 vec_tab (i3) * (+ux*vx-px*vx-px*ux) / 6.d0 1233 endif 1234 enddo 1235 1236 deallocate (q) 1237 deallocate ( gk ) 1238 if (spline_ps) deallocate(xdata) 1239 return 1240end subroutine gen_us_djl 1241! 1242!---------------------------------------------------------------------- 1243subroutine gen_us_vkb0 (ik,npw,vkb0,size_tab,vec_tab, spline_ps, vec_tab_d2y) 1244 !---------------------------------------------------------------------- 1245 ! 1246 ! Calculates the kleinman-bylander pseudopotentials with the 1247 ! derivative of the spherical harmonics projected on vector u 1248 ! 1249 USE kinds, ONLY : DP 1250 USE io_global, ONLY : stdout 1251 USE constants, ONLY : tpi 1252 USE cell_base, ONLY : tpiba 1253 USE klist, ONLY : xk, igk_k 1254 USE gvect, ONLY : g 1255 USE us, ONLY : nqx, dq 1256 USE splinelib, ONLY : splint 1257 USE uspp_param, ONLY : upf 1258 ! 1259 implicit none 1260 ! 1261 integer, intent(in) :: ik, npw 1262 real(DP), intent(inout) ::vkb0(1:npw) 1263 integer, intent(in) :: size_tab 1264 real(DP), intent(in) :: vec_tab(1:size_tab) 1265 real(DP), intent(in) :: vec_tab_d2y(1:size_tab) 1266 logical :: spline_ps 1267 ! 1268 integer :: na, nt, nb, ikb,i0, i1, i2, & 1269 i3, ig 1270 real(DP), allocatable :: gk(:,:), q (:) 1271 real(DP) :: px, ux, vx, wx 1272 complex(DP), allocatable :: sk (:) 1273 1274 integer :: iq 1275 real(DP), allocatable :: xdata(:) 1276 1277 allocate ( gk(3,npw) ) 1278 allocate ( q(npw) ) 1279 1280 do ig = 1, npw 1281 gk (1, ig) = xk (1, ik) + g (1, igk_k (ig,ik) ) 1282 gk (2, ig) = xk (2, ik) + g (2, igk_k (ig,ik) ) 1283 gk (3, ig) = xk (3, ik) + g (3, igk_k (ig,ik) ) 1284 q (ig) = gk(1, ig)**2 + gk(2, ig)**2 + gk(3, ig)**2 1285 enddo 1286 1287 do ig = 1, npw 1288 q (ig) = sqrt ( q(ig) ) * tpiba 1289 end do 1290 1291 if (spline_ps) then 1292 allocate(xdata(nqx)) 1293 do iq = 1, nqx 1294 xdata(iq) = (iq - 1) * dq 1295 enddo 1296 endif 1297 1298 ! calculate beta in G-space using an interpolation table 1299 do ig = 1, npw 1300 if (spline_ps) then 1301 vkb0(ig) = splint(xdata, vec_tab(:), & 1302 vec_tab_d2y(:), q(ig)) 1303 else 1304 px = q (ig) / dq - int (q (ig) / dq) 1305 ux = 1.d0 - px 1306 vx = 2.d0 - px 1307 wx = 3.d0 - px 1308 i0 = q (ig) / dq + 1 1309 i1 = i0 + 1 1310 i2 = i0 + 2 1311 i3 = i0 + 3 1312 vkb0 (ig) = vec_tab (i0) * ux * vx * wx / 6.d0 + & 1313 vec_tab (i1) * px * vx * wx / 2.d0 - & 1314 vec_tab (i2) * px * ux * wx / 2.d0 + & 1315 vec_tab (i3) * px * ux * vx / 6.d0 1316 endif 1317 enddo 1318 1319 deallocate (q) 1320 deallocate ( gk ) 1321 if (spline_ps) deallocate(xdata) 1322 return 1323end subroutine gen_us_vkb0 1324