1! 2! Copyright (C) 2001-2019 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!---------------------------------------------------------------------- 9SUBROUTINE lr_dvpsi_e(ik,ipol,dvpsi) 10 !---------------------------------------------------------------------- 11 ! 12 ! On output: dvpsi contains P_c^+ x | psi_ik > in crystal axis 13 ! (projected on at(*,ipol) ) 14 ! 15 ! dvpsi is computed here (vkb and evc must be set) 16 ! and it is written on file in the routine lr_solve_e. 17 ! 18 ! See Ref.[1] : J. Tobik and A. Dal Corso, JCP 120, 9934 (2004) 19 ! for the details of the theory implemented in this routine. 20 ! 21 ! Modified by Osman Baris Malcioglu (2009) 22 ! Rebased wrt PHONON routines. S J Binnie (2011) 23 ! Modified by Iurii Timrov (2016) 24 ! 25 USE kinds, ONLY : DP 26 USE cell_base, ONLY : tpiba2, at 27 USE ions_base, ONLY : ntyp => nsp 28 USE io_global, ONLY : stdout 29 USE klist, ONLY : xk, ngk, igk_k 30 USE wvfct, ONLY : npwx, nbnd, g2kin, et 31 USE wavefunctions, ONLY : evc 32 USE gvect, ONLY : g 33 USE noncollin_module, ONLY : npol 34 USE becmod, ONLY : allocate_bec_type, calbec, becp, & 35 & deallocate_bec_type, bec_type 36 USE uspp, ONLY : okvan, nkb, vkb 37 USE uspp_param, ONLY : nh, nhm 38 USE control_flags, ONLY : gamma_only 39 USE control_lr, ONLY : nbnd_occ 40 USE lr_variables, ONLY : lr_verbosity, sevc0 41 USE io_global, ONLY : stdout 42 USE lrus, ONLY : dpqq 43 ! 44 IMPLICIT NONE 45 ! 46 INTEGER, INTENT(in) :: ipol, ik 47 COMPLEX(kind=dp), INTENT(out) :: dvpsi(npwx*npol,nbnd) 48 ! 49 ! Local variables 50 ! 51 INTEGER :: ig, ibnd, lter, npw 52 ! counters 53 ! npw: number of plane-waves at point ik 54 REAL(kind=dp) :: atnorm 55 COMPLEX(kind=dp),ALLOCATABLE :: d0psi(:,:) 56 REAL(DP), ALLOCATABLE :: h_diag (:,:) 57 ! diagonal part of h_scf 58 real(DP) :: anorm 59 ! preconditioning cut-off 60 REAL(DP), PARAMETER :: thresh = 1.0e-5_DP 61 ! the desired convergence of linter 62 LOGICAL :: conv_root 63 ! true if convergence has been achieved 64 COMPLEX(DP), ALLOCATABLE :: spsi(:,:) 65 ! 66 TYPE(bec_type) :: becp1 67 TYPE(bec_type) :: becp2 68 ! 69 EXTERNAL ch_psi_all, cg_psi 70 ! 71 CALL start_clock ('lr_dvpsi_e') 72 ! 73 conv_root = .TRUE. 74 ! 75 ALLOCATE ( d0psi(npwx*npol,nbnd) ) 76 d0psi = (0.d0, 0.d0) 77 dvpsi = (0.d0, 0.d0) 78 ! 79 npw = ngk(ik) 80 ! 81 CALL allocate_bec_type ( nkb, nbnd, becp1 ) 82 ! 83 CALL calbec ( npw, vkb, evc, becp1 ) 84 ! 85 CALL allocate_bec_type ( nkb, nbnd, becp2 ) 86 ! 87 CALL commutator_Hx_psi (ik, nbnd_occ(ik), becp1, becp2, ipol, d0psi ) 88 ! 89 IF (okvan) CALL calbec ( npw, vkb, evc, becp, nbnd) 90 ! 91 ! Orthogonalize d0psi to the valence subspace. Apply P_c^+ 92 ! 93 CALL orthogonalize(d0psi, evc, ik, ik, sevc0(:,:,ik), npw, .true.) 94 d0psi = -d0psi 95 ! 96 ! Calculate the kinetic energy g2kin: (k+G)^2 97 ! 98 CALL g2_kin(ik) 99 ! 100 ! Calculate the preconditioning matrix h_diag used by cgsolve_all 101 ! 102 ALLOCATE ( h_diag(npwx*npol, nbnd) ) 103 CALL h_prec(ik, evc, h_diag) 104 ! 105 ! d0psi contains P^+_c [H-eS,x] psi_v for the polarization direction ipol 106 ! Now solve the linear systems (H+Q-e_vS)*P_c(x*psi_v)=P_c^+ [H-e_vS,x]*psi_v 107 ! See Eq.(9) in Ref. [1] 108 ! 109 CALL cgsolve_all (ch_psi_all, cg_psi, et (1, ik), d0psi, dvpsi, & 110 h_diag, npwx, npw, thresh, ik, lter, conv_root, anorm, nbnd_occ(ik), 1) 111 ! 112 IF (.not.conv_root) WRITE( stdout, '(5x,"ik",i4," ibnd",i4, & 113 & " lr_dvpsi_e: root not converged ",e10.3)') & 114 ik, ibnd, anorm 115 ! 116 FLUSH( stdout ) 117 DEALLOCATE (h_diag) 118 ! 119 ! In the US case we obtain P_c x |psi>, but we need P_c^+ x | psi>, 120 ! therefore we apply S again, and then subtract the additional term 121 ! furthermore we add the term due to dipole of the augmentation charges. 122 ! See Eq.(10) in Ref. [1] 123 ! 124 IF (okvan) THEN 125 ALLOCATE (spsi ( npwx*npol, nbnd)) 126 CALL calbec (npw, vkb, dvpsi, becp ) 127 CALL s_psi(npwx,npw,nbnd,dvpsi,spsi) 128 CALL DCOPY(2*npwx*npol*nbnd,spsi,1,dvpsi,1) 129 DEALLOCATE (spsi) 130 ALLOCATE (dpqq( nhm, nhm, 3, ntyp)) 131 CALL compute_qdipol(dpqq) 132 CALL qdipol_cryst() 133 CALL adddvepsi_us(becp1,becp2,ipol,ik,dvpsi) 134 DEALLOCATE (dpqq) 135 ENDIF 136 ! 137 IF (okvan) CALL calbec ( npw, vkb, evc, becp, nbnd) 138 ! 139 ! Orthogonalize dvpsi to the valence subspace. Apply P_c^+ 140 ! 141 CALL orthogonalize(dvpsi, evc, ik, ik, sevc0(:,:,ik), npw, .true.) 142 dvpsi = -dvpsi 143 ! 144 DEALLOCATE (d0psi) 145 ! 146 ! US case: apply the S^{-1} operator 147 ! 148 IF (okvan) THEN 149 ALLOCATE (spsi ( npwx*npol, nbnd)) 150 CALL lr_sm1_initialize() 151 CALL lr_sm1_psi(ik,npwx,ngk(ik),nbnd,dvpsi,spsi) 152 dvpsi(:,:) = spsi(:,:) 153 DEALLOCATE(spsi) 154 ENDIF 155 ! 156 ! For some ibrav the crystal axes are not normalized 157 ! Here we include the correct normalization 158 ! for Lanczos initial wfcs 159 ! 160 atnorm = DSQRT(at(1,ipol)**2 + at(2,ipol)**2 + at(3,ipol)**2) 161 ! 162 dvpsi(:,:) = dvpsi(:,:) / atnorm 163 ! 164 CALL deallocate_bec_type ( becp1 ) 165 IF (nkb > 0) CALL deallocate_bec_type ( becp2 ) 166 ! 167 CALL stop_clock ('lr_dvpsi_e') 168 ! 169 RETURN 170 ! 171END SUBROUTINE lr_dvpsi_e 172 173 174