! This file is copied and modified from QUANTUM ESPRESSO ! Kun Cao, Henry Lambert, Feliciano Giustino ! Copyright (C) 2001-2004 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 . ! ! !---------------------------------------------------------------------- subroutine cgsolve_all (h_psi, cg_psi, e, d0psi, dpsi, h_diag, & ndmx, ndim, ethr, ik, kter, conv_root, anorm, nbnd, npol) !---------------------------------------------------------------------- ! ! iterative solution of the linear system: ! ! ( h - e + Q ) * dpsi = d0psi (1) ! ! where h is a complex hermitean matrix, e is a real sca ! dpsi and d0psi are complex vectors ! ! on input: ! h_psi EXTERNAL name of a subroutine: ! h_psi(ndim,psi,psip) ! Calculates H*psi products. ! Vectors psi and psip should be dimensined ! (ndmx,nvec). nvec=1 is used! ! ! cg_psi EXTERNAL name of a subroutine: ! g_psi(ndmx,ndim,notcnv,psi,e) ! which calculates (h-e)^-1 * psi, with ! some approximation, e.g. (diag(h)-e) ! ! e real unperturbed eigenvalue. ! ! dpsi contains an estimate of the solution ! vector. ! ! d0psi contains the right hand side vector ! of the system. ! ! ndmx integer row dimension of dpsi, ecc. ! ! ndim integer actual row dimension of dpsi ! ! ethr real convergence threshold. solution ! improvement is stopped when the error in ! eq (1), defined as l.h.s. - r.h.s., becomes ! less than ethr in norm. ! ! on output: dpsi contains the refined estimate of the ! solution vector. ! ! d0psi is corrupted on exit ! ! revised (extensively) 6 Apr 1997 by A. Dal Corso & F. Mauri ! revised (to reduce memory) 29 May 2004 by S. de Gironcoli ! USE kinds, ONLY : DP USE mp_global, ONLY : intra_pool_comm USE mp, ONLY : mp_sum USE control_flags, ONLY : gamma_only USE gvect, ONLY : gstart USE io_global, ONLY : stdout, ionode implicit none ! ! first the I/O variables ! integer :: ndmx, & ! input: the maximum dimension of the vectors ndim, & ! input: the actual dimension of the vectors kter, & ! output: counter on iterations nbnd, & ! input: the number of bands npol, & ! input: number of components of the wavefunctions ik ! input: the k point real(DP) :: & e(nbnd), & ! input: the actual eigenvalue anorm, & ! output: the norm of the error in the solution h_diag(ndmx*npol,nbnd), & ! input: an estimate of ( H - \epsilon ) ethr ! input: the required precision complex(DP) :: & dpsi (ndmx*npol, nbnd), & ! output: the solution of the linear syst d0psi (ndmx*npol, nbnd) ! input: the known term logical :: conv_root ! output: if true the root is converged external h_psi ! input: the routine computing h_psi external cg_psi ! input: the routine computing cg_psi ! ! here the local variables ! integer, parameter :: maxiter = 200 ! the maximum number of iterations integer :: iter, ibnd, lbnd ! counters on iteration, bands integer , allocatable :: conv (:) ! if 1 the root is converged complex(DP), allocatable :: g (:,:), t (:,:), h (:,:), hold (:,:) ! the gradient of psi ! the preconditioned gradient ! the delta gradient ! the conjugate gradient ! work space complex(DP) :: dcgamma, dclambda ! the ratio between rho ! step length complex(DP), external :: zdotc REAL(kind=dp), EXTERNAL :: ddot ! the scalar product real(DP), allocatable :: rho (:), rhoold (:), eu (:), a(:), c(:) ! the residue ! auxiliary for h_diag real(DP) :: kter_eff ! account the number of iterations with b ! coefficient of quadratic form ! call start_clock ('cgsolve') allocate ( g(ndmx*npol,nbnd), t(ndmx*npol,nbnd), h(ndmx*npol,nbnd), & hold(ndmx*npol ,nbnd) ) allocate (a(nbnd), c(nbnd)) allocate (conv ( nbnd)) allocate (rho(nbnd),rhoold(nbnd)) allocate (eu ( nbnd)) ! WRITE( stdout,*) g,t,h,hold kter_eff = 0.d0 do ibnd = 1, nbnd conv (ibnd) = 0 enddo g=(0.d0,0.d0) t=(0.d0,0.d0) h=(0.d0,0.d0) hold=(0.d0,0.d0) do iter = 1, maxiter ! ! compute the gradient. can reuse information from previous step ! if (iter == 1) then call h_psi (ndim, dpsi, g, e, ik, nbnd) ! (H-e+\alpha |p>