quantum-espresso/TDDFPT/src/lr_readin.f90

565 lines
18 KiB
Fortran
Executable File

!
! Copyright (C) 2001-2018 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 lr_readin
!-----------------------------------------------------------------------
!
! This routine reads the control variables from standard input (unit 5).
! A second routine read_file reads the variables saved to file
! by PW scf (and PW nscf for EELS).
!
USE lr_variables
USE lr_dav_variables
USE kinds, ONLY : DP
USE io_files, ONLY : tmp_dir, prefix, wfc_dir, create_directory
USE lsda_mod, ONLY : current_spin, nspin, isk, lsda
USE control_flags, ONLY : use_para_diag, tqr, gamma_only,&
& do_makov_payne
USE scf, ONLY : vltot, v, vrs, vnew, &
& destroy_scf_type, rho
USE fft_base, ONLY : dfftp, dffts
USE gvecs, ONLY : doublegrid
USE wvfct, ONLY : nbnd, et, wg, current_k
USE lsda_mod, ONLY : isk
USE ener, ONLY : ef
USE io_global, ONLY : ionode, ionode_id, stdout
USE klist, ONLY : nks, wk, nelec, lgauss, ltetra
USE fixed_occ, ONLY : tfixed_occ
USE input_parameters, ONLY : degauss, nosym, wfcdir, outdir
USE check_stop, ONLY : max_seconds
USE realus, ONLY : real_space, init_realspace_vars, generate_qpointlist, &
betapointlist
USE funct, ONLY : dft_is_meta
USE charg_resp, ONLY : w_T_prefix, omeg, w_T_npol, epsil
USE mp, ONLY : mp_bcast
USE mp_global, ONLY : my_pool_id, intra_image_comm, &
& intra_bgrp_comm, nproc_image, &
& nproc_pool, nproc_pool_file, &
& nproc_image_file, nproc_bgrp, &
& nproc_bgrp_file, my_image_id
USE DFUNCT, ONLY : newd
USE vlocal, ONLY : strf
USE martyna_tuckerman, ONLY : do_comp_mt
USE esm, ONLY : do_comp_esm
USE qpoint, ONLY : xq
USE io_rho_xml, ONLY : write_scf
USE noncollin_module, ONLY : noncolin
USE mp_bands, ONLY : ntask_groups
USE constants, ONLY : eps4
USE control_lr, ONLY : lrpa
IMPLICIT NONE
!
CHARACTER(LEN=256), EXTERNAL :: trimcheck
!
CHARACTER(LEN=256) :: beta_gamma_z_prefix
! Fine control of beta_gamma_z file
CHARACTER(LEN=80) :: disk_io
! Specify the amount of I/O activities
CHARACTER(LEN=6) :: int_to_char
INTEGER :: ios, iunout, ierr, ipol
LOGICAL, EXTERNAL :: check_para_diag
!
NAMELIST / lr_input / restart, restart_step ,lr_verbosity, prefix, outdir, &
& test_case_no, wfcdir, disk_io, max_seconds
NAMELIST / lr_control / itermax, ipol, ltammd, lrpa, &
& charge_response, no_hxc, n_ipol, project, &
& scissor, pseudo_hermitian, d0psi_rs, lshift_d0psi, &
& q1, q2, q3, approximation
NAMELIST / lr_post / omeg, beta_gamma_z_prefix, w_T_npol, plot_type, epsil, itermax_int,sum_rule
namelist / lr_dav / num_eign, num_init, num_basis_max, residue_conv_thr, precondition, &
& dav_debug, reference,single_pole, sort_contr, diag_of_h, close_pre, &
& broadening,print_spectrum,start,finish,step,if_check_orth, if_random_init, &
& if_check_her,p_nbnd_occ,p_nbnd_virt,poor_of_ram,poor_of_ram2,max_iter, &
& conv_assistant,if_dft_spectrum,no_hxc,d0psi_rs,lshift_d0psi, &
& lplot_drho, vccouple_shift, ltammd
!
#if defined(__MPI)
IF (ionode) THEN
#endif
!
! Checking for the path to the output directory.
!
CALL get_environment_variable( 'ESPRESSO_TMPDIR', outdir )
IF ( trim( outdir ) == ' ' ) outdir = './'
!
! Set default values for variables in namelist.
!
itermax = 500
restart = .FALSE.
restart_step = itermax+1
lr_verbosity = 1
prefix = 'pwscf'
disk_io = 'default'
ltammd = .FALSE.
d0psi_rs=.false.
lshift_d0psi = .true.
pseudo_hermitian=.true.
ipol = 1
n_ipol = 1
no_hxc = .FALSE.
lrpa = .false.
charge_response = 0
sum_rule = -99
test_case_no = 0
beta_gamma_z_prefix = 'undefined'
omeg= 0.0_DP
epsil = 0.0_DP
w_T_npol = 1
plot_type = 1
project = .FALSE.
max_seconds = 1.0E+7_DP
scissor = 0.d0
!
! For EELS
!
q1 = 1.0d0
q2 = 1.0d0
q3 = 1.0d0
approximation = 'TDDFT'
!
! For lr_dav (Davidson program)
!
num_eign=1
num_init=2
num_basis_max=20
broadening=0.005d0
residue_conv_thr=1.0E-4
close_pre=1.0E-5
turn2planb=1.0E-3
precondition=.true.
dav_debug=.false.
reference=0.0d0
vccouple_shift=0.0d0
single_pole=.false.
sort_contr=.true.
print_spectrum=.true.
start=0.0d0
finish=1.0d0
step=0.001d0
if_check_orth=.false.
if_check_her=.false.
if_random_init=.false.
p_nbnd_occ=10
p_nbnd_virt=10
poor_of_ram=.false.
poor_of_ram2=.false.
max_iter=100
conv_assistant=.false.
if_dft_spectrum=.false.
lplot_drho=.false.
!
! Reading possible plugin arguments
!
CALL plugin_arguments()
!
! Reading the namelist lr_input
!
CALL input_from_file( )
!
READ (5, lr_input, err = 200, iostat = ios)
200 CALL errore ('lr_readin', 'reading lr_input namelist', ABS (ios) )
!
! Reading the namelist lr_dav or lr_control
!
IF (davidson) THEN
READ (5, lr_dav, err = 201, iostat = ios)
201 CALL errore ('lr_readin', 'reading lr_dav namelist', ABS (ios) )
ELSE
READ (5, lr_control, err = 202, iostat = ios)
202 CALL errore ('lr_readin', 'reading lr_control namelist', ABS (ios) )
ENDIF
!
! Reading the namelist lr_post (only for optical case)
!
IF (charge_response == 1 .AND. .NOT.eels) THEN
!
READ (5, lr_post, err = 203, iostat = ios)
203 CALL errore ('lr_readin', 'reading lr_post namelist', ABS (ios) )
!
bgz_suffix = TRIM ( "-stage2.beta_gamma_z." )
WRITE(stdout,'(/5x,"Prefix of current run is appended by -stage2")')
!
IF ( beta_gamma_z_prefix == 'undefined' ) THEN
beta_gamma_z_prefix = TRIM(prefix)
ENDIF
!
ELSE
bgz_suffix = TRIM ( ".beta_gamma_z." )
ENDIF
!
! The status of the real space flags should be read manually
!
! Set-up all the dir and suffix variables.
!
outdir = trimcheck(outdir)
tmp_dir = outdir
!
!IF ( .NOT. TRIM( wfcdir ) == 'undefined' ) THEN
! wfc_dir = trimcheck ( wfcdir )
!ENDIF
!
IF (.NOT.eels) THEN
w_T_prefix = TRIM( tmp_dir ) // &
& TRIM( beta_gamma_z_prefix ) // ".beta_gamma_z."
ENDIF
!
ierr = 0
!
IF (eels) THEN
!
IF (n_ipol /= 1) THEN
WRITE(stdout,'(5X,"n_ipol /= 1 is not allowed for EELS.")')
WRITE(stdout,'(5X,"Setting n_ipol = 1 ...")')
n_ipol = 1
ENDIF
IF (ipol /= 1) THEN
WRITE(stdout,'(5X,"ipol /= 1 is not allowed for EELS.")')
WRITE(stdout,'(5X,"Setting ipol = 1 ...")')
ipol = 1
ENDIF
LR_polarization = 1
!
ELSE
!
! Optics: set up polarization direction(s)
!
IF ( ipol==4 ) THEN
!
n_ipol = 3
LR_polarization = 1
!
ELSE
!
LR_polarization = ipol
!
ENDIF
!
ENDIF
!
IF (itermax_int < itermax) itermax_int = itermax
!
! Limited disk_io support: currently only one setting is supported
!
SELECT CASE( TRIM( disk_io ) )
!
CASE ( 'reduced' )
!
lr_io_level = -1
restart = .FALSE.
!
CASE DEFAULT
!
lr_io_level = 1
!
END SELECT
!
IF (eels) THEN
!
! Level of approximation in turboEELS.
!
SELECT CASE( trim(approximation) )
!
CASE ( 'TDDFT' )
!
no_hxc = .FALSE.
lrpa = .FALSE.
!
CASE ( 'IPA' )
!
no_hxc = .TRUE.
lrpa = .TRUE.
!
CASE ( 'RPA_with_CLFE' )
!
no_hxc = .FALSE.
lrpa = .TRUE.
!
CASE DEFAULT
!
CALL errore( 'lr_readin', 'Approximation ' // &
& trim( approximation ) // ' not implemented', 1 )
!
END SELECT
!
! We do this trick because xq is used in LR_Modules/dv_of_drho.f90
! in the Hartree term ~1/|xq+k|^2
!
xq(1) = q1
xq(2) = q2
xq(3) = q3
!
IF ( (q1.lt.eps4) .AND. (q2.lt.eps4) .AND. (q3.lt.eps4) ) &
CALL errore( 'lr_readin', 'The transferred momentum |q| is too small, the limit is not implemented.', 1 )
!
ENDIF
!
#if defined(__MPI)
ENDIF
!
CALL bcast_lr_input
#endif
!
! Required for restart runs as this never gets initialized.
!
current_k = 1
!
outdir = TRIM( tmp_dir ) // TRIM( prefix ) // '.save'
!
! EELS: Create a temporary directory for nscf files, and for
! writing of the turboEELS restart files.
!
IF (eels) THEN
tmp_dir_lr = TRIM (tmp_dir) // 'tmp_eels/'
CALL create_directory(tmp_dir_lr)
ENDIF
!
! Now PWSCF XML file will be read, and various initialisations will be done.
! Allocate space for PW scf variables, read and check them.
! Optical case: the variables igk_k and ngk are set up through this path:
! read_file -> init_igk.
! EELS: the variables igk_k and ngk will be re-set up later (because there
! will be not only poins k but also points k+q) through the path:
! lr_run_nscf -> init_run -> hinit0 -> init_igk
!
CALL read_file()
!
IF (.NOT.eels .AND. (tqr .OR. real_space)) &
WRITE(stdout,'(/5x,"Status of real space flags: TQR=", L5 , &
& " REAL_SPACE=", L5)') tqr, real_space
!
! Set wfc_dir - this is done here because read_file sets wfc_dir = tmp_dir
! FIXME:,if wfcdir is not present in input, wfc_dir is set to "undefined"
! instead of tmp_dir, because of the logic used in the rest of TDDFPT
!
wfc_dir = trimcheck ( wfcdir )
!
IF (eels) THEN
!
! Specify the temporary directory.
!
tmp_dir = tmp_dir_lr
!
! Copy the scf-charge-density to the tmp_dir (PH/check_initial_status.f90).
! Needed for the nscf calculation.
!
IF (.NOT.restart) CALL write_scf( rho, nspin )
!
ENDIF
!
! Make sure all the features used in the PWscf calculation
! are actually supported by TDDFPT.
!
CALL input_sanity()
!
! Compute and add plugin contributions to SCF potential
!
CALL plugin_read_input("TD")
!
CALL plugin_initbase()
!
CALL plugin_init_ions()
!
CALL plugin_init_cell()
!
CALL plugin_init_potential( v%of_r(:,1) )
!
CALL plugin_scf_potential( rho, .FALSE., -1.D0, v%of_r(:,1) )
!
CALL plugin_clean( 'TD', .FALSE. )
!
! Deallocate some variables created with read_file but not used in TDDFPT
!
DEALLOCATE( strf )
CALL destroy_scf_type(vnew)
!
! Re-initialize all needed quantities from the scf run
! I. Timrov: this was already done in read_file.
current_spin = 1
!
! I. Timrov: The routine init_us_1 was already called in read_file above.
CALL init_us_1 ( )
!
! I. Timrov: The routine newd was already called in read_file above.
!
CALL newd() !OBM: this is for the ground-state charge density
!
IF (tqr .AND. .NOT.eels) CALL generate_qpointlist()
!
IF ( real_space .AND. .NOT.eels) THEN
!
WRITE(stdout,'(/5x,"Real space implementation V.1 D190908",1x)')
! OBM - correct parellism issues
CALL init_realspace_vars()
CALL betapointlist()
WRITE(stdout,'(5X,"Real space initialisation completed")')
ENDIF
!
! Now put the potential calculated in read_file into the correct place
! and deallocate the redundant associated variables.
! Set the total local potential vrs on the smooth mesh
! adding the scf (Hartree + XC) part and the sum of
! all the local pseudopotential contributions.
! vrs = vltot + v%of_r
!
CALL set_vrs ( vrs, vltot, v%of_r, 0, 0, dfftp%nnr, nspin, doublegrid )
!
DEALLOCATE( vltot )
CALL destroy_scf_type(v)
!
! Recalculate the weights of the Kohn-Sham orbitals.
! (Should this not be a call to weights() to make this
! less insulator specific?)
!
IF (.NOT.eels) CALL iweights( nks, wk, nbnd, nelec, et, ef, wg, 0, isk)
!
IF ( charge_response == 2 .AND. .NOT.eels) CALL lr_set_boxes_density()
!
! Scalapack related stuff.
!
use_para_diag = check_para_diag( nbnd )
!
RETURN
!
CONTAINS
!
SUBROUTINE input_sanity()
!--------------------------------------------------------------------------
!
! This subroutine aims to gather all of the input sanity checks
! (features enabled in PWscf which are unsupported in TDDFPT).
! Written by Simone Binnie 2011
! Modified by Iurii Timrov 2015
!
USE paw_variables, ONLY : okpaw
USE uspp, ONLY : okvan
USE funct, ONLY : dft_is_hybrid
USE ldaU, ONLY : lda_plus_u
IMPLICIT NONE
!
! Charge response mode 1 is the "do Lanczos chains twice, conserve memory" scheme.
!
IF (.NOT.eels) THEN
!
IF (charge_response == 1 .AND. ( omeg == 0.D0 .AND. sum_rule == -99 ) ) &
& CALL errore ('lr_readin', &
& 'omeg must be defined for charge response mode 1', 1 )
!
IF ( project .AND. charge_response /= 1) &
& CALL errore ('lr_readin', &
& 'projection is possible only in charge response mode 1', 1 )
!
ENDIF
!
! Meta-DFT currently not supported by TDDFPT
!
IF (dft_is_meta()) CALL errore( 'lr_readin', 'Meta DFT is not implemented yet', 1 )
!
! Hubbard U is not supported
!
IF (lda_plus_u) CALL errore('lr_readin', 'TDDFPT with Hubbard U is not implemented',1)
!
! Tetrahedron method and fixed occupations are not implemented.
!
IF (ltetra) CALL errore( 'lr_readin', 'ltetra is not implemented', 1 )
IF (tfixed_occ) CALL errore( 'lr_readin', 'tfixed_occ is not implemented', 1 )
!
! Some limitations of turboTDDFT (and not of turboEELS).
!
IF (.NOT. eels) THEN
!
! Non-insulating systems currently not supported by turboTDDFPT, but
! supported by turboEELS.
!
IF (lgauss) CALL errore( 'lr_readin', 'turboTDDFT is not extended to metals', 1 )
!
! Symmetry is not supported.
!
IF (.NOT.nosym ) CALL errore( 'lr_readin', 'Linear response calculation' // &
& 'is not implemented with symmetry', 1 )
!
! K-points are implemented but still unsupported (use at your own risk!)
!
IF (.NOT. gamma_only ) CALL errore('lr_readin', 'k-point algorithm is not tested yet',1)
!
ENDIF
!
! No taskgroups and EXX.
!
IF (dffts%has_task_groups .AND. dft_is_hybrid()) &
& CALL errore( 'lr_readin', ' Linear response calculation ' // &
& 'not implemented for EXX+Task groups', 1 )
!
! Experimental task groups warning.
!
IF (dffts%has_task_groups) &
& CALL infomsg( 'lr_readin','Usage of task &
& groups with TDDFPT is still experimental. Use at your own risk.' )
! & CALL errore( 'lr_readin', ' Linear response calculation ' // &
! & 'not implemented for task groups', 1 )
!
! No PAW support.
!
IF (okpaw) &
& CALL errore( 'lr_readin', ' Linear response calculation ' // &
& 'not implemented for PAW', 1 )
!
! No USPP+EXX support.
!
IF (okvan .AND. dft_is_hybrid()) &
& CALL errore( 'lr_readin', ' Linear response calculation ' // &
& 'not implemented for EXX+Ultrasoft', 1 )
!
! Spin-polarised case is not implemented, but partially accounted in
! some routines.
!
IF (lsda) CALL errore( 'lr_readin', 'LSDA is not implemented', 1 )
!
IF (real_space) THEN
IF (eels) THEN
CALL errore( 'lr_readin', 'Option real_space=.true. is not implemented', 1 )
ELSE
CALL errore( 'lr_readin', 'Option real_space=.true. is not tested', 1 )
ENDIF
ENDIF
!
! EELS-related restrictions
!
IF (eels) THEN
!
IF (okvan .AND. noncolin) CALL errore( 'lr_readin', 'Ultrasoft PP + noncolin is not fully implemented', 1 )
IF (gamma_only) CALL errore( 'lr_readin', 'gamma_only is not supported', 1 )
!
! Tamm-Dancoff approximation is not recommended to be used with EELS, and
! thus it was not implemented.
!
IF (ltammd) CALL errore( 'lr_readin', 'EELS + Tamm-Dancoff approximation is not supported', 1 )
!
IF (project) CALL errore( 'lr_readin', 'project is not allowed', 1 )
IF (tqr) CALL errore( 'lr_readin', 'tqr is not supported', 1 )
IF (charge_response /= 0) CALL errore( 'lr_readin', 'charge_response /= 0 is not allowed', 1 )
IF (dft_is_hybrid()) CALL errore( 'lr_readin', 'EXX is not supported', 1 )
IF (do_comp_mt) CALL errore( 'lr_readin', 'Martyna-Tuckerman PBC is not supported.', 1 )
IF (d0psi_rs) CALL errore( 'lr_readin', 'd0psi_rs is not allowed', 1 )
!
! Note, all variables of the turboDavidson code cannot be used by turboEELS.
!
! EELS + plugins is not supported.
!
CALL plugin_check('lr_readin')
!
ENDIF
!
RETURN
!
END SUBROUTINE input_sanity
!
END SUBROUTINE lr_readin