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 10!this routine calculate the terms \psi_i(r)\_psi_v(sc)(r) 11!and write them on disk on global G grid 12 13 subroutine semicore(n_semicore, num_nbnds,ispin) 14!NOT_TO_BE_INCLUDED_START 15 USE io_global, ONLY : stdout, ionode,ionode_id 16 USE io_files, ONLY : diropn,prefix,tmp_dir 17 use pwcom 18 USE wavefunctions, ONLY : evc 19 USE kinds, ONLY : DP 20 USE gvect, ONLY : ngm_g, ig_l2g, gstart 21 USE mp, ONLY : mp_sum, mp_barrier, mp_bcast 22 USE mp_wave, ONLY : mergewf,splitwf 23 USE mp_pools, ONLY : intra_pool_comm, inter_pool_comm, intra_pool_comm 24 USE mp_world, ONLY : world_comm, mpime, nproc 25 USE fft_base, ONLY : dfftp, dffts 26 USE fft_interfaces, ONLY : fwfft, invfft 27 USE wavefunctions, ONLY : psic 28 USE wvfct, ONLY : et 29 30 implicit none 31 32 INTEGER, EXTERNAL :: find_free_unit 33 INTEGER, INTENT(in) :: n_semicore!number of semicore states 34 INTEGER, INTENT(in) :: num_nbnds!total KS states considered 35 INTEGER, INTENT(in) :: ispin!spin channel 36 37 REAL(kind=DP), ALLOCATABLE :: psi_sc(:,:) 38 COMPLEX(kind=DP), ALLOCATABLE :: prod(:), prod_g(:,:), prod_g_tot(:) 39 INTEGER :: iv, iun, ii 40 INTEGER :: npwx_g 41 42 43!fft trasform semicore states to R space 44 allocate(psi_sc(dfftp%nnr,n_semicore)) 45 allocate(prod(dfftp%nnr), prod_g(npw,2),prod_g_tot(ngm_g)) 46 47 do iv=1,n_semicore,2 48 psic(:)=(0.d0,0.d0) 49 if(iv<n_semicore) then 50 psic(dffts%nl(igk_k(1:npw,1))) = evc(1:npw,iv) + & 51 ( 0.D0, 1.D0 ) * evc(1:npw,iv+1) 52 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,iv) - & 53 ( 0.D0, 1.D0 ) * evc(1:npw,iv+1) ) 54 55 else 56 psic(dffts%nl(igk_k(1:npw,1))) = evc(1:npw,iv) 57 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,iv) ) 58 endif 59 CALL invfft ('Wave', psic, dffts) 60 psi_sc(1:dfftp%nnr,iv)=dble(psic(1:dfftp%nnr)) 61 if(iv< n_semicore) psi_sc(1:dfftp%nnr,iv+1)=dimag(psic(1:dfftp%nnr)) 62 63 enddo 64 65 66!write header of file with KS energies of all states (just in case..) 67 68 npwx_g=npwx 69 call mp_sum(npwx_g, world_comm) 70 71 if(ionode) then 72 iun = find_free_unit() 73 if(ispin==1) then 74 open( unit= iun, file=trim(tmp_dir)//trim(prefix)//'.sc_states', status='unknown',form='unformatted') 75 else 76 open( unit= iun, file=trim(tmp_dir)//trim(prefix)//'.sc_states2', status='unknown',form='unformatted') 77 endif 78 write(iun) num_nbnds 79 write(iun) n_semicore 80 write(iun) npwx_g 81 write(iun) et(1:num_nbnds,ispin) 82 endif 83 84 85!loop (double) on KS states 86 do ii=n_semicore+1,num_nbnds 87 88!fft 89 90 psic(:)=(0.d0,0.d0) 91 psic(dffts%nl(igk_k(1:npw,1))) = evc(1:npw,ii) 92 psic(dffts%nlm(igk_k(1:npw,1))) = CONJG( evc(1:npw,ii) ) 93 94 CALL invfft ('Wave', psic, dffts) 95 96!calculate product 97 do iv=1,n_semicore,2 98 if(iv<n_semicore) then 99 prod(1:dfftp%nnr)=dcmplx(dble(psic(1:dfftp%nnr))*psi_sc(1:dfftp%nnr,iv),& 100 &dble(psic(1:dfftp%nnr))*psi_sc(1:dfftp%nnr,iv+1)) 101 else 102 prod(1:dfftp%nnr)=dble(psic(1:dfftp%nnr)) 103 prod(1:dfftp%nnr)=prod(1:dfftp%nnr)*psi_sc(1:dfftp%nnr,iv) 104 endif 105!fft back 106 107 CALL fwfft ('Rho', prod, dfftp) 108 if(iv==n_semicore) then 109 prod_g(1:npw,1)=prod(dfftp%nl(1:npw)) 110 if(gstart==2) then 111 write(stdout,*) 'Putting to zero:', iv,ii, prod_g(1,1) 112 prod_g(1,1)=(0.d0,0.d0) 113 endif 114 else 115 prod_g(1:npw, 1)= 0.5d0*(prod(dfftp%nl(1:npw))+conjg( prod(dfftp%nlm(1:npw)))) 116 prod_g(1:npw, 2)= (0.d0,-0.5d0)*(prod(dfftp%nl(1:npw)) - conjg(prod(dfftp%nlm(1:npw)))) 117 if(gstart==2) then 118 write(stdout,*)'Putting to zero:', iv,ii, prod_g(1,1) 119 write(stdout,*)'Putting to zero:', iv+1,ii, prod_g(1,2) 120 prod_g(1,1:2)=(0.d0,0.d0) 121 endif 122 endif 123 124!merge 125 126 prod_g_tot(:)=(0.d0,0.d0) 127 call mergewf(prod_g(:,1),prod_g_tot,npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm) 128!put it on disk 129 if(ionode) then 130 write(iun) prod_g_tot(1:npwx_g) 131 endif 132 if(iv<n_semicore) then 133 134!merge 135 prod_g_tot(:)=(0.d0,0.d0) 136 call mergewf(prod_g(:,2),prod_g_tot,npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm) 137!put it on disk 138 if(ionode) then 139 write(iun) prod_g_tot(1:npwx_g) 140 endif 141 142 143 endif 144 enddo 145 146 147 148 149 enddo 150 151 152!now put on disk KS states on global G 153 do ii=n_semicore+1,num_nbnds 154 prod_g_tot(:)=(0.d0,0.d0) 155 call mergewf(evc(:,ii),prod_g_tot,npw,ig_l2g,mpime,nproc,ionode_id,intra_pool_comm) 156!put it on disk 157 if(ionode) then 158 write(iun) prod_g_tot(1:npwx_g) 159 endif 160 enddo 161 162 163 deallocate(psi_sc,prod_g_tot) 164 close(iun) 165 return 166!NOT_TO_BE_INCLUDED_END 167 end subroutine semicore 168