TABLE OF CONTENTS
ABINIT/mkresi [ Functions ]
NAME
mkresi
FUNCTION
Make residuals from knowledge of wf in G space and application of Hamiltonian.
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
cg(2,mcg)=<G|Cnk>=Fourier coefficients of wavefunction gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k icg=shift to be applied on the location of data in the array cg ikpt=index of k-point isppol=index of spin mcg=second dimension of the cg array mpi_enreg=information about MPI parallelization nband=number of bands involved in subspace matrix. npw=number of planewaves in basis sphere at this k point. prtvol=control print volume and debugging output usepaw= 0 for non paw calculation; =1 for paw calculation
OUTPUT
eig_k(nband)$= \langle C_n \mid H \mid C_n \rangle $ for each band. resid_k(nband)=residual for each band $= \langle C_n \mid H H \mid C_n \rangle- \langle C_n \mid H \mid C_n \rangle^2 $.
PARENTS
energy
CHILDREN
dotprod_g,getghc,prep_getghc,sqnorm_g,timab
SOURCE
41 #if defined HAVE_CONFIG_H 42 #include "config.h" 43 #endif 44 45 #include "abi_common.h" 46 47 subroutine mkresi(cg,eig_k,gs_hamk,icg,ikpt,isppol,mcg,mpi_enreg,nband,prtvol,resid_k) 48 49 use defs_basis 50 use defs_abitypes 51 use m_profiling_abi 52 use m_errors 53 use m_xmpi 54 use m_cgtools 55 56 use m_pawcprj, only : pawcprj_type 57 use m_hamiltonian, only : gs_hamiltonian_type 58 59 !This section has been created automatically by the script Abilint (TD). 60 !Do not modify the following lines by hand. 61 #undef ABI_FUNC 62 #define ABI_FUNC 'mkresi' 63 use interfaces_18_timing 64 use interfaces_32_util 65 use interfaces_66_wfs 66 !End of the abilint section 67 68 implicit none 69 70 !Arguments ------------------------------------ 71 !scalars 72 integer,intent(in) :: icg,ikpt,isppol,mcg,nband,prtvol 73 type(MPI_type),intent(inout) :: mpi_enreg 74 type(gs_hamiltonian_type),intent(inout) :: gs_hamk 75 !arrays 76 real(dp),intent(in) :: cg(2,mcg) 77 real(dp),intent(out) :: eig_k(nband),resid_k(nband) 78 79 !Local variables------------------------------- 80 !scalars 81 integer,parameter :: tim_getghc=3 82 integer :: blocksize,cpopt,iband,iband_last,iblock,iblocksize,ipw,ipw_shift 83 integer :: my_nspinor,nblockbd,npw_k 84 real(dp) :: doti,dotr 85 !arrays 86 real(dp) :: tsec(2) 87 real(dp),allocatable,target :: cwavef(:,:),ghc(:,:),gsc(:,:),gvnlc(:,:) 88 real(dp), ABI_CONTIGUOUS pointer :: cwavef_ptr(:,:),ghc_ptr(:,:),gsc_ptr(:,:) 89 type(pawcprj_type) :: cwaveprj(0,0) 90 91 ! ************************************************************************* 92 93 !Keep track of total time spent in mkresi 94 call timab(13,1,tsec) 95 96 !Parallelism over FFT and/or bands: define sizes and tabs 97 my_nspinor=max(1,gs_hamk%nspinor/mpi_enreg%nproc_spinor) 98 if (mpi_enreg%paral_kgb==1) then 99 nblockbd=nband/(mpi_enreg%nproc_band*mpi_enreg%bandpp) 100 else 101 nblockbd=nband/mpi_enreg%nproc_fft 102 if (nband/=nblockbd*mpi_enreg%nproc_fft) nblockbd=nblockbd+1 103 end if 104 blocksize=nband/nblockbd 105 106 npw_k=gs_hamk%npw_k 107 ABI_ALLOCATE(cwavef,(2,npw_k*my_nspinor)) 108 ABI_ALLOCATE(ghc,(2,npw_k*my_nspinor)) 109 ABI_ALLOCATE(gvnlc,(2,npw_k*my_nspinor)) 110 if (gs_hamk%usepaw==1) then 111 ABI_ALLOCATE(gsc,(2,npw_k*my_nspinor)) 112 else 113 ABI_ALLOCATE(gsc,(0,0)) 114 end if 115 116 !Loop over (blocks of) bands 117 do iblock=1,nblockbd 118 iband=(iblock-1)*blocksize+1;iband_last=min(iband+blocksize-1,nband) 119 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband_last,isppol,mpi_enreg%me_kpt)) cycle 120 121 ! Load |Cn> 122 ipw_shift=(iblock-1)*npw_k*my_nspinor*blocksize+icg 123 !$OMP PARALLEL DO 124 do ipw=1,npw_k*my_nspinor*blocksize 125 cwavef(1,ipw)=cg(1,ipw+ipw_shift) 126 cwavef(2,ipw)=cg(2,ipw+ipw_shift) 127 end do 128 129 ! Compute H|Cn> 130 cpopt=-1 131 if (mpi_enreg%paral_kgb==0) then 132 call getghc(cpopt,cwavef,cwaveprj,ghc,gsc,gs_hamk,gvnlc,zero,mpi_enreg,1,& 133 & prtvol,gs_hamk%usepaw,tim_getghc,0) 134 else 135 call prep_getghc(cwavef,gs_hamk,gvnlc,ghc,gsc,zero,nband,mpi_enreg,& 136 & prtvol,gs_hamk%usepaw,cpopt,cwaveprj,& 137 & already_transposed=.false.) 138 end if 139 140 ! Compute the residual, <Cn|(H-<Cn|H|Cn>)**2|Cn>: 141 do iblocksize=1,blocksize 142 iband=(iblock-1)*blocksize+iblocksize 143 ipw_shift=(iblocksize-1)*npw_k*my_nspinor 144 cwavef_ptr => cwavef(:,1+ipw_shift:npw_k*my_nspinor+ipw_shift) 145 ghc_ptr => ghc (:,1+ipw_shift:npw_k*my_nspinor+ipw_shift) 146 147 ! First get eigenvalue <Cn|H|Cn>: 148 call dotprod_g(dotr,doti,gs_hamk%istwf_k,npw_k*my_nspinor,1,cwavef_ptr,ghc_ptr,& 149 & mpi_enreg%me_g0,mpi_enreg%comm_spinorfft) 150 eig_k(iband)=dotr 151 152 ! Next need <G|(H-S<Cn|H|Cn>)|Cn> (in ghc): 153 if (gs_hamk%usepaw==0) then 154 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(cwavef_ptr,ghc_ptr,eig_k,iband,npw_k,my_nspinor) 155 do ipw=1,npw_k*my_nspinor 156 ghc_ptr(1,ipw)=ghc_ptr(1,ipw)-eig_k(iband)*cwavef_ptr(1,ipw) 157 ghc_ptr(2,ipw)=ghc_ptr(2,ipw)-eig_k(iband)*cwavef_ptr(2,ipw) 158 end do 159 else 160 gsc_ptr => gsc(:,1+ipw_shift:npw_k*my_nspinor+ipw_shift) 161 !$OMP PARALLEL DO PRIVATE(ipw) SHARED(gsc_ptr,ghc_ptr,eig_k,iband,npw_k,my_nspinor) 162 do ipw=1,npw_k*my_nspinor 163 ghc_ptr(1,ipw)=ghc_ptr(1,ipw)-eig_k(iband)*gsc_ptr(1,ipw) 164 ghc_ptr(2,ipw)=ghc_ptr(2,ipw)-eig_k(iband)*gsc_ptr(2,ipw) 165 end do 166 end if 167 168 ! Then simply square the result: 169 call sqnorm_g(dotr,gs_hamk%istwf_k,npw_k*my_nspinor,ghc_ptr,& 170 & mpi_enreg%me_g0,mpi_enreg%comm_fft) 171 resid_k(iband)=dotr 172 173 end do ! iblocksize 174 175 end do ! iblock 176 177 ABI_DEALLOCATE(cwavef) 178 ABI_DEALLOCATE(ghc) 179 ABI_DEALLOCATE(gvnlc) 180 ABI_DEALLOCATE(gsc) 181 182 call timab(13,2,tsec) 183 184 end subroutine mkresi