! ! 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 . ! !----------------------------------------------------------------------- subroutine ccgdiagg (nmax, n, nbnd, psi, e, precondition, eps, & maxter, reorder, notconv, avg_iter, g2kin, vr) !----------------------------------------------------------------------- ! ! FG - slightly modified to be consistent with the modules, ! to remove preprocessing stuff, and the clock, and add g2kin and vr for h_psi ! ! "poor man" iterative diagonalization of a complex hermitian matrix ! through preconditioned conjugate gradient algorithm ! Band-by-band algorithm with minimal use of memory ! Calls h_1psi and s_1psi to calculate H|psi> and S|psi> ! Works for generalized eigenvalue problem (US pseudopotentials) as well ! !@@@ use parameters, only : DP use gspace, only : ngm, nr !@@@ implicit none ! logical :: reorder integer :: nmax, n, nbnd, notconv, maxter complex(kind=DP) :: psi(nmax,nbnd) real(kind=DP) :: e (nbnd), precondition (n), eps, avg_iter ! integer :: i, j, m, iter, moved complex(kind=DP), allocatable :: spsi (:), lagrange (:), & hpsi (:), g (:), cg (:), scg (:), ppsi (:), g0 (:) complex(kind=DP) :: ZDOTC real(kind=DP) :: norma, a0, b0, gg0, gamma, gg, gg1, theta, cg0, e0, & es (2), DDOT external ZDOTC, DDOT real(kind=8), parameter:: pi = 3.14159265358979d0 ! !@@@ real(kind=DP) :: g2kin(ngm) complex(kind=DP) :: vr(nr) !@@@ allocate (spsi(nmax)) allocate (scg( nmax)) allocate (hpsi(nmax)) allocate (g( nmax)) allocate (cg( nmax)) allocate (g0( nmax)) allocate (ppsi(nmax)) allocate (lagrange( nbnd)) avg_iter = 0.d0 notconv = 0 moved = 0 ! every eigenfunction is calculated separately do m = 1, nbnd ! ! calculate S|psi> ! call s_1psi (nmax, n, psi (1, m), spsi) ! ! orthogonalize starting eigenfunction to those already calculated ! do j = 1, m lagrange (j) = ZDOTC (n, psi (1, j), 1, spsi, 1) enddo norma = DREAL (lagrange (m) ) do j = 1, m - 1 call ZAXPY (n, - lagrange (j), psi (1, j), 1, psi (1, m), 1) norma = norma - (DREAL (lagrange (j) ) **2 + DIMAG (lagrange (j) ) & **2) enddo norma = sqrt (norma) call DSCAL (2 * n, 1.d0 / norma, psi (1, m), 1) ! ! calculate starting gradient (|hpsi> = H|psi>)... ! ! call h_1psi (nmax, n, psi (1, m), hpsi, spsi) call h_1psi (nmax, n, psi (1, m), hpsi, spsi, g2kin, vr) ! ! ...and starting eigenvalue (e = = ) ! NB: DDOT(2*n,a,1,b,1) = DREAL(ZDOTC(n,a,1,b,1)) ! e (m) = DDOT (2 * n, psi (1, m), 1, hpsi, 1) ! ! start iteration for this band ! do iter = 1, maxter ! ! calculate P (PHP)|y> (P=preconditioning matrix, assumed diagonal) ! do i = 1, n g (i) = hpsi (i) / precondition (i) ppsi (i) = spsi (i) / precondition (i) enddo ! ! ppsi is now S P(P^2)|y> = S P^2|psi>) ! es (1) = DDOT (2 * n, spsi, 1, g, 1) es (2) = DDOT (2 * n, spsi, 1, ppsi, 1) es (1) = es (1) / es (2) call DAXPY (2 * n, - es (1), ppsi, 1, g, 1) ! ! e1 = / ensures that = 0 ! orthogonalize to lowest eigenfunctions (already calculated) ! ! scg is used as workspace ! call s_1psi (nmax, n, g, scg) do j = 1, m - 1 lagrange (j) = ZDOTC (n, psi (1, j), 1, scg, 1) enddo do j = 1, m - 1 call ZAXPY (n, - lagrange (j), psi (1, j), 1, g, 1) call ZAXPY (n, - lagrange (j), psi (1, j), 1, scg, 1) enddo if (iter.ne.1) then ! ! gg1 is (used in Polak-Ribiere formula) ! gg1 = DDOT (2 * n, g, 1, g0, 1) endif ! ! gg is ! call ZCOPY (n, scg, 1, g0, 1) do i = 1, n g0 (i) = g0 (i) * precondition (i) enddo gg = DDOT (2 * n, g, 1, g0, 1) ! Here one can test on the norm of the gradient ! if (sqrt(gg).lt. eps.or. iter.eq.maxter) go to 10 ! WRITE( stdout,*) iter, gg if (iter.eq.1) then ! ! starting iteration, the conjugate gradient |cg> = |g> ! gg0 = gg call ZCOPY (n, g, 1, cg, 1) else ! ! following iterations: |cg(n+1)> = |g(n+1)> + gamma(n) |cg(n)> ! This is Fletcher-Reeves formula for gamma ! gamma = gg/gg0 ! and this is Polak-Ribiere formula ! gamma = (gg - gg1) / gg0 gg0 = gg call DSCAL (2 * n, gamma, cg, 1) call DAXPY (2 * n, 1.d0, g, 1, cg, 1) ! ! The following is needed because is not 0 ! In fact, = sin(theta) !** norma = DDOT(2*n,psi(1,m),1,cg,1) !** call reduce(1,norma) ! norma = gamma * cg0 * sin (theta) call DAXPY (2 * n, - norma, psi (1, m), 1, cg, 1) endif ! ! |cg> contains now the conjugate gradient ! ! |scg> is S|cg> ! ! call h_1psi (nmax, n, cg, ppsi, scg) call h_1psi (nmax, n, cg, ppsi, scg, g2kin, vr) cg0 = DDOT (2 * n, cg, 1, scg, 1) cg0 = sqrt (cg0) ! ! |ppsi> contains now HP|cg> ! minimize , where |y(t)>=cos(t)|y> + sin(t)/cg0 |cg> ! Note that = 1 , = 0 , = cg0^2 ! so that the result is correctly normalized: = 1 ! a0 = 2.d0 * DDOT (2 * n, psi (1, m), 1, ppsi, 1) / cg0 b0 = DDOT (2 * n, cg, 1, ppsi, 1) / cg0**2 e0 = e (m) theta = atan (a0 / (e0 - b0) ) / 2.d0 es (1) = ( (e0 - b0) * cos (2.d0 * theta) + a0 * sin (2.d0 * & theta) + e0 + b0) / 2.d0 es (2) = ( - (e0 - b0) * cos (2.d0 * theta) - a0 * sin (2.d0 * & theta) + e0 + b0) / 2.d0 ! ! there are two possible solutions, choose the minimum ! if (es (2) .lt.es (1) ) theta = theta + pi / 2.d0 ! ! new estimate of the eigenvalue ! e (m) = min (es (1), es (2) ) ! ! upgrade |psi> ... ! call DSCAL (2 * n, cos (theta), psi (1, m), 1) call DAXPY (2 * n, sin (theta) / cg0, cg, 1, psi (1, m), 1) ! ! here one could test convergence on the energy ! if (abs (e (m) - e0) .lt.eps) goto 10 ! ! .... S|psi>.... ! call DSCAL (2 * n, cos (theta), spsi, 1) call DAXPY (2 * n, sin (theta) / cg0, scg, 1, spsi, 1) ! ! ...and H|psi> ! call DSCAL (2 * n, cos (theta), hpsi, 1) call DAXPY (2 * n, sin (theta) / cg0, ppsi, 1, hpsi, 1) enddo notconv = notconv + 1 10 avg_iter = avg_iter + iter + 1 ! reorder eigenvalues if they are not in the right order ! (this CAN and WILL happen in not-so-special cases) if (m.gt.1.and.reorder) then if (e (m) - e (m - 1) .lt. - 2.d0 * eps) then ! if the last calculated eigenvalue is not the largest... do i = m - 2, 1, - 1 if (e (m) - e (i) .gt.2.d0 * eps) goto 20 enddo 20 i = i + 1 moved = moved+1 ! last calculated eigenvalue should be in the i-th position: reorder e0 = e (m) call ZCOPY (n, psi (1, m), 1, ppsi, 1) do j = m, i + 1, - 1 e (j) = e (j - 1) call ZCOPY (n, psi (1, j - 1), 1, psi (1, j), 1) enddo e (i) = e0 call ZCOPY (n, ppsi, 1, psi (1, i), 1) ! this procedure should be good if only a few inversions occur, ! extremely inefficient if eigenvectors are often in bad order ! (but this should not happen) endif endif enddo avg_iter = avg_iter / nbnd deallocate (lagrange) deallocate (ppsi) deallocate (g0) deallocate (cg) deallocate (g) deallocate (hpsi) deallocate (scg) deallocate (spsi) return end subroutine ccgdiagg ! ! 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 s_1psi( lda, n, psi, spsi ) !----------------------------------------------------------------------- ! ! fake subroutine ... FG ! use parameters, only : DP IMPLICIT NONE ! INTEGER :: lda, n COMPLEX(KIND=DP) :: psi(n), spsi(n) ! ! CALL ZCOPY( n, psi, 1, spsi, 1 ) ! RETURN ! END SUBROUTINE s_1psi ! ! 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 . ! !----------------------------------------------------------------------- subroutine h_1psi (lda, n, psic, hpsic, spsi, g2kin, vr) !----------------------------------------------------------------------- ! ! FG - wrapper for h_psi ! h_psi is for real vectors, while the original subs were for complex ! here makes the conversion... ! !@@@ use parameters, only : DP use gspace, only : ngm, nr !@@@ implicit none ! integer :: lda, n complex(kind=DP) :: psic (n), hpsic (n), spsi(n) !@@@ real(kind=DP) :: g2kin(ngm) complex(kind=DP) :: vr(nr) real(kind=DP) :: psi (n), hpsi (n) !@@@ psi = dreal (psic) hpsi = dreal (hpsic) ! call h_psi (psi, hpsi, g2kin, vr) call s_1psi (lda, n, psic, spsi) ! psic = dcmplx (psi, 0.d0) hpsic = dcmplx (hpsi, 0.d0) ! return end subroutine h_1psi