2005-03-05 01:46:02 +08:00
|
|
|
!
|
|
|
|
! Copyright (C) 2001 PWSCF 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 .
|
|
|
|
!
|
2005-05-25 10:58:35 +08:00
|
|
|
#include "f_defs.h"
|
2005-03-05 01:46:02 +08:00
|
|
|
!
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
subroutine solve_e2
|
|
|
|
!-----------------------------------------------------------------------
|
|
|
|
!
|
|
|
|
! Self consistent cycle to compute the second order derivatives
|
|
|
|
! of the wavefunctions with respect to electric fields
|
|
|
|
!
|
|
|
|
use kinds, only : DP
|
2008-08-24 01:55:06 +08:00
|
|
|
USE io_global, ONLY : stdout
|
2005-03-05 01:46:02 +08:00
|
|
|
use pwcom
|
|
|
|
use becmod
|
2008-08-24 01:55:06 +08:00
|
|
|
USE io_files, ONLY: prefix, iunigk
|
2005-03-05 01:46:02 +08:00
|
|
|
USE ions_base, ONLY: nat
|
2008-08-24 01:55:06 +08:00
|
|
|
USE uspp, ONLY: okvan, nkb, vkb
|
|
|
|
USE uspp_param,ONLY : nhm
|
2005-03-05 01:46:02 +08:00
|
|
|
USE wavefunctions_module, ONLY: evc
|
|
|
|
USE phcom
|
|
|
|
USE ramanm
|
2006-11-16 07:43:43 +08:00
|
|
|
USE check_stop, ONLY: check_stop_now
|
2008-01-25 07:43:05 +08:00
|
|
|
USE mp_global, ONLY : inter_pool_comm, intra_pool_comm
|
|
|
|
USE mp, ONLY : mp_sum
|
2005-03-05 01:46:02 +08:00
|
|
|
implicit none
|
|
|
|
|
2005-08-28 22:09:42 +08:00
|
|
|
real(DP) :: thresh, weight, avg_iter, dr2
|
2005-03-05 01:46:02 +08:00
|
|
|
! convergence threshold for the solution of the
|
|
|
|
! linear system
|
|
|
|
! used for summation over k points
|
|
|
|
! average number of iterations
|
|
|
|
! convergence limit
|
|
|
|
|
2005-08-28 22:09:42 +08:00
|
|
|
complex(DP) , pointer :: dvscfin (:,:,:), dvscfins (:,:,:)
|
2005-03-05 01:46:02 +08:00
|
|
|
! change of the scf potential (input)
|
|
|
|
! change of the scf potential (smooth)
|
|
|
|
|
2005-08-28 22:09:42 +08:00
|
|
|
complex(DP) , allocatable :: dvscfout (:,:,:), dbecsum (:,:), &
|
2005-04-15 23:43:33 +08:00
|
|
|
aux1 (:)
|
2005-03-05 01:46:02 +08:00
|
|
|
! change of the scf potential (output)
|
|
|
|
! auxiliary space
|
|
|
|
! auxiliary space
|
|
|
|
|
2005-08-28 22:09:42 +08:00
|
|
|
complex(DP) :: ZDOTC
|
2005-03-05 01:46:02 +08:00
|
|
|
! the scalar product function
|
|
|
|
|
|
|
|
logical :: exst
|
|
|
|
! used to open the recover file
|
|
|
|
|
2005-04-21 20:41:18 +08:00
|
|
|
integer :: kter, iter0, ipol, ibnd, iter, ik, is, ig, iig, irr, ir, nrec, ios
|
2005-03-05 01:46:02 +08:00
|
|
|
! counter on iterations
|
|
|
|
! counter on perturbations
|
|
|
|
! counter on bands
|
|
|
|
! counter on iterations
|
|
|
|
! counter on k points
|
|
|
|
! counter on G vectors
|
|
|
|
! counter on g vectors
|
|
|
|
! counter on mesh points
|
|
|
|
! the record number
|
|
|
|
! integer variable for I/O control
|
|
|
|
|
2005-05-20 05:17:42 +08:00
|
|
|
character (len=256) :: flmixdpot
|
2005-03-05 01:46:02 +08:00
|
|
|
! the name of the file with the
|
|
|
|
! mixing potential
|
|
|
|
|
|
|
|
external ch_psi_all, cg_psi
|
|
|
|
|
|
|
|
if (lsda) call errore ('solve_e2', ' LSDA not implemented', 1)
|
2006-12-12 00:52:18 +08:00
|
|
|
if (okvan) call errore ('solve_e2', ' Ultrasoft PP not implemented', 1)
|
2005-03-05 01:46:02 +08:00
|
|
|
|
2005-04-21 22:37:21 +08:00
|
|
|
call start_clock('solve_e2')
|
2005-03-05 01:46:02 +08:00
|
|
|
allocate (dvscfin( nrxx, nspin, 6))
|
|
|
|
if (doublegrid) then
|
|
|
|
allocate (dvscfins( nrxxs, nspin, 6))
|
|
|
|
else
|
|
|
|
dvscfins => dvscfin
|
|
|
|
endif
|
|
|
|
allocate (dvscfout( nrxx , nspin, 6))
|
|
|
|
allocate (dbecsum( nhm*(nhm+1)/2, nat))
|
|
|
|
allocate (aux1( nrxxs))
|
2008-07-23 16:46:48 +08:00
|
|
|
if (rec_code == -10) then
|
2005-04-21 20:41:18 +08:00
|
|
|
! restarting in Raman
|
2006-12-11 23:14:02 +08:00
|
|
|
read (iunrec) iter0, dr2
|
2008-09-16 22:40:58 +08:00
|
|
|
read (iunrec) this_pcxpsi_is_on_file
|
|
|
|
read (iunrec) zstareu0, zstarue0
|
2005-04-21 20:41:18 +08:00
|
|
|
read (iunrec) dvscfin
|
2006-12-12 00:52:18 +08:00
|
|
|
if (okvan) read (iunrec) int1, int2, int3
|
2005-03-05 01:46:02 +08:00
|
|
|
close (unit = iunrec, status = 'keep')
|
|
|
|
if (doublegrid) then
|
|
|
|
do is = 1, nspin
|
|
|
|
do ipol = 1, 6
|
|
|
|
call cinterpolate (dvscfin (1, is, ipol), &
|
|
|
|
dvscfins (1, is, ipol), -1)
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
end if
|
2005-04-21 20:41:18 +08:00
|
|
|
else
|
|
|
|
iter0 = 0
|
2005-03-05 01:46:02 +08:00
|
|
|
end if
|
2005-04-21 20:41:18 +08:00
|
|
|
!
|
2006-12-19 02:28:42 +08:00
|
|
|
if (lgauss.or..not.lgamma) &
|
2005-03-05 01:46:02 +08:00
|
|
|
call errore ('solve_e2', 'called in the wrong case', 1)
|
|
|
|
!
|
|
|
|
! The outside loop is over the iterations
|
|
|
|
!
|
|
|
|
if (reduce_io) then
|
|
|
|
flmixdpot = ' '
|
|
|
|
else
|
2005-06-17 21:27:38 +08:00
|
|
|
flmixdpot = 'mixd'
|
2005-03-05 01:46:02 +08:00
|
|
|
endif
|
|
|
|
|
|
|
|
do kter = 1, niter_ph
|
|
|
|
|
|
|
|
iter = kter + iter0
|
|
|
|
avg_iter = 0.d0
|
|
|
|
|
|
|
|
dvscfout (:,:,:) = (0.d0, 0.d0)
|
|
|
|
dbecsum (:,:) = (0.d0, 0.d0)
|
|
|
|
if (nksq.gt.1) rewind (unit = iunigk)
|
|
|
|
|
|
|
|
do ik = 1, nksq
|
|
|
|
if (nksq.gt.1) then
|
|
|
|
read (iunigk, err = 100, iostat = ios) npw, igk
|
|
|
|
100 call errore ('solve_e', 'reading igk', abs (ios) )
|
|
|
|
endif
|
|
|
|
!
|
|
|
|
! reads unperturbed wavefuctions psi_k in G_space, for all bands
|
|
|
|
!
|
|
|
|
if (nksq.gt.1) call davcio (evc, lrwfc, iuwfc, ik, -1)
|
|
|
|
npwq = npw
|
|
|
|
call init_us_2 (npw, igk, xk (1, ik), vkb)
|
|
|
|
!
|
|
|
|
! compute the kinetic energy
|
|
|
|
!
|
|
|
|
do ig = 1, npwq
|
|
|
|
iig = igkq (ig)
|
|
|
|
g2kin (ig) = ( (xk (1, ik) + g (1, iig) ) **2 + &
|
|
|
|
(xk (2, ik) + g (2, iig) ) **2 + &
|
|
|
|
(xk (3, ik) + g (3, iig) ) **2 ) * tpiba2
|
|
|
|
enddo
|
|
|
|
!
|
|
|
|
! The counter on the polarizations runs only on the 6 inequivalent
|
|
|
|
! indexes --see the comment on raman.F--
|
|
|
|
!
|
|
|
|
do ipol = 1, 6
|
|
|
|
nrec = (ipol - 1) * nksq + ik
|
|
|
|
|
|
|
|
if (kter.eq.1) then
|
|
|
|
dpsi (:,:) = (0.d0, 0.d0)
|
|
|
|
else
|
|
|
|
call davcio (dpsi, lrd2w, iud2w, nrec, -1)
|
|
|
|
endif
|
|
|
|
|
|
|
|
if (iter.eq.1) then
|
|
|
|
dvscfin (:,:,:) = (0.d0, 0.d0)
|
|
|
|
call davcio (dvpsi, lrba2, iuba2, nrec, -1)
|
|
|
|
thresh = 1.0d-2
|
|
|
|
else
|
|
|
|
call davcio (dvpsi, lrba2, iuba2, nrec, -1)
|
|
|
|
do ibnd = 1, nbnd_occ (ik)
|
|
|
|
call cft_wave (evc (1, ibnd), aux1, +1)
|
|
|
|
do ir = 1, nrxxs
|
|
|
|
aux1 (ir) = aux1 (ir) * dvscfins (ir, 1, ipol)
|
|
|
|
enddo
|
|
|
|
call cft_wave (dvpsi (1, ibnd), aux1, -1)
|
|
|
|
enddo
|
|
|
|
thresh = min (0.1d0 * sqrt(dr2), 1.0d-2)
|
|
|
|
endif
|
|
|
|
|
2005-04-15 23:43:33 +08:00
|
|
|
call pcgreen (avg_iter, thresh, ik, et (1, ik) )
|
2005-03-05 01:46:02 +08:00
|
|
|
call davcio ( dpsi, lrd2w, iud2w, nrec, +1)
|
|
|
|
!
|
|
|
|
! calculates dvscf, sum over k => dvscf_q_ipert
|
|
|
|
!
|
|
|
|
weight = wk (ik)
|
|
|
|
call incdrhoscf (dvscfout (1,1,ipol), weight, ik, &
|
2007-02-08 21:57:01 +08:00
|
|
|
dbecsum (1, 1))
|
2005-03-05 01:46:02 +08:00
|
|
|
enddo ! on perturbations
|
|
|
|
enddo ! on k points
|
|
|
|
#ifdef __PARA
|
2008-04-19 18:14:45 +08:00
|
|
|
call mp_sum ( dbecsum, intra_pool_comm )
|
2005-03-05 01:46:02 +08:00
|
|
|
#endif
|
|
|
|
if (doublegrid) then
|
|
|
|
do is = 1, nspin
|
|
|
|
do ipol = 1, 6
|
|
|
|
call cinterpolate (dvscfout (1, is, ipol), &
|
|
|
|
dvscfout (1, is, ipol), 1)
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
endif
|
|
|
|
|
|
|
|
! call addusddense (dvscfout, dbecsum)
|
|
|
|
|
|
|
|
!
|
|
|
|
! After the loop over the perturbations we have the change of the pote
|
|
|
|
! for all the modes, and we symmetrize this potential
|
|
|
|
!
|
|
|
|
#ifdef __PARA
|
2008-01-25 07:43:05 +08:00
|
|
|
call mp_sum ( dvscfout, inter_pool_comm )
|
2005-03-05 01:46:02 +08:00
|
|
|
#endif
|
|
|
|
do ipol = 1, 6
|
|
|
|
call dv_of_drho (0, dvscfout (1, 1, ipol), .false.)
|
|
|
|
enddo
|
|
|
|
|
|
|
|
#ifdef __PARA
|
|
|
|
call psyme2(dvscfout)
|
|
|
|
#else
|
|
|
|
call syme2(dvscfout)
|
|
|
|
#endif
|
|
|
|
!
|
|
|
|
! Mixing with the old potential
|
|
|
|
!
|
|
|
|
call mix_potential (2 * 6 * nrxx * nspin, dvscfout, dvscfin, &
|
|
|
|
alpha_mix (kter), dr2, 6 * tr2_ph, kter, &
|
|
|
|
nmix_ph, flmixdpot, convt)
|
|
|
|
|
|
|
|
if (doublegrid) then
|
|
|
|
do is = 1, nspin
|
|
|
|
do ipol = 1, 6
|
|
|
|
call cinterpolate (dvscfin (1, is, ipol), &
|
|
|
|
dvscfins (1, is, ipol), -1)
|
|
|
|
enddo
|
|
|
|
enddo
|
|
|
|
end if
|
|
|
|
|
|
|
|
write (6, "(//,5x,' iter # ',i3, &
|
|
|
|
& ' av.it.: ',f5.1)") iter, avg_iter / (6.d0 * nksq)
|
|
|
|
dr2 = dr2 / 6
|
|
|
|
write (6, "(5x,' thresh=',e10.3, ' alpha_mix = ',f6.3, &
|
|
|
|
& ' |ddv_scf|^2 = ',e10.3 )") thresh, alpha_mix (kter), dr2
|
2005-05-25 10:58:35 +08:00
|
|
|
!
|
|
|
|
CALL flush_unit( stdout )
|
|
|
|
!
|
2008-07-23 16:46:48 +08:00
|
|
|
! rec_code: state of the calculation
|
|
|
|
! rec_code=-10 to -19 Raman
|
|
|
|
rec_code=-10
|
|
|
|
CALL write_rec('solve_e2..', irr, dr2, iter, convt, dvscfin, 6)
|
2005-04-21 20:41:18 +08:00
|
|
|
|
2008-07-23 16:46:48 +08:00
|
|
|
if ( check_stop_now() ) call stop_ph (.false.)
|
2006-11-16 07:43:43 +08:00
|
|
|
if ( convt ) goto 155
|
2005-03-05 01:46:02 +08:00
|
|
|
|
|
|
|
enddo
|
2005-04-21 20:41:18 +08:00
|
|
|
155 continue
|
2005-03-05 01:46:02 +08:00
|
|
|
deallocate (dvscfin )
|
|
|
|
if (doublegrid) deallocate (dvscfins )
|
|
|
|
deallocate (dvscfout )
|
|
|
|
deallocate (dbecsum )
|
|
|
|
deallocate (aux1 )
|
|
|
|
|
2005-04-21 22:37:21 +08:00
|
|
|
call stop_clock('solve_e2')
|
2005-03-05 01:46:02 +08:00
|
|
|
|
|
|
|
return
|
|
|
|
end subroutine solve_e2
|