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