TABLE OF CONTENTS
ABINIT/wvl_initro [ Functions ]
NAME
wvl_initro
FUNCTION
FIXME: add description.
COPYRIGHT
Copyright (C) 2012-2018 ABINIT group (TRangel) 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
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
PARENTS
gstate
CHILDREN
ext_buffers,ind_positions,metric,mkdenpos,pawrad_free,pawrad_init sort_dp,splint,wrtout,xred2xcart
SOURCE
34 #if defined HAVE_CONFIG_H 35 #include "config.h" 36 #endif 37 38 39 #include "abi_common.h" 40 41 subroutine wvl_initro(& 42 & atindx1,geocode,h,me,& 43 & natom,nattyp,nfft,nspden,ntypat,& 44 & n1,n1i,n2,n2i,n3,& 45 & pawrad,pawtab,psppar,& 46 & rhor,rprimd,spinat,wvl_den,xc_denpos,xred,zion) 47 48 use defs_basis 49 use m_profiling_abi 50 use m_splines 51 use m_errors 52 use defs_wvltypes 53 use m_sort 54 55 use m_geometry, only : xred2xcart, metric 56 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 57 use m_pawtab, only : pawtab_type 58 #if defined HAVE_BIGDFT 59 use BigDFT_API, only : ELECTRONIC_DENSITY, ext_buffers, ind_positions 60 #endif 61 62 !This section has been created automatically by the script Abilint (TD). 63 !Do not modify the following lines by hand. 64 #undef ABI_FUNC 65 #define ABI_FUNC 'wvl_initro' 66 use interfaces_14_hidewrite 67 use interfaces_41_xc_lowlevel 68 !End of the abilint section 69 70 implicit none 71 72 !Arguments ------------------------------------ 73 integer,intent(in) :: me,natom,ntypat,nfft,nspden 74 integer,intent(in)::n1,n2,n1i,n2i,n3 75 real(dp),intent(in) :: h(3) 76 type(pawrad_type),intent(in) :: pawrad(ntypat) 77 type(pawtab_type),intent(in) :: pawtab(ntypat) 78 real(dp),intent(in) :: spinat(3,natom),zion(ntypat) 79 real(dp),intent(inout) :: rhor(nfft,nspden) 80 real(dp),intent(in) :: xc_denpos 81 character(1),intent(in)::geocode 82 type(wvl_denspot_type), intent(inout) :: wvl_den 83 !arrays 84 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 85 real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3) 86 real(dp),intent(inout)::xred(3,natom) 87 88 !Local variables------------------------------- 89 #if defined HAVE_BIGDFT 90 integer :: ia1,ia2 91 integer :: iat,iatm,iatom,iatom_tot,iex,iey,iez,ii,ind 92 integer :: ifft,ispden 93 integer :: isx,isy,isz,itypat,i1,i2,i3,iwarn,i3s 94 integer :: j1,j2,j3,msz 95 integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3 96 integer :: ncmax,nfgd,nspden_updn,n3pi,shift 97 real(dp) :: cutoff,fact,fact0 98 real(dp) :: rloc,rr2,rx,ry,rz 99 real(dp) :: rshp,r2shp 100 real(dp) :: ucvol,xx,yy,zz 101 type(pawrad_type)::vale_mesh 102 !arrays 103 logical :: perx,pery,perz,gox,goy,goz 104 real(dp) :: hh(3) !fine grid spacing for wavelets 105 real(dp) :: gmet(3,3),gprimd(3,3),rcart(3),rmet(3,3),xcart(3,natom) 106 character(len=500) :: message ! to be uncommented, if needed 107 !allocatable arrays 108 integer,allocatable :: ifftsph_tmp(:),iindex(:) 109 real(dp),allocatable:: raux(:),raux2(:) 110 real(dp),allocatable:: rr(:)!,rred(:,:) 111 #endif 112 113 ! ************************************************************************* 114 115 DBG_ENTER("COLL") 116 117 #if defined HAVE_BIGDFT 118 119 !PENDING: PARALLELIZATION OVER ATOMS 120 121 write(message,'(a,a)') ch10,& 122 & ' wvl_initro: Initialize valence density from atomic data by splines' 123 call wrtout(std_out,message,'COLL') 124 125 !initialize 126 rhor(:,:)=zero 127 128 if(nspden==4)then 129 write(std_out,*)' initro : might work yet for nspden=4 (not checked)' 130 write(std_out,*)'spinat',spinat(1:3,1:natom) 131 ! stop 132 end if 133 134 !Check whether the values of spinat are acceptable 135 if(nspden==2)then 136 do itypat=1,ntypat 137 do iat=1,nattyp(itypat) 138 iatm=iatm+1;iatom=atindx1(iatm) 139 iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom) 140 141 if( sqrt(spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2) & 142 & > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then 143 write(message, '(a,a,a,a,i4,a,a,3es11.4,a,a,a,es11.4)' ) ch10,& 144 & ' initro : WARNING - ',ch10,& 145 & ' For atom number ',iatom,ch10,& 146 & ' input spinat=',spinat(:,iatom),' is larger, in magnitude,',ch10,& 147 & ' than zion(ia)=',zion(itypat) 148 call wrtout(std_out,message,'COLL') 149 call wrtout(ab_out,message,'COLL') 150 end if 151 end do 152 ia1=ia2+1 153 end do 154 end if 155 156 !Fine grid 157 hh(:)=0.5d0*h(:) 158 159 !mpi: 160 !Obtain n3pi, BigDFT quantity: 161 n3pi=wvl_den%denspot%dpbox%n3pi 162 i3s=wvl_den%denspot%dpbox%nscatterarr(me,3)+1-wvl_den%denspot%dpbox%nscatterarr(me,4) 163 shift=n1i*n2i*wvl_den%denspot%dpbox%nscatterarr(me,4) 164 165 !Compute xcart from xred 166 call xred2xcart(natom,rprimd,xcart,xred) 167 168 !Compute metric tensors and ucvol from rprimd 169 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 170 171 !Conditions for periodicity in the three directions 172 perx=(geocode /= 'F') 173 pery=(geocode == 'P') 174 perz=(geocode /= 'F') 175 176 !Compute values of external buffers 177 call ext_buffers(perx,nbl1,nbr1) 178 call ext_buffers(pery,nbl2,nbr2) 179 call ext_buffers(perz,nbl3,nbr3) 180 181 iatm=0 182 !Big loop on atom types 183 do itypat=1,ntypat 184 ! 185 rloc=psppar(0,0,itypat) 186 cutoff=10.d0*rloc 187 188 ! Create mesh_core object 189 ! since tnvale_mesh_size can be bigger than pawrad%mesh_size, 190 msz=pawtab(itypat)%tnvale_mesh_size 191 call pawrad_init(vale_mesh,mesh_size=msz,mesh_type=pawrad(itypat)%mesh_type,& 192 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 193 ! 194 ! Set radius size: 195 rshp=vale_mesh%rmax 196 r2shp=1.0000001_dp*rshp**2 197 198 ! allocate arrays 199 if (n3pi > 0) then 200 ! sphere: cycle i1,i2,i3 201 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 202 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 203 ! 1+int(1.1* factors are included just for cautioness 204 ! circle: cycle only i1 and i2 205 ! ncmax=1+int(1.1d0*((rshp/hh(1))*(rshp/hh(2))*pi)) 206 ! line: 207 ncmax=1+int(1.1_dp*rshp/hh(1)*2.d0) 208 else 209 ncmax=1 210 end if 211 ! 212 ABI_ALLOCATE(ifftsph_tmp,(ncmax)) 213 ABI_ALLOCATE(iindex,(ncmax)) 214 ABI_ALLOCATE(rr,(ncmax)) 215 ABI_ALLOCATE(raux,(ncmax)) 216 if(nspden==2) then 217 ABI_ALLOCATE(raux2,(ncmax)) 218 end if 219 220 ! Big loop on atoms 221 do iat=1,nattyp(itypat) 222 iatm=iatm+1;iatom=atindx1(iatm) 223 iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom) 224 225 ! Spin 226 if(nspden==2) then 227 fact0=half/zion(itypat) 228 fact=fact0*(zion(itypat)+spinat(3,iatom)) 229 end if 230 231 ! 232 ! Define a "box" around each atom 233 rx=xcart(1,iatom_tot) 234 ry=xcart(2,iatom_tot) 235 rz=xcart(3,iatom_tot) 236 ! 237 isx=floor((rx-cutoff)/hh(1)) 238 isy=floor((ry-cutoff)/hh(2)) 239 isz=floor((rz-cutoff)/hh(3)) 240 241 iex=ceiling((rx+cutoff)/hh(1)) 242 iey=ceiling((ry+cutoff)/hh(2)) 243 iez=ceiling((rz+cutoff)/hh(3)) 244 ! 245 do i3=isz,iez 246 zz=real(i3,kind=8)*hh(3)-rz 247 call ind_positions(perz,i3,n3,j3,goz) 248 j3=j3+nbl3+1 249 ! 250 do i2=isy,iey 251 yy=real(i2,kind=8)*hh(2)-ry 252 call ind_positions(pery,i2,n2,j2,goy) 253 ! 254 ! Initialize counters 255 nfgd=0 256 ! nfgd_r0=0 257 ! 258 do i1=isx,iex 259 xx=real(i1,kind=8)*hh(1)-rx 260 call ind_positions(perx,i1,n1,j1,gox) 261 rr2=xx**2+yy**2+zz**2 262 if (j3 >= i3s .and. j3 <= i3s+n3pi-1 .and. goy .and. gox ) then 263 ! 264 if(rr2<=r2shp) then 265 if(rr2>tol5) then 266 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i 267 nfgd=nfgd+1 268 rcart=[xx,yy,zz] 269 rr(nfgd)=(rr2)**0.5 270 ifftsph_tmp(nfgd)=shift+ind 271 ! DEBUG 272 ! write(itmp,'(i10,3(f13.7,x))')ind,xx+rx,yy+ry,zz+rz 273 ! write(itmp,'(6(f13.7,x))')rcart,rred(:,nfgd) 274 ! ENDDEBUG 275 ! else 276 ! ! We save r=0 vectors 277 ! ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i 278 ! ! We reuse the same variable "ifftshp_tmp", 279 ! ! but we start from the higher index 280 ! nfgd_r0=nfgd_r0+1 281 ! ifftsph_tmp(ncmax-nfgd_r0+1)=shift+ind 282 end if !rr2>tol5 283 end if !rr2<r2shp 284 end if !j3.. 285 end do !i1 286 287 ! All of the following could be done inside or outside the loops (i2,i1,i3) 288 ! Outside the loops: the memory consuption increases. 289 ! Inside the inner loop: the time of calculation increases. 290 291 if(nfgd==0) cycle 292 293 ! Evaluate spline fit of 1st der of core charge density 294 ! from tcoredens(:,2) and tcoredens(:,4) 295 do ii=1,nfgd 296 iindex(ii)=ii 297 end do 298 ! write(600,'(i4,x,9999f14.7)')nfgd, rr(1:nfgd) 299 call sort_dp(nfgd,rr(1:nfgd),iindex(1:nfgd),tol16) 300 call splint(msz,vale_mesh%rad,& 301 & pawtab(itypat)%tvalespl(:,1),pawtab(itypat)%tvalespl(:,2),& 302 & nfgd,rr(1:nfgd),raux(1:nfgd)) 303 304 305 ! Accumulate contributions to valence density on the entire cell 306 rhor(ifftsph_tmp(1:nfgd),1)=rhor(ifftsph_tmp(1:nfgd),1)+raux(iindex(1:nfgd)) 307 308 if(nspden==2) then 309 raux2(1:nfgd)=raux(iindex(1:nfgd))*fact 310 rhor(ifftsph_tmp(1:nfgd),2)=rhor(ifftsph_tmp(1:nfgd),1)+raux2(1:nfgd) 311 end if 312 ! DEBUG 313 ! do ii=1,msz 314 ! write(itmp,'(2(f15.7,1x))')vale_mesh%rad(ii),pawtab(itypat)%tvalespl(ii,1) 315 ! end do 316 ! do ii=1,nfgd 317 ! write(itmp,'(2i10)')ii,iindex(ii) 318 ! write(itmp,'(2(f15.7,1x))')rr(iindex(ii)),rhor(ifftsph_tmp(ii),1)!,raux(iindex(ii)) 319 ! end do 320 ! END DEBUG 321 322 end do !i2 323 end do !i1 324 end do !iat 325 326 ! Deallocate 327 call pawrad_free(vale_mesh) 328 ABI_DEALLOCATE(ifftsph_tmp) 329 ABI_DEALLOCATE(iindex) 330 ABI_DEALLOCATE(rr) 331 ABI_DEALLOCATE(raux) 332 if(nspden==2) then 333 ABI_DEALLOCATE(raux2) 334 end if 335 336 end do !itypat 337 338 !nspden_updn: 1 for non-polarized, 2 for polarized 339 nspden_updn=min(nspden,2) 340 341 !Make the density positive everywhere 342 call mkdenpos(iwarn,nfft,nspden_updn,1,rhor(:,1:nspden_updn),xc_denpos) 343 344 !There seems to be a bug in the intel11 compiler 345 !rhor = reshape(wvl_den%denspot%rhov, shape(rhor)) 346 do ispden=1,nspden 347 do ifft=1,nfft 348 ii=ifft+nfft*(ispden-1) 349 ! rhor(ifft,ispden)=wvl_den%denspot%rhov(ii) 350 wvl_den%denspot%rhov(ii)=rhor(ifft,ispden) 351 end do 352 end do 353 wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY 354 write(message, '(a,a,a,a)' ) ch10, ' wvl_initro : but why are you copying me :..o(' 355 call wrtout(std_out,message,'COLL') 356 357 358 #else 359 BIGDFT_NOTENABLED_ERROR() 360 if (.false.) write(std_out,*) me,natom,ntypat,nfft,nspden,n1,n2,n1i,n2i,n3,h(1),& 361 & pawrad(1)%mesh_size,pawtab(1)%mesh_size,spinat(1,1),zion(1),rhor(1,1),xc_denpos,& 362 & geocode,wvl_den%symObj,atindx1(1),nattyp(1),psppar(1,1,1),rprimd(1,1),xred(1,1) 363 #endif 364 365 DBG_EXIT("COLL") 366 367 end subroutine wvl_initro