mirror of https://gitlab.com/QEF/q-e.git
464 lines
15 KiB
Fortran
464 lines
15 KiB
Fortran
!
|
|
! 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 phq_setup
|
|
!-----------------------------------------------------------------------
|
|
!! This subroutine prepares several variables which are needed in the
|
|
!! \(\texttt{phonon}\) program:
|
|
!! 1) computes the total local potential (external+scf) on the smooth
|
|
!! grid to be used in \(\texttt{h_psi}\) and similia;
|
|
!! 2) computes the local magnetization (if necessary);
|
|
!! 3) computes dmuxc (with GC if needed);
|
|
!! 4) set the inverse of every matrix invs;
|
|
!! 5) for metals sets the occupied bands;
|
|
!! 6) computes \(\text{alpha_pv}\);
|
|
!! 7) computes the variables needed to pass to the pattern representation:
|
|
!! \(\text{u}\): the patterns;
|
|
!! \(\text{t}\): the matrices of the small group of q on the pattern basis;
|
|
!! \(\text{tmq}\): the matrix of the symmetry which sends \(q\rightarrow -q+G\);
|
|
!! \(\text{gi}\): the G associated to each symmetry operation;
|
|
!! \(\text{gimq}\): the G of the \(q\rightarrow -q+G\) symmetry;
|
|
!! \(\text{nsymq}\): the order of the small group of \(q\);
|
|
!! \(\text{irotmq}\): the index of the \(q\rightarrow -q+G\) symmetry;
|
|
!! \(\text{nirr}\): the number of irreducible representation;
|
|
!! \(\text{npert}\): the dimension of each irreducible representation;
|
|
!! \(\text{nmodes}\): the number of modes;
|
|
!! \(\text{minus_q}\): true if there is a symmetry sending \(q\rightarrow -q+G\);
|
|
!! 8) for testing purposes it sets \(\text{ubar}\);
|
|
!! 9) set the variables needed to deal with \(\text{nlcc}\);
|
|
!! 10) set the variables needed for the partial computation
|
|
! of the dynamical matrix.
|
|
!
|
|
!! IMPORTANT NOTE ABOUT SYMMETRIES:
|
|
!
|
|
!! * \(\text{nrot}\) is the number of sym.ops. of the Bravais lattice,
|
|
!! read from data file, only used in \(\texttt{set_default_pw}\);
|
|
!! * \(\text{nsym}\) is the number of sym.ops. of the crystal symmetry group,
|
|
!! read from data file, should never be changed;
|
|
!! * \(\text{nsymq}\) is the number of sym.ops. of the small group of q,
|
|
!! it is calculated in \(\texttt{set_defaults_pw}\) for each q;
|
|
!! * The matrices "s" of sym.ops are ordered as follows:
|
|
!! first the nsymq sym.ops. of the small group of q
|
|
!! (the ordering is done in subroutine \(\texttt{copy_sym} in
|
|
!! \(\texttt{set_defaults_pw}\)),
|
|
!! followed by the remaining nsym-nsymq sym.ops. of the crystal group,
|
|
!! followed by the remaining nrot-nsym sym.ops. of the Bravais group.
|
|
!
|
|
!
|
|
USE kinds, ONLY : DP
|
|
USE ions_base, ONLY : tau, nat, ntyp => nsp, ityp
|
|
USE cell_base, ONLY : at, bg
|
|
USE io_global, ONLY : ionode, stdout
|
|
USE io_files, ONLY : tmp_dir
|
|
USE klist, ONLY : xk, nks, nkstot
|
|
USE lsda_mod, ONLY : nspin, starting_magnetization
|
|
USE scf, ONLY : v, vrs, vltot, kedtau, rho
|
|
USE dfunct, ONLY : newd
|
|
USE fft_base, ONLY : dfftp
|
|
USE gvect, ONLY : ngm
|
|
USE gvecs, ONLY : doublegrid
|
|
USE symm_base, ONLY : nrot, nsym, s, irt, t_rev, time_reversal, &
|
|
sr, invs, inverse_s, d1, d2, d3, check_grid_sym
|
|
USE uspp_param, ONLY : upf
|
|
USE uspp, ONLY : nlcc_any, deeq_nc, okvan
|
|
USE noncollin_module, ONLY : noncolin, domag, m_loc, angle1, angle2, ux
|
|
USE nlcc_ph, ONLY : drc
|
|
USE control_ph, ONLY : search_sym, start_irr, &
|
|
last_irr, all_done, trans, epsil, recover, &
|
|
done_epsil, zeu, done_zeu, current_iq, u_from_file
|
|
USE el_phon, ONLY : elph, comp_elph, done_elph, elph_nbnd_min, elph_nbnd_max
|
|
USE output, ONLY : fildrho
|
|
USE modes, ONLY : u, npertx, npert, nirr, t, tmq, nmodes, num_rap_mode
|
|
USE dynmat, ONLY : dyn, dyn_rec, dyn00
|
|
USE efield_mod, ONLY : epsilon, zstareu
|
|
USE partial, ONLY : comp_irr, atomo, nat_todo, all_comp, &
|
|
done_irr
|
|
USE gamma_gamma, ONLY : has_equivalent, asr, nasr, n_diff_sites, &
|
|
equiv_atoms, n_equiv_atoms, with_symmetry
|
|
USE ph_restart, ONLY : ph_writefile, ph_readfile
|
|
USE control_flags, ONLY : modenum, noinv
|
|
USE grid_irr_iq, ONLY : comp_irr_iq
|
|
USE xc_lib, ONLY : xclib_dft_is
|
|
USE ramanm, ONLY : lraman, elop, ramtns, eloptns, done_lraman, &
|
|
done_elop
|
|
USE mp_pools, ONLY : inter_pool_comm, npool
|
|
USE acfdtest, ONLY : acfdt_is_active, acfdt_num_der
|
|
USE elph_tetra_mod, ONLY : elph_tetra
|
|
USE wvfct, ONLY : nbnd, et
|
|
USE ener, ONLY : ef
|
|
USE mp, ONLY : mp_max, mp_min
|
|
USE lr_symm_base, ONLY : gi, gimq, irotmq, minus_q, invsymq, nsymq, rtau
|
|
USE qpoint, ONLY : xq, xk_col
|
|
USE lr_nc_mag, ONLY : deeq_nc_save
|
|
USE control_lr, ONLY : lgamma, lgamma_gamma, niter_ph, alpha_mix, flmixdpot, &
|
|
reduce_io, rec_code, rec_code_read, where_rec
|
|
USE ldaU, ONLY : lda_plus_u, Hubbard_U, Hubbard_J0
|
|
USE ldaU_lr, ONLY : effU
|
|
USE constants, ONLY : rytoev
|
|
USE dvscf_interpolate, ONLY : ldvscf_interpolate, dvscf_interpol_setup
|
|
USE ahc, ONLY : elph_ahc, elph_ahc_setup
|
|
|
|
USE el_phon, ONLY : elph_mat
|
|
implicit none
|
|
|
|
real(DP) :: sr_is(3,3,48)
|
|
|
|
integer :: isym, jsym, irot, ik, ibnd, ipol, nah, &
|
|
mu, nu, imode0, irr, ipert, na, it, nt, nsym_is, last_irr_eff
|
|
! counters
|
|
|
|
logical :: sym (48), magnetic_sym
|
|
LOGICAL :: symmorphic_or_nzb
|
|
! the symmetry operations
|
|
integer, allocatable :: ifat(:)
|
|
integer :: ierr
|
|
|
|
call start_clock ('phq_setup')
|
|
! 0) A few checks
|
|
!
|
|
IF (xclib_dft_is('gradient').and.(lraman.or.elop)) call errore('phq_setup', &
|
|
'third order derivatives not implemented with GGA', 1)
|
|
|
|
IF (nsymq==0) CALL errore('phq_setup', &
|
|
'The small group of q is no more calculated in phq_setup',1)
|
|
!
|
|
! read the displacement patterns
|
|
!
|
|
IF (u_from_file) THEN
|
|
CALL ph_readfile('data_u',current_iq,0,ierr)
|
|
IF (ierr /= 0) CALL errore('phq_setup', 'problem with modes file',1)
|
|
ENDIF
|
|
!
|
|
! 1) Computes the total local potential (external+scf) on the smooth grid
|
|
!
|
|
!************************ ACFDT TEST ***********************
|
|
! write(*,*) " acfdt_is_active ",acfdt_is_active, " acfdt_num_der ", acfdt_num_der
|
|
IF (acfdt_is_active) THEN
|
|
! discard set_vrs for numerical derivatives
|
|
if (.not.acfdt_num_der) then
|
|
call set_vrs (vrs, vltot, v%of_r, kedtau, v%kin_r, dfftp%nnr, nspin, doublegrid)
|
|
end if
|
|
ELSE
|
|
call set_vrs (vrs, vltot, v%of_r, kedtau, v%kin_r, dfftp%nnr, nspin, doublegrid)
|
|
ENDIF
|
|
!************************END OF ACFDT TEST ******************
|
|
!
|
|
! Set non linear core correction stuff
|
|
!
|
|
nlcc_any = ANY ( upf(1:ntyp)%nlcc )
|
|
if (nlcc_any) allocate (drc( ngm, ntyp))
|
|
!
|
|
! 2) If necessary calculate the local magnetization. This information is
|
|
! needed in find_sym
|
|
!
|
|
IF (.not.ALLOCATED(m_loc)) ALLOCATE( m_loc( 3, nat ) )
|
|
IF (noncolin.and.domag) THEN
|
|
DO na = 1, nat
|
|
!
|
|
m_loc(1,na) = starting_magnetization(ityp(na)) * &
|
|
SIN( angle1(ityp(na)) ) * COS( angle2(ityp(na)) )
|
|
m_loc(2,na) = starting_magnetization(ityp(na)) * &
|
|
SIN( angle1(ityp(na)) ) * SIN( angle2(ityp(na)) )
|
|
m_loc(3,na) = starting_magnetization(ityp(na)) * &
|
|
COS( angle1(ityp(na)) )
|
|
END DO
|
|
ux=0.0_DP
|
|
if (xclib_dft_is('gradient')) call compute_ux(m_loc,ux,nat)
|
|
IF (okvan) THEN
|
|
!
|
|
! Change the sign of the magnetic field in the screened US coefficients
|
|
! and save also the coefficients computed with -B_xc.
|
|
!
|
|
deeq_nc_save(:,:,:,:,1)=deeq_nc(:,:,:,:)
|
|
v%of_r(:,2:4)=-v%of_r(:,2:4)
|
|
CALL newd()
|
|
v%of_r(:,2:4)=-v%of_r(:,2:4)
|
|
deeq_nc_save(:,:,:,:,2)=deeq_nc(:,:,:,:)
|
|
deeq_nc(:,:,:,:)=deeq_nc_save(:,:,:,:,1)
|
|
!$acc update device(deeq_nc)
|
|
ENDIF
|
|
ENDIF
|
|
!
|
|
! 3) Computes the derivative of the XC potential
|
|
!
|
|
call setup_dmuxc()
|
|
!
|
|
! Setup all gradient correction stuff
|
|
!
|
|
call setup_dgc()
|
|
!
|
|
! 4) Computes the number of occupied bands for each k point
|
|
!
|
|
if(.not.elph_mat)&
|
|
call setup_nbnd_occ()
|
|
!
|
|
! 5) Computes alpha_pv
|
|
!
|
|
if(.not.elph_mat)&
|
|
call setup_alpha_pv()
|
|
!
|
|
! 6) Set all symmetries and variables needed to use the pattern representation
|
|
!
|
|
call inverse_s()
|
|
IF ( .NOT. check_grid_sym (dfftp%nr1,dfftp%nr2,dfftp%nr3) ) &
|
|
CALL errore('phq_setup','FFT grid incompatible with symmetry',1)
|
|
!
|
|
magnetic_sym = noncolin .AND. domag
|
|
time_reversal = .NOT. noinv .AND. .NOT. magnetic_sym
|
|
|
|
nmodes = 3 * nat
|
|
|
|
! FIXME: xk_col is used only to write the output
|
|
IF (npool > 1) THEN
|
|
CALL poolcollect( 3, nks, xk, nkstot, xk_col )
|
|
ELSE
|
|
xk_col(:,1:nks) = xk(:,1:nks)
|
|
ENDIF
|
|
!
|
|
! The small group of q may be known. At a given q it is calculated
|
|
! by set_nscf, at gamma it coincides with the point group and we
|
|
! take nsymq=nsym
|
|
!
|
|
IF (lgamma.AND.modenum==0) THEN
|
|
nsymq=nsym
|
|
minus_q=.TRUE.
|
|
ENDIF
|
|
!
|
|
! If the code arrives here and nsymq is still 0 the small group of q has
|
|
! not been calculated by set_nscf because this is a recover run.
|
|
! We recalculate here the small group of q.
|
|
!
|
|
IF (nsymq==0) CALL set_small_group_of_q(nsymq, invsymq, minus_q)
|
|
IF ( .NOT. time_reversal ) minus_q = .FALSE.
|
|
!
|
|
!
|
|
IF (modenum > 0) THEN
|
|
search_sym=.FALSE.
|
|
minus_q = .FALSE.
|
|
ENDIF
|
|
!
|
|
! allocate and calculate rtau, the bravais lattice vector associated
|
|
! to a rotation
|
|
!
|
|
call sgam_lr (at, bg, nsym, s, irt, tau, rtau, nat)
|
|
!
|
|
! and calculate the vectors G associated to the symmetry Sq = q + G
|
|
! if minus_q is true calculate also irotmq and the G associated to Sq=-g+G
|
|
!
|
|
CALL set_giq (xq,s,nsymq,nsym,irotmq,minus_q,gi,gimq)
|
|
|
|
search_sym = search_sym .AND. symmorphic_or_nzb()
|
|
|
|
num_rap_mode=-1
|
|
IF (search_sym) CALL prepare_sym_analysis(nsymq,sr,t_rev,magnetic_sym)
|
|
|
|
IF (.NOT.u_from_file) THEN
|
|
CALL find_irrep()
|
|
CALL ph_writefile('data_u',current_iq,0,ierr)
|
|
ENDIF
|
|
CALL find_irrep_sym()
|
|
|
|
|
|
IF (lgamma_gamma) THEN
|
|
ALLOCATE(has_equivalent(nat))
|
|
ALLOCATE(with_symmetry(3*nat))
|
|
ALLOCATE(n_equiv_atoms(nat))
|
|
ALLOCATE(equiv_atoms(nat,nat))
|
|
CALL find_equiv_sites (nat,nsym,irt,has_equivalent,n_diff_sites, &
|
|
n_equiv_atoms,equiv_atoms)
|
|
|
|
IF (n_diff_sites .LE. 0 .OR. n_diff_sites .GT. nat) &
|
|
& CALL errore('phq_setup','problem with n_diff_sites',1)
|
|
!
|
|
! look if ASR can be exploited to reduce the number of calculations
|
|
! we need to locate an independent atom with no equivalent atoms
|
|
nasr=0
|
|
IF (asr.AND.n_diff_sites.GT.1) THEN
|
|
DO na = 1, n_diff_sites
|
|
IF (n_equiv_atoms(na).EQ.1 ) THEN
|
|
nasr = equiv_atoms(na, 1)
|
|
GO TO 1
|
|
END IF
|
|
END DO
|
|
1 CONTINUE
|
|
END IF
|
|
END IF
|
|
|
|
|
|
if (fildrho.ne.' '.and.ionode) call io_pattern (nat,fildrho,nirr,npert,u,xq,tmp_dir,+1)
|
|
|
|
last_irr_eff=last_irr
|
|
if (last_irr > nirr.or.last_irr<0) last_irr_eff=nirr
|
|
!
|
|
! set the alpha_mix parameter
|
|
!
|
|
do it = 2, niter_ph
|
|
if (alpha_mix (it) .eq.0.d0) alpha_mix (it) = alpha_mix (it - 1)
|
|
enddo
|
|
!
|
|
! Set flmixdpot
|
|
!
|
|
if (reduce_io) then
|
|
flmixdpot = ' '
|
|
else
|
|
flmixdpot = 'mixd'
|
|
endif
|
|
!
|
|
!
|
|
! 8) Set the ubar
|
|
!
|
|
! ubar removed on 16/02/2012, used only for debugging
|
|
!
|
|
! 9) set the variables needed for the partial computation:
|
|
! nat_todo, atomo, comp_irr
|
|
|
|
comp_irr(0:nirr) = comp_irr_iq(0:nirr, current_iq)
|
|
IF (elph) comp_elph(1:nirr) = comp_irr_iq(1:nirr, current_iq)
|
|
!
|
|
! The gamma_gamma case needs a different treatment
|
|
!
|
|
if (lgamma_gamma) then
|
|
with_symmetry=1
|
|
comp_irr = .FALSE.
|
|
comp_irr(0)=.TRUE.
|
|
do na=1,nat
|
|
if (has_equivalent(na)==0) then
|
|
do ipol=1,3
|
|
comp_irr(3*(na-1)+ipol)=.TRUE.
|
|
with_symmetry(3*(na-1)+ipol)=0
|
|
enddo
|
|
endif
|
|
enddo
|
|
if (nasr>0) then
|
|
do ipol=1,3
|
|
comp_irr(3*(nasr-1)+ipol)=.FALSE.
|
|
with_symmetry(3*(nasr-1)+ipol)=0
|
|
enddo
|
|
endif
|
|
IF (start_irr <= last_irr_eff) THEN
|
|
DO irr=1,start_irr-1
|
|
comp_irr(irr) = .FALSE.
|
|
ENDDO
|
|
DO irr=last_irr_eff+1,3*nat
|
|
comp_irr(irr) = .FALSE.
|
|
ENDDO
|
|
ENDIF
|
|
endif
|
|
!
|
|
! Compute how many atoms moves and set the list atomo
|
|
!
|
|
ALLOCATE(ifat(nat))
|
|
ifat = 0
|
|
imode0 = 0
|
|
DO irr = 1, nirr
|
|
if (comp_irr (irr)) then
|
|
do ipert = 1, npert (irr)
|
|
do na = 1, nat
|
|
do ipol = 1, 3
|
|
mu = 3 * (na - 1) + ipol
|
|
if (abs (u (mu, imode0+ipert) ) > 1.d-12) ifat (na) = 1
|
|
enddo
|
|
enddo
|
|
enddo
|
|
endif
|
|
imode0 = imode0 + npert (irr)
|
|
ENDDO
|
|
nat_todo = 0
|
|
DO na = 1, nat
|
|
IF (ifat (na) == 1) THEN
|
|
nat_todo = nat_todo + 1
|
|
atomo (nat_todo) = na
|
|
ENDIF
|
|
ENDDO
|
|
|
|
DEALLOCATE(ifat)
|
|
!
|
|
! Initialize done_irr, find max dimension of the irreps
|
|
!
|
|
all_comp=.true.
|
|
DO irr=1,nirr
|
|
IF (.NOT.comp_irr(irr)) all_comp=.false.
|
|
ENDDO
|
|
all_comp = all_comp.OR.lgamma_gamma
|
|
all_done = .FALSE.
|
|
npertx = 0
|
|
done_irr = .FALSE.
|
|
IF (elph) done_elph = .FALSE.
|
|
DO irr = 1, nirr
|
|
npertx = max (npertx, npert (irr) )
|
|
ENDDO
|
|
!
|
|
! set to zero the variable written on file
|
|
!
|
|
dyn=(0.0_DP,0.0_DP)
|
|
dyn00=(0.0_DP,0.0_DP)
|
|
dyn_rec=(0.0_DP,0.0_DP)
|
|
IF (epsil.and..not.done_epsil) epsilon=0.0_DP
|
|
IF (zeu.and..not.done_zeu) zstareu=0.0_DP
|
|
IF (lraman.and..not.done_lraman) ramtns=0.0_DP
|
|
IF (elop.and..not.done_elop) eloptns=0.0_DP
|
|
|
|
where_rec='phq_setup.'
|
|
rec_code=-40
|
|
CALL ph_writefile('status_ph',current_iq,0,ierr)
|
|
!
|
|
! Bands for the electron-phonon calculation
|
|
!
|
|
IF(elph_tetra == 1 .OR. elph_tetra == 2) THEN
|
|
!
|
|
elph_nbnd_min = nbnd
|
|
elph_nbnd_max = 1
|
|
!
|
|
DO ibnd = 1, nbnd
|
|
IF(MAXVAL(et(ibnd, 1:nks)) > ef .AND. ibnd < elph_nbnd_min) elph_nbnd_min = ibnd
|
|
IF(MINVAL(et(ibnd, 1:nks)) < ef .AND. ibnd > elph_nbnd_max) elph_nbnd_max = ibnd
|
|
END DO
|
|
!
|
|
call mp_min(elph_nbnd_min, inter_pool_comm)
|
|
call mp_max(elph_nbnd_max, inter_pool_comm)
|
|
!
|
|
END IF
|
|
!
|
|
! DFPT+U
|
|
!
|
|
IF (lda_plus_u) THEN
|
|
!
|
|
! Define effU
|
|
!
|
|
effU = 0.d0
|
|
DO nah = 1, nat
|
|
nt = ityp(nah)
|
|
! For U only calculations
|
|
effU(nt) = Hubbard_U(nt)
|
|
! When there is also Hubbard_J0/=0
|
|
IF (Hubbard_J0(nt).NE.0.d0) &
|
|
effU(nt) = Hubbard_U(nt) - Hubbard_J0(nt)
|
|
ENDDO
|
|
!
|
|
! Initialize d1, d2, d3 to rotate the spherical harmonics
|
|
!
|
|
CALL d_matrix (d1, d2, d3)
|
|
!
|
|
ENDIF
|
|
!
|
|
! dVscf Fourier interpolation
|
|
!
|
|
IF (ldvscf_interpolate) THEN
|
|
CALL dvscf_interpol_setup()
|
|
ENDIF
|
|
!
|
|
! AHC e-ph coupling
|
|
!
|
|
IF (elph_ahc) CALL elph_ahc_setup()
|
|
!
|
|
CALL stop_clock ('phq_setup')
|
|
!
|
|
RETURN
|
|
!
|
|
END SUBROUTINE phq_setup
|