2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
! Copyright (C) 2003-2004 PWSCF group
|
2003-12-10 22:57:07 +08:00
|
|
|
! 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 .
|
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
#define ZERO ( 0.D0, 0.D0 )
|
|
|
|
#define ONE ( 1.D0, 0.D0 )
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
#include "f_defs.h"
|
|
|
|
!
|
|
|
|
!----------------------------------------------------------------------------
|
2005-03-30 22:37:55 +08:00
|
|
|
SUBROUTINE regterg( ndim, ndmx, nvec, nvecx, evc, ethr, overlap, &
|
|
|
|
gstart, e, btype, notcnv, lrot, dav_iter )
|
2004-07-05 14:53:39 +08:00
|
|
|
!----------------------------------------------------------------------------
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
! ... iterative solution of the eigenvalue problem:
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
! ... ( H - e S ) * evc = 0
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
! ... where H is an hermitean operator, e is a real scalar,
|
|
|
|
! ... S is an overlap matrix, evc is a complex vector
|
|
|
|
! ... (real wavefunctions with only half plane waves stored)
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
|
|
|
USE io_global, ONLY : stdout
|
2004-07-05 14:53:39 +08:00
|
|
|
USE kinds, ONLY : DP
|
|
|
|
!
|
|
|
|
IMPLICIT NONE
|
|
|
|
!
|
|
|
|
! ... on INPUT
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
INTEGER, INTENT(IN) :: ndim, ndmx, nvec, nvecx, gstart
|
2004-07-05 14:53:39 +08:00
|
|
|
! dimension of the matrix to be diagonalized
|
|
|
|
! leading dimension of matrix evc, as declared in the calling pgm unit
|
|
|
|
! integer number of searched low-lying roots
|
|
|
|
! maximum dimension of the reduced basis set
|
|
|
|
! (the basis set is refreshed when its dimension would exceed nvecx)
|
2005-03-30 22:37:55 +08:00
|
|
|
COMPLEX (KIND=DP), INTENT(INOUT) :: evc(ndmx,nvec)
|
2004-07-05 14:53:39 +08:00
|
|
|
! evc contains the refined estimates of the eigenvectors
|
2005-03-30 22:37:55 +08:00
|
|
|
REAL (KIND=DP), INTENT(IN) :: ethr
|
2004-07-05 14:53:39 +08:00
|
|
|
! energy threshold for convergence: root improvement is stopped,
|
|
|
|
! when two consecutive estimates of the root differ by less than ethr.
|
2005-03-30 22:37:55 +08:00
|
|
|
LOGICAL, INTENT(IN) :: overlap
|
2004-07-05 14:53:39 +08:00
|
|
|
! if .FALSE. : S|psi> not needed
|
|
|
|
INTEGER, INTENT(IN) :: btype(nvec)
|
|
|
|
! band type ( 1 = occupied, 0 = empty )
|
2005-03-30 22:37:55 +08:00
|
|
|
INTEGER, INTENT(IN) :: lrot
|
|
|
|
! .TRUE. if the wfc have already be rotated
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... on OUTPUT
|
|
|
|
!
|
|
|
|
REAL (KIND=DP) :: e(nvec)
|
|
|
|
! contains the estimated roots.
|
2005-03-30 22:37:55 +08:00
|
|
|
INTEGER :: dav_iter, notcnv
|
2004-07-05 14:53:39 +08:00
|
|
|
! integer number of iterations performed
|
|
|
|
! number of unconverged roots
|
|
|
|
!
|
|
|
|
! ... LOCAL variables
|
|
|
|
!
|
|
|
|
INTEGER, PARAMETER :: maxter = 20
|
|
|
|
! maximum number of iterations
|
|
|
|
!
|
|
|
|
INTEGER :: kter, nbase, np, n, m
|
|
|
|
! counter on iterations
|
|
|
|
! dimension of the reduced basis
|
|
|
|
! counter on the reduced basis vectors
|
|
|
|
! do-loop counters
|
2004-08-20 00:22:51 +08:00
|
|
|
REAL (KIND=DP), ALLOCATABLE :: hr(:,:), sr(:,:), vr(:,:), ew(:)
|
2004-07-05 14:53:39 +08:00
|
|
|
! Hamiltonian on the reduced basis
|
|
|
|
! S matrix on the reduced basis
|
|
|
|
! eigenvectors of the Hamiltonian
|
|
|
|
! eigenvalues of the reduced hamiltonian
|
|
|
|
COMPLEX (KIND=DP), ALLOCATABLE :: psi(:,:), hpsi(:,:), spsi(:,:)
|
|
|
|
! work space, contains psi
|
|
|
|
! the product of H and psi
|
|
|
|
! the product of S and psi
|
|
|
|
LOGICAL, ALLOCATABLE :: conv(:)
|
|
|
|
! true if the root is converged
|
|
|
|
REAL (KIND=DP) :: empty_ethr
|
|
|
|
! threshold for empty bands
|
2004-08-20 00:22:51 +08:00
|
|
|
INTEGER :: ndim2, ndmx2
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... Called routines:
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
REAL (KIND=DP), EXTERNAL :: DDOT
|
|
|
|
!
|
2004-08-27 16:29:48 +08:00
|
|
|
EXTERNAL h_psi, s_psi, g_psi
|
2004-07-05 14:53:39 +08:00
|
|
|
! h_psi(ndmx,ndim,nvec,psi,hpsi)
|
|
|
|
! calculates H|psi>
|
|
|
|
! s_psi(ndmx,ndim,nvec,psi,spsi)
|
|
|
|
! calculates S|psi> (if needed)
|
|
|
|
! Vectors psi,hpsi,spsi are dimensioned (ndmx,nvec)
|
|
|
|
! g_psi(ndmx,ndim,notcnv,psi,e)
|
|
|
|
! calculates (diag(h)-e)^-1 * psi, diagonal approx. to (h-e)^-1*psi
|
|
|
|
! the first nvec columns contain the trial eigenvectors
|
|
|
|
!
|
|
|
|
!
|
|
|
|
CALL start_clock( 'cegterg' )
|
|
|
|
!
|
|
|
|
! ... ALLOCATE the work arrays
|
|
|
|
!
|
|
|
|
ALLOCATE( psi( ndmx, nvecx ) )
|
|
|
|
ALLOCATE( hpsi( ndmx, nvecx ) )
|
|
|
|
IF ( overlap ) ALLOCATE( spsi( ndmx, nvecx ) )
|
|
|
|
ALLOCATE( sr( nvecx, nvecx ) )
|
|
|
|
ALLOCATE( hr( nvecx, nvecx ) )
|
|
|
|
ALLOCATE( vr( nvecx, nvecx ) )
|
|
|
|
ALLOCATE( ew( nvecx ) )
|
|
|
|
ALLOCATE( conv( nvec ) )
|
|
|
|
!
|
|
|
|
IF ( nvec > nvecx / 2 ) CALL errore( 'regter', 'nvecx is too small', 1 )
|
|
|
|
!
|
|
|
|
! ... threshold for empty bands
|
|
|
|
!
|
|
|
|
empty_ethr = MAX( ( ethr * 5.D0 ), 1.D-5 )
|
|
|
|
!
|
|
|
|
! ... prepare the hamiltonian for the first iteration
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
ndim2 = 2 * ndim
|
|
|
|
ndmx2 = 2 * ndmx
|
2003-12-10 22:57:07 +08:00
|
|
|
notcnv = nvec
|
2004-07-05 14:53:39 +08:00
|
|
|
nbase = nvec
|
2004-08-20 00:22:51 +08:00
|
|
|
conv = .FALSE.
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( overlap ) spsi = ZERO
|
|
|
|
psi = ZERO
|
|
|
|
hpsi = ZERO
|
|
|
|
psi(:,1:nvec) = evc(:,1:nvec)
|
|
|
|
!
|
|
|
|
! ... hpsi contains h times the basis vectors
|
|
|
|
!
|
|
|
|
CALL h_psi( ndmx, ndim, nvec, psi, hpsi )
|
|
|
|
!
|
|
|
|
IF ( overlap ) CALL s_psi( ndmx, ndim, nvec, psi, spsi )
|
|
|
|
!
|
|
|
|
! ... hr contains the projection of the hamiltonian onto the reduced space
|
|
|
|
! ... vr contains the eigenvectors of hr
|
|
|
|
!
|
|
|
|
hr(:,:) = 0.D0
|
2005-03-30 22:37:55 +08:00
|
|
|
sr(:,:) = 0.D0
|
2004-07-05 14:53:39 +08:00
|
|
|
vr(:,:) = 0.D0
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'T', 'N', nbase, nbase, ndim2, 2.D0 , &
|
|
|
|
psi, ndmx2, hpsi, ndmx2, 0.D0, hr, nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) &
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGER( nbase, nbase, -1.D0, psi, ndmx2, hpsi, ndmx2, hr, nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL reduce( nbase * nvecx, hr )
|
|
|
|
!
|
|
|
|
IF ( overlap ) THEN
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'T', 'N', nbase, nbase, ndim2, 2.D0, &
|
|
|
|
psi, ndmx2, spsi, ndmx2, 0.D0, sr, nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) &
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGER( nbase, nbase, -1.D0, psi, ndmx2, spsi, ndmx2, sr, nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
ELSE
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'T', 'N', nbase, nbase, ndim2, 2.D0, &
|
|
|
|
psi, ndmx2, psi, ndmx2, 0.D0, sr, nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) &
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGER( nbase, nbase, -1.D0, psi, ndmx2, psi, ndmx2, sr, nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
|
|
|
CALL reduce( nbase * nvecx, sr )
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
IF ( lrot ) THEN
|
2004-08-20 00:22:51 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
DO n = 1, nbase
|
|
|
|
!
|
|
|
|
e(n) = hr(n,n)
|
|
|
|
!
|
|
|
|
vr(n,n) = 1.D0
|
|
|
|
!
|
|
|
|
END DO
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
ELSE
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
! ... diagonalize the reduced hamiltonian
|
|
|
|
!
|
|
|
|
CALL rdiaghg( nbase, nvec, hr, sr, nvecx, ew, vr )
|
|
|
|
!
|
|
|
|
e(1:nvec) = ew(1:nvec)
|
|
|
|
!
|
|
|
|
END IF
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... iterate
|
|
|
|
!
|
|
|
|
iterate: DO kter = 1, maxter
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
dav_iter = kter
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL start_clock( 'update' )
|
|
|
|
!
|
2003-12-10 22:57:07 +08:00
|
|
|
np = 0
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
DO n = 1, nvec
|
|
|
|
!
|
|
|
|
IF ( .NOT. conv(n) ) THEN
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
! ... this root not yet converged ...
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
|
|
|
np = np + 1
|
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
! ... reorder eigenvectors so that coefficients for unconverged
|
|
|
|
! ... roots come first. This allows to use quick matrix-matrix
|
|
|
|
! ... multiplications to set a new basis vector (see below)
|
|
|
|
!
|
|
|
|
IF ( np /= n ) vr(:,np) = vr(:,n)
|
|
|
|
!
|
|
|
|
! ... for use in g_psi
|
|
|
|
!
|
|
|
|
ew(nbase+np) = e(n)
|
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
|
|
|
END DO
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
! ... expand the basis set with new basis vectors ( H - e*S )|psi> ...
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( overlap ) THEN
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'N', 'N', ndim2, notcnv, nbase, 1.D0, spsi, &
|
|
|
|
ndmx2, vr, nvecx, 0.D0, psi(1,nbase+1), ndmx2 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
ELSE
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'N', 'N', ndim2, notcnv, nbase, 1.D0, psi, &
|
|
|
|
ndmx2, vr, nvecx, 0.D0, psi(1,nbase+1), ndmx2 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
2004-10-25 23:26:30 +08:00
|
|
|
DO np = 1, notcnv
|
2005-03-30 22:37:55 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
psi(:,nbase+np) = - ew(nbase+np) * psi(:,nbase+np)
|
2005-03-30 22:37:55 +08:00
|
|
|
!
|
2004-10-25 23:26:30 +08:00
|
|
|
END DO
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'N', 'N', ndim2, notcnv, nbase, 1.D0, hpsi, &
|
|
|
|
ndmx2, vr, nvecx, 1.D0, psi(1,nbase+1), ndmx2 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL stop_clock( 'update' )
|
|
|
|
!
|
|
|
|
! ... approximate inverse iteration
|
|
|
|
!
|
|
|
|
CALL g_psi( ndmx, ndim, notcnv, psi(1,nbase+1), ew(nbase+1) )
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
! ... "normalize" correction vectors psi(:,nbase+1:nbase+notcnv) in order
|
2005-03-30 22:37:55 +08:00
|
|
|
! ... to improve numerical stability of subspace diagonalization cdiaghg
|
|
|
|
! ... ew is used as work array :
|
|
|
|
!
|
|
|
|
! ... ew = <psi_i|psi_i>, i = nbase + 1, nbase + notcnv
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
DO n = 1, notcnv
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
ew(n) = 2.D0 * DDOT( ndim2, psi(1,nbase+n), 1, psi(1,nbase+n), 1 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) ew(n) = ew(n) - psi(1,nbase+n) * psi(1,nbase+n)
|
|
|
|
!
|
|
|
|
END DO
|
|
|
|
!
|
|
|
|
CALL reduce( notcnv, ew )
|
|
|
|
!
|
2004-10-25 23:26:30 +08:00
|
|
|
DO n = 1, notcnv
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2004-07-19 21:36:17 +08:00
|
|
|
psi(:,nbase+n) = psi(:,nbase+n) / SQRT( ew(n) )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2004-10-25 23:26:30 +08:00
|
|
|
END DO
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... here compute the hpsi and spsi of the new functions
|
|
|
|
!
|
|
|
|
CALL h_psi( ndmx, ndim, notcnv, psi(1,nbase+1), hpsi(1,nbase+1) )
|
|
|
|
!
|
|
|
|
IF ( overlap ) &
|
|
|
|
CALL s_psi( ndmx, ndim, notcnv, psi(1,nbase+1), spsi(1,nbase+1) )
|
|
|
|
!
|
|
|
|
! ... update the reduced hamiltonian
|
|
|
|
!
|
|
|
|
CALL start_clock( 'overlap' )
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'T', 'N', nbase+notcnv, notcnv, ndim2, 2.D0, psi, &
|
|
|
|
ndmx2, hpsi(1,nbase+1), ndmx2, 0.D0, hr(1,nbase+1), nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) &
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGER( nbase+notcnv, notcnv, -1.D0, psi, ndmx2, &
|
|
|
|
hpsi(1,nbase+1), ndmx2, hr(1,nbase+1), nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL reduce( nvecx * notcnv, hr(1,nbase+1) )
|
|
|
|
!
|
|
|
|
IF ( overlap ) THEN
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'T', 'N', nbase+notcnv, notcnv, ndim2, 2.D0, psi, ndmx2, &
|
|
|
|
spsi(1,nbase+1), ndmx2, 0.D0, sr(1,nbase+1), nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) &
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGER( nbase+notcnv, notcnv, -1.D0, psi, ndmx2, &
|
|
|
|
spsi(1,nbase+1), ndmx2, sr(1,nbase+1), nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
ELSE
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'T', 'N', nbase+notcnv, notcnv, ndim2, 2.D0, psi, ndmx2, &
|
|
|
|
psi(1,nbase+1), ndmx2, 0.D0, sr(1,nbase+1) , nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( gstart == 2 ) &
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGER( nbase+notcnv, notcnv, -1.D0, psi, ndmx2, &
|
|
|
|
psi(1,nbase+1), ndmx2, sr(1,nbase+1), nvecx )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
|
|
|
CALL reduce( nvecx * notcnv, sr(1,nbase+1) )
|
|
|
|
!
|
|
|
|
CALL stop_clock( 'overlap' )
|
|
|
|
!
|
|
|
|
nbase = nbase + notcnv
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
DO n = 1, nbase
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
DO m = n + 1, nbase
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
hr(m,n) = hr(n,m)
|
|
|
|
sr(m,n) = sr(n,m)
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
END DO
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
END DO
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... diagonalize the reduced hamiltonian
|
|
|
|
!
|
|
|
|
CALL rdiaghg( nbase, nvec, hr, sr, nvecx, ew, vr )
|
|
|
|
!
|
|
|
|
! ... test for convergence
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
WHERE( btype(1:nvec) == 1 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
conv(1:nvec) = ( ( ABS( ew(1:nvec) - e(1:nvec) ) < ethr ) )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
ELSEWHERE
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
conv(1:nvec) = ( ( ABS( ew(1:nvec) - e(1:nvec) ) < empty_ethr ) )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
END WHERE
|
|
|
|
!
|
|
|
|
notcnv = COUNT( .NOT. conv(:) )
|
|
|
|
!
|
|
|
|
e(1:nvec) = ew(1:nvec)
|
|
|
|
!
|
|
|
|
! ... if overall convergence has been achieved, OR
|
|
|
|
! ... the dimension of the reduced basis set is becoming too large, OR
|
|
|
|
! ... in any case if we are at the last iteration
|
|
|
|
! ... refresh the basis set. i.e. replace the first nvec elements
|
|
|
|
! ... with the current estimate of the eigenvectors;
|
|
|
|
! ... set the basis dimension to nvec.
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
IF ( notcnv == 0 .OR. nbase+notcnv > nvecx .OR. dav_iter == maxter ) THEN
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL start_clock( 'last' )
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'N', 'N', ndim2, nvec, nbase, 1.D0, &
|
|
|
|
psi, ndmx2, vr, nvecx, 0.D0, evc, ndmx2 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
IF ( notcnv == 0 ) THEN
|
|
|
|
!
|
|
|
|
! ... all roots converged: return
|
|
|
|
!
|
|
|
|
CALL stop_clock( 'last' )
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
2004-07-05 14:53:39 +08:00
|
|
|
EXIT iterate
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
ELSE IF ( dav_iter == maxter ) THEN
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... last iteration, some roots not converged: return
|
|
|
|
!
|
|
|
|
WRITE( UNIT = stdout, &
|
2005-03-30 22:37:55 +08:00
|
|
|
FMT = '(5X,"WARNING: ",I5," eigenvalues not converged")' ) &
|
2003-12-10 22:57:07 +08:00
|
|
|
notcnv
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL stop_clock( 'last' )
|
|
|
|
!
|
|
|
|
EXIT iterate
|
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
|
|
|
! ... refresh psi, H*psi and S*psi
|
|
|
|
!
|
|
|
|
psi(:,1:nvec) = evc(:,1:nvec)
|
|
|
|
!
|
|
|
|
IF ( overlap ) THEN
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'N', 'N', ndim2, nvec, nbase, 1.D0, spsi, &
|
|
|
|
ndmx2, vr, nvecx, 0.D0, psi(1,nvec+1), ndmx2 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
spsi(:,1:nvec) = psi(:,nvec+1:nvec+nvec)
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
2004-08-20 00:22:51 +08:00
|
|
|
CALL DGEMM( 'N', 'N', ndim2, nvec, nbase, 1.D0, hpsi, &
|
|
|
|
ndmx2, vr, nvecx, 0.D0, psi(1,nvec+1), ndmx2 )
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
hpsi(:,1:nvec) = psi(:,nvec+1:nvec+nvec)
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
! ... refresh the reduced hamiltonian
|
2003-12-10 22:57:07 +08:00
|
|
|
!
|
|
|
|
nbase = nvec
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
hr(:,1:nbase) = 0.D0
|
|
|
|
sr(:,1:nbase) = 0.D0
|
|
|
|
vr(:,1:nbase) = 0.D0
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
DO n = 1, nbase
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
hr(n,n) = e(n)
|
|
|
|
sr(n,n) = 1.D0
|
|
|
|
vr(n,n) = 1.D0
|
|
|
|
!
|
2005-03-30 22:37:55 +08:00
|
|
|
END DO
|
2004-07-05 14:53:39 +08:00
|
|
|
!
|
|
|
|
CALL stop_clock( 'last' )
|
|
|
|
!
|
|
|
|
END IF
|
|
|
|
!
|
|
|
|
END DO iterate
|
|
|
|
!
|
|
|
|
DEALLOCATE( conv )
|
|
|
|
DEALLOCATE( ew )
|
|
|
|
DEALLOCATE( vr )
|
|
|
|
DEALLOCATE( hr )
|
|
|
|
DEALLOCATE( sr )
|
|
|
|
IF ( overlap ) DEALLOCATE( spsi )
|
|
|
|
DEALLOCATE( hpsi )
|
|
|
|
DEALLOCATE( psi )
|
|
|
|
!
|
|
|
|
CALL stop_clock( 'cegterg' )
|
|
|
|
!
|
|
|
|
RETURN
|
|
|
|
!
|
|
|
|
END SUBROUTINE regterg
|