! ! Copyright (C) 2001-2020 Quantum ESPRESSO group ! This file is distributed under the terms of the ! GNU General Public License. See the file `License' ! in the root directory of the present distribution, ! or http://www.gnu.org/copyleft/gpl.txt . ! ! !---------------------------------------------------------------------------- SUBROUTINE phq_readin() !---------------------------------------------------------------------------- !! This routine reads the control variables for the program \(\texttt{phononq}\). !! A second routine, \(\texttt{read_file_ph}\), reads the variables saved to file !! by the self-consistent program. ! ! USE kinds, ONLY : DP USE ions_base, ONLY : nat, ntyp => nsp USE mp, ONLY : mp_bcast USE mp_world, ONLY : world_comm USE ions_base, ONLY : amass, atm USE check_stop, ONLY : max_seconds USE start_k, ONLY : reset_grid USE klist, ONLY : xk, nks, nkstot, lgauss, two_fermi_energies, ltetra USE control_flags, ONLY : gamma_only, tqr, restart, io_level, & ts_vdw, ldftd3, lxdm, isolve, dfpt_hub USE xc_lib, ONLY : xclib_dft_is USE uspp, ONLY : okvan USE fixed_occ, ONLY : tfixed_occ USE lsda_mod, ONLY : lsda, nspin USE cellmd, ONLY : lmovecell USE run_info, ONLY : title USE control_ph, ONLY : maxter, alpha_mix, lgamma_gamma, epsil, & zue, zeu, xmldyn, newgrid, & trans, reduce_io, tr2_ph, niter_ph, & nmix_ph, ldisp, recover, lnoloc, start_irr, & last_irr, start_q, last_q, current_iq, tmp_dir_ph, & ext_recover, ext_restart, u_from_file, ldiag, & search_sym, lqdir, electron_phonon, tmp_dir_phq, & rec_code_read, qplot, only_init, only_wfc, & low_directory_check, nk1, nk2, nk3, k1, k2, k3, & dftd3_hess USE save_ph, ONLY : tmp_dir_save, save_ph_input_variables USE gamma_gamma, ONLY : asr USE partial, ONLY : atomo, nat_todo, nat_todo_input USE output, ONLY : fildyn, fildvscf, fildrho USE disp, ONLY : nq1, nq2, nq3, x_q, wq, nqs, lgamma_iq USE io_files, ONLY : tmp_dir, prefix, postfix, create_directory, & check_tempdir, xmlpun_schema USE noncollin_module, ONLY : domag, i_cons, noncolin, lspinorb USE control_flags, ONLY : iverbosity, modenum USE io_global, ONLY : meta_ionode, meta_ionode_id, ionode, ionode_id, & qestdin, stdout USE mp_images, ONLY : nimage, my_image_id, intra_image_comm, & me_image, nproc_image USE mp_pools, ONLY : npool USE paw_variables, ONLY : okpaw USE ramanm, ONLY : eth_rps, eth_ns, lraman, elop, dek USE freq_ph, ONLY : fpol, fiu, nfs USE cryst_ph, ONLY : magnetic_sym USE ph_restart, ONLY : ph_readfile USE el_phon, ONLY : elph,elph_mat,elph_simple,elph_epa,elph_nbnd_min, elph_nbnd_max, & el_ph_sigma, el_ph_nsigma, el_ph_ngauss,auxdvscf USE dfile_star, ONLY : drho_star, dvscf_star USE qpoint, ONLY : nksq, xq USE control_lr, ONLY : lgamma, lrpa ! YAMBO > USE YAMBO, ONLY : elph_yambo,dvscf_yambo ! YAMBO < USE elph_tetra_mod,ONLY : elph_tetra, lshift_q, in_alpha2f USE ktetra, ONLY : tetra_type USE ldaU, ONLY : lda_plus_u, Hubbard_projectors, lda_plus_u_kind, & is_hubbard_back USE ldaU_ph, ONLY : read_dns_bare, d2ns_type USE dvscf_interpolate, ONLY : ldvscf_interpolate, do_long_range, & do_charge_neutral, wpot_dir USE ahc, ONLY : elph_ahc, ahc_dir, ahc_nbnd, ahc_nbndskip, & skip_upperfan USE read_namelists_module, ONLY : check_namelist_read USE open_close_input_file, ONLY : open_input_file, close_input_file USE el_phon, ONLY : kx, ky, kz, elph_print ! IMPLICIT NONE ! CHARACTER(LEN=256), EXTERNAL :: trimcheck ! INTEGER :: ios, ipol, iter, na, it, ierr, ierr1 ! integer variable for I/O control ! counter on polarizations ! counter on iterations ! counter on atoms ! counter on types REAL(DP), ALLOCATABLE :: amass_input(:) ! save masses read from input here CHARACTER (LEN=256) :: outdir, filename CHARACTER (LEN=8) :: verbosity CHARACTER(LEN=80) :: card CHARACTER(LEN=6) :: int_to_char INTEGER :: i LOGICAL :: nogg LOGICAL :: q2d, q_in_band_form INTEGER, EXTERNAL :: atomic_number REAL(DP), EXTERNAL :: atom_weight LOGICAL, EXTERNAL :: imatches LOGICAL, EXTERNAL :: has_xml LOGICAL :: exst, parallelfs REAL(DP), ALLOCATABLE :: xqaux(:,:) INTEGER, ALLOCATABLE :: wqaux(:) INTEGER :: nqaux, iq CHARACTER(len=80) :: diagonalization='david' LOGICAL :: needwf_ph=.TRUE. ! NAMELIST / INPUTPH / tr2_ph, amass, alpha_mix, niter_ph, nmix_ph, & nat_todo, verbosity, iverbosity, outdir, epsil, & trans, zue, zeu, max_seconds, reduce_io, & modenum, prefix, fildyn, fildvscf, fildrho, & ldisp, nq1, nq2, nq3, & eth_rps, eth_ns, lraman, elop, dek, recover, & fpol, asr, lrpa, lnoloc, start_irr, last_irr, & start_q, last_q, nogg, ldiag, search_sym, lqdir, & nk1, nk2, nk3, k1, k2, k3, & drho_star, dvscf_star, only_init, only_wfc, & elph_nbnd_min, elph_nbnd_max, el_ph_ngauss, & el_ph_nsigma, el_ph_sigma, electron_phonon, & q_in_band_form, q2d, qplot, low_directory_check, & lshift_q, read_dns_bare, d2ns_type, diagonalization, & ldvscf_interpolate, do_long_range, do_charge_neutral, & wpot_dir, ahc_dir, ahc_nbnd, ahc_nbndskip, & skip_upperfan, dftd3_hess, kx, ky, kz ! tr2_ph : convergence threshold ! amass : atomic masses ! alpha_mix : the mixing parameter ! niter_ph : maximum number of iterations ! nmix_ph : number of previous iterations used in mixing ! nat_todo : number of atom to be displaced ! verbosity : verbosity control (iverbosity is obsolete) ! outdir : directory where input, output, temporary files reside ! epsil : if true calculate dielectric constant ! trans : if true calculate phonon ! electron-phonon : select the kind of electron-phonon calculation ! elph : if true calculate electron-phonon coefficients ! elph_mat : if true eph coefficients for wannier ! zue : if .true. calculate effective charges ( d force / dE ) ! zeu : if .true. calculate effective charges ( d P / du ) ! lraman : if true calculate raman tensor ! elop : if true calculate electro-optic tensor ! max_seconds : maximum cputime for this run ! reduce_io : reduce I/O to the strict minimum ! modenum : single mode calculation ! prefix : the prefix of files produced by pwscf ! fildyn : output file for the dynamical matrix ! fildvscf : output file containing deltavsc ! fildrho : output file containing deltarho ! fildrho_dir : directory where fildrho files will be stored (default: outdir or ESPRESSO_FILDRHO_DIR variable) ! eth_rps : threshold for calculation of Pc R |psi> (Raman) ! eth_ns : threshold for non-scf wavefunction calculation (Raman) ! dek : delta_xk used for wavefunctions derivation (Raman) ! recover : recover=.true. to restart from an interrupted run ! asr : in the gamma_gamma case apply acoustic sum rule ! start_q : in q list does the q points from start_q to last_q ! last_q : ! start_irr : does the irred. representation from start_irr to last_irr ! last_irr : ! nogg : if .true. lgamma_gamma tricks are not used ! ldiag : if .true. force diagonalization of the dyn mat ! lqdir : if .true. each q writes in its own directory ! search_sym : if .true. analyze symmetry if possible ! nk1,nk2,nk3, k1,k2,k3 : ! when specified in input, the phonon run uses a different ! k-point mesh from that used for the charge density. ! ! dvscf_star%open : if .true. write in dvscf_star%dir the dvscf_q ! 'for all q' in the star of q with suffix dvscf_star%ext. ! The dvscf_q' is written in the basis dvscf_star%basis; ! if dvscf_star%pat is .true. also save a pattern file. ! dvscf_star%dir, dvscf_star%ext, dvscf_star%basis : see dvscf_star%open ! drho_star%open : like dvscf_star%open but for drho_q ! drho_star%dir, drho_star%ext, drho_star%basis : see drho_star%open ! ! elph_nbnd_min, ! elph_nbnd_max: if (elph_mat=.true.) it dumps the eph matrix element from elph_nbnd_min ! to elph_nbnd_max ! el_ph_ngauss, ! el_ph_nsigma, ! el_ph_sigma : if (elph_mat=.true.) it defines the kind and the val-ue of the ! smearing to be used in the eph coupling calculation. ! qplot, : if true a list of q points is given in input ! q_in_band_form: if true the input list of q points defines paths ! q2d, : if .true. the q list define a mesh in a square. ! low_directory_check : if .true. only the requested representations ! are searched on file ! ! read_dns_bare : If .true. the code tries to read three files in DFPT+U calculations: ! dnsorth, dnsbare, d2nsbare ! d2ns_type : DFPT+U - the 2nd bare derivative of occupation matrices ns ! (d2ns_bare matrix). Experimental! This is why it is not documented in Doc. ! d2ns_type='full': matrix calculated with no approximation. ! d2ns_type='fmmp': assume a m <=> m' symmetry. ! d2ns_type='diag': if okvan=.true. the matrix is calculated retaining only ! for <\beta_J|\phi_I> products where for J==I. ! d2ns_type='dmmp': same as 'diag', but also assuming a m <=> m'. ! diagonalization : diagonalization method used in the nscf calc ! ldvscf_interpolate: if .true., use Fourier interpolation of phonon potential ! do_long_range: if .true., add the long-range part of the potential to dvscf ! do_charge_neutral: if .true., impose the neutrality condition on Born effective charges ! wpot_dir: folder where w_pot binary files are located ! ahc_dir: Directory where the output binary files for AHC e-ph coupling are written ! ahc_nbnd: Number of bands where the electron self-energy is to be computed. ! ahc_nbndskip: Number of bands to exclude when computing the self-energy. ! skip_upperfan: If .true., skip the calculation of upper Fan self-energy. ! ! dftd3_hess: file from where the dftd3 hessian is read ! ! Note: meta_ionode is a single processor that reads the input ! (ionode is also a single processor but per image) ! Data read from input is subsequently broadcast to all processors ! from ionode_id (using the default communicator world_comm) ! IF (meta_ionode) THEN ! ! ... Input from file (ios=0) or standard input (ios=-1) on unit "qestdin" ! ios = open_input_file ( ) ! ! ... Read the first line of the input file ! IF ( ios <= 0 ) READ( qestdin, '(A)', IOSTAT = ios ) title ! ENDIF ! CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore( 'phq_readin', 'reading input file ', ABS( ios ) ) CALL mp_bcast(title, meta_ionode_id, world_comm ) ! ! Rewind the input if the title is actually the beginning of inputph namelist ! IF( imatches("&inputph", title) ) THEN WRITE(stdout,'(6x,a)') "Title line not specified: using 'default'." title='default' IF (meta_ionode) REWIND(qestdin, iostat=ios) CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore('phq_readin', 'Title line missing from input.', abs(ios)) ENDIF ! ! ... set default values for variables in namelist ! tr2_ph = 1.D-12 eth_rps = 1.D-9 eth_ns = 1.D-12 amass(:) = 0.D0 alpha_mix(:) = 0.D0 alpha_mix(1) = 0.7D0 niter_ph = maxter nmix_ph = 4 nat_todo = 0 modenum = 0 iverbosity = 1234567 verbosity = 'default' trans = .TRUE. lrpa = .FALSE. lnoloc = .FALSE. epsil = .FALSE. zeu = .TRUE. zue = .FALSE. fpol = .FALSE. electron_phonon=' ' elph_nbnd_min = 1 elph_nbnd_max = 0 el_ph_sigma = 0.02 el_ph_nsigma = 10 el_ph_ngauss = 1 lraman = .FALSE. elop = .FALSE. max_seconds = 1.E+7_DP reduce_io = .FALSE. CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir ) IF ( TRIM( outdir ) == ' ' ) outdir = './' prefix = 'pwscf' fildyn = 'matdyn' fildrho = ' ' fildvscf = ' ' dftd3_hess = ' ' ldisp = .FALSE. nq1 = 0 nq2 = 0 nq3 = 0 dek = 1.0d-3 nogg = .FALSE. recover = .FALSE. asr = .FALSE. start_irr = 0 last_irr = -1000 start_q = 1 last_q =-1000 ldiag =.FALSE. lqdir =.FALSE. qplot =.FALSE. q_in_band_form=.FALSE. q2d = .FALSE. only_init = .FALSE. only_wfc = .FALSE. low_directory_check=.FALSE. search_sym =.TRUE. nk1 = 0 nk2 = 0 nk3 = 0 k1 = 0 k2 = 0 k3 = 0 kx = 0.D0 ky = 0.D0 kz = 0.D0 ! ! dvscf_interpolate ldvscf_interpolate = .FALSE. do_charge_neutral = .FALSE. do_long_range = .FALSE. wpot_dir = ' ' ! ! electron_phonon == 'ahc' ahc_dir = ' ' ahc_nbnd = 0 ahc_nbndskip = 0 skip_upperfan = .FALSE. ! drho_star%open = .FALSE. drho_star%basis = 'modes' drho_star%pat = .TRUE. drho_star%ext = 'drho' CALL get_environment_variable( 'ESPRESSO_FILDRHO_DIR', drho_star%dir) IF ( TRIM( drho_star%dir ) == ' ' ) & drho_star%dir = TRIM(outdir)//"/Rotated_DRHO/" ! dvscf_star%open = .FALSE. dvscf_star%basis = 'modes' dvscf_star%pat = .FALSE. dvscf_star%ext = 'dvscf' CALL get_environment_variable( 'ESPRESSO_FILDVSCF_DIR', dvscf_star%dir) IF ( TRIM( dvscf_star%dir ) == ' ' ) & dvscf_star%dir = TRIM(outdir)//"/Rotated_DVSCF/" ! lshift_q = .false. read_dns_bare =.false. d2ns_type = 'full' ! ! ... reading the namelist inputph ! IF (meta_ionode) THEN READ( qestdin, INPUTPH, IOSTAT = ios ) ! ! ... iverbosity/verbosity hack ! IF ( iverbosity == 1234567 ) THEN SELECT CASE( trim( verbosity ) ) CASE( 'debug', 'high', 'medium' ) iverbosity = 1 CASE( 'low', 'default', 'minimal' ) iverbosity = 0 CASE DEFAULT iverbosity = 0 END SELECT ELSE ios = 1234567 END IF END IF CALL check_namelist_read(ios, qestdin, "inputph") 30 CONTINUE CALL mp_bcast(ios, meta_ionode_id, world_comm ) IF ( ios == 1234567 ) THEN CALL infomsg( 'phq_readin' , & 'iverbosity is obsolete, use "verbosity" instead' ) ELSE IF ( ABS(ios) /= 0 ) THEN CALL errore( 'phq_readin', 'reading inputph namelist', ABS( ios ) ) END IF ! diagonalization option SELECT CASE(TRIM(diagonalization)) CASE ('david','davidson') isolve = 0 CASE ('cg') isolve = 1 CASE DEFAULT CALL errore('phq_readin','diagonalization '//trim(diagonalization)//' not implemented',1) END SELECT ! ! ... broadcast all input variables ! tmp_dir = trimcheck (outdir) CALL bcast_ph_input ( ) CALL mp_bcast(nogg, meta_ionode_id, world_comm ) CALL mp_bcast(q2d, meta_ionode_id, world_comm ) CALL mp_bcast(q_in_band_form, meta_ionode_id, world_comm ) CALL mp_bcast(kx, meta_ionode_id, world_comm ) CALL mp_bcast(ky, meta_ionode_id, world_comm ) CALL mp_bcast(kz, meta_ionode_id, world_comm ) ! drho_star%dir=trimcheck(drho_star%dir) dvscf_star%dir=trimcheck(dvscf_star%dir) ! filename for the star must always be automatically generated: IF(drho_star%ext(1:5)/='auto:') drho_star%ext = 'auto:'//drho_star%ext IF(dvscf_star%ext(1:5)/='auto:') dvscf_star%ext = 'auto:'//dvscf_star%ext ! ! ... Check all namelist variables ! IF (tr2_ph <= 0.D0) CALL errore (' phq_readin', ' Wrong tr2_ph ', 1) IF (eth_rps<= 0.D0) CALL errore ( 'phq_readin', ' Wrong eth_rps', 1) IF (eth_ns <= 0.D0) CALL errore ( 'phq_readin', ' Wrong eth_ns ', 1) DO iter = 1, maxter IF (alpha_mix (iter) .LT.0.D0.OR.alpha_mix (iter) .GT.1.D0) CALL & errore ('phq_readin', ' Wrong alpha_mix ', iter) ENDDO IF (niter_ph.LT.1.OR.niter_ph.GT.maxter) CALL errore ('phq_readin', & ' Wrong niter_ph ', 1) IF (iverbosity.NE.0.AND.iverbosity.NE.1) CALL errore ('phq_readin', & &' Wrong iverbosity ', 1) IF (fildyn.EQ.' ') CALL errore ('phq_readin', ' Wrong fildyn ', 1) IF (max_seconds.LT.0.1D0) CALL errore ('phq_readin', ' Wrong max_seconds', 1) IF (only_init.AND.only_wfc) CALL errore('phq_readin', & 'only_init or only_wfc can be .true.', 1) IF (modenum < 0) CALL errore ('phq_readin', ' Wrong modenum ', 1) IF (dek <= 0.d0) CALL errore ( 'phq_readin', ' Wrong dek ', 1) ! ! elph_simple= .FALSE. elph_mat = .FALSE. elph_ahc = .FALSE. elph_tetra = 0 elph_epa = .FALSE. SELECT CASE( trim( electron_phonon ) ) CASE( 'simple' ) elph=.true. elph_simple=.true. CASE( 'epa' ) elph=.true. elph_epa=.true. CASE( 'Wannier' ) elph=.true. elph_mat=.true. auxdvscf=trim(fildvscf) CASE( 'interpolated' ) elph=.true. CASE( 'yambo' ) elph=.true. elph_yambo=.true. nogg=.true. auxdvscf=trim(fildvscf) ! also for Yambo CASE( 'dvscf' ) elph=.false. dvscf_yambo=.true. nogg=.true. auxdvscf=trim(fildvscf) CASE( 'lambda_tetra' ) elph=.true. trans = .false. elph_tetra = 1 CASE( 'gamma_tetra' ) elph=.true. trans = .false. elph_tetra = 2 CASE( 'scdft_input' ) elph=.true. trans = .false. elph_tetra = 3 CASE( 'ahc' ) elph = .true. elph_ahc = .true. trans = .false. nogg = .true. CASE( 'prt' ) elph = .true. elph_print=.true. CASE DEFAULT elph=.false. END SELECT ! YAMBO > IF (.not.elph_yambo.and..not.dvscf_yambo) then ! YAMBO < IF (qplot.AND..NOT.ldisp) CALL errore('phq_readin','qplot requires ldisp=.true.',1) ! ENDIF IF (ldisp.AND.only_init.AND.(.NOT.lqdir)) & CALL errore('phq_readin', & 'only_init=.TRUE. requires lqdir=.TRUE. or data are lost',1) epsil = epsil .OR. lraman .OR. elop IF (modenum /= 0) search_sym=.FALSE. IF (elph_mat .OR. elph_ahc) THEN trans = .FALSE. ELSEIF (.NOT. elph) THEN trans = trans .OR. ldisp ENDIF ! ! Set default value for fildrho and fildvscf if they are required IF ( (lraman.OR.elop.OR.drho_star%open) .AND. fildrho == ' ') fildrho = 'drho' IF ( (elph_mat.OR.dvscf_star%open) .AND. fildvscf == ' ') fildvscf = 'dvscf' ! ! We can calculate dielectric, raman or elop tensors and no Born effective ! charges dF/dE, but we cannot calculate Born effective charges dF/dE ! without epsil. ! IF (zeu) zeu = epsil ! ! reads the q point (just if ldisp = .false.) ! IF (meta_ionode) THEN ios = 0 IF (qplot) THEN READ (qestdin, *, iostat = ios) nqaux ELSE IF (.NOT. ldisp) READ (qestdin, *, iostat = ios) (xq (ipol), ipol=1,3) ENDIF END IF CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore ('phq_readin', 'reading xq', ABS (ios) ) IF (qplot) THEN CALL mp_bcast(nqaux, meta_ionode_id, world_comm ) ALLOCATE(xqaux(3,nqaux)) ALLOCATE(wqaux(nqaux)) IF (meta_ionode) THEN DO iq=1, nqaux READ (qestdin, *, iostat = ios) (xqaux (ipol,iq), ipol=1,3), wqaux(iq) ENDDO ENDIF CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore ('phq_readin', 'reading xq', ABS (ios) ) CALL mp_bcast(xqaux, meta_ionode_id, world_comm ) CALL mp_bcast(wqaux, meta_ionode_id, world_comm ) ELSE CALL mp_bcast(xq, meta_ionode_id, world_comm ) ENDIF IF (.NOT.ldisp) THEN lgamma = xq (1) .EQ.0.D0.AND.xq (2) .EQ.0.D0.AND.xq (3) .EQ.0.D0 IF ( (epsil.OR.zue.or.lraman.or.elop) .AND..NOT.lgamma) & CALL errore ('phq_readin', 'gamma is needed for elec.field', 1) ENDIF IF (zue.AND..NOT.trans) CALL errore ('phq_readin', 'trans must be & &.t. for Zue calc.', 1) IF (trans.AND.(lrpa.OR.lnoloc)) CALL errore('phq_readin', & 'only dielectric constant with lrpa or lnoloc',1) IF (lrpa.or.lnoloc) THEN zeu=.FALSE. lraman=.FALSE. elop = .FALSE. ENDIF ! ! dvscf_interpolate ! IF (ldvscf_interpolate) THEN ! IF (wpot_dir == ' ') wpot_dir = TRIM(tmp_dir) // "w_pot/" wpot_dir = trimcheck(wpot_dir) ! IF (do_charge_neutral .AND. (.NOT. do_long_range)) THEN WRITE(stdout, '(5x,a)') 'charge neutrality for dvscf_interpolate is & & meaningful only if do_long_range is true.' WRITE(stdout, '(5x,a)') 'Set do_charge_neutral = .false.' do_charge_neutral = .FALSE. ENDIF ! ENDIF ! IF (trans .AND. ldvscf_interpolate) CALL errore ('phq_readin', & 'ldvscf_interpolate should be used only when trans = .false.', 1) IF (domag .AND. ldvscf_interpolate) CALL errore ('phq_readin', & 'ldvscf_interpolate and magnetism not implemented', 1) IF (okpaw .AND. ldvscf_interpolate) CALL errore ('phq_readin', & 'PAW and ldvscf_interpolate not tested.', 1) ! ! AHC e-ph coupling ! IF (elph_ahc) THEN ! IF (ahc_nbnd <= 0) CALL errore('phq_readin', & 'ahc_nbnd must be specified as a positive integer', 1) IF (ahc_nbndskip < 0) CALL errore('phq_readin', & 'ahc_nbndskip cannot be negative', 1) ! IF (ahc_dir == ' ') ahc_dir = TRIM(tmp_dir) // "ahc_dir/" ahc_dir = trimcheck(ahc_dir) ! ENDIF ! elph_ahc ! IF (elph_ahc .AND. trans) CALL errore ('phq_readin', & 'elph_ahc can be used only when trans = .false.', 1) ! ! reads the frequencies ( just if fpol = .true. ) ! IF ( fpol ) THEN IF ( .NOT. epsil) CALL errore ('phq_readin', & 'fpol=.TRUE. needs epsil=.TRUE.', 1 ) nfs=0 IF (meta_ionode) THEN READ (qestdin, *, iostat = ios) card IF ( TRIM(card)=='FREQUENCIES'.OR. & TRIM(card)=='frequencies'.OR. & TRIM(card)=='Frequencies') THEN READ (qestdin, *, iostat = ios) nfs ENDIF ENDIF CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore ('phq_readin', 'reading number of FREQUENCIES', ABS(ios) ) CALL mp_bcast(nfs, meta_ionode_id, world_comm ) if (nfs < 1) call errore('phq_readin','Too few frequencies',1) ALLOCATE(fiu(nfs)) IF (meta_ionode) THEN IF ( TRIM(card) == 'FREQUENCIES' .OR. & TRIM(card) == 'frequencies' .OR. & TRIM(card) == 'Frequencies' ) THEN DO i = 1, nfs READ (qestdin, *, iostat = ios) fiu(i) END DO END IF END IF CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore ('phq_readin', 'reading FREQUENCIES card', ABS(ios) ) CALL mp_bcast(fiu, meta_ionode_id, world_comm ) ELSE nfs=1 ALLOCATE(fiu(1)) fiu=0.0_DP END IF ! ! ! Here we finished the reading of the input file. ! Now allocate space for pwscf variables, read and check them. ! ! amass will be read again from file: ! save its content in auxiliary variables ! ALLOCATE ( amass_input( SIZE(amass) ) ) amass_input(:)= amass(:) ! tmp_dir_save=tmp_dir tmp_dir_ph= trimcheck( TRIM (tmp_dir) // '_ph' // int_to_char(my_image_id) ) CALL check_tempdir ( tmp_dir_ph, exst, parallelfs ) tmp_dir_phq=tmp_dir_ph ext_restart=.FALSE. ext_recover=.FALSE. rec_code_read=-1000 IF (recover) THEN ! ! With a recover run we read here the mesh of q points, the current iq, ! and the current frequency ! CALL ph_readfile('init', 0, 0, ierr) CALL ph_readfile('status_ph', 0, 0, ierr1) ! ! If some error occured here, we cannot recover the run ! IF (ierr /= 0 .OR. ierr1 /= 0 ) THEN write(stdout,'(5x,"Run is not recoverable starting from scratch")') recover=.FALSE. goto 1001 ENDIF ! ! We check if the bands and the information on the pw run are in the directory ! written by the phonon code for the current q point. If the file exists ! we read from there, otherwise use the information in tmp_dir. ! IF (lqdir) THEN tmp_dir_phq= trimcheck ( TRIM(tmp_dir_ph) // TRIM(prefix) // & & '.q_' // int_to_char(current_iq) ) CALL check_restart_recover(ext_recover, ext_restart) IF (.NOT.ext_recover.AND..NOT.ext_restart) tmp_dir_phq=tmp_dir_ph ENDIF ! filename=TRIM(tmp_dir_phq)//TRIM(prefix)//postfix//xmlpun_schema IF (ionode) inquire (file =TRIM(filename), exist = exst) ! CALL mp_bcast( exst, ionode_id, intra_image_comm ) ! ! If this file exist we use it to recover the pw.x informations ! IF (exst) tmp_dir=tmp_dir_phq u_from_file=.true. ENDIF 1001 CONTINUE IF (qplot.AND..NOT.recover) THEN IF (q2d) THEN nqs=wqaux(2)*wqaux(3) ALLOCATE(x_q(3,nqs)) ALLOCATE(wq(nqs)) CALL generate_k_in_plane(nqaux, xqaux, wqaux, x_q, wq, nqs) ELSEIF (q_in_band_form) THEN nqs=SUM(wqaux(1:nqaux-1))+1 DO i=1,nqaux-1 IF (wqaux(i)==0) nqs=nqs+1 ENDDO ALLOCATE(x_q(3,nqs)) ALLOCATE(wq(nqs)) CALL generate_k_along_lines(nqaux, xqaux, wqaux, x_q, wq, nqs) ELSE nqs=nqaux ALLOCATE(x_q(3,nqs)) ALLOCATE(wq(nqs)) wq(:)=wqaux(:) x_q(:,1:nqs)=xqaux(:,1:nqs) ENDIF DEALLOCATE(xqaux) DEALLOCATE(wqaux) ALLOCATE(lgamma_iq(nqs)) DO iq=1, nqs lgamma_iq(iq)= ( ABS(x_q(1,iq)) .LT. 1.0e-10_dp ) .AND. & ( ABS(x_q(2,iq)) .LT. 1.0e-10_dp ) .AND. & ( ABS(x_q(3,iq)) .LT. 1.0e-10_dp ) ENDDO WRITE(stdout, '(//5x,"Dynamical matrices for q-points given in input")') WRITE(stdout, '(5x,"(",i4,"q-points):")') nqs WRITE(stdout, '(5x," N xq(1) xq(2) xq(3) " )') DO iq = 1, nqs WRITE(stdout, '(5x,i3, 3f14.9)') iq, x_q(1,iq), x_q(2,iq), x_q(3,iq) END DO ENDIF ! IF (reduce_io) io_level=0 ! ! Here all needed data from the scf calculation are read ! ! IO of wfcs in collected format is not needed for recover run needwf_ph = .NOT. recover ! CALL read_file_ph ( needwf_ph ) ! magnetic_sym=noncolin .AND. domag ! ! init_start_grid returns .true. if a new k-point grid is set from values ! read from input (this happens if nk1*nk2*nk3 > 0; otherwise reset_grid ! returns .false., leaves the current values, read in read_file, unchanged) ! newgrid = reset_grid (nk1, nk2, nk3, k1, k2, k3) if(newgrid.and.elph_mat)then WRITE(stdout, '(//5x,"WARNING: Wannier elph do not use explicit new grid: nk1 nk2 nk3 ignored")') newgrid=.false. end if ! tmp_dir=tmp_dir_save ! IF (modenum > 3*nat) CALL errore ('phq_readin', ' Wrong modenum ', 2) IF (gamma_only) CALL errore('phq_readin',& 'cannot start from pw.x data file using Gamma-point tricks',1) IF (lda_plus_u) THEN ! WRITE(stdout,'(/5x,a)') "Phonon calculation with DFPT+U; please cite" WRITE(stdout,'(5x,a)') "A. Floris et al., Phys. Rev. B 84, 161102(R) (2011)" WRITE(stdout,'(5x,a)') "A. Floris et al., Phys. Rev. B 101, 064305 (2020)" WRITE(stdout,'(5x,a)') "in publications or presentations arising from this work." ! IF (Hubbard_projectors.NE."atomic") CALL errore("phq_readin", & " The phonon code for this Hubbard projectors type is not implemented",1) IF (lda_plus_u_kind.NE.0) CALL errore("phq_readin", & " The phonon code for this lda_plus_u_kind is not implemented",1) IF (ANY(is_hubbard_back(:))) CALL errore ("phq_readin", & " Two (or more) Hubbard channels per atomic type is not implemented", 1) IF (elph) CALL errore("phq_readin", & " Electron-phonon with Hubbard U is not supported",1) IF (lraman) CALL errore("phq_readin", & " The phonon code with Raman and Hubbard U is not implemented",1) ! ENDIF ! checks IF (elph_ahc .AND. domag) CALL errore ('phq_readin', & 'elph_ahc and magnetism not implemented', 1) IF (elph_ahc .AND. okpaw) CALL errore ('phq_readin', & 'elph_ahc and PAW not tested.', 1) IF (elph_ahc .AND. okvan) CALL errore ('phq_readin', & 'elph_ahc and PAW not tested.', 1) IF (elph_ahc .AND. lda_plus_u) CALL errore ('phq_readin', & 'elph_ahc and lda_plus_u not tested.', 1) IF (ts_vdw) CALL errore('phq_readin',& 'The phonon code with TS-VdW is not yet available',1) !IF (ldftd3) CALL errore('phq_readin',& ! 'The phonon code with Grimme''s DFT-D3 is not yet available',1) IF ( xclib_dft_is('meta') ) CALL errore('phq_readin',& 'The phonon code with meta-GGA functionals is not yet available',1) IF ( xclib_dft_is('hybrid') ) CALL errore('phq_readin',& 'The phonon code with hybrid functionals is not yet available',1) IF (magnetic_sym) THEN WRITE(stdout,'(/5x,a)') "Phonon calculation in the non-collinear magnetic case;" WRITE(stdout,'(5x,a)') "please cite A. Urru and A. Dal Corso, Phys. Rev. B 100," WRITE(stdout,'(5x,a)') "045115 (2019) for the theoretical background." IF (okpaw) CALL errore('phq_readin',& 'The phonon code with paw and magnetism is not yet available',1) ENDIF IF (okpaw.and.(lraman.or.elop)) CALL errore('phq_readin',& 'The phonon code with paw and raman or elop is not yet available',1) IF (okvan.and.(lraman.or.elop)) CALL errore('phq_readin',& 'The phonon code with US-PP and raman or elop not yet available',1) IF (noncolin.and.(lraman.or.elop)) CALL errore('phq_readin', & 'lraman, elop, and noncolin not programmed',1) IF ( domag .and. lspinorb .and. elph ) CALL errore('phq_readin', & 'el-ph coefficient calculation disabled in magnetic spinorbit case',1) IF (lmovecell) CALL errore('phq_readin', & 'The phonon code is not working after vc-relax',1) IF(elph_mat) THEN if(fildvscf.eq.' ') call errore('phq_readin',& 'el-ph with wannier requires fildvscf',1) IF(npool.ne.1) call errore('phq_readin',& 'el-ph with wannier : pools not implemented',1) IF ((elph_nbnd_min .eq. 1).and.(elph_nbnd_max .eq. 0))& call errore('phq_readin',& 'el-ph with wannier : specify bands range with elph_nbnd_min,elph_nbnd_max',1) END IF IF (elph .AND. okpaw) CALL errore('phq_readin',& 'Electron-phonon calculations with PAW not tested',1) IF(elph.and.nimage>1) call errore('phq_readin',& 'el-ph with images not implemented',1) IF (elph.OR.fildvscf /= ' ') lqdir=.TRUE. IF(dvscf_star%open.and.nimage>1) CALL errore('phq_readin',& 'dvscf_star with image parallelization is not yet available',1) IF(drho_star%open.and.nimage>1) CALL errore('phq_readin',& 'drho_star with image parallelization is not yet available',1) IF (lda_plus_u .AND. read_dns_bare .AND. ldisp) lqdir=.TRUE. IF (.NOT.ldisp) lqdir=.FALSE. IF (i_cons /= 0) & CALL errore('phq_readin',& 'The phonon code with constrained magnetization is not yet available',1) IF (two_fermi_energies .AND. (ltetra .OR. lgauss)) & CALL errore('phq_readin',& 'The phonon code with two fermi energies is not available for metals',1) IF (tqr) CALL errore('phq_readin',& 'The phonon code with Q in real space not available',1) IF (start_irr < 0 ) CALL errore('phq_readin', 'wrong start_irr',1) ! IF (start_q <= 0 ) CALL errore('phq_readin', 'wrong start_q',1) ! ! the dynamical matrix is written in xml format if fildyn ends in ! .xml or in the noncollinear case. ! xmldyn=has_xml(fildyn) !IF (noncolin) xmldyn=.TRUE. ! ! If a band structure calculation needs to be done do not open a file ! for k point ! restart = recover ! ! set masses to values read from input, if available; ! leave values read from file otherwise ! DO it = 1, ntyp IF (amass_input(it) < 0.0_DP) amass_input(it)= & atom_weight(atomic_number(TRIM(atm(it)))) IF (amass_input(it) > 0.D0) amass(it) = amass_input(it) IF (amass(it) <= 0.D0) CALL errore ('phq_readin', 'Wrong masses', it) ENDDO DEALLOCATE (amass_input) lgamma_gamma=.FALSE. IF (.NOT.ldisp) THEN IF (nkstot==1.OR.(nkstot==2.AND.nspin==2)) THEN lgamma_gamma=(lgamma.AND.(ABS(xk(1,1))<1.D-12) & .AND.(ABS(xk(2,1))<1.D-12) & .AND.(ABS(xk(3,1))<1.D-12) ) ENDIF IF (nogg) lgamma_gamma=.FALSE. IF ((nat_todo /= 0) .and. lgamma_gamma) CALL errore( & 'phq_readin', 'gamma_gamma tricks with nat_todo & & not available. Use nogg=.true.', 1) IF (lda_plus_u .AND. lgamma_gamma) THEN WRITE(stdout,'(5x,a)') "DFPT+U does not support k=gamma and q=gamma tricks: disabling them..." lgamma_gamma=.FALSE. ENDIF ! IF (nimage > 1 .AND. lgamma_gamma) CALL errore( & 'phq_readin','gamma_gamma tricks with images not implemented',1) IF (lgamma) THEN nksq = nks ELSE nksq = nks / 2 ENDIF ENDIF IF (lgamma_gamma.AND.ldiag) CALL errore('phq_readin','incompatible flags',1) IF (lgamma_gamma.AND.elph ) CALL errore('phq_readin','lgamma_gamma and electron-phonon are incompatible',1) ! IF (tfixed_occ) & CALL errore('phq_readin','phonon with arbitrary occupations not tested',1) ! !YAMBO > IF (elph .AND. .NOT.(lgauss .OR. ltetra) & .AND. .NOT. (elph_yambo .OR. elph_ahc .OR. elph_print).and..not.elph_mat) & CALL errore ('phq_readin', 'Electron-phonon only for metals', 1) !YAMBO < IF (elph .AND. fildvscf.EQ.' ' .AND. .NOT. ldvscf_interpolate) & CALL errore ('phq_readin', 'El-ph needs a DeltaVscf file', 1) ! There might be other variables in the input file which describe ! partial computation of the dynamical matrix. Read them here ! CALL allocate_part ( nat ) ! IF ( nat_todo < 0 .OR. nat_todo > nat ) & CALL errore ('phq_readin', 'nat_todo is wrong', 1) IF (nat_todo.NE.0) THEN IF (meta_ionode) READ (qestdin, *, iostat = ios) (atomo (na), na=1,nat_todo) CALL mp_bcast(ios, meta_ionode_id, world_comm ) CALL errore ('phq_readin', 'reading atoms', ABS (ios) ) CALL mp_bcast(atomo, meta_ionode_id, world_comm ) ENDIF nat_todo_input=nat_todo ! ! end of reading, close unit qestdin, remove temporary input file if existing ! FIXME: closing input file here breaks alpha2F.x that reads what follows ! IF (meta_ionode) ios = close_input_file () IF (epsil.AND.(lgauss .OR. ltetra)) & CALL errore ('phq_readin', 'no elec. field with metals', 1) IF (modenum > 0) THEN IF ( ldisp ) & CALL errore('phq_readin','Dispersion calculation and & & single mode calculation not possibile !',1) nat_todo = 0 ENDIF ! ! Tetrahedron method for DFPT and El-Ph ! IF(ltetra .AND. tetra_type == 0) & & CALL errore ('phq_readin', 'DFPT with the Blochl correction is not implemented', 1) ! IF(.NOT. ltetra .AND. elph_tetra /= 0) & & CALL errore ('phq_readin', '"electron_phonon" and "occupation" are inconsistent.', 1) ! IF (modenum > 0 .OR. lraman ) lgamma_gamma=.FALSE. IF (.NOT.lgamma_gamma) asr=.FALSE. ! IF ((ldisp.AND..NOT.qplot) .AND. & (nq1 .LE. 0 .OR. nq2 .LE. 0 .OR. nq3 .LE. 0)) & CALL errore('phq_readin','nq1, nq2, and nq3 must be greater than 0',1) ! ! If dftd3_hess is not specified, use a default name set from prefix ! IF ( dftd3_hess == ' ' ) dftd3_hess = trim(prefix)//'.hess' dftd3_hess = TRIM(tmp_dir)//TRIM(dftd3_hess) CALL save_ph_input_variables() ! RETURN ! END SUBROUTINE phq_readin