1!=============================================================================== 2! 3! Routines: 4! 5! (1) inread() Originally By ? Last Modified 7/8/2008 (JRD) 6! 7! Reads parameters from sigma.inp for the current job. 8! 9!=============================================================================== 10 11#include "f_defs.h" 12 13subroutine inread(sig) 14 15 use global_m 16 use scissors_m 17 use inread_common_m 18 implicit none 19 20 type (siginfo), intent(out) :: sig 21 22 character*256 :: blockword,keyword,line,errmsg,tmpstr 23 integer :: ii,jj,kk,nkpt_read,nqpt_read,ndiag_read,noff_read 24 integer :: nbnmin,nbnmax,itestq,iostat 25 26 real(DP) :: kpt_read(3,MAX_KPTS),qpt_read(3,MAX_KPTS),div,max_freq_eval 27 integer :: off1_read(MAX_BANDS),off2_read(MAX_BANDS),off3_read(MAX_BANDS) 28 integer :: band_occ(MAX_BANDS),diag(MAX_BANDS),spinmin,spinmax 29 logical :: occ_set, q0vec_read, VXC_exists, found, old_grid_set 30 31 integer :: ik ! ZL: for loop generating k_phonq 32 real(DP), allocatable :: k_phonq_read(:) 33 34 PUSH_SUB(inread) 35 36#ifdef MPI 37 ! Non-root nodes should wait for root to read the whole file. 38 ! That way, we can be sure root gets a chance to write errors before 39 ! any die call is issued by another node. Root calls MPI_Barrier below. 40 if(peinf%inode /= 0) call MPI_Barrier(MPI_COMM_WORLD, mpierr) 41#endif 42 43!---------------------- 44! Initialize sig%igamma 45 46 sig%igamma = 0 47 48!---------------------- 49! Initialize wcoul0 50 sig%wcoul0 = ZERO 51 52 53!--------------------------------- 54! Set default values 55 56 occ_set = .false. 57 q0vec_read = .false. 58 sig%invalid_gpp_mode = -1 59 nkpt_read=0 60 nqpt_read=0 61 noff_read=0 62 band_occ=0 63 nbnmin=0 64 nbnmax=0 65 spinmin=1 66 spinmax=1 67 diag=0 68#ifdef HDF5 69 sig%use_hdf5 = .true. 70#else 71 sig%use_hdf5 = .false. 72#endif 73 sig%wfn_hdf5 = .false. 74 sig%loff=0 75 sig%toff=0 76 sig%freq_dep=1 77 sig%freq_dep_method=2 78 sig%nFreq=1 79 sig%exact_ch=-1 80 sig%freq_grid_shift=2 81 sig%freqevalmin=0d0 82 old_grid_set=.false. 83 sig%freqevalstep=0.2d0 84 max_freq_eval=2d0 85 sig%nfreqeval=1 86 sig%iuseremainder=0 87 sig%nkn=0 88 sig%nq=0 89 sig%nq0=0 90 sig%nq1=0 91 sig%subsample=.false. 92 sig%elph=.false. 93 sig%ep_bracket=0 94 sig%nphonq = 0 95 sig%ndiag=0 96 sig%noffdiag=0 97 sig%nspin=0 98 sig%nfreq_imag=0 99 sig%cd_int_method=0 100 call scissors_zero(sig%scis) 101 call scissors_zero(sig%scis_outer) 102 sig%spl_tck%n=0 103 sig%spl_tck_outer%n=0 104 sig%avgpot=0.0d0 105 sig%avgpot_outer=0.0d0 106 sig%tol = TOL_Degeneracy 107 sig%use_xdat=.false. 108 sig%use_vxcdat=.TRUE. 109 sig%use_vxc2dat=.TRUE. 110 sig%is_EXX=.FALSE. 111 sig%ecuts=0.d0 112 sig%ecutb=0.d0 113 sig%ntband=0 ! ZL: number of total bands used in sigma 114 sig%ncore_excl=0 115 sig%nvband=0 116 sig%fdf=2 117 sig%dw=1.0d0 118 sig%xfrac=1.0d0 119 sig%gamma=0.5d0 120 sig%sexcut=4.0d0 121 sig%icutv=TRUNC_NONE 122 sig%iscreen=SCREEN_SEMICOND 123 sig%iwritecoul=0 124 sig%truncval(:)=0.0d0 125 sig%ncrit=0 126 sig%efermi_input=0.0d0 127 sig%rfermi=.true. 128 sig%fullConvLog=0 129 sig%qgrid(1:3)=0 130 sig%avgcut=-1 !FHJ: <0 means "auto", see code below 131 sig%freplacebz=.false. 132 sig%fwritebz=.false. 133 sig%degeneracy_check_override=.false. 134 sig%offdiagsym=.true. 135 sig%qgridsym=.true. 136 sig%die_outside_sphere=.true. 137 sig%averagew=.true. 138 peinf%npools=0 139 sig%eqp_corrections=.false. 140 sig%eqp_outer_corrections=.false. 141 sig%coulomb_mod%short_range_frac_fock=1.0d0 142 sig%coulomb_mod%long_range_frac_fock=1.0d0 143 sig%coulomb_mod%screening_length=0.0d0 144 sig%coul_mod_flag=.false. 145 sig%sigma_correction=.false. 146 sig%symmetrize=.true. 147 sig%do_sigma_subspace=.false. 148 sig%sub_collective_redistr=.false. 149 150!--------------------------------- 151! Never-ending loop... 152 153 do while(0.eq.0) 154 155! Actually the loop ends when the end of the file is reached 156 157 read(55,'(a256)',iostat=iostat) line 158 if(iostat < 0) exit 159 160! Skip comment lines 161 162 if(len_trim(line).eq.0) cycle 163 if(line(1:1).eq.'#') cycle 164 165! Determine keyword 166 167 keyword=line(1:scan(line," ")-1) 168 line=adjustl(line(scan(line," ")+1:256)) 169 170 if (trim(keyword).eq.'begin') then 171 blockword=line(1:scan(line," ")-1) 172 ii=0 173 do while(trim(line).ne.'end') 174 read(55,'(a256)',iostat=iostat) line 175 if(iostat /= 0) then 176 write(errmsg,'(3a)') 'The end of the file was reached while reading elements of the ', & 177 trim(blockword),' block.' 178 call die(errmsg, only_root_writes = .true.) 179 endif 180 181 if(trim(line).ne.'end') then 182 ii=ii+1 183 if(trim(blockword).eq.'kpoints') then 184 call check_bounds_nkq(ii, 'k', 'begin kpoints') 185 read(line,*,err=112) (kpt_read(jj,ii),jj=1,3),div 186 kpt_read(1:3,ii)=kpt_read(1:3,ii)/div 187 elseif(trim(blockword).eq.'qpoints') then 188 call check_bounds_nkq(ii, 'q', 'begin qpoints') 189 read(line,*,err=112) (qpt_read(jj,ii),jj=1,3),div,itestq 190 qpt_read(1:3,ii)=qpt_read(1:3,ii)/div 191 ! for Hartree-Fock 192 if (itestq.eq.1) then 193 sig%nq0 = sig%nq0 + 1 194 if (.not.q0vec_read) then 195 sig%q0vec(1:3)=qpt_read(1:3,ii) 196 q0vec_read = .true. 197 endif 198 endif 199 elseif(trim(blockword).eq.'diag') then 200 call check_bounds_nbands(ii, 'begin diag') 201 read(line,*,err=112) diag(ii) 202 elseif(trim(blockword).eq.'offdiag') then 203 call check_bounds_nbands(ii, 'begin offdiag') 204 read(line,*,err=112) off1_read(ii),off2_read(ii),off3_read(ii) 205 else 206 write(errmsg,'(3a)') 'Unexpected blockword ', trim(blockword), ' was found in sigma.inp.' 207 call die(errmsg, only_root_writes = .true.) 208 end if 209 end if 210 end do 211 if(trim(blockword).eq.'kpoints') then 212 nkpt_read=ii 213 elseif(trim(blockword).eq.'qpoints') then 214 nqpt_read=ii 215 elseif(trim(blockword).eq.'diag') then 216 ndiag_read=ii 217 elseif(trim(blockword).eq.'offdiag') then 218 noff_read=ii 219 endif 220 221! Spline scissors 222 elseif(trim(keyword).eq.'spline_scissors') then 223 if (peinf%inode==0) then 224 write(6,*) 'Reading spline coefficients from `spline_scissors.dat`' 225 write(6,*) 226 ! read number of pts, knots, coefficients, degree 227 call open_file(20,file='spline_scissors.dat',form='formatted',status='old') 228 read(20,*) sig%spl_tck%n 229 SAFE_ALLOCATE(sig%spl_tck%t, (sig%spl_tck%n)) 230 SAFE_ALLOCATE(sig%spl_tck%c, (sig%spl_tck%n)) 231 read(20,*) (sig%spl_tck%t(ii), ii=1,sig%spl_tck%n) 232 read(20,*) (sig%spl_tck%c(ii), ii=1,sig%spl_tck%n) 233 read(20,*) sig%spl_tck%k 234 call close_file(20) 235 endif 236 237! Spline scissors, outer 238 elseif(trim(keyword).eq.'spline_scissors_outer') then 239 if (peinf%inode==0) then 240 write(6,*) 'Reading outer spline coefficients from `spline_scissors_outer.dat`' 241 write(6,*) 242 ! read number of pts, knots, coefficients, degree 243 call open_file(20,file='spline_scissors_outer.dat',form='formatted',status='old') 244 read(20,*) sig%spl_tck_outer%n 245 SAFE_ALLOCATE(sig%spl_tck_outer%t, (sig%spl_tck_outer%n)) 246 SAFE_ALLOCATE(sig%spl_tck_outer%c, (sig%spl_tck_outer%n)) 247 read(20,*) (sig%spl_tck_outer%t(ii), ii=1,sig%spl_tck_outer%n) 248 read(20,*) (sig%spl_tck_outer%c(ii), ii=1,sig%spl_tck_outer%n) 249 read(20,*) sig%spl_tck_outer%k 250 call close_file(20) 251 endif 252 253 elseif(trim(keyword).eq.'verbosity') then 254 read(line,*,err=110) peinf%verbosity 255! The average potential on the faces of the unit cell in the non-periodic directions for the bands in WFN_inner 256 elseif(trim(keyword).eq.'avgpot') then 257 read(line,*,err=110) sig%avgpot 258! The average potential on the faces of the unit cell in the non-periodic directions for the bands in WFN_outer 259 elseif(trim(keyword).eq.'avgpot_outer') then 260 read(line,*,err=110) sig%avgpot_outer 261! Frequency dependence of the inverse dielectric matrix 262 elseif(trim(keyword).eq.'frequency_dependence') then 263 read(line,*,err=110) sig%freq_dep 264! Is calculation only of bare exchange for EXX calculation? If so, don`t need VXC or vxc.dat 265 elseif(trim(keyword).eq.'for_EXX') then 266 sig%is_EXX = .TRUE. 267 elseif(trim(keyword).eq.'frequency_dependence_method') then 268 read(line,*,err=110) sig%freq_dep_method 269 elseif(trim(keyword).eq.'cd_integration_method') then 270 read(line,*,err=110) sig%cd_int_method 271! Frequency grid for numerical integration of full-frequency sigma (may be different from Epsilon) 272 elseif(trim(keyword).eq.'init_frequency_eval') then 273 read(line,*,err=110) sig%freqevalmin 274 old_grid_set = .true. 275 elseif(trim(keyword).eq.'delta_frequency_eval') then 276 read(line,*,err=110) sig%freqevalstep 277 elseif(trim(keyword).eq.'number_frequency_eval') then 278 read(line,*,err=110) sig%nfreqeval 279 old_grid_set = .true. 280 elseif(trim(keyword).eq.'max_frequency_eval') then 281 read(line,*,err=110) max_freq_eval 282 elseif(trim(keyword).eq.'frequency_grid_shift') then 283 read(line,*,err=110) sig%freq_grid_shift 284! Use full frequency tail for better energies 285 elseif(trim(keyword).eq.'use_epsilon_remainder') then 286 sig%iuseremainder = 1 287! Grid of Qs in Epsilon 288 elseif(trim(keyword).eq.'qgrid') then 289 read(line,*,err=110) sig%qgrid(1), sig%qgrid(2), sig%qgrid(3) 290! Compute the exact static CH 291 elseif(trim(keyword).eq.'exact_static_ch') then 292 read(line,*,err=110) sig%exact_ch 293! Full CH convergence logging of all calculated bands 294 elseif(trim(keyword).eq.'full_ch_conv_log') then 295 read(line,*,err=110) sig%fullConvLog 296! Use use_xdat to skip the computation of bare exchange 297 elseif(trim(keyword).eq.'use_xdat') then 298 sig%use_xdat = .true. 299! Use dont_use_vxcdat to compute exchange-correlation matrix elements from VXC 300 elseif(trim(keyword).eq.'dont_use_vxcdat') then 301 sig%use_vxcdat = .false. 302 elseif(trim(keyword).eq.'dont_use_vxc2dat') then 303 sig%use_vxc2dat = .false. 304 elseif(trim(keyword).eq.'dont_use_hdf5') then 305 sig%use_hdf5 = .false. 306 elseif(trim(keyword).eq.'number_kpoints') then 307 read(line,*,err=110) sig%nkn ! FHJ: deprecated 308 elseif(trim(keyword).eq.'number_qpoints') then 309 read(line,*,err=110) sig%nq ! FHJ: deprecated 310 elseif(trim(keyword).eq.'dont_symmetrize') then 311 sig%symmetrize = .false. 312 elseif(trim(keyword).eq.'subsample') then 313 sig%subsample = .true. 314 elseif(trim(keyword).eq.'number_sigma_pools') then 315 read(line,*,err=110) peinf%npools 316 elseif(trim(keyword).eq.'bare_coulomb_cutoff') then 317 read(line,*,err=110) sig%ecutb 318 elseif(trim(keyword).eq.'cell_average_cutoff') then 319 read(line,*,err=110) sig%avgcut 320 elseif(trim(keyword).eq.'screened_coulomb_cutoff') then 321 read(line,*,err=110) sig%ecuts 322 elseif(trim(keyword).eq.'number_bands') then 323 read(line,*,err=110) sig%ntband 324 call check_bounds_nbands(sig%ntband, 'number_bands') 325 elseif(trim(keyword).eq.'number_core_excluded') then 326 read(line,*,err=110) sig%ncore_excl 327 elseif(trim(keyword).eq.'band_index_min') then 328 read(line,*,err=110) nbnmin 329 elseif(trim(keyword).eq.'band_index_max') then 330 read(line,*,err=110) nbnmax 331 elseif(trim(keyword).eq.'spin_index_min') then 332 read(line,*,err=110) spinmin 333 elseif(trim(keyword).eq.'spin_index_max') then 334 read(line,*,err=110) spinmax 335 elseif(trim(keyword).eq.'finite_difference_form') then 336 read(line,*,err=110) sig%fdf 337 elseif(trim(keyword).eq.'finite_difference_spacing') then 338 read(line,*,err=110) sig%dw 339 elseif(trim(keyword).eq.'bare_exchange_fraction') then 340 read(line,*,err=110) sig%xfrac 341 elseif(trim(keyword).eq.'short_range_frac_fock') then 342 read(line,*,err=110) sig%coulomb_mod%short_range_frac_fock 343 elseif(trim(keyword).eq.'long_range_frac_fock') then 344 read(line,*,err=110) sig%coulomb_mod%long_range_frac_fock 345 elseif(trim(keyword).eq.'screening_length') then 346 read(line,*,err=110) sig%coulomb_mod%screening_length 347 elseif(trim(keyword).eq.'write_vcoul') then 348 sig%iwritecoul=1 349 elseif(trim(keyword).eq.'tol_degeneracy') then 350 read(line,*,err=110) sig%tol 351 elseif(trim(keyword).eq.'gpp_broadening') then 352 read(line,*,err=110) sig%gamma 353 elseif(trim(keyword).eq.'gpp_sexcutoff') then 354 read(line,*,err=110) sig%sexcut 355 elseif(trim(keyword).eq.'number_diag') then 356 read(line,*,err=110) sig%ndiag 357 elseif(trim(keyword).eq.'number_offdiag') then 358 read(line,*,err=110) sig%noffdiag 359 elseif(trim(keyword).eq.'sigma_matrix') then 360 read(line,*,err=110) sig%loff, sig%toff 361 elseif(trim(keyword).eq.'band_occupation') then 362 read(line,*,err=110) (band_occ(ii),ii=1,sig%ntband) 363 occ_set = .true. 364 elseif(trim(keyword).eq.'number_partial_occup') then 365 read(line,*,err=110) sig%ncrit 366 occ_set = .true. 367 elseif(trim(keyword).eq.'fermi_level') then 368 read(line,*,err=110) sig%efermi_input 369 elseif(trim(keyword).eq.'fermi_level_absolute') then 370 sig%rfermi=.false. 371 elseif(trim(keyword).eq.'fermi_level_relative') then 372 sig%rfermi=.true. 373 elseif(trim(keyword).eq.'fullbz_replace') then 374 sig%freplacebz=.true. 375 elseif(trim(keyword).eq.'fullbz_write') then 376 sig%fwritebz=.true. 377 elseif(trim(keyword).eq.'degeneracy_check_override') then 378 sig%degeneracy_check_override=.true. 379 elseif(trim(keyword).eq.'no_symmetries_offdiagonals') then 380 sig%offdiagsym=.false. 381 elseif(trim(keyword).eq.'no_symmetries_q_grid') then 382 sig%qgridsym=.false. 383 elseif(trim(keyword).eq.'use_symmetries_q_grid') then 384 sig%qgridsym=.true. 385 elseif(trim(keyword).eq.'die_outside_sphere') then 386 sig%die_outside_sphere=.true. 387 elseif(trim(keyword).eq.'ignore_outside_sphere') then 388 sig%die_outside_sphere=.false. 389 elseif(trim(keyword).eq.'eqp_corrections') then 390 sig%eqp_corrections=.true. 391 elseif(trim(keyword).eq.'eqp_outer_corrections') then 392 sig%eqp_outer_corrections=.true. 393 elseif(trim(keyword).eq.'skip_averagew') then 394 sig%averagew = .false. 395 elseif(trim(keyword).eq.'do_sigma_subspace') then 396 sig%do_sigma_subspace = .true. 397 elseif(trim(keyword).eq.'sub_collective_redistr') then 398 sig%sub_collective_redistr = .true. 399 elseif(trim(keyword).eq.'invalid_gpp_mode') then 400 read(line,*,err=110) sig%invalid_gpp_mode 401 elseif(try_inread_truncation(trim(keyword), trim(line), sig%icutv, sig%truncval(1))) then 402 ! subroutine already does the job 403 elseif(try_inread_screening(trim(keyword), trim(line), sig%iscreen)) then 404 ! subroutine already does the job 405 else 406 call scissors_inread(keyword, line, sig%scis, found) 407 if(.not. found) call scissors_inread(keyword, line, sig%scis_outer, found, "_outer") 408 if(.not. found) then 409 write(errmsg,'(3a)') 'Unexpected keyword ', trim(keyword), ' was found in sigma.inp.' 410 call die(errmsg, only_root_writes = .true.) 411 endif 412 end if 413 enddo 414 415! for the moment subspace works only in combination with MPI and SCALAPACK 416#if !defined MPI || !defined USESCALAPACK 417 if(sig%do_sigma_subspace) then 418 if (peinf%inode==0) then 419 write(0,*) 420 write(0,*) 'WARNING: Static Subspace method only works with MPI and SCALAPACK.' 421 write(0,*) 'Subspace method turned off.' 422 write(0,*) 423 endif 424 sig%do_sigma_subspace = .false. 425 end if 426#endif 427 428 ! entered in Ryd, stored in eV since kp%el and kp%elda are soon converted to eV. 429 sig%tol = sig%tol * RYD 430 431 call peinfo_set_verbosity() 432 if (peinf%inode==0 .and. sig%nkn>0) then 433 write(0,'(/,a)') 'WARNING: the `number_kpoints` flag is deprecated. The code now' 434 write(0,'(a,/)') 'automatically figures out the number of q-points from the input.' 435 endif 436 sig%nkn = nkpt_read 437 if (peinf%inode==0 .and. sig%nq>0) then 438 write(0,'(/,a)') 'WARNING: the `number_qpoints` flag is deprecated. The code now' 439 write(0,'(a,/)') 'automatically figures out the number of k-points from the input.' 440 endif 441 sig%nq = nqpt_read 442 if(sig%freq_dep==-1 .and. sig%nq<1) then 443 write(0,'(/,a)') 'ERROR: Hartree-Fock calculations require a list of q-points' 444 call die('The `begin qpoints` block could not be found.', only_root_writes=.true.) 445 endif 446 sig%nq1 = sig%nq - sig%nq0 447 448 if(any(sig%qgrid(1:3) == 0).and.sig%freq_dep.eq.-1) then 449 call die('qgrid must be specified for Hartree-Fock.', only_root_writes = .true.) 450 endif 451 452 if(.not. q0vec_read .and. sig%freq_dep .eq. -1) then 453 call die('No q0 specified in qpoints block.', only_root_writes = .true.) 454 endif 455 456 SAFE_ALLOCATE(sig%kpt, (3,sig%nkn)) 457 sig%kpt(1:3,1:sig%nkn)=kpt_read(1:3,1:sig%nkn) 458 459 if (sig%nq>0) then 460 if (sig%freq_dep/=-1) then 461 if (peinf%inode==0) & 462 write(0,'(/,a,/)') 'WARNING: ignoring the qpoints block for calculations other than HF.' 463 else 464 SAFE_ALLOCATE(sig%qpt, (3,sig%nq)) 465 sig%qpt(1:3,1:sig%nq)=qpt_read(1:3,1:sig%nq) 466 endif 467 endif 468 469 if(sig%ndiag.ne.0) then 470 SAFE_ALLOCATE(sig%diag, (sig%ndiag)) 471 sig%diag(1:sig%ndiag)=diag(1:sig%ndiag) 472 endif 473 474 if(sig%ndiag.eq.0.and.nbnmin*nbnmax.ne.0) then 475 sig%ndiag= nbnmax - nbnmin + 1 476 ndiag_read = sig%ndiag 477 SAFE_ALLOCATE(sig%diag, (sig%ndiag)) 478 do ii=1,sig%ndiag 479 sig%diag(ii)= ii + nbnmin - 1 480 enddo 481 endif 482 483 if(sig%ndiag.eq.0.and.nbnmin*nbnmax.eq.0) then 484 write(errmsg,'(a,3i6)') 'Incomprehensible list of energy bands: ', sig%ndiag, nbnmin, nbnmax 485 call die(errmsg, only_root_writes = .true.) 486 endif 487 488 if(ndiag_read.lt.sig%ndiag) then 489 if(peinf%inode.eq.0) then 490 write(0,*) 'The number of diagonal elements found in the diag block (',ndiag_read,')' 491 write(0,*) ' is smaller than the one specified by the keyword number_diag (',sig%ndiag,').' 492 endif 493 call die("ndiag too small", only_root_writes = .true.) 494 endif 495 496 ! no finite-difference evaluation unless GPP 497 if(sig%freq_dep /= 1 .and. sig%freq_dep /= 3) sig%fdf=-2 498 499 if(ndiag_read.gt.sig%ndiag) then 500 if(peinf%inode.eq.0) then 501 write(0,887) ndiag_read, sig%ndiag, sig%ndiag 502 endif 503 endif 504887 format(1x,"WARNING: The number of diag elements in the diag block (",i4,") is larger",/,& 505 3x,"than the one specified by the keyword number_diag (",i4,").",/,& 506 3x,"The program will continue using only",i4,1x,"diagonal elements.",/) 507 508 if(nbnmin*nbnmax.eq.0) then 509 nbnmin=max_bands 510 nbnmax=-max_bands 511 do ii=1,sig%ndiag 512 if (sig%diag(ii).lt.nbnmin) nbnmin=sig%diag(ii) 513 if (sig%diag(ii).gt.nbnmax) nbnmax=sig%diag(ii) 514 enddo 515 endif 516 517 sig%bmin=nbnmin 518 sig%bmax=nbnmax 519 520 if(sig%loff.lt.-2.or.(sig%loff.gt.0.and.sig%loff.lt.nbnmin).or. & 521 sig%loff.gt.nbnmax.or.sig%toff.lt.-1.or.sig%toff.gt.1) then 522 call die("sigma_matrix parameters out of range", only_root_writes = .true.) 523 endif 524 525 if(sig%loff.ne.0) then 526 if(noff_read.ne.0.or.sig%noffdiag.gt.0) then 527 if(peinf%inode.eq.0) then 528 write(0,994) 529 endif 530 endif 531 kk=0 532 do ii=nbnmin,nbnmax 533 do jj=nbnmin,nbnmax 534 if (sig%toff.eq.-1.and.ii.lt.jj) cycle 535 if (sig%toff.eq.1.and.ii.gt.jj) cycle 536 kk=kk+1 537 off1_read(kk)=ii 538 off2_read(kk)=jj 539 if (sig%loff.eq.-1) then 540 off3_read(kk)=ii 541 else if (sig%loff.eq.-2) then 542 off3_read(kk)=jj 543 else 544 off3_read(kk)=sig%loff 545 endif 546 enddo 547 enddo 548 sig%noffdiag=kk 549 noff_read=sig%noffdiag 550 endif 551994 format(1x,"WARNING: both offdiag and sigma_matrix found in sigma.inp.",/,& 552 3x,"The latter overrides the former.",/) 553 554 if(noff_read.lt.sig%noffdiag) then 555 if(peinf%inode.eq.0) then 556 write(0,997) noff_read, sig%noffdiag 557 endif 558 call die("offdiag error", only_root_writes = .true.) 559 endif 560997 format(1x,"The number of offdiag elements in the offdiag block (",i4,") is",/,& 561 3x,"smaller than the one specified by the keyword number_offdiag (",i4,").") 562 563 if(noff_read.gt.sig%noffdiag) then 564 if(peinf%inode.eq.0) then 565 write(0,996) noff_read, sig%noffdiag, sig%noffdiag 566 endif 567 endif 568996 format(1x,"WARNING: The number of offdiag elements in the offdiag block (",i4,") is",/,& 569 3x,"larger than the one specified by the keyword number_offdiag (",i4,").",/,& 570 3x,"The program will continue using only",i4,1x,"offdiag elements.",/) 571 572 ! allocate even if noffdiag = 0 since sunf90 -xcheck will complain when passed unallocated to read_matrix_elements 573 SAFE_ALLOCATE(sig%off1, (sig%noffdiag)) 574 SAFE_ALLOCATE(sig%off2, (sig%noffdiag)) 575 SAFE_ALLOCATE(sig%off3, (sig%noffdiag)) 576 if(sig%noffdiag.gt.0) then 577 sig%off1(1:sig%noffdiag)=off1_read(1:sig%noffdiag) 578 sig%off2(1:sig%noffdiag)=off2_read(1:sig%noffdiag) 579 sig%off3(1:sig%noffdiag)=off3_read(1:sig%noffdiag) 580 endif 581 582 if (peinf%inode==0 .and. (sig%ecuts>TOL_ZERO .or. sig%ecutb>TOL_ZERO)) then 583 write(6,'(/1x,a/)') 'NOTE: `screened_coulomb_cutoff` and `bare_coulomb_cutoff` are now optional flags.' 584 endif 585 586 if(spinmin.eq.1.and.spinmax.eq.1) then 587 sig%nspin=1 588 sig%spin_index(1)=1 589 sig%spin_index(2)=2 ! only used for spinor case 590 elseif(spinmin.eq.2.and.spinmax.eq.2) then 591 sig%nspin=1 592 sig%spin_index(1)=2 593 elseif(spinmin.eq.1.and.spinmax.eq.2) then 594 sig%nspin=2 595 sig%spin_index(1)=1 596 sig%spin_index(2)=2 597 else 598 write(errmsg,'(a,i2,a,i2,a,i2)') 'Illegal range of spin indices from ', spinmin, ' to ', spinmax 599 call die(errmsg, only_root_writes = .true.) 600 endif 601 602!------------------------------ 603! Build the map sig%offmap 604! sig%off1(ii) = sig%diag(sig%offmap(ii,1)) 605! sig%off2(ii) = sig%diag(sig%offmap(ii,2)) 606! sig%off3(ii) = sig%diag(sig%offmap(ii,3)) 607! This is done only if sig%noffdiag .ne. 0 608 609 if(sig%noffdiag.gt.0) then 610 SAFE_ALLOCATE(sig%offmap, (sig%noffdiag,3)) 611 sig%offmap(1:sig%noffdiag, 1:3) = 0 612 do ii=1,sig%noffdiag 613 do jj=1,sig%ndiag 614 if(sig%diag(jj) .eq. sig%off1(ii)) sig%offmap(ii, 1) = jj 615 if(sig%diag(jj) .eq. sig%off2(ii)) sig%offmap(ii, 2) = jj 616 if(sig%diag(jj) .eq. sig%off3(ii)) sig%offmap(ii, 3) = jj 617 enddo 618 619 if(sig%offmap(ii, 1) == 0) then 620 write(errmsg,'(a,i6,a)') 'Off-diagonal matrix el. requested for band ', & 621 sig%off1(ii),' but no corresponding diagonal matrix el. is requested' 622 call die(errmsg, only_root_writes = .true.) 623 endif 624 625 if(sig%offmap(ii, 2) == 0) then 626 write(errmsg,'(a,i6,a)') 'Off-diagonal matrix el. requested for band ', & 627 sig%off2(ii),' but no corresponding diagonal matrix el. is requested' 628 call die(errmsg, only_root_writes = .true.) 629 endif 630 631 if(sig%offmap(ii, 3) == 0) then 632 write(errmsg,'(a,i6,a)') 'Off-diagonal matrix el. requested at energy of band ', & 633 sig%off3(ii),' but no corresponding diagonal matrix el. is requested' 634 call die(errmsg, only_root_writes = .true.) 635 endif 636 enddo 637 endif ! sig%noffdiag.gt.0 638 639 if (occ_set) then 640 if (peinf%inode==0) then 641 write(0,'(/1x,a)') 'WARNING: keywords `number_partial_occup` and `band_occupations` are deprecated.' 642 write(0,'(1x,a/)') 'BerkeleyGW now figures out these parameters automatically.' 643 endif 644 645 sig%nvband = count(band_occ==1) 646 647 if ((sig%nvband+sig%ncrit)==0) & 648 call die("There are no occupied or partially occupied bands.", only_root_writes = .true.) 649 650 if (any(band_occ/=0 .and. band_occ/=1) .and. & 651 any(band_occ(2:sig%ntband)>band_occ(1:sig%ntband-1))) then 652 ! FHJ: non-equilibrium occ completely disabled. Go change your mean-field! 653 call die("Non-equilibrium occupations not supported.", only_root_writes = .true.) 654 endif 655 endif 656 657 if(sig%fullConvLog<0.or.sig%fullConvLog>1) then 658 call die('Invalid full_ch_conv_log', only_root_writes = .true.) 659 endif 660 661! gsm: What frequency dependence we are using? 662 663 if(peinf%inode.eq.0) then 664 if(sig%freq_dep.eq.-1) then 665 write(6,699) 666 elseif(sig%freq_dep.eq.0) then 667 write(6,700) 668 elseif(sig%freq_dep.eq.1) then 669 write(6,701) 670 elseif(sig%freq_dep.eq.2) then 671 write(6,702) 672 elseif(sig%freq_dep.eq.3) then 673 write(6,703) 674 else 675 call die('Need to specify frequency dependence') 676 endif 677 endif 678699 format(1x,'Using the Hartree-Fock approximation',/) 679700 format(1x,'Using the static COHSEX approximation',/) 680701 format(1x,'Using the Generalized Plasmon Pole model',/) 681702 format(1x,'Using the full frequency-dependent inverse dielectric matrix',/) 682703 format(1x,'Using the Generalized Plasmon Pole model (GN flavor)',/) 683 684 if(peinf%inode == 0) then 685 if(peinf%npes > 1) then 686 write(6,803) 687 else 688 write(6,805) 689 endif 690 endif 691803 format(1x,'We are communicating via MPI',/) 692805 format(1x,'We are not communicating',/) 693 694! gsm: Do we compute the exact static CH? 695 696 if(sig%exact_ch.eq.-1) then ! set default according to freq_dep 697 if(sig%freq_dep .eq. 0) then 698 sig%exact_ch = 1 699 else 700 sig%exact_ch = 0 701 endif 702 endif 703 if(sig%freq_dep.ne.-1) then ! HF has no CH sum so don`t write any comments about it 704 if(sig%exact_ch.eq.0) then 705 if(peinf%inode.eq.0) write(6,750) 706 elseif(sig%exact_ch.eq.1) then 707 if(sig%freq_dep.eq.0) then 708 if(peinf%inode.eq.0) write(6,751) 709 else 710 if(peinf%inode.eq.0) write(6,752) 711 endif 712 else 713 call die('Illegal value for exact_static_ch') 714 endif 715 endif 716750 format(1x,'Computing CH as a partial sum over empty bands',/) 717751 format(1x,'Computing the exact static CH',/) 718752 format(1x,'Computing CH as a partial sum over empty bands with the static remainder',/) 719 720 ! JRD: What screening is present? 721 if(peinf%inode.eq.0) then 722 if(sig%ncrit < 0) then 723 call die("number_partial_occup < 0") 724 endif 725 726 select case (sig%iscreen) 727 case (SCREEN_SEMICOND) 728 if(sig%ncrit > 0) then 729 write(0,'(a)') "WARNING: Semiconductor screening is inappropriate for number_partial_occup > 0." 730 write(0,'(a)') "If a band really crosses the Fermi level, graphene or metal screening must be used." 731 ! this only makes sense if ncrit was used to handle two spins, or 732 ! to tune the number of conduction and valence bands for optimal parallelization. 733 endif 734 write(6,'(1x,a/)') 'Running with semiconductor screening' 735 case (SCREEN_GRAPHENE) 736 if(sig%ncrit == 0) then 737 write(0,*) "WARNING: Graphene screening usually should have number_partial_occup > 0." 738 endif 739 write(6,'(1x,a/)') 'Running with graphene screening' 740 case (SCREEN_METAL) 741 if(sig%ncrit == 0) then 742 write(0,*) "WARNING: Metal screening usually should have number_partial_occup > 0." 743 endif 744 write(6,'(1x,a/)') 'Running with metal screening' 745 case default 746 call die('Unknown screening type', only_root_writes=.true.) 747 endselect 748 endif 749 750 call print_truncation_summary(sig%icutv, sig%truncval(1)) 751 752 ! FHJ: set cell_average_cutoff, if not set by user 753 if (sig%avgcut<0) then 754 sig%avgcut = TOL_ZERO 755 if (sig%icutv==TRUNC_NONE .and. sig%iscreen==SCREEN_SEMICOND) then 756 sig%avgcut = 1d0/TOL_ZERO 757 endif 758 endif 759 if (peinf%inode==0) then 760 write(6,'(1x,a,es10.3e3,a/)') & 761 'Cutoff for Monte-Carlo average of Coulomb potential: ', sig%avgcut, ' Ry' 762 endif 763 764 if(peinf%inode.eq.0) then 765 if(sig%ncrit.gt.0) then 766 write(6,951)sig%ncrit 767951 format(1x,'We have partially occupied bands',/, & 768 1x,'number_partial_occup (ncrit) =',i3,/) 769 endif 770 endif 771 772 if(sig%qgridsym .and. .not. sig%symmetrize) then 773 write(errmsg,'(a,i6,a)') 'Must use no_symmetries_q_grid flag with dont_symmetrize flag' 774 call die(errmsg, only_root_writes = .true.) 775 endif 776 777 778 if(sig%use_xdat) then 779 inquire(file='x.dat', exist=sig%use_xdat) 780 if(sig%use_xdat .and. peinf%inode == 0) then 781 write(6,899) 782899 format(1x,"Reading bare exchange matrix elements from x.dat file."/) 783 endif 784 endif 785 if(sig%use_vxcdat) then 786 inquire(file='vxc.dat', exist=sig%use_vxcdat) 787 if(sig%use_vxcdat .and. peinf%inode == 0) then 788 write(6,898) 789898 format(1x,"Reading exchange-correlation matrix elements from vxc.dat file",/) 790 endif 791 endif 792 if(.not.sig%use_vxcdat .and. .not.sig%sigma_correction .and. .not. sig%is_EXX) then 793 inquire(file='VXC', exist=VXC_exists) 794 if(.not. VXC_exists) call die("VXC file is missing.", only_root_writes = .true.) 795 if(peinf%inode == 0) then 796 write(6,897) 797897 format(1x,"Computing exchange-correlation matrix elements from VXC file",/) 798 endif 799 endif 800 ! This is for hybrid functional like calculation (one shot) 801 if (sig%freq_dep.eq.-1 .and. ((1.0d0 - sig%coulomb_mod%long_range_frac_fock > TOL_SMALL) .or. & 802 (1.0d0 - sig%coulomb_mod%short_range_frac_fock > TOL_SMALL))) then 803 804 sig%coul_mod_flag = .true. 805 if(sig%use_vxc2dat) then 806 inquire(file='vxc2.dat', exist=sig%use_vxc2dat) 807 if(sig%use_vxc2dat .and. peinf%inode == 0) then 808 write(6,896) 809896 format(1x,"Reading exchange-correlation matrix elements from vxc2.dat file",/) 810 endif 811 endif 812 if(.not. sig%use_vxc2dat) then 813 inquire(file='VXC2', exist=VXC_exists) 814 if(.not. VXC_exists) call die("VXC2 file is missing.", only_root_writes = .true.) 815 if(peinf%inode == 0) then 816 write(6,895) 817895 format(1x,"Computing exchange-correlation matrix elements from VXC2 file",/) 818 endif 819 endif 820 endif 821 822 !FHJ: report on the frequency dependence method 823 if (peinf%inode==0) then 824 select case (sig%freq_dep) 825 case (-1) 826 tmpstr = 'Hartree-Fock or hybrid functional approximation' 827 case (0) 828 tmpstr = 'Static COHSEX approximation' 829 case (1) 830 tmpstr = 'Hybertsen-Louie Generalized Plasmon Pole model' 831 case (3) 832 tmpstr = 'Godby-Needs Generalized Plasmon Pole model' 833 case (2) 834 select case (sig%freq_dep_method) 835 case (0) 836 tmpstr = 'full frequency Real Axis Integration method' 837 case (2) 838 tmpstr = 'full frequency Contour Deformation method' 839 case default 840 write(0,'(a,i0)') 'ERROR: Got sig%freq_dep_method = ', sig%freq_dep_method 841 call die('Invalid option for `frequency_dependence_method` flag.') 842 endselect 843 case default 844 write(0,'(a,i0)') 'ERROR: Got sig%freq_dep = ', sig%freq_dep 845 call die('Invalid option for `frequency_dependence` flag.') 846 endselect 847 write(6,'(1x,a)') 'Treating W within the '//trim(tmpstr) 848 849 if (sig%freq_dep==2 .and. sig%freq_dep_method==2) then 850 if (sig%cd_int_method/=0 .and. sig%cd_int_method/=2 .and. sig%cd_int_method/=3) then 851 write(0,'(a,i0)') 'ERROR: Got sig%cd_int_method = ', sig%cd_int_method 852 call die('Invalid option for `cd_integration_method` flag.') 853 endif 854 write(6,'(1x,a,i0,a/)') 'We`ll use an integration scheme of order ', & 855 sig%cd_int_method,' in the Contour Deformation method' 856 else 857 write(6,*) 858 endif 859 endif 860 sig%need_advanced = SCALARSIZE==2 .and. sig%freq_dep==2 .and. sig%freq_dep_method/=2 861 862 !FHJ: report on the kind of grid we are using 863 if (sig%freq_dep==2 .or. (sig%freq_dep==1 .and. sig%fdf==-3)) then 864 if (peinf%inode==0) write(6,'(1x,a)', advance='no') 'Frequency sampling: ' 865 if (sig%freq_grid_shift==2) then 866 sig%nfreqeval = 2*int((max_freq_eval+TOL_ZERO)/sig%freqevalstep) + 1 867 if (peinf%inode==0) write(6,'(a)', advance='no') 'using a QP-centered ' 868 elseif (sig%freq_grid_shift==1) then 869 if (peinf%inode==0) write(6,'(a)', advance='no') 'using a Fermi-energy-shifted ' 870 elseif (sig%freq_grid_shift==0) then 871 if (peinf%inode==0) write(6,'(a)', advance='no') 'using an absolute energy ' 872 else 873 call die('Invalid value for frequency_grid_shift.') 874 endif 875 if (peinf%inode==0) write(6,'(a,i0,a/)') 'grid with ', sig%nfreqeval, ' points.' 876 endif 877 !FHJ: Make sure the user is not using incompatible flags. 878 !if ((old_grid_set .and. sig%freq_grid_shift==2).and.sig%freq_dep==2) then 879 if (old_grid_set .and. sig%freq_grid_shift==2) then 880 if (peinf%inode==0) then 881 write(0,'(/a)') 'ERROR: flags `init_frequency_eval` and/or `number_frequency_eval` were set, but' 882 write(0,'(a)') 'you didn`t change the `frequency_grid_shift`. Either:' 883 write(0,'(a)') ' - set `frequency_grid_shift` to 0 or 1; or' 884 write(0,'(a/)') ' - don`t set the `init_frequency_eval` and `number_frequency_eval` flags.' 885 endif 886 call die('Inconsistent flags with frequency_grid_shift') 887 endif 888 889#ifdef MPI 890 ! root lets the others go after it is done reading (see beginning of function) 891 if(peinf%inode == 0) call MPI_Barrier(MPI_COMM_WORLD, mpierr) 892 893 ! FHJ: Broadcast spline scissors (this must be called after the barrier) 894 ! FHJ: inner splines 895 call MPI_BCAST(sig%spl_tck%n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 896 if ( sig%spl_tck%n > 0 ) then 897 if(peinf%inode/=0) then 898 SAFE_ALLOCATE(sig%spl_tck%t, (sig%spl_tck%n)) 899 SAFE_ALLOCATE(sig%spl_tck%c, (sig%spl_tck%n)) 900 endif 901 call MPI_BCAST(sig%spl_tck%t, sig%spl_tck%n, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr) 902 call MPI_BCAST(sig%spl_tck%c, sig%spl_tck%n, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr) 903 call MPI_BCAST(sig%spl_tck%k, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 904 endif 905 ! FHJ: outer splines 906 call MPI_BCAST(sig%spl_tck_outer%n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 907 if ( sig%spl_tck_outer%n > 0 ) then 908 if(peinf%inode/=0) then 909 SAFE_ALLOCATE(sig%spl_tck_outer%t, (sig%spl_tck_outer%n)) 910 SAFE_ALLOCATE(sig%spl_tck_outer%c, (sig%spl_tck_outer%n)) 911 endif 912 call MPI_BCAST(sig%spl_tck_outer%t, sig%spl_tck_outer%n, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr) 913 call MPI_BCAST(sig%spl_tck_outer%c, sig%spl_tck_outer%n, MPI_REAL_DP, 0, MPI_COMM_WORLD, mpierr) 914 call MPI_BCAST(sig%spl_tck_outer%k, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr) 915 endif 916#endif 917 918 POP_SUB(inread) 919 920 return 921 922110 write(errmsg,'(3a)') 'Unexpected characters were found while reading the value for the keyword ', & 923 trim(keyword), '. ' 924 call die(errmsg, only_root_writes = .true.) 925 926112 write(errmsg,'(3a)') 'Unexpected characters were found while reading elements of the ', & 927 trim(blockword),' block.' 928 call die(errmsg, only_root_writes = .true.) 929 930end subroutine inread 931 932 933! ZL: fold back kpoint into [0, 1) range 934! kpoint MUST be in fractional coordinate 935subroutine fold_back_zero_one(ka, kb, kc) 936 use global_m 937 implicit none 938 939 real(DP), intent(inout) :: ka, kb, kc 940 real(DP) :: infsmall = 1.0d-5 ! ZL: k/q/phonq-points should have higher accuracy than this 941 integer :: i 942 943 PUSH_SUB(fold_back_zero_one) 944 ! ZL: TODO: get rid of this function by using macros UNIT_RANGE() and TOL_SMALL 945 946 ! ka 947 do while (ka .gt. 1.0-infsmall) 948 ka = ka - 1.0 949 enddo 950 do while (ka .lt. 0.0-infsmall) 951 ka = ka + 1.0 952 enddo 953 954 ! kb 955 do while (kb .gt. 1.0-infsmall) 956 kb = kb - 1.0 957 enddo 958 do while (kb .lt. 0.0-infsmall) 959 kb = kb + 1.0 960 enddo 961 962 ! kc 963 do while (kc .gt. 1.0-infsmall) 964 kc = kc - 1.0 965 enddo 966 do while (kc .lt. 0.0-infsmall) 967 kc = kc + 1.0 968 enddo 969 970 POP_SUB(fold_back_zero_one) 971 972 return 973end subroutine fold_back_zero_one 974 975