1! 2! Copyright (C) 2001-2018 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! 8!---------------------------------------------------------------------------- 9SUBROUTINE lr_readin 10 !----------------------------------------------------------------------- 11 ! 12 ! This routine reads the control variables from standard input (unit 5). 13 ! A second routine read_file reads the variables saved to file 14 ! by PW scf (and PW nscf for EELS). 15 ! 16 USE lr_variables 17 USE lr_dav_variables 18 USE kinds, ONLY : DP 19 USE io_files, ONLY : tmp_dir, prefix, wfc_dir, create_directory 20 USE lsda_mod, ONLY : current_spin, nspin, isk, lsda 21 USE control_flags, ONLY : use_para_diag, tqr, gamma_only,& 22 & do_makov_payne 23 USE scf, ONLY : vltot, v, vrs, vnew, & 24 & destroy_scf_type, rho 25 USE fft_base, ONLY : dfftp, dffts 26 USE gvecs, ONLY : doublegrid 27 USE wvfct, ONLY : nbnd, et, wg, current_k 28 USE lsda_mod, ONLY : isk 29 USE ener, ONLY : ef 30 USE io_global, ONLY : ionode, ionode_id, stdout, meta_ionode, meta_ionode_id 31 USE klist, ONLY : nks, wk, nelec, lgauss, ltetra 32 USE fixed_occ, ONLY : tfixed_occ 33 USE symm_base, ONLY : nosym 34 USE check_stop, ONLY : max_seconds 35 USE realus, ONLY : real_space, init_realspace_vars, generate_qpointlist, & 36 betapointlist 37 USE funct, ONLY : dft_is_meta 38 USE charg_resp, ONLY : w_T_prefix, omeg, w_T_npol, epsil 39 USE mp, ONLY : mp_bcast 40 USE mp_global, ONLY : my_pool_id, intra_image_comm, & 41 & intra_bgrp_comm, nproc_image, & 42 & nproc_pool, nproc_pool_file, & 43 & nproc_image_file, nproc_bgrp, & 44 & nproc_bgrp_file, my_image_id 45 USE DFUNCT, ONLY : newd 46 USE vlocal, ONLY : strf 47 USE martyna_tuckerman, ONLY : do_comp_mt 48 USE esm, ONLY : do_comp_esm 49 USE qpoint, ONLY : xq 50 USE io_rho_xml, ONLY : write_scf 51 USE mp_bands, ONLY : ntask_groups 52 USE constants, ONLY : eps4, rytoev 53 USE control_lr, ONLY : lrpa, alpha_mix 54 USE mp_world, ONLY : world_comm 55 56 IMPLICIT NONE 57 ! 58 CHARACTER(LEN=256) :: wfcdir = 'undefined', outdir 59 CHARACTER(LEN=256), EXTERNAL :: trimcheck 60 ! 61 CHARACTER(LEN=256) :: beta_gamma_z_prefix 62 ! Fine control of beta_gamma_z file 63 CHARACTER(LEN=80) :: disk_io 64 ! Specify the amount of I/O activities 65 CHARACTER(LEN=6) :: int_to_char 66 INTEGER :: ios, iunout, ierr, ipol 67 LOGICAL, EXTERNAL :: check_para_diag 68 ! 69 CHARACTER(LEN=80) :: card 70 INTEGER :: i 71 ! 72 NAMELIST / lr_input / restart, restart_step ,lr_verbosity, prefix, outdir, & 73 & test_case_no, wfcdir, disk_io, max_seconds 74 NAMELIST / lr_control / itermax, ipol, ltammd, lrpa, & 75 & charge_response, no_hxc, n_ipol, project, & 76 & scissor, pseudo_hermitian, d0psi_rs, lshift_d0psi, & 77 & q1, q2, q3, approximation, calculator, alpha_mix, start, & 78 & end, increment, epsil, units 79 NAMELIST / lr_post / omeg, beta_gamma_z_prefix, w_T_npol, plot_type, epsil, itermax_int,sum_rule 80 namelist / lr_dav / num_eign, num_init, num_basis_max, residue_conv_thr, precondition, & 81 & dav_debug, reference,single_pole, sort_contr, diag_of_h, close_pre, & 82 & broadening,print_spectrum,start,finish,step,if_check_orth, if_random_init, & 83 & if_check_her,p_nbnd_occ,p_nbnd_virt,poor_of_ram,poor_of_ram2,max_iter, & 84 & conv_assistant,if_dft_spectrum,no_hxc,d0psi_rs,lshift_d0psi, & 85 & lplot_drho, vccouple_shift, ltammd 86 ! 87#if defined(__MPI) 88 IF (ionode) THEN 89#endif 90 ! 91 ! Checking for the path to the output directory. 92 ! 93 CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) 94 IF ( trim( outdir ) == ' ' ) outdir = './' 95 ! 96 ! Set default values for variables in namelist. 97 ! 98 itermax = 500 99 restart = .FALSE. 100 restart_step = itermax+1 101 lr_verbosity = 1 102 prefix = 'pwscf' 103 disk_io = 'default' 104 ltammd = .FALSE. 105 d0psi_rs=.false. 106 lshift_d0psi = .true. 107 pseudo_hermitian=.true. 108 ipol = 1 109 n_ipol = 1 110 no_hxc = .FALSE. 111 lrpa = .false. 112 charge_response = 0 113 sum_rule = -99 114 test_case_no = 0 115 beta_gamma_z_prefix = 'undefined' 116 omeg= 0.0_DP 117 w_T_npol = 1 118 plot_type = 1 119 project = .FALSE. 120 max_seconds = 1.0E+7_DP 121 scissor = 0.d0 122 ! 123 ! For EELS 124 ! 125 q1 = 1.0d0 126 q2 = 1.0d0 127 q3 = 1.0d0 128 approximation = 'TDDFT' 129 calculator = 'lanczos' 130 ! 131 ! Sternheimer 132 ! 133 start_freq=1 134 last_freq=0 135 alpha_mix(:) = 0.0D0 136 alpha_mix(1) = 0.7D0 137 units = 0 138 end = 2.5D0 139 epsil = 0.02D0 140 increment = 0.001D0 141 ! 142 ! For lr_dav (Davidson program) 143 ! 144 num_eign=1 145 num_init=2 146 num_basis_max=20 147 broadening=0.005d0 148 residue_conv_thr=1.0E-4 149 close_pre=1.0E-5 150 turn2planb=1.0E-3 151 precondition=.true. 152 dav_debug=.false. 153 reference=0.0d0 154 vccouple_shift=0.0d0 155 single_pole=.false. 156 sort_contr=.true. 157 print_spectrum=.true. 158 start=0.0d0 159 finish=1.0d0 160 step=0.001d0 161 if_check_orth=.false. 162 if_check_her=.false. 163 if_random_init=.false. 164 p_nbnd_occ=10 165 p_nbnd_virt=10 166 poor_of_ram=.false. 167 poor_of_ram2=.false. 168 max_iter=100 169 conv_assistant=.false. 170 if_dft_spectrum=.false. 171 lplot_drho=.false. 172 ! 173 ! Reading possible plugin arguments 174 ! 175 CALL plugin_arguments() 176 ! 177 ! Reading the namelist lr_input 178 ! 179 CALL input_from_file( ) 180 ! 181 READ (5, lr_input, err = 200, iostat = ios) 182200 CALL errore ('lr_readin', 'reading lr_input namelist', ABS (ios) ) 183 ! 184 ! Reading the namelist lr_dav or lr_control 185 ! 186 IF (davidson) THEN 187 READ (5, lr_dav, err = 201, iostat = ios) 188201 CALL errore ('lr_readin', 'reading lr_dav namelist', ABS (ios) ) 189 ELSE 190 READ (5, lr_control, err = 202, iostat = ios) 191202 CALL errore ('lr_readin', 'reading lr_control namelist', ABS (ios) ) 192 ENDIF 193 ! 194 ! Reading the namelist lr_post (only for optical case) 195 ! 196 IF (charge_response == 1 .AND. .NOT.eels) THEN 197 ! 198 READ (5, lr_post, err = 203, iostat = ios) 199203 CALL errore ('lr_readin', 'reading lr_post namelist', ABS (ios) ) 200 ! 201 bgz_suffix = TRIM ( "-stage2.beta_gamma_z." ) 202 WRITE(stdout,'(/5x,"Prefix of current run is appended by -stage2")') 203 ! 204 IF ( beta_gamma_z_prefix == 'undefined' ) THEN 205 beta_gamma_z_prefix = TRIM(prefix) 206 ENDIF 207 ! 208 ELSE 209 bgz_suffix = TRIM ( ".beta_gamma_z." ) 210 ENDIF 211 ! 212 ! The status of the real space flags should be read manually 213 ! 214 ! Set-up all the dir and suffix variables. 215 ! 216 outdir = trimcheck(outdir) 217 tmp_dir = outdir 218 ! 219 !IF ( .NOT. TRIM( wfcdir ) == 'undefined' ) THEN 220 ! wfc_dir = trimcheck ( wfcdir ) 221 !ENDIF 222 ! 223 IF (.NOT.eels) THEN 224 w_T_prefix = TRIM( tmp_dir ) // & 225 & TRIM( beta_gamma_z_prefix ) // ".beta_gamma_z." 226 ENDIF 227 ! 228 ierr = 0 229 ! 230 IF (eels) THEN 231 ! 232 IF (n_ipol /= 1) THEN 233 WRITE(stdout,'(5X,"n_ipol /= 1 is not allowed for EELS.")') 234 WRITE(stdout,'(5X,"Setting n_ipol = 1 ...")') 235 n_ipol = 1 236 ENDIF 237 IF (ipol /= 1) THEN 238 WRITE(stdout,'(5X,"ipol /= 1 is not allowed for EELS.")') 239 WRITE(stdout,'(5X,"Setting ipol = 1 ...")') 240 ipol = 1 241 ENDIF 242 LR_polarization = 1 243 ! 244 ELSE 245 ! 246 ! Optics: set up polarization direction(s) 247 ! 248 IF ( ipol==4 ) THEN 249 ! 250 n_ipol = 3 251 LR_polarization = 1 252 ! 253 ELSE 254 ! 255 LR_polarization = ipol 256 ! 257 ENDIF 258 ! 259 ENDIF 260 ! 261 IF (itermax_int < itermax) itermax_int = itermax 262 ! 263 ! Limited disk_io support: currently only one setting is supported 264 ! 265 SELECT CASE( TRIM( disk_io ) ) 266 ! 267 CASE ( 'reduced' ) 268 ! 269 lr_io_level = -1 270 restart = .FALSE. 271 ! 272 CASE DEFAULT 273 ! 274 lr_io_level = 1 275 ! 276 END SELECT 277 ! 278 IF (eels) THEN 279 ! 280 ! Level of approximation in turboEELS. 281 ! 282 SELECT CASE( trim(approximation) ) 283 ! 284 CASE ( 'TDDFT' ) 285 ! 286 no_hxc = .FALSE. 287 lrpa = .FALSE. 288 ! 289 CASE ( 'IPA' ) 290 ! 291 no_hxc = .TRUE. 292 lrpa = .TRUE. 293 ! 294 CASE ( 'RPA_with_CLFE' ) 295 ! 296 no_hxc = .FALSE. 297 lrpa = .TRUE. 298 ! 299 CASE DEFAULT 300 ! 301 CALL errore( 'lr_readin', 'Approximation ' // & 302 & trim( approximation ) // ' not implemented', 1 ) 303 ! 304 END SELECT 305 ! 306 ! We do this trick because xq is used in LR_Modules/dv_of_drho.f90 307 ! in the Hartree term ~1/|xq+k|^2 308 ! 309 xq(1) = q1 310 xq(2) = q2 311 xq(3) = q3 312 ! 313 IF ( (q1.lt.eps4) .AND. (q2.lt.eps4) .AND. (q3.lt.eps4) ) & 314 CALL errore( 'lr_readin', 'The transferred momentum |q| is too small, the limit is not implemented.', 1 ) 315 ! 316 ENDIF 317 ! 318#if defined(__MPI) 319 ENDIF 320 ! 321 CALL bcast_lr_input 322 ! 323#endif 324 ! 325 IF ( trim(calculator)=='sternheimer' ) THEN 326 nfs=0 327 nfs = ((end-start) / increment) + 1 328 if (nfs < 1) call errore('lr_readin','Too few frequencies',1) 329 ALLOCATE(fiu(nfs)) 330 ALLOCATE(fru(nfs)) 331 ALLOCATE(comp_f(nfs)) 332 comp_f=.TRUE. 333 IF (units == 0) THEN 334 fru(1) =start 335 fru(nfs) = end 336 deltaf = increment 337 ELSEIF (units == 1) THEN 338 fru(1) =start/rytoev 339 fru(nfs) = end/rytoev 340 deltaf = increment/rytoev 341 ENDIF 342 fiu(:) = epsil 343 DO i=2,nfs-1 344 fru(i)=fru(1) + (i-1) * deltaf 345 ENDDO 346 IF (start_freq<=0) start_freq=1 347 IF (last_freq<=0.OR.last_freq>nfs) last_freq=nfs 348 ELSE 349 nfs=1 350 ALLOCATE(fru(1)) 351 ALLOCATE(fiu(1)) 352 ALLOCATE(comp_f(1)) 353 fru=0.0_DP 354 fiu=0.0_DP 355 comp_f=.TRUE. 356 END IF 357 ! 358 ! Required for restart runs as this never gets initialized. 359 ! 360 current_k = 1 361 ! 362 outdir = TRIM( tmp_dir ) // TRIM( prefix ) // '.save' 363 ! 364 ! EELS: Create a temporary directory for nscf files, and for 365 ! writing of the turboEELS restart files. 366 ! 367 IF (eels) THEN 368 tmp_dir_lr = TRIM (tmp_dir) // 'tmp_eels/' 369 CALL create_directory(tmp_dir_lr) 370 ENDIF 371 ! 372 ! Now PWSCF XML file will be read, and various initialisations will be done. 373 ! Allocate space for PW scf variables, read and check them. 374 ! Optical case: the variables igk_k and ngk are set up through this path: 375 ! read_file -> init_igk. 376 ! EELS: the variables igk_k and ngk will be re-set up later (because there 377 ! will be not only poins k but also points k+q) through the path: 378 ! lr_run_nscf -> init_run -> allocate_wfc_k -> init_igk 379 ! 380 CALL read_file() 381 ! 382 IF (.NOT.eels .AND. (tqr .OR. real_space)) & 383 WRITE(stdout,'(/5x,"Status of real space flags: TQR=", L5 , & 384 & " REAL_SPACE=", L5)') tqr, real_space 385 ! 386 ! Set wfc_dir - this is done here because read_file sets wfc_dir = tmp_dir 387 ! FIXME:,if wfcdir is not present in input, wfc_dir is set to "undefined" 388 ! instead of tmp_dir, because of the logic used in the rest of TDDFPT 389 ! 390 wfc_dir = trimcheck ( wfcdir ) 391 ! 392 IF (eels) THEN 393 ! 394 ! Specify the temporary directory. 395 ! 396 tmp_dir = tmp_dir_lr 397 ! 398 ! Copy the scf-charge-density to the tmp_dir (PH/check_initial_status.f90). 399 ! Needed for the nscf calculation. 400 ! 401 IF (.NOT.restart) CALL write_scf( rho, nspin ) 402 ! 403 ENDIF 404 ! 405 ! Make sure all the features used in the PWscf calculation 406 ! are actually supported by TDDFPT. 407 ! 408 CALL input_sanity() 409 ! 410 ! Compute and add plugin contributions to SCF potential 411 ! 412 CALL plugin_read_input("TD") 413 ! 414 CALL plugin_initbase() 415 ! 416 CALL plugin_init_ions() 417 ! 418 CALL plugin_init_cell() 419 ! 420 CALL plugin_init_potential( v%of_r(:,1) ) 421 ! 422 CALL plugin_scf_potential( rho, .FALSE., -1.D0, v%of_r(:,1) ) 423 ! 424 CALL plugin_clean( 'TD', .FALSE. ) 425 ! 426 ! Deallocate some variables created with read_file but not used in TDDFPT 427 ! 428 DEALLOCATE( strf ) 429 CALL destroy_scf_type(vnew) 430 ! 431 ! Re-initialize all needed quantities from the scf run 432 ! I. Timrov: this was already done in read_file. 433 current_spin = 1 434 ! 435 ! Now put the potential calculated in read_file into the correct place 436 ! and deallocate the redundant associated variables. 437 ! Set the total local potential vrs on the smooth mesh 438 ! adding the scf (Hartree + XC) part and the sum of 439 ! all the local pseudopotential contributions. 440 ! vrs = vltot + v%of_r 441 ! 442 CALL set_vrs ( vrs, vltot, v%of_r, 0, 0, dfftp%nnr, nspin, doublegrid ) 443 ! 444 DEALLOCATE( vltot ) 445 CALL destroy_scf_type(v) 446 ! 447 ! Recalculate the weights of the Kohn-Sham orbitals. 448 ! (Should this not be a call to weights() to make this 449 ! less insulator specific?) 450 ! 451 IF (.NOT.eels) CALL iweights( nks, wk, nbnd, nelec, et, ef, wg, 0, isk) 452 ! 453 IF ( charge_response == 2 .AND. .NOT.eels) CALL lr_set_boxes_density() 454 ! 455 ! Scalapack related stuff. 456 ! 457 use_para_diag = check_para_diag( nbnd ) 458 ! 459 RETURN 460 ! 461CONTAINS 462 ! 463 SUBROUTINE input_sanity() 464 !-------------------------------------------------------------------------- 465 ! 466 ! This subroutine aims to gather all of the input sanity checks 467 ! (features enabled in PWscf which are unsupported in TDDFPT). 468 ! Written by Simone Binnie 2011 469 ! Modified by Iurii Timrov 2015 470 ! 471 USE paw_variables, ONLY : okpaw 472 USE uspp, ONLY : okvan 473 USE funct, ONLY : dft_is_hybrid 474 USE ldaU, ONLY : lda_plus_u 475 476 IMPLICIT NONE 477 ! 478 ! Charge response mode 1 is the "do Lanczos chains twice, conserve memory" scheme. 479 ! 480 IF (.NOT.eels) THEN 481 ! 482 IF (charge_response == 1 .AND. ( omeg == 0.D0 .AND. sum_rule == -99 ) ) & 483 & CALL errore ('lr_readin', & 484 & 'omeg must be defined for charge response mode 1', 1 ) 485 ! 486 IF ( project .AND. charge_response /= 1) & 487 & CALL errore ('lr_readin', & 488 & 'projection is possible only in charge response mode 1', 1 ) 489 ! 490 IF (gamma_only) THEN 491 nosym=.true. 492 WRITE(stdout,*) "Symmetries are disabled for the gamma_only case" 493 ENDIF 494 ! 495 ENDIF 496 ! 497 ! Meta-DFT currently not supported by TDDFPT 498 ! 499 IF (dft_is_meta()) CALL errore( 'lr_readin', 'Meta DFT is not implemented yet', 1 ) 500 ! 501 ! Hubbard U is not supported 502 ! 503 IF (lda_plus_u) CALL errore('lr_readin', 'TDDFPT with Hubbard U is not implemented',1) 504 ! 505 ! Tetrahedron method and fixed occupations are not implemented. 506 ! 507 IF (ltetra) CALL errore( 'lr_readin', 'ltetra is not implemented', 1 ) 508 IF (tfixed_occ) CALL errore( 'lr_readin', 'tfixed_occ is not implemented', 1 ) 509 ! 510 ! Some limitations of turboTDDFT (and not of turboEELS). 511 ! 512 IF (.NOT. eels) THEN 513 ! 514 ! Non-insulating systems currently not supported by turboTDDFPT, but 515 ! supported by turboEELS. 516 ! 517 IF (lgauss) CALL errore( 'lr_readin', 'turboTDDFT is not extended to metals', 1 ) 518 ! 519 ! Symmetry is not supported. 520 ! 521 IF (.NOT.nosym ) CALL errore( 'lr_readin', 'Linear response calculation' // & 522 & 'is not implemented with symmetry', 1 ) 523 ! 524 ! K-points are implemented but still unsupported (use at your own risk!) 525 ! 526 IF (.NOT. gamma_only ) CALL errore('lr_readin', 'k-point algorithm is not tested yet',1) 527 ! 528 ENDIF 529 ! 530 ! No taskgroups and EXX. 531 ! 532 IF (dffts%has_task_groups .AND. dft_is_hybrid()) & 533 & CALL errore( 'lr_readin', ' Linear response calculation ' // & 534 & 'not implemented for EXX+Task groups', 1 ) 535 ! 536 ! Experimental task groups warning. 537 ! 538 IF (dffts%has_task_groups) & 539 & CALL infomsg( 'lr_readin','Usage of task & 540 & groups with TDDFPT is still experimental. Use at your own risk.' ) 541 ! & CALL errore( 'lr_readin', ' Linear response calculation ' // & 542 ! & 'not implemented for task groups', 1 ) 543 ! 544 ! No PAW support. 545 ! 546 IF (okpaw) & 547 & CALL errore( 'lr_readin', ' Linear response calculation ' // & 548 & 'not implemented for PAW', 1 ) 549 ! 550 ! No USPP+EXX support. 551 ! 552 IF (okvan .AND. dft_is_hybrid()) & 553 & CALL errore( 'lr_readin', ' Linear response calculation ' // & 554 & 'not implemented for EXX+Ultrasoft', 1 ) 555 ! 556 ! Spin-polarised case is not implemented, but partially accounted in 557 ! some routines. 558 ! 559 IF (lsda) CALL errore( 'lr_readin', 'LSDA is not implemented', 1 ) 560 ! 561 IF (real_space) THEN 562 IF (eels) THEN 563 CALL errore( 'lr_readin', 'Option real_space=.true. is not implemented', 1 ) 564 ELSE 565 CALL errore( 'lr_readin', 'Option real_space=.true. is not tested', 1 ) 566 ENDIF 567 ENDIF 568 ! 569 ! EELS-related restrictions 570 ! 571 IF (eels) THEN 572 ! 573 IF (gamma_only) CALL errore( 'lr_readin', 'gamma_only is not supported', 1 ) 574 ! 575 ! Tamm-Dancoff approximation is not recommended to be used with EELS, and 576 ! thus it was not implemented. 577 ! 578 IF (ltammd) CALL errore( 'lr_readin', 'EELS + Tamm-Dancoff approximation is not supported', 1 ) 579 ! 580 IF (project) CALL errore( 'lr_readin', 'project is not allowed', 1 ) 581 IF (tqr) CALL errore( 'lr_readin', 'tqr is not supported', 1 ) 582 IF (charge_response /= 0) CALL errore( 'lr_readin', 'charge_response /= 0 is not allowed', 1 ) 583 IF (dft_is_hybrid()) CALL errore( 'lr_readin', 'EXX is not supported', 1 ) 584 IF (do_comp_mt) CALL errore( 'lr_readin', 'Martyna-Tuckerman PBC is not supported.', 1 ) 585 IF (d0psi_rs) CALL errore( 'lr_readin', 'd0psi_rs is not allowed', 1 ) 586 ! 587 ! Note, all variables of the turboDavidson code cannot be used by turboEELS. 588 ! 589 ! EELS + plugins is not supported. 590 ! 591 CALL plugin_check('lr_readin') 592 ! 593 ENDIF 594 ! 595 RETURN 596 ! 597 END SUBROUTINE input_sanity 598 ! 599END SUBROUTINE lr_readin 600