1! 2! Copyright (C) 2001-2020 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! 9!---------------------------------------------------------------------------- 10SUBROUTINE phq_readin() 11 !---------------------------------------------------------------------------- 12 ! 13 ! This routine reads the control variables for the program phononq. 14 ! from standard input (unit 5). 15 ! A second routine readfile reads the variables saved on a file 16 ! by the self-consistent program. 17 ! 18 ! 19 USE kinds, ONLY : DP 20 USE ions_base, ONLY : nat, ntyp => nsp 21 USE mp, ONLY : mp_bcast 22 USE mp_world, ONLY : world_comm 23 USE ions_base, ONLY : amass, atm 24 USE check_stop, ONLY : max_seconds 25 USE start_k, ONLY : reset_grid 26 USE klist, ONLY : xk, nks, nkstot, lgauss, two_fermi_energies, ltetra 27 USE control_flags, ONLY : gamma_only, tqr, restart, io_level, & 28 ts_vdw, ldftd3, lxdm, isolve 29 USE funct, ONLY : dft_is_meta, dft_is_hybrid 30 USE uspp, ONLY : okvan 31 USE fixed_occ, ONLY : tfixed_occ 32 USE lsda_mod, ONLY : lsda, nspin 33 USE fft_base, ONLY : dffts 34 USE spin_orb, ONLY : domag 35 USE cellmd, ONLY : lmovecell 36 USE run_info, ONLY : title 37 USE control_ph, ONLY : maxter, alpha_mix, lgamma_gamma, epsil, & 38 zue, zeu, xmldyn, newgrid, & 39 trans, reduce_io, tr2_ph, niter_ph, & 40 nmix_ph, ldisp, recover, lnoloc, start_irr, & 41 last_irr, start_q, last_q, current_iq, tmp_dir_ph, & 42 ext_recover, ext_restart, u_from_file, ldiag, & 43 search_sym, lqdir, electron_phonon, tmp_dir_phq, & 44 rec_code_read, qplot, only_init, only_wfc, & 45 low_directory_check, nk1, nk2, nk3, k1, k2, k3 46 USE save_ph, ONLY : tmp_dir_save, save_ph_input_variables 47 USE gamma_gamma, ONLY : asr 48 USE partial, ONLY : atomo, nat_todo, nat_todo_input 49 USE output, ONLY : fildyn, fildvscf, fildrho 50 USE disp, ONLY : nq1, nq2, nq3, x_q, wq, nqs, lgamma_iq 51 USE io_files, ONLY : tmp_dir, prefix, postfix, create_directory, & 52 check_tempdir, xmlpun_schema 53 USE noncollin_module, ONLY : i_cons, noncolin 54 USE control_flags, ONLY : iverbosity, modenum 55 USE io_global, ONLY : meta_ionode, meta_ionode_id, ionode, ionode_id, stdout 56 USE mp_images, ONLY : nimage, my_image_id, intra_image_comm, & 57 me_image, nproc_image 58 USE mp_pools, ONLY : npool 59 USE paw_variables, ONLY : okpaw 60 USE ramanm, ONLY : eth_rps, eth_ns, lraman, elop, dek 61 USE freq_ph, ONLY : fpol, fiu, nfs 62 USE cryst_ph, ONLY : magnetic_sym 63 USE ph_restart, ONLY : ph_readfile 64 USE el_phon, ONLY : elph,elph_mat,elph_simple,elph_epa,elph_nbnd_min, elph_nbnd_max, & 65 el_ph_sigma, el_ph_nsigma, el_ph_ngauss,auxdvscf 66 USE dfile_star, ONLY : drho_star, dvscf_star 67 68 USE qpoint, ONLY : nksq, xq 69 USE control_lr, ONLY : lgamma, lrpa 70 ! YAMBO > 71 USE YAMBO, ONLY : elph_yambo,dvscf_yambo 72 ! YAMBO < 73 USE elph_tetra_mod,ONLY : elph_tetra, lshift_q, in_alpha2f 74 USE ktetra, ONLY : tetra_type 75 USE ldaU, ONLY : lda_plus_u, U_projection, lda_plus_u_kind 76 USE ldaU_ph, ONLY : read_dns_bare, d2ns_type 77 USE dvscf_interpolate, ONLY : ldvscf_interpolate, do_long_range, & 78 do_charge_neutral, wpot_dir 79 USE ahc, ONLY : elph_ahc, ahc_dir, ahc_nbnd, ahc_nbndskip, & 80 skip_upperfan 81 USE wannier_gw, ONLY : l_head, omega_gauss, n_gauss, grid_type, nsteps_lanczos,second_grid_n,second_grid_i,& 82 &l_scissor,scissor, len_head_block_freq, len_head_block_wfc 83 ! 84 IMPLICIT NONE 85 ! 86 CHARACTER(LEN=256), EXTERNAL :: trimcheck 87 ! 88 INTEGER :: ios, ipol, iter, na, it, ierr, ierr1 89 ! integer variable for I/O control 90 ! counter on polarizations 91 ! counter on iterations 92 ! counter on atoms 93 ! counter on types 94 REAL(DP), ALLOCATABLE :: amass_input(:) 95 ! save masses read from input here 96 CHARACTER (LEN=256) :: outdir, filename 97 CHARACTER (LEN=8) :: verbosity 98 CHARACTER(LEN=80) :: card 99 CHARACTER(LEN=6) :: int_to_char 100 INTEGER :: i 101 LOGICAL :: nogg 102 LOGICAL :: q2d, q_in_band_form 103 INTEGER, EXTERNAL :: atomic_number 104 REAL(DP), EXTERNAL :: atom_weight 105 LOGICAL, EXTERNAL :: imatches 106 LOGICAL, EXTERNAL :: has_xml 107 LOGICAL :: exst, parallelfs 108 REAL(DP), ALLOCATABLE :: xqaux(:,:) 109 INTEGER, ALLOCATABLE :: wqaux(:) 110 INTEGER :: nqaux, iq 111 CHARACTER(len=80) :: diagonalization='david' 112 ! 113 NAMELIST / INPUTPH / tr2_ph, amass, alpha_mix, niter_ph, nmix_ph, & 114 nat_todo, verbosity, iverbosity, outdir, epsil, & 115 trans, zue, zeu, max_seconds, reduce_io, & 116 modenum, prefix, fildyn, fildvscf, fildrho, & 117 ldisp, nq1, nq2, nq3, & 118 eth_rps, eth_ns, lraman, elop, dek, recover, & 119 fpol, asr, lrpa, lnoloc, start_irr, last_irr, & 120 start_q, last_q, nogg, ldiag, search_sym, lqdir, & 121 nk1, nk2, nk3, k1, k2, k3, & 122 drho_star, dvscf_star, only_init, only_wfc, & 123 elph_nbnd_min, elph_nbnd_max, el_ph_ngauss, & 124 el_ph_nsigma, el_ph_sigma, electron_phonon, & 125 q_in_band_form, q2d, qplot, low_directory_check, & 126 lshift_q, read_dns_bare, d2ns_type, diagonalization, & 127 ldvscf_interpolate, do_long_range, do_charge_neutral, & 128 wpot_dir, ahc_dir, ahc_nbnd, ahc_nbndskip, & 129 skip_upperfan, & 130 l_head, omega_gauss, n_gauss, grid_type,nsteps_lanczos,l_scissor,scissor,& 131 second_grid_n,second_grid_i,len_head_block_wfc,len_head_block_freq 132 133 ! tr2_ph : convergence threshold 134 ! amass : atomic masses 135 ! alpha_mix : the mixing parameter 136 ! niter_ph : maximum number of iterations 137 ! nmix_ph : number of previous iterations used in mixing 138 ! nat_todo : number of atom to be displaced 139 ! verbosity : verbosity control (iverbosity is obsolete) 140 ! outdir : directory where input, output, temporary files reside 141 ! epsil : if true calculate dielectric constant 142 ! trans : if true calculate phonon 143 ! electron-phonon : select the kind of electron-phonon calculation 144 ! elph : if true calculate electron-phonon coefficients 145 ! elph_mat : if true eph coefficients for wannier 146 ! zue : if .true. calculate effective charges ( d force / dE ) 147 ! zeu : if .true. calculate effective charges ( d P / du ) 148 ! lraman : if true calculate raman tensor 149 ! elop : if true calculate electro-optic tensor 150 ! max_seconds : maximum cputime for this run 151 ! reduce_io : reduce I/O to the strict minimum 152 ! modenum : single mode calculation 153 ! prefix : the prefix of files produced by pwscf 154 ! fildyn : output file for the dynamical matrix 155 ! fildvscf : output file containing deltavsc 156 ! fildrho : output file containing deltarho 157 ! fildrho_dir : directory where fildrho files will be stored (default: outdir or ESPRESSO_FILDRHO_DIR variable) 158 ! eth_rps : threshold for calculation of Pc R |psi> (Raman) 159 ! eth_ns : threshold for non-scf wavefunction calculation (Raman) 160 ! dek : delta_xk used for wavefunctions derivation (Raman) 161 ! recover : recover=.true. to restart from an interrupted run 162 ! asr : in the gamma_gamma case apply acoustic sum rule 163 ! start_q : in q list does the q points from start_q to last_q 164 ! last_q : 165 ! start_irr : does the irred. representation from start_irr to last_irr 166 ! last_irr : 167 ! nogg : if .true. lgamma_gamma tricks are not used 168 ! ldiag : if .true. force diagonalization of the dyn mat 169 ! lqdir : if .true. each q writes in its own directory 170 ! search_sym : if .true. analyze symmetry if possible 171 ! nk1,nk2,nk3, k1,k2,k3 : 172 ! when specified in input, the phonon run uses a different 173 ! k-point mesh from that used for the charge density. 174 ! 175 ! dvscf_star%open : if .true. write in dvscf_star%dir the dvscf_q 176 ! 'for all q' in the star of q with suffix dvscf_star%ext. 177 ! The dvscf_q' is written in the basis dvscf_star%basis; 178 ! if dvscf_star%pat is .true. also save a pattern file. 179 ! dvscf_star%dir, dvscf_star%ext, dvscf_star%basis : see dvscf_star%open 180 ! drho_star%open : like dvscf_star%open but for drho_q 181 ! drho_star%dir, drho_star%ext, drho_star%basis : see drho_star%open 182 ! 183 ! elph_nbnd_min, 184 ! elph_nbnd_max: if (elph_mat=.true.) it dumps the eph matrix element from elph_nbnd_min 185 ! to elph_nbnd_max 186 ! el_ph_ngauss, 187 ! el_ph_nsigma, 188 ! el_ph_sigma : if (elph_mat=.true.) it defines the kind and the val-ue of the 189 ! smearing to be used in the eph coupling calculation. 190 ! qplot, : if true a list of q points is given in input 191 ! q_in_band_form: if true the input list of q points defines paths 192 ! q2d, : if .true. the q list define a mesh in a square. 193 ! low_directory_check : if .true. only the requested representations 194 ! are searched on file 195 ! 196 ! read_dns_bare : If .true. the code tries to read three files in DFPT+U calculations: 197 ! dnsorth, dnsbare, d2nsbare 198 ! d2ns_type : DFPT+U - the 2nd bare derivative of occupation matrices ns 199 ! (d2ns_bare matrix). Experimental! This is why it is not documented in Doc. 200 ! d2ns_type='full': matrix calculated with no approximation. 201 ! d2ns_type='fmmp': assume a m <=> m' symmetry. 202 ! d2ns_type='diag': if okvan=.true. the matrix is calculated retaining only 203 ! for <\beta_J|\phi_I> products where for J==I. 204 ! d2ns_type='dmmp': same as 'diag', but also assuming a m <=> m'. 205 ! diagonalization : diagonalization method used in the nscf calc 206 ! ldvscf_interpolate: if .true., use Fourier interpolation of phonon potential 207 ! do_long_range: if .true., add the long-range part of the potential to dvscf 208 ! do_charge_neutral: if .true., impose the neutrality condition on Born effective charges 209 ! wpot_dir: folder where w_pot binary files are located 210 ! ahc_dir: Directory where the output binary files for AHC e-ph coupling are written 211 ! ahc_nbnd: Number of bands where the electron self-energy is to be computed. 212 ! ahc_nbndskip: Number of bands to exclude when computing the self-energy. 213 ! skip_upperfan: If .true., skip the calculation of upper Fan self-energy. 214 ! 215 ! Note: meta_ionode is a single processor that reads the input 216 ! (ionode is also a single processor but per image) 217 ! Data read from input is subsequently broadcast to all processors 218 ! from ionode_id (using the default communicator world_comm) 219 ! 220 ! l_head : if true calculates the head of the symmetrized dielectric matrix -1 221 ! n_gauss : number of frequency steps for head calculation 222 ! omega_gauss : period for frequency calculation 223 ! grid_type : 0 GL -T,T 2 GL 0 T 224 225 IF (meta_ionode) THEN 226 ! 227 ! ... Input from file ? 228 ! 229 CALL input_from_file ( ) 230 ! 231 ! ... Read the first line of the input file 232 ! 233 READ( 5, '(A)', IOSTAT = ios ) title 234 ! 235 ENDIF 236 ! 237 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 238 CALL errore( 'phq_readin', 'reading title ', ABS( ios ) ) 239 CALL mp_bcast(title, meta_ionode_id, world_comm ) 240 ! 241 ! Rewind the input if the title is actually the beginning of inputph namelist 242 ! 243 IF( imatches("&inputph", title) ) THEN 244 WRITE(stdout,'(6x,a)') "Title line not specified: using 'default'." 245 title='default' 246 IF (meta_ionode) REWIND(5, iostat=ios) 247 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 248 CALL errore('phq_readin', 'Title line missing from input.', abs(ios)) 249 ENDIF 250 ! 251 ! ... set default values for variables in namelist 252 ! 253 tr2_ph = 1.D-12 254 eth_rps = 1.D-9 255 eth_ns = 1.D-12 256 amass(:) = 0.D0 257 alpha_mix(:) = 0.D0 258 alpha_mix(1) = 0.7D0 259 niter_ph = maxter 260 nmix_ph = 4 261 nat_todo = 0 262 modenum = 0 263 iverbosity = 1234567 264 verbosity = 'default' 265 trans = .TRUE. 266 lrpa = .FALSE. 267 lnoloc = .FALSE. 268 epsil = .FALSE. 269 zeu = .TRUE. 270 zue = .FALSE. 271 fpol = .FALSE. 272 electron_phonon=' ' 273 elph_nbnd_min = 1 274 elph_nbnd_max = 0 275 el_ph_sigma = 0.02 276 el_ph_nsigma = 10 277 el_ph_ngauss = 1 278 lraman = .FALSE. 279 elop = .FALSE. 280 max_seconds = 1.E+7_DP 281 reduce_io = .FALSE. 282 CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) 283 IF ( TRIM( outdir ) == ' ' ) outdir = './' 284 prefix = 'pwscf' 285 fildyn = 'matdyn' 286 fildrho = ' ' 287 fildvscf = ' ' 288 ldisp = .FALSE. 289 nq1 = 0 290 nq2 = 0 291 nq3 = 0 292 dek = 1.0d-3 293 nogg = .FALSE. 294 recover = .FALSE. 295 asr = .FALSE. 296 start_irr = 0 297 last_irr = -1000 298 start_q = 1 299 last_q =-1000 300 ldiag =.FALSE. 301 lqdir =.FALSE. 302 qplot =.FALSE. 303 q_in_band_form=.FALSE. 304 q2d = .FALSE. 305 only_init = .FALSE. 306 only_wfc = .FALSE. 307 low_directory_check=.FALSE. 308 search_sym =.TRUE. 309 nk1 = 0 310 nk2 = 0 311 nk3 = 0 312 k1 = 0 313 k2 = 0 314 k3 = 0 315 ! 316 ! dvscf_interpolate 317 ldvscf_interpolate = .FALSE. 318 do_charge_neutral = .FALSE. 319 do_long_range = .FALSE. 320 wpot_dir = ' ' 321 ! 322 ! electron_phonon == 'ahc' 323 ahc_dir = ' ' 324 ahc_nbnd = 0 325 ahc_nbndskip = 0 326 skip_upperfan = .FALSE. 327 elph_ahc = .FALSE. 328 ! 329 drho_star%open = .FALSE. 330 drho_star%basis = 'modes' 331 drho_star%pat = .TRUE. 332 drho_star%ext = 'drho' 333 334 CALL get_environment_variable( 'ESPRESSO_FILDRHO_DIR', drho_star%dir) 335 IF ( TRIM( drho_star%dir ) == ' ' ) & 336 drho_star%dir = TRIM(outdir)//"/Rotated_DRHO/" 337 ! 338 dvscf_star%open = .FALSE. 339 dvscf_star%basis = 'modes' 340 dvscf_star%pat = .FALSE. 341 dvscf_star%ext = 'dvscf' 342 CALL get_environment_variable( 'ESPRESSO_FILDVSCF_DIR', dvscf_star%dir) 343 IF ( TRIM( dvscf_star%dir ) == ' ' ) & 344 dvscf_star%dir = TRIM(outdir)//"/Rotated_DVSCF/" 345 ! 346 lshift_q = .false. 347 read_dns_bare =.false. 348 d2ns_type = 'full' 349 len_head_block_freq=0 350 len_head_block_wfc=0 351 ! 352 ! ... reading the namelist inputph 353 ! 354 IF (meta_ionode) THEN 355 READ( 5, INPUTPH, ERR=30, IOSTAT = ios ) 356 ! 357 ! ... iverbosity/verbosity hack 358 ! 359 IF ( iverbosity == 1234567 ) THEN 360 SELECT CASE( trim( verbosity ) ) 361 CASE( 'debug', 'high', 'medium' ) 362 iverbosity = 1 363 CASE( 'low', 'default', 'minimal' ) 364 iverbosity = 0 365 CASE DEFAULT 366 iverbosity = 0 367 END SELECT 368 ELSE 369 ios = 1234567 370 END IF 371 END IF 37230 CONTINUE 373 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 374 IF ( ios == 1234567 ) THEN 375 CALL infomsg( 'phq_readin' , & 376 'iverbosity is obsolete, use "verbosity" instead' ) 377 ELSE IF ( ABS(ios) /= 0 ) THEN 378 CALL errore( 'phq_readin', 'reading inputph namelist', ABS( ios ) ) 379 END IF 380 ! diagonalization option 381 SELECT CASE(TRIM(diagonalization)) 382 CASE ('david','davidson') 383 isolve = 0 384 CASE ('cg') 385 isolve = 1 386 CASE DEFAULT 387 CALL errore('phq_readin','diagonalization '//trim(diagonalization)//' not implemented',1) 388 END SELECT 389 390 ! 391 ! ... broadcast all input variables 392 ! 393 tmp_dir = trimcheck (outdir) 394 CALL bcast_ph_input ( ) 395 CALL mp_bcast(nogg, meta_ionode_id, world_comm ) 396 CALL mp_bcast(q2d, meta_ionode_id, world_comm ) 397 CALL mp_bcast(q_in_band_form, meta_ionode_id, world_comm ) 398 ! 399 drho_star%dir=trimcheck(drho_star%dir) 400 dvscf_star%dir=trimcheck(dvscf_star%dir) 401 ! filename for the star must always be automatically generated: 402 IF(drho_star%ext(1:5)/='auto:') drho_star%ext = 'auto:'//drho_star%ext 403 IF(dvscf_star%ext(1:5)/='auto:') dvscf_star%ext = 'auto:'//dvscf_star%ext 404 ! 405 ! ... Check all namelist variables 406 ! 407 IF (tr2_ph <= 0.D0) CALL errore (' phq_readin', ' Wrong tr2_ph ', 1) 408 IF (eth_rps<= 0.D0) CALL errore ( 'phq_readin', ' Wrong eth_rps', 1) 409 IF (eth_ns <= 0.D0) CALL errore ( 'phq_readin', ' Wrong eth_ns ', 1) 410 411 DO iter = 1, maxter 412 IF (alpha_mix (iter) .LT.0.D0.OR.alpha_mix (iter) .GT.1.D0) CALL & 413 errore ('phq_readin', ' Wrong alpha_mix ', iter) 414 ENDDO 415 IF (niter_ph.LT.1.OR.niter_ph.GT.maxter) CALL errore ('phq_readin', & 416 ' Wrong niter_ph ', 1) 417 IF (nmix_ph.LT.1.OR.nmix_ph.GT.5) CALL errore ('phq_readin', ' Wrong & 418 &nmix_ph ', 1) 419 IF (iverbosity.NE.0.AND.iverbosity.NE.1) CALL errore ('phq_readin', & 420 &' Wrong iverbosity ', 1) 421 IF (fildyn.EQ.' ') CALL errore ('phq_readin', ' Wrong fildyn ', 1) 422 IF (max_seconds.LT.0.1D0) CALL errore ('phq_readin', ' Wrong max_seconds', 1) 423 424 IF (only_init.AND.only_wfc) CALL errore('phq_readin', & 425 'only_init or only_wfc can be .true.', 1) 426 427 IF (modenum < 0) CALL errore ('phq_readin', ' Wrong modenum ', 1) 428 IF (dek <= 0.d0) CALL errore ( 'phq_readin', ' Wrong dek ', 1) 429 ! 430 ! 431 elph_tetra = 0 432 SELECT CASE( trim( electron_phonon ) ) 433 CASE( 'simple' ) 434 elph=.true. 435 elph_mat=.false. 436 elph_simple=.true. 437 elph_epa=.false. 438 CASE( 'epa' ) 439 elph=.true. 440 elph_mat=.false. 441 elph_simple=.false. 442 elph_epa=.true. 443 CASE( 'Wannier' ) 444 elph=.true. 445 elph_mat=.true. 446 elph_simple=.false. 447 elph_epa=.false. 448 auxdvscf=trim(fildvscf) 449 CASE( 'interpolated' ) 450 elph=.true. 451 elph_mat=.false. 452 elph_simple=.false. 453 elph_epa=.false. 454 ! YAMBO > 455 CASE( 'yambo' ) 456 elph=.true. 457 elph_mat=.false. 458 elph_simple=.false. 459 elph_epa=.false. 460 elph_yambo=.true. 461 nogg=.true. 462 auxdvscf=trim(fildvscf) 463 CASE( 'dvscf' ) 464 elph=.false. 465 elph_mat=.false. 466 elph_simple=.false. 467 elph_epa=.false. 468 elph_yambo=.false. 469 dvscf_yambo=.true. 470 nogg=.true. 471 auxdvscf=trim(fildvscf) 472 ! YAMBO < 473 CASE( 'lambda_tetra' ) 474 elph=.true. 475 elph_mat=.false. 476 elph_simple=.false. 477 trans = .false. 478 elph_tetra = 1 479 CASE( 'gamma_tetra' ) 480 elph=.true. 481 elph_mat=.false. 482 elph_simple=.false. 483 trans = .false. 484 elph_tetra = 2 485 CASE( 'scdft_input' ) 486 elph=.true. 487 elph_mat=.false. 488 elph_simple=.false. 489 trans = .false. 490 elph_tetra = 3 491 CASE( 'ahc' ) 492 elph = .true. 493 elph_ahc = .true. 494 elph_mat = .false. 495 elph_simple = .false. 496 elph_epa = .false. 497 trans = .false. 498 nogg = .true. 499 CASE DEFAULT 500 elph=.false. 501 elph_mat=.false. 502 elph_simple=.false. 503 elph_epa=.false. 504 END SELECT 505 506 ! YAMBO > 507 IF (.not.elph_yambo.and..not.dvscf_yambo) then 508 ! YAMBO < 509 IF (qplot.AND..NOT.ldisp) CALL errore('phq_readin','qplot requires ldisp=.true.',1) 510 ! 511 ENDIF 512 513 IF (ldisp.AND.only_init.AND.(.NOT.lqdir)) & 514 CALL errore('phq_readin', & 515 'only_init=.TRUE. requires lqdir=.TRUE. or data are lost',1) 516 517 epsil = epsil .OR. lraman .OR. elop 518 519 IF (modenum /= 0) search_sym=.FALSE. 520 521 IF (elph_mat .OR. elph_ahc) THEN 522 trans = .FALSE. 523 ELSEIF (.NOT. elph) THEN 524 trans = trans .OR. ldisp 525 ENDIF 526 ! 527 ! Set default value for fildrho and fildvscf if they are required 528 IF ( (lraman.OR.elop.OR.drho_star%open) .AND. fildrho == ' ') fildrho = 'drho' 529 IF ( (elph_mat.OR.dvscf_star%open) .AND. fildvscf == ' ') fildvscf = 'dvscf' 530 ! 531 ! We can calculate dielectric, raman or elop tensors and no Born effective 532 ! charges dF/dE, but we cannot calculate Born effective charges dF/dE 533 ! without epsil. 534 ! 535 IF (zeu) zeu = epsil 536 ! 537 ! reads the q point (just if ldisp = .false.) 538 ! 539 IF (meta_ionode) THEN 540 ios = 0 541 IF (qplot) THEN 542 READ (5, *, iostat = ios) nqaux 543 ELSE 544 IF (.NOT. ldisp) READ (5, *, iostat = ios) (xq (ipol), ipol = 1, 3) 545 ENDIF 546 END IF 547 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 548 CALL errore ('phq_readin', 'reading xq', ABS (ios) ) 549 IF (qplot) THEN 550 CALL mp_bcast(nqaux, meta_ionode_id, world_comm ) 551 ALLOCATE(xqaux(3,nqaux)) 552 ALLOCATE(wqaux(nqaux)) 553 IF (meta_ionode) THEN 554 DO iq=1, nqaux 555 READ (5, *, iostat = ios) (xqaux (ipol,iq), ipol = 1, 3), wqaux(iq) 556 ENDDO 557 ENDIF 558 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 559 CALL errore ('phq_readin', 'reading xq', ABS (ios) ) 560 CALL mp_bcast(xqaux, meta_ionode_id, world_comm ) 561 CALL mp_bcast(wqaux, meta_ionode_id, world_comm ) 562 ELSE 563 CALL mp_bcast(xq, meta_ionode_id, world_comm ) 564 ENDIF 565 566 IF (.NOT.ldisp) THEN 567 lgamma = xq (1) .EQ.0.D0.AND.xq (2) .EQ.0.D0.AND.xq (3) .EQ.0.D0 568 IF ( (epsil.OR.zue.or.lraman.or.elop) .AND..NOT.lgamma) & 569 CALL errore ('phq_readin', 'gamma is needed for elec.field', 1) 570 ENDIF 571 IF (zue.AND..NOT.trans) CALL errore ('phq_readin', 'trans must be & 572 &.t. for Zue calc.', 1) 573 574 IF (trans.AND.(lrpa.OR.lnoloc)) CALL errore('phq_readin', & 575 'only dielectric constant with lrpa or lnoloc',1) 576 IF (lrpa.or.lnoloc) THEN 577 zeu=.FALSE. 578 lraman=.FALSE. 579 elop = .FALSE. 580 ENDIF 581 ! 582 ! dvscf_interpolate 583 ! 584 IF (ldvscf_interpolate) THEN 585 ! 586 IF (wpot_dir == ' ') wpot_dir = TRIM(tmp_dir) // "w_pot/" 587 wpot_dir = trimcheck(wpot_dir) 588 ! 589 IF (do_charge_neutral .AND. (.NOT. do_long_range)) THEN 590 WRITE(stdout, '(5x,a)') 'charge neutrality for dvscf_interpolate is & 591 & meaningful only if do_long_range is true.' 592 WRITE(stdout, '(5x,a)') 'Set do_charge_neutral = .false.' 593 do_charge_neutral = .FALSE. 594 ENDIF 595 ! 596 ENDIF 597 ! 598 IF (trans .AND. ldvscf_interpolate) CALL errore ('phq_readin', & 599 'ldvscf_interpolate should be used only when trans = .false.', 1) 600 IF (domag .AND. ldvscf_interpolate) CALL errore ('phq_readin', & 601 'ldvscf_interpolate and magnetism not implemented', 1) 602 IF (okpaw .AND. ldvscf_interpolate) CALL errore ('phq_readin', & 603 'PAW and ldvscf_interpolate not tested.', 1) 604 ! 605 ! AHC e-ph coupling 606 ! 607 IF (elph_ahc) THEN 608 ! 609 IF (ahc_nbnd <= 0) CALL errore('phq_readin', & 610 'ahc_nbnd must be specified as a positive integer') 611 IF (ahc_nbndskip < 0) CALL errore('phq_readin', & 612 'ahc_nbndskip cannot be negative') 613 ! 614 IF (ahc_dir == ' ') ahc_dir = TRIM(tmp_dir) // "ahc_dir/" 615 ahc_dir = trimcheck(ahc_dir) 616 ! 617 ENDIF ! elph_ahc 618 ! 619 IF (elph_ahc .AND. trans) CALL errore ('phq_readin', & 620 'elph_ahc can be used only when trans = .false.', 1) 621 ! 622 ! reads the frequencies ( just if fpol = .true. ) 623 ! 624 IF ( fpol ) THEN 625 IF ( .NOT. epsil) CALL errore ('phq_readin', & 626 'fpol=.TRUE. needs epsil=.TRUE.', 1 ) 627 nfs=0 628 IF (meta_ionode) THEN 629 READ (5, *, iostat = ios) card 630 IF ( TRIM(card)=='FREQUENCIES'.OR. & 631 TRIM(card)=='frequencies'.OR. & 632 TRIM(card)=='Frequencies') THEN 633 READ (5, *, iostat = ios) nfs 634 ENDIF 635 ENDIF 636 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 637 CALL errore ('phq_readin', 'reading number of FREQUENCIES', ABS(ios) ) 638 CALL mp_bcast(nfs, meta_ionode_id, world_comm ) 639 if (nfs < 1) call errore('phq_readin','Too few frequencies',1) 640 ALLOCATE(fiu(nfs)) 641 IF (meta_ionode) THEN 642 IF ( TRIM(card) == 'FREQUENCIES' .OR. & 643 TRIM(card) == 'frequencies' .OR. & 644 TRIM(card) == 'Frequencies' ) THEN 645 DO i = 1, nfs 646 READ (5, *, iostat = ios) fiu(i) 647 END DO 648 END IF 649 END IF 650 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 651 CALL errore ('phq_readin', 'reading FREQUENCIES card', ABS(ios) ) 652 CALL mp_bcast(fiu, meta_ionode_id, world_comm ) 653 ELSE 654 nfs=1 655 ALLOCATE(fiu(1)) 656 fiu=0.0_DP 657 END IF 658 ! 659 ! 660 ! Here we finished the reading of the input file. 661 ! Now allocate space for pwscf variables, read and check them. 662 ! 663 ! amass will be read again from file: 664 ! save its content in auxiliary variables 665 ! 666 ALLOCATE ( amass_input( SIZE(amass) ) ) 667 amass_input(:)= amass(:) 668 ! 669 tmp_dir_save=tmp_dir 670 tmp_dir_ph= trimcheck( TRIM (tmp_dir) // '_ph' // int_to_char(my_image_id) ) 671 CALL check_tempdir ( tmp_dir_ph, exst, parallelfs ) 672 tmp_dir_phq=tmp_dir_ph 673 674 ext_restart=.FALSE. 675 ext_recover=.FALSE. 676 rec_code_read=-1000 677 IF (recover) THEN 678! 679! With a recover run we read here the mesh of q points, the current iq, 680! and the current frequency 681! 682 CALL ph_readfile('init', 0, 0, ierr) 683 CALL ph_readfile('status_ph', 0, 0, ierr1) 684! 685! If some error occured here, we cannot recover the run 686! 687 IF (ierr /= 0 .OR. ierr1 /= 0 ) THEN 688 write(stdout,'(5x,"Run is not recoverable starting from scratch")') 689 recover=.FALSE. 690 goto 1001 691 ENDIF 692! 693! We check if the bands and the information on the pw run are in the directory 694! written by the phonon code for the current q point. If the file exists 695! we read from there, otherwise use the information in tmp_dir. 696! 697 IF (lqdir) THEN 698 tmp_dir_phq= trimcheck ( TRIM(tmp_dir_ph) // TRIM(prefix) // & 699 & '.q_' // int_to_char(current_iq) ) 700 CALL check_restart_recover(ext_recover, ext_restart) 701 IF (.NOT.ext_recover.AND..NOT.ext_restart) tmp_dir_phq=tmp_dir_ph 702 ENDIF 703 ! 704 filename=TRIM(tmp_dir_phq)//TRIM(prefix)//postfix//xmlpun_schema 705 IF (ionode) inquire (file =TRIM(filename), exist = exst) 706 ! 707 CALL mp_bcast( exst, ionode_id, intra_image_comm ) 708 ! 709 ! If this file exist we use it to recover the pw.x informations 710 ! 711 IF (exst) tmp_dir=tmp_dir_phq 712 u_from_file=.true. 713 ENDIF 7141001 CONTINUE 715 716 IF (qplot.AND..NOT.recover) THEN 717 IF (q2d) THEN 718 nqs=wqaux(2)*wqaux(3) 719 ALLOCATE(x_q(3,nqs)) 720 ALLOCATE(wq(nqs)) 721 CALL generate_k_in_plane(nqaux, xqaux, wqaux, x_q, wq, nqs) 722 ELSEIF (q_in_band_form) THEN 723 nqs=SUM(wqaux(1:nqaux-1))+1 724 DO i=1,nqaux-1 725 IF (wqaux(i)==0) nqs=nqs+1 726 ENDDO 727 ALLOCATE(x_q(3,nqs)) 728 ALLOCATE(wq(nqs)) 729 CALL generate_k_along_lines(nqaux, xqaux, wqaux, x_q, wq, nqs) 730 ELSE 731 nqs=nqaux 732 ALLOCATE(x_q(3,nqs)) 733 ALLOCATE(wq(nqs)) 734 wq(:)=wqaux(:) 735 x_q(:,1:nqs)=xqaux(:,1:nqs) 736 ENDIF 737 DEALLOCATE(xqaux) 738 DEALLOCATE(wqaux) 739 ALLOCATE(lgamma_iq(nqs)) 740 DO iq=1, nqs 741 lgamma_iq(iq)= ( ABS(x_q(1,iq)) .LT. 1.0e-10_dp ) .AND. & 742 ( ABS(x_q(2,iq)) .LT. 1.0e-10_dp ) .AND. & 743 ( ABS(x_q(3,iq)) .LT. 1.0e-10_dp ) 744 ENDDO 745 WRITE(stdout, '(//5x,"Dynamical matrices for q-points given in input")') 746 WRITE(stdout, '(5x,"(",i4,"q-points):")') nqs 747 WRITE(stdout, '(5x," N xq(1) xq(2) xq(3) " )') 748 DO iq = 1, nqs 749 WRITE(stdout, '(5x,i3, 3f14.9)') iq, x_q(1,iq), x_q(2,iq), x_q(3,iq) 750 END DO 751 ENDIF 752 ! 753 ! DFPT+U: the occupation matrix ns is read via read_file 754 ! 755 CALL read_file ( ) 756 ! 757 magnetic_sym=noncolin .AND. domag 758 ! 759 ! init_start_grid returns .true. if a new k-point grid is set from values 760 ! read from input (this happens if nk1*nk2*nk3 > 0; otherwise reset_grid 761 ! returns .false., leaves the current values, read in read_file, unchanged) 762 ! 763 newgrid = reset_grid (nk1, nk2, nk3, k1, k2, k3) 764 ! 765 tmp_dir=tmp_dir_save 766 ! 767 IF (modenum > 3*nat) CALL errore ('phq_readin', ' Wrong modenum ', 2) 768 769 IF (gamma_only) CALL errore('phq_readin',& 770 'cannot start from pw.x data file using Gamma-point tricks',1) 771 772 IF (lda_plus_u) THEN 773 ! 774 WRITE(stdout,'(/5x,a)') "Phonon calculation with DFPT+U; please cite" 775 WRITE(stdout,'(5x,a)') "A. Floris et al., Phys. Rev. B 84, 161102(R) (2011)" 776 WRITE(stdout,'(5x,a)') "A. Floris et al., Phys. Rev. B 101, 064305 (2020)" 777 WRITE(stdout,'(5x,a)') "in publications or presentations arising from this work." 778 ! 779 IF (U_projection.NE."atomic") CALL errore("phq_readin", & 780 " The phonon code for this U_projection_type is not implemented",1) 781 IF (lda_plus_u_kind.NE.0) CALL errore("phq_readin", & 782 " The phonon code for this lda_plus_u_kind is not implemented",1) 783 IF (elph) CALL errore("phq_readin", & 784 " Electron-phonon with Hubbard U is not supported",1) 785 IF (lraman) CALL errore("phq_readin", & 786 " The phonon code with Raman and Hubbard U is not implemented",1) 787 ! 788 ENDIF 789 ! checks 790 IF (elph_ahc .AND. domag) CALL errore ('phq_readin', & 791 'elph_ahc and magnetism not implemented', 1) 792 IF (elph_ahc .AND. okpaw) CALL errore ('phq_readin', & 793 'elph_ahc and PAW not tested.', 1) 794 IF (elph_ahc .AND. okvan) CALL errore ('phq_readin', & 795 'elph_ahc and PAW not tested.', 1) 796 IF (elph_ahc .AND. lda_plus_u) CALL errore ('phq_readin', & 797 'elph_ahc and lda_plus_u not tested.', 1) 798 799 IF (ts_vdw) CALL errore('phq_readin',& 800 'The phonon code with TS-VdW is not yet available',1) 801 802 IF (ldftd3) CALL errore('phq_readin',& 803 'The phonon code with Grimme''s DFT-D3 is not yet available',1) 804 805 IF ( dft_is_meta() ) CALL errore('phq_readin',& 806 'The phonon code with meta-GGA functionals is not yet available',1) 807 808 IF ( dft_is_hybrid() ) CALL errore('phq_readin',& 809 'The phonon code with hybrid functionals is not yet available',1) 810 811 IF (okpaw.and.(lraman.or.elop)) CALL errore('phq_readin',& 812 'The phonon code with paw and raman or elop is not yet available',1) 813 814 IF (magnetic_sym) THEN 815 816 WRITE(stdout,'(/5x,a)') "Phonon calculation in the non-collinear magnetic case;" 817 WRITE(stdout,'(5x,a)') "please cite A. Urru and A. Dal Corso, Phys. Rev. B 100," 818 WRITE(stdout,'(5x,a)') "045115 (2019) for the theoretical background." 819 820 IF (epsil) CALL errore('phq_readin',& 821 'The calculation of Born effective charges in the non collinear & 822 magnetic case does not work yet and is temporarily disabled',1) 823 824 IF (okpaw) CALL errore('phq_readin',& 825 'The phonon code with paw and domag is not available yet',1) 826 ENDIF 827 828 IF (okvan.and.(lraman.or.elop)) CALL errore('phq_readin',& 829 'The phonon code with US-PP and raman or elop not yet available',1) 830 831 IF (noncolin.and.(lraman.or.elop)) CALL errore('phq_readin', & 832 'lraman, elop, and noncolin not programed',1) 833 834 IF (lmovecell) CALL errore('phq_readin', & 835 'The phonon code is not working after vc-relax',1) 836 837 IF (reduce_io) io_level=1 838 839 if(elph_mat.and.fildvscf.eq.' ') call errore('phq_readin',& 840 'el-ph with wannier requires fildvscf',1) 841 842 IF(elph_mat.and.npool.ne.1) call errore('phq_readin',& 843 'el-ph with wannier : pools not implemented',1) 844 845 IF(elph.and.nimage>1) call errore('phq_readin',& 846 'el-ph with images not implemented',1) 847 848 IF (elph.OR.fildvscf /= ' ') lqdir=.TRUE. 849 850 IF(dvscf_star%open.and.nimage>1) CALL errore('phq_readin',& 851 'dvscf_star with image parallelization is not yet available',1) 852 IF(drho_star%open.and.nimage>1) CALL errore('phq_readin',& 853 'drho_star with image parallelization is not yet available',1) 854 855 IF (lda_plus_u .AND. read_dns_bare .AND. ldisp) lqdir=.TRUE. 856 857 IF (.NOT.ldisp) lqdir=.FALSE. 858 859 IF (i_cons /= 0) & 860 CALL errore('phq_readin',& 861 'The phonon code with constrained magnetization is not yet available',1) 862 863 IF (two_fermi_energies .AND. (ltetra .OR. lgauss)) & 864 CALL errore('phq_readin',& 865 'The phonon code with two fermi energies is not available for metals',1) 866 867 IF (tqr) CALL errore('phq_readin',& 868 'The phonon code with Q in real space not available',1) 869 870 IF (start_irr < 0 ) CALL errore('phq_readin', 'wrong start_irr',1) 871 ! 872 IF (start_q <= 0 ) CALL errore('phq_readin', 'wrong start_q',1) 873 ! 874 ! the dynamical matrix is written in xml format if fildyn ends in 875 ! .xml or in the noncollinear case. 876 ! 877 xmldyn=has_xml(fildyn) 878 !IF (noncolin) xmldyn=.TRUE. 879 ! 880 ! If a band structure calculation needs to be done do not open a file 881 ! for k point 882 ! 883 restart = recover 884 ! 885 ! set masses to values read from input, if available; 886 ! leave values read from file otherwise 887 ! 888 DO it = 1, ntyp 889 IF (amass_input(it) < 0.0_DP) amass_input(it)= & 890 atom_weight(atomic_number(TRIM(atm(it)))) 891 IF (amass_input(it) > 0.D0) amass(it) = amass_input(it) 892 IF (amass(it) <= 0.D0) CALL errore ('phq_readin', 'Wrong masses', it) 893 ENDDO 894 DEALLOCATE (amass_input) 895 lgamma_gamma=.FALSE. 896 IF (.NOT.ldisp) THEN 897 IF (nkstot==1.OR.(nkstot==2.AND.nspin==2)) THEN 898 lgamma_gamma=(lgamma.AND.(ABS(xk(1,1))<1.D-12) & 899 .AND.(ABS(xk(2,1))<1.D-12) & 900 .AND.(ABS(xk(3,1))<1.D-12) ) 901 ENDIF 902 IF (nogg) lgamma_gamma=.FALSE. 903 IF ((nat_todo /= 0) .and. lgamma_gamma) CALL errore( & 904 'phq_readin', 'gamma_gamma tricks with nat_todo & 905 & not available. Use nogg=.true.', 1) 906 ! 907 IF (nimage > 1 .AND. lgamma_gamma) CALL errore( & 908 'phq_readin','gamma_gamma tricks with images not implemented',1) 909 IF (lgamma) THEN 910 nksq = nks 911 ELSE 912 nksq = nks / 2 913 ENDIF 914 ENDIF 915 IF (lgamma_gamma.AND.ldiag) CALL errore('phq_readin','incompatible flags',1) 916 IF (lgamma_gamma.AND.elph ) CALL errore('phq_readin','lgamma_gamma and electron-phonon are incompatible',1) 917 ! 918 IF (tfixed_occ) & 919 CALL errore('phq_readin','phonon with arbitrary occupations not tested',1) 920 ! 921 !YAMBO > 922 IF (elph .AND. .NOT.(lgauss .OR. ltetra) & 923 .AND. .NOT. (elph_yambo .OR. elph_ahc)) & 924 CALL errore ('phq_readin', 'Electron-phonon only for metals', 1) 925 !YAMBO < 926 IF (elph .AND. fildvscf.EQ.' ' .AND. .NOT. ldvscf_interpolate) & 927 CALL errore ('phq_readin', 'El-ph needs a DeltaVscf file', 1) 928 ! There might be other variables in the input file which describe 929 ! partial computation of the dynamical matrix. Read them here 930 ! 931 CALL allocate_part ( nat ) 932 ! 933 IF ( nat_todo < 0 .OR. nat_todo > nat ) & 934 CALL errore ('phq_readin', 'nat_todo is wrong', 1) 935 IF (nat_todo.NE.0) THEN 936 IF (meta_ionode) READ (5, *, iostat = ios) (atomo (na), na = 1, nat_todo) 937 CALL mp_bcast(ios, meta_ionode_id, world_comm ) 938 CALL errore ('phq_readin', 'reading atoms', ABS (ios) ) 939 CALL mp_bcast(atomo, meta_ionode_id, world_comm ) 940 ENDIF 941 nat_todo_input=nat_todo 942 943 IF (epsil.AND.(lgauss .OR. ltetra)) & 944 CALL errore ('phq_readin', 'no elec. field with metals', 1) 945 IF (modenum > 0) THEN 946 IF ( ldisp ) & 947 CALL errore('phq_readin','Dispersion calculation and & 948 & single mode calculation not possibile !',1) 949 nat_todo = 0 950 ENDIF 951 ! 952 ! Tetrahedron method for DFPT and El-Ph 953 ! 954 IF(ltetra .AND. tetra_type == 0) & 955 & CALL errore ('phq_readin', 'DFPT with the Blochl correction is not implemented', 1) 956 ! 957 IF(.NOT. ltetra .AND. elph_tetra /= 0) & 958 & CALL errore ('phq_readin', '"electron_phonon" and "occupation" are inconsistent.', 1) 959 ! 960 IF (modenum > 0 .OR. lraman ) lgamma_gamma=.FALSE. 961 IF (.NOT.lgamma_gamma) asr=.FALSE. 962 ! 963 IF ((ldisp.AND..NOT.qplot) .AND. & 964 (nq1 .LE. 0 .OR. nq2 .LE. 0 .OR. nq3 .LE. 0)) & 965 CALL errore('phq_readin','nq1, nq2, and nq3 must be greater than 0',1) 966 967 CALL save_ph_input_variables() 968 ! 969 RETURN 970 ! 971END SUBROUTINE phq_readin 972