TABLE OF CONTENTS
ABINIT/calc_ucrpa [ Functions ]
NAME
calc_ucrpa
FUNCTION
Calculate the screening interaction in the correlated orbitals
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (TApplencourt,BA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
npwe : number of plane wave for the dielectric constant npw : number of plane wave nomega : number of frequencis bandinf,bandsup : kohn sham band optimisation : string for the optimisation Wfd:: MPI communicator mesh <kmesh_t> %nbz=Number of points in the BZ %nibz=Number of points in IBZ %kibz(,nibz)=k-point coordinates, irreducible Brillouin zone %kbz(3,nbz)=k-point coordinates, full Brillouin zone %ktab(nbz)= table giving for each k-point in the BZ (kBZ), the corresponding %ktabi(nbz)= for each k-point in the BZ defines whether inversion has to be considered %ktabp(nbz)= phase factor associated to tnons M1_q_m(bandinf:bandsup,bandinf:bandsup,npw,Qmesh%nibz): Oscillator strengh in Wannier basis rhot1_q_m(bandinf:bandsup,bandinf:bandsup,npw,Qmesh%nibz): Oscillator strengh multiplied by coulomb potential in Wannier basis
OUTPUT
NOTES
PARENTS
sigma
CHILDREN
affichage,checkk,cpu_time,get_bz_item,read_screening,wrtout xmpi_barrier,xmpi_sum
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 subroutine calc_ucrpa(itypatcor,cryst,Kmesh,lpawu,M1_q_m,Qmesh,npwe,& 54 & npw,nsym,rhot1_q_m,nomega,omegamin,omegamax,bandinf,bandsup,optimisation,ucvol,Wfd,fname) 55 56 use defs_basis 57 use m_profiling_abi 58 use m_xmpi 59 use m_errors 60 61 use m_io_tools, only : open_file 62 use m_wfd, only : wfd_t 63 use m_io_screening, only : read_screening, em1_ncname 64 use m_bz_mesh, only : kmesh_t, get_BZ_item 65 use m_crystal, only : crystal_t 66 67 !This section has been created automatically by the script Abilint (TD). 68 !Do not modify the following lines by hand. 69 #undef ABI_FUNC 70 #define ABI_FUNC 'calc_ucrpa' 71 use interfaces_14_hidewrite 72 !End of the abilint section 73 74 implicit none 75 ! _____ _ 76 ! |_ _| | | 77 ! | | _ __ _ __ _ _| |_ 78 ! | | | '_ \| '_ \| | | | __| 79 ! _| |_| | | | |_) | |_| | |_ 80 ! |_____|_| |_| .__/ \__,_|\__| 81 ! | | 82 ! |_| 83 84 !Arguments ------------------------------------ 85 integer, intent(in) :: itypatcor,lpawu,npw,npwe,nsym 86 integer, intent(in) :: nomega 87 integer, intent(in) :: bandinf 88 integer, intent(in) :: bandsup 89 character(len=fnlen), intent(in) :: fname 90 character(len=*), intent(in) :: optimisation 91 real(dp), intent(in) :: ucvol,omegamin,omegamax 92 93 type(wfd_t),intent(inout) :: Wfd 94 type(kmesh_t),intent(in) :: Kmesh,Qmesh 95 type(crystal_t),intent(in) :: Cryst 96 complex(dpc), intent(in) :: rhot1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,npw,Qmesh%nibz) 97 complex(dpc), intent(in) :: M1_q_m(cryst%nattyp(itypatcor),Wfd%nspinor,Wfd%nspinor,2*lpawu+1,2*lpawu+1,npw,Qmesh%nibz) 98 99 !Local variables ------------------------------ 100 !scalars 101 real(dp) :: x 102 real(dp) :: t1,t2 103 real(dp):: tol,eVnorme 104 complex(dpc) :: uu,jj 105 106 complex :: nC,ualter 107 108 integer :: iat,im_paral,iqalloc,ib1,ib2,m1,m2,m3,m4 109 integer :: ierr,ik_bz,ik_ibz,iq_ibz,i,iG1,iG2,iG,iiG,iomega,iomega1,ispinor1,ispinor2,ispinor3,ispinor4 110 integer :: lpawu_read,nkibz,nbband,nkbz,nprocs,nqalloc,nqibz,ms1,ms2,ms3,ms4,mbband,nspinor 111 integer :: isym_kgw,iik,unt 112 complex(dpc) ::ph_mkt 113 114 logical :: wannier=.TRUE. 115 logical :: verbose=.FALSE. 116 logical :: bug=.FALSE. 117 logical :: lscr_one 118 119 character(len=500) :: message 120 121 !arrays 122 complex(dpc), allocatable :: V_m(:,:,:,:) 123 complex(dpc), allocatable :: U_m(:,:,:,:) 124 125 ! complex(dpc), allocatable :: coeffW_BZ(:,:,:),coeffW_IBZ(:,:,:) 126 complex(dpc), allocatable :: rhot_q_m(:,:,:,:,:,:) 127 complex(dpc), allocatable :: rhot_q_m_npwe(:,:,:,:,:,:),trrho(:,:),sumrhorhoeps(:) 128 complex(gwpc), allocatable :: scr(:,:,:,:) 129 130 real(dp),allocatable :: k_coord(:,:)!,k_coordIBZ(:,:) 131 real(dp),allocatable :: q_coord(:,:) 132 real(dp),allocatable:: normG(:) 133 complex(dpc),allocatable:: uomega(:),jomega(:) 134 real(dp),allocatable:: omega(:) 135 136 integer,allocatable:: ikmq_bz_t(:,:) 137 138 logical,allocatable :: bijection(:) 139 !************************************************************************ 140 141 write(message,*) ch10, '==== Calculation of the screened interaction ====' 142 call wrtout(std_out,message,'COLL') 143 call wrtout(ab_out,message,'COLL') 144 write(message,*) "" 145 call wrtout(std_out,message,'COLL') 146 call wrtout(ab_out,message,'COLL') 147 nkbz = Kmesh%nbz 148 nqibz= Qmesh%nibz 149 nspinor=Wfd%nspinor 150 151 nbband=1+bandsup-bandinf 152 ! _ __ ____ 153 ! | |/ / ___ / __ \ 154 ! | ' / ( _ ) | | | | 155 ! | < / _ \/\ | | | | 156 ! | . \ | (_> < | |__| | 157 ! |_|\_\ \___/\/ \___\_\ 158 159 write(message,*) "Read K and Q mesh" 160 call wrtout(std_out,message,'COLL') 161 call wrtout(ab_out,message,'COLL') 162 163 ABI_ALLOCATE(k_coord,(nkbz,3)) 164 ABI_ALLOCATE(q_coord,(nqibz,4)) 165 166 167 !==Read k and q==! 168 !open(unit=2012,file='ikbz_COORD',form='formatted',status='unknown') 169 !read(2012,*) (ik_bz,k_coord(ik_bz,:),i=1,nkbz) 170 !close(2012) 171 172 do ik_bz=1,nkbz 173 call get_BZ_item(Kmesh,ik_bz,k_coord(ik_bz,:),ik_ibz,isym_kgw,iik,ph_mkt) 174 end do 175 176 ! open(unit=2012,file='iqbz_COORD',form='formatted',status='unknown') 177 ! read(2012,*) 178 do i=1,nqibz 179 ! read(2012,*) iq_ibz,q_coord(iq_ibz,:) 180 q_coord(i,1)=Qmesh%ibz(1,i) 181 q_coord(i,2)=Qmesh%ibz(2,i) 182 q_coord(i,3)=Qmesh%ibz(3,i) 183 q_coord(i,4)=Qmesh%wt(i) 184 ! if (iq_ibz > nqibz) then 185 ! write(message,*) iq_ibz,nqibz," Error on line",i,"Are you in iBZ ?" 186 ! call wrtout(std_out,message,'COLL') 187 ! end if 188 end do 189 ! close(2012) 190 191 !==Bijection and array for k-q==! 192 ABI_ALLOCATE(bijection,(nkbz)) 193 ABI_ALLOCATE(ikmq_bz_t,(nkbz,nqibz)) 194 bijection(:)=.FALSE. 195 do ik_bz=1,nkbz 196 do iq_ibz=1,nqibz 197 ikmq_bz_t(ik_bz,iq_ibz)=findkmq(ik_bz,k_coord,q_coord(iq_ibz,:),nkbz) 198 if (ikmq_bz_t(ik_bz,iq_ibz)>nkbz.and.nsym==1) then 199 BUG=.TRUE. 200 write(message,*) "No K-Q for K/Q =",ik_bz,iq_ibz 201 MSG_ERROR(message) 202 end if 203 bijection(ikmq_bz_t(ik_bz,iq_ibz))=.TRUE. 204 end do 205 206 if (count(bijection).NE.nqibz.and.nsym==1) then 207 BUG=.TRUE. 208 write(message,*) 'No bijection ',ik_bz 209 MSG_ERROR(message) 210 end if 211 212 bijection(:)=.FALSE. 213 end do 214 215 if (.NOT.BUG.and.nsym==1) then 216 write(message,*) "Bijection Ok." 217 call wrtout(std_out,message,'COLL') 218 end if 219 ! _____ 220 ! / ____| 221 ! _ __ ___ _ __ _ __ ___ ___ | | __ 222 ! | '_ \ / _ \| '__| '_ ` _ \ / _ \ | | |_ | 223 ! | | | | (_) | | | | | | | | __/ | |__| | 224 ! |_| |_|\___/|_| |_| |_| |_|\___| \_____| 225 226 ABI_ALLOCATE(normG,(npw)) 227 if (verbose) then 228 write(message,*) 'Read the potential and G norm' 229 call wrtout(std_out,message,'COLL') 230 if (open_file('normeG',message,newunit=unt,form='formatted',status='unknown') /= 0) then 231 MSG_ERROR(message) 232 end if 233 read(unt,*) (iiG,x,normG(iiG),iG=1,npw) 234 close(unt) 235 !!False norme for G=0 idd G is the inverse of the potential inverse du potentiel (q=0) 236 normG(1)=0 237 end if 238 239 !======================================================================== 240 !------------------------------------------------------------------------ 241 242 ! FIRST PART OF THE ROUTINE: USE M_G^(nn')(q,k) to do formal checks. 243 ! USE rhot_q_n to do compute bare interaction 244 245 !------------------------------------------------------------------------ 246 !======================================================================== 247 248 ! _ 249 ! ( ) 250 ! _ __ ___ _ __ _ __ |/ 251 ! | '_ ` _ \ | '_ \| '_ \ 252 ! | | | | | | | | | | | | | 253 ! |_| |_| |_| |_| |_|_| |_| 254 255 !========================================================== 256 !========================================================== 257 ! tol=1E-1 258 ! tolerance for the normalization of wfc: should be around 0.01. 259 tol = 1 ! very large for test. 260 write(message,*) 'Check the norm of M' 261 call wrtout(std_out,message,'COLL') 262 write(message,*) 'Tolerance :',tol 263 call wrtout(std_out,message,'COLL') 264 ! __ __ 265 ! \ \ / / 266 ! \ \ / / _ __ 267 ! \ \/ / | '_ \ 268 ! \ / | | | | 269 ! \/ |_| |_| 270 !========================================================================== 271 !==Compute V_{n,n'}: bare interaction in the KS basis <- rhot_q_n -> V_n 272 !========================================================================== 273 !========================================================== 274 !==Compute V_{n,n'} 275 !========================================================== 276 if(verbose) then 277 write(message,*) "" 278 call wrtout(std_out,message,'COLL') 279 call wrtout(ab_out,message,'COLL') 280 write(message,*) "==Calcul of the bare kohn-sham interaction V n==" 281 call wrtout(std_out,message,'COLL') 282 endif 283 tol=1E+1 284 285 if (.NOT.wannier) RETURN 286 287 !======================================================================== 288 !------------------------------------------------------------------------ 289 290 ! SECOND PART OF THE ROUTINE: Read Wannier projections 291 292 !------------------------------------------------------------------------ 293 !======================================================================== 294 ! 295 ! \ \ / / (_) 296 ! \ \ /\ / /_ _ _ __ _ __ _ ___ _ __ 297 ! \ \/ \/ / _` | '_ \| '_ \| |/ _ \ '__| 298 ! \ /\ / (_| | | | | | | | | __/ | 299 ! \/ \/ \__,_|_| |_|_| |_|_|\___|_| 300 301 !========================================================== 302 !========================================================== 303 !== Read Wannier coefficient in forlb.ovlp 304 !========================================================== 305 write(message,*) "" 306 call wrtout(std_out,message,'COLL') 307 call wrtout(ab_out,message,'COLL') 308 write(message,*) "Read wannier in iBZ" 309 call wrtout(ab_out,message,'COLL') 310 call wrtout(std_out,message,'COLL') 311 nkibz=Kmesh%nibz 312 313 !Read "l" 314 if (open_file('forlb.ovlp',message,newunit=unt,form='formatted',status='unknown') /= 0) then 315 MSG_ERROR(message) 316 end if 317 rewind(unt) 318 read(unt,*) message, lpawu_read 319 read(unt,*) message, ib1, ib2 320 close(unt) 321 if(ib1/=bandinf.and.ib2/=bandsup) then 322 write(message,*) "Error with bands",ib1,bandinf,ib2,bandsup 323 MSG_ERROR(message) 324 endif 325 !!Read the bandinf, bandinf redondance information 326 mbband=2*lpawu_read+1 327 328 !******************************************************* 329 !if (3==4) then 330 !! USELESS START