1!
2! Copyright (C) 2001-2008 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!
9!-----------------------------------------------------------------------
10subroutine solve_e2
11  !-----------------------------------------------------------------------
12  !
13  !   Self consistent cycle to compute the second order derivatives
14  !   of the wavefunctions with respect to electric fields
15  !
16  USe kinds,                 ONLY : DP
17  USE io_global,             ONLY : stdout
18  USE cell_base,             ONLY : tpiba2
19  USE klist,                 ONLY : ltetra, lgauss, wk, xk, ngk, igk_k
20  USE lsda_mod,              ONLY : lsda, nspin
21  USE gvect,                 ONLY : g
22  USE gvecs,                 ONLY : doublegrid
23  USE fft_base,              ONLY : dfftp, dffts
24  USE fft_interfaces,        ONLY : fft_interpolate
25  USE wvfct,                 ONLY : npwx, nbnd, et
26  USE buffers,   ONLY: get_buffer
27  USE ions_base, ONLY: nat
28  USE uspp,      ONLY: okvan, nkb, vkb
29  USE uspp_param,ONLY : nhm
30  USE wavefunctions,  ONLY: evc
31  USE control_ph, ONLY : convt, nmix_ph, alpha_mix, tr2_ph, &
32                         niter_ph, rec_code, flmixdpot, rec_code_read
33  USE units_lr,   ONLY : lrwfc, iuwfc
34  USE ramanm,     ONLY : lrba2, iuba2, lrd2w, iud2w
35  USE recover_mod, ONLY : read_rec, write_rec
36
37  USE check_stop, ONLY: check_stop_now
38  USE mp_pools,   ONLY : inter_pool_comm
39  USE mp_bands,   ONLY : intra_bgrp_comm
40  USE mp,         ONLY : mp_sum
41
42  USE eqv,       ONLY : dpsi, dvpsi
43  USE qpoint,    ONLY : nksq, ikks, ikqs
44  USE control_lr, ONLY : nbnd_occ, lgamma
45  USE dv_of_drho_lr
46
47  implicit none
48
49  real(DP) ::  thresh, weight, avg_iter, dr2
50  ! convergence threshold for the solution of the
51  ! linear system
52  ! used for summation over k points
53  ! average number of iterations
54  ! convergence limit
55
56  complex(DP) , pointer :: dvscfin (:,:,:), dvscfins (:,:,:)
57  ! change of the scf potential (input)
58  ! change of the scf potential (smooth)
59
60  complex(DP) , allocatable :: dvscfout (:,:,:), dbecsum (:,:), &
61                                    aux1 (:)
62  ! change of the scf potential (output)
63  ! auxiliary space
64  ! auxiliary space
65
66  logical :: exst
67  ! used to open the recover file
68
69  integer :: npw, npwq, ikk, ikq
70  integer :: kter, iter0, ipol, ibnd, iter, ik, is, ig, iig, irr, ir, nrec, ios
71  ! counter on iterations
72  ! counter on perturbations
73  ! counter on bands
74  ! counter on iterations
75  ! counter on k points
76  ! counter on G vectors
77  ! counter on g vectors
78  ! counter on mesh points
79  ! the record number
80  ! integer variable for I/O control
81
82  if (lsda) call errore ('solve_e2', ' LSDA not implemented', 1)
83  if (okvan) call errore ('solve_e2', ' Ultrasoft PP not implemented', 1)
84
85  call start_clock('solve_e2')
86  allocate (dvscfin( dfftp%nnr, nspin, 6))
87  if (doublegrid) then
88     allocate (dvscfins(dffts%nnr, nspin, 6))
89  else
90     dvscfins => dvscfin
91  endif
92  allocate (dvscfout( dfftp%nnr , nspin, 6))
93  allocate (dbecsum( nhm*(nhm+1)/2, nat))
94  allocate (aux1(dffts%nnr))
95  convt=.FALSE.
96  if (rec_code_read == -10) then
97     ! restarting in Raman
98     CALL read_rec(dr2, iter0, 6, dvscfin, dvscfins)
99  else
100     iter0 = 0
101  end if
102  if (convt) goto 155
103  !
104  if ( (lgauss .or. ltetra) .or..not.lgamma) &
105        call errore ('solve_e2', 'called in the wrong case', 1)
106  !
107  !   The outside loop is over the iterations
108  !
109
110  do kter = 1, niter_ph
111
112     iter = kter + iter0
113     avg_iter = 0.d0
114
115     dvscfout (:,:,:) = (0.d0, 0.d0)
116     dbecsum (:,:) = (0.d0, 0.d0)
117
118     do ik = 1, nksq
119        ! in this routine, ikk=ikq=ik always (q=0)
120        ikk = ikks(ik)
121        ikq = ikqs(ik)
122        !
123        ! reads unperturbed wavefunctions psi_k in G_space, for all bands
124        !
125        if (nksq.gt.1) call get_buffer(evc, lrwfc, iuwfc, ikk)
126        npw = ngk(ikk)
127        npwq= ngk(ikq)
128        !
129        ! compute beta functions and kinetic energy for k-point ikk
130        ! needed by h_psi, called by pcgreen
131        !
132        CALL init_us_2 (npw, igk_k(1,ikk), xk (1, ikk), vkb)
133        CALL g2_kin(ikk)
134        !
135        ! The counter on the polarizations runs only on the 6 inequivalent
136        ! indexes --see the comment on raman.F--
137        !
138        do ipol = 1, 6
139           nrec = (ipol - 1) * nksq + ik
140
141           if (kter.eq.1) then
142              dpsi (:,:) = (0.d0, 0.d0)
143           else
144              call davcio (dpsi, lrd2w, iud2w, nrec, -1)
145           endif
146
147           if (iter.eq.1) then
148              dvscfin (:,:,:) = (0.d0, 0.d0)
149              call davcio (dvpsi, lrba2, iuba2, nrec, -1)
150              thresh = 1.0d-2
151           else
152              call davcio (dvpsi, lrba2, iuba2, nrec, -1)
153              do ibnd = 1, nbnd_occ (ik)
154                 call cft_wave (ik, evc (1, ibnd), aux1, +1)
155                 do ir = 1, dffts%nnr
156                    aux1 (ir) = aux1 (ir) * dvscfins (ir, 1, ipol)
157                 enddo
158                 call cft_wave (ik, dvpsi (1, ibnd), aux1, -1)
159              enddo
160              thresh = min (0.1d0 * sqrt(dr2), 1.0d-2)
161           endif
162
163           call pcgreen (avg_iter, thresh, ik, et (1, ik) )
164           call davcio ( dpsi, lrd2w, iud2w, nrec, +1)
165           !
166           ! calculates dvscf, sum over k => dvscf_q_ipert
167           !
168           weight = wk (ik)
169           call incdrhoscf (dvscfout (1,1,ipol), weight, ik, &
170                            dbecsum (1, 1), dpsi)
171           enddo   ! on perturbations
172        enddo      ! on k points
173     call mp_sum ( dbecsum, intra_bgrp_comm )
174     if (doublegrid) then
175        do is = 1, nspin
176           do ipol = 1, 6
177              call fft_interpolate (dffts, dvscfout (:, is, ipol), dfftp, dvscfout (:, is, ipol))
178           enddo
179        enddo
180     endif
181     !
182     !   After the loop over the perturbations we have the change of the pote
183     !   for all the modes, and we symmetrize this potential
184     !
185     call mp_sum ( dvscfout, inter_pool_comm )
186     do ipol = 1, 6
187        call dv_of_drho (dvscfout (1, 1, ipol), .false.)
188     enddo
189
190     call psyme2(dvscfout)
191     !
192     ! Mixing with the old potential
193     !
194     call mix_potential (2 * 6 * dfftp%nnr* nspin, dvscfout, dvscfin,  &
195                         alpha_mix (kter), dr2, 6 * tr2_ph, iter,  &
196                         nmix_ph, flmixdpot, convt)
197
198     if (doublegrid) then
199        do is = 1, nspin
200           do ipol = 1, 6
201              call fft_interpolate (dfftp,dvscfin (:, is, ipol), dffts, dvscfins (:, is, ipol))
202           enddo
203        enddo
204     end if
205
206     write (6, "(//,5x,' iter # ',i3, &
207          &      '   av.it.: ',f5.1)") iter, avg_iter / (6.d0 * nksq)
208     dr2 = dr2 / 6
209     write (6, "(5x,' thresh=',es10.3, ' alpha_mix = ',f6.3, &
210          &      ' |ddv_scf|^2 = ',es10.3 )") thresh, alpha_mix (kter), dr2
211     !
212     FLUSH( stdout )
213     !
214     ! rec_code: state of the calculation
215     ! rec_code=-10  to -19 Raman
216     rec_code=-10
217     CALL write_rec('solve_e2..', irr, dr2, iter, convt, 6, dvscfin)
218
219     if ( check_stop_now() ) call stop_smoothly_ph (.false.)
220     if ( convt ) goto 155
221
222  enddo
223155 continue
224  deallocate (dvscfin )
225  if (doublegrid) deallocate (dvscfins )
226  deallocate (dvscfout )
227  deallocate (dbecsum )
228  deallocate (aux1 )
229
230  call stop_clock('solve_e2')
231
232  return
233end subroutine solve_e2
234