1! 2! Copyright (C) 2001-2013 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 o_rinitcgg( npwx, npw, nstart, numwp, psi, o_evc, e, numv, v_states,hdiag,ptype,fcw_number,fcw_state,fcw_mat,ethr) 11 !---------------------------------------------------------------------------- 12 ! 13 ! ... Operator O diagonalization in the subspace spanned 14 ! ... by nstart states psi (atomic or random wavefunctions). 15 ! ... Produces on output numwp eigenvectors (numwp <= nstart) in o_evc. 16 ! ... Minimal memory use - o_evc and psi may overlap 17 ! ... Calls o_1psi to calculate O|psi> 18 ! 19 USE kinds, ONLY : DP 20 USE gvect, ONLY : gstart 21 USE io_global, ONLY : stdout 22 USE mp, ONLY : mp_sum 23 USE mp_world, ONLY : world_comm 24 USE fft_base, ONLY : dffts 25 USE mp_bands, ONLY : me_bgrp, root_bgrp, intra_bgrp_comm 26 ! 27 IMPLICIT NONE 28 ! 29 include 'laxlib.fh' 30 ! 31 INTEGER :: npw, npwx, nstart, numwp 32 ! dimension of the matrix to be diagonalized 33 ! leading dimension of matrix psi, as declared in the calling pgm unit 34 ! input number of states 35 ! output number of states 36 COMPLEX (DP) :: psi(npwx,nstart), o_evc(npwx,numwp) 37 ! input and output eigenvectors (may overlap) 38 REAL(DP) :: e(numwp) 39 40 INTEGER, INTENT(in) :: numv!number of valence states 41 REAL(kind=DP), INTENT(in) :: v_states(dffts%nnr,numv) !valence states in real space 42 REAL(kind=DP), INTENT(in) :: hdiag(npw)!inverse of estimation of diagonal part of hamiltonian 43 INTEGER, INTENT(in) :: ptype!type of approximation for O operator 44 INTEGER, INTENT(in) :: fcw_number!number of "fake conduction" states for O matrix method 45 COMPLEX(kind=DP) :: fcw_state(npw,fcw_number)! "fake conduction" states for O matrix method 46 REAL(kind=DP) :: fcw_mat(fcw_number,fcw_number)! "fake conduction" matrix 47 REAL(kind=DP), INTENT(in) :: ethr!threshold for o_1psi_gamma 48 49 50 ! eigenvalues 51 ! 52 !... local variables 53 ! 54 INTEGER :: m, i, j, npw2, npwx2 55 REAL (DP) :: rtmp(2) 56 COMPLEX (DP), ALLOCATABLE :: aux(:,:) 57 COMPLEX (DP), ALLOCATABLE :: ctmp(:) 58 REAL (DP), ALLOCATABLE :: hr(:,:,:), sr(:,:) 59 REAL (DP), ALLOCATABLE :: en(:) 60 ! 61 ! 62 CALL start_clock( 'wfcrot1' ) 63 ! 64 npw2 = 2 * npw 65 npwx2 = 2 * npwx 66 ! 67 ALLOCATE( aux( npwx, 2 ) ) 68 ALLOCATE( ctmp( numwp ) ) 69 ALLOCATE( hr( nstart, nstart, 2 ) ) 70 ALLOCATE( sr( nstart, nstart ) ) 71 ALLOCATE( en( nstart ) ) 72 ! 73 ! ... Set up the Hamiltonian and Overlap matrix 74 ! 75 DO m = 1, nstart 76 ! 77 !CALL hs_1psi( npwx, npw, psi(1,m), aux(1,1), aux(1,2) ) 78 79 write(stdout,*) 'Call o_1psi_gamma',m,nstart 80 FLUSH(stdout) 81 call o_1psi_gamma( numv, v_states, psi(1,m), aux(1,1),.false.,hdiag,ptype,fcw_number,fcw_state,fcw_mat,ethr) 82 write(stdout,*) 'Done' 83 FLUSH(stdout) 84 85 !call o_1psi_gamma_real( numv, v_states, psi(1,m), aux(1,1)) 86 aux(:,2)=psi(:,m) 87 88 ! 89 CALL DGEMV( 'T', npw2, 2, 2.D0, aux, npwx2, psi(1,m), 1, 0.D0, rtmp, 1 ) 90 ! 91 IF ( gstart == 2 ) rtmp(:) = rtmp(:) - psi(1,m) * aux(1,:) 92 ! 93 hr(m,m,1) = rtmp(1) 94 sr(m,m) = rtmp(2) 95 ! 96 DO j = m + 1, nstart 97 ! 98 CALL DGEMV( 'T', npw2, 2, 2.D0, aux, npwx2, psi(1,j), 1, 0.D0, rtmp, 1 ) 99 ! 100 IF ( gstart == 2 ) rtmp(:) = rtmp(:) - psi(1,j) * aux(1,:) 101 ! 102 hr(j,m,1) = rtmp(1) 103 sr(j,m) = rtmp(2) 104 ! 105 hr(m,j,1) = rtmp(1) 106 sr(m,j) = rtmp(2) 107 ! 108 END DO 109 ! 110 END DO 111 ! 112 !CALL reduce( nstart * nstart, hr(1,1,1) ) 113 call mp_sum(hr(:,:,1),world_comm) 114 !CALL reduce( nstart * nstart, sr(1,1) ) 115 CALL mp_sum(sr(:,:),world_comm) 116 ! 117 ! ... diagonalize 118 ! 119 write(stdout,*) 'Call rdiaghg' 120 FLUSH(stdout) 121 122 CALL diaghg( nstart, numwp, hr(:,:,1), sr, nstart, en, hr(:,:,2), me_bgrp, root_bgrp, intra_bgrp_comm ) 123 write(stdout,*) 'Done' 124 FLUSH(stdout) 125 126 ! 127 e(1:numwp) = en(1:numwp) 128 ! 129 ! ... update the basis set 130 ! 131 DO i = 1, npw 132 ! 133 DO m = 1, numwp 134 ! 135 ctmp(m) = SUM( hr(:,m,2) * psi(i,:) ) 136 ! 137 END DO 138 ! 139 o_evc(i,1:numwp) = ctmp(1:numwp) 140 ! 141 END DO 142 ! 143 DEALLOCATE( en ) 144 DEALLOCATE( sr ) 145 DEALLOCATE( hr ) 146 DEALLOCATE( ctmp ) 147 DEALLOCATE( aux ) 148 ! 149 CALL stop_clock( 'wfcrot1' ) 150 ! 151 RETURN 152 ! 153END SUBROUTINE o_rinitcgg 154