TABLE OF CONTENTS
ABINIT/m_wvl_rho [ Modules ]
NAME
m_wvl_rho
FUNCTION
COPYRIGHT
Copyright (C) 2012-2024 ABINIT group (TRangel, DC) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_wvl_rho 23 24 use defs_basis 25 use m_abicore 26 use m_splines 27 use m_errors 28 use defs_wvltypes 29 use m_sort 30 use m_abi2big 31 use m_xmpi 32 use m_dtset 33 34 use defs_abitypes, only : MPI_type 35 use m_geometry, only : xred2xcart, metric 36 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_free 37 use m_pawtab, only : pawtab_type 38 use m_pawrhoij, only : pawrhoij_type 39 use m_drivexc, only : mkdenpos 40 41 implicit none 42 43 private
ABINIT/wvl_initro [ Functions ]
NAME
wvl_initro
FUNCTION
FIXME: add description.
INPUTS
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
SOURCE
73 subroutine wvl_initro(& 74 & atindx1,geocode,h,me,& 75 & natom,nattyp,nfft,nspden,ntypat,& 76 & n1,n1i,n2,n2i,n3,& 77 & pawrad,pawtab,psppar,& 78 & rhor,rprimd,spinat,wvl_den,xc_denpos,xred,zion) 79 80 #if defined HAVE_BIGDFT 81 use BigDFT_API, only : ELECTRONIC_DENSITY, ext_buffers, ind_positions 82 #endif 83 84 !Arguments ------------------------------------ 85 integer,intent(in) :: me,natom,ntypat,nfft,nspden 86 integer,intent(in)::n1,n2,n1i,n2i,n3 87 real(dp),intent(in) :: h(3) 88 type(pawrad_type),intent(in) :: pawrad(ntypat) 89 type(pawtab_type),intent(in) :: pawtab(ntypat) 90 real(dp),intent(in) :: spinat(3,natom),zion(ntypat) 91 real(dp),intent(inout) :: rhor(nfft,nspden) 92 real(dp),intent(in) :: xc_denpos 93 character(1),intent(in)::geocode 94 type(wvl_denspot_type), intent(inout) :: wvl_den 95 !arrays 96 integer,intent(in) :: atindx1(natom),nattyp(ntypat) 97 real(dp),intent(in) :: psppar(0:4,0:6,ntypat),rprimd(3,3) 98 real(dp),intent(inout)::xred(3,natom) 99 100 !Local variables------------------------------- 101 #if defined HAVE_BIGDFT 102 integer :: ia1,ia2 103 integer :: iat,iatm,iatom,iatom_tot,iex,iey,iez,ii,ind 104 integer :: ifft,ispden 105 integer :: isx,isy,isz,itypat,i1,i2,i3,iwarn,i3s 106 integer :: j1,j2,j3,msz 107 integer :: nbl1,nbr1,nbl2,nbr2,nbl3,nbr3 108 integer :: ncmax,nfgd,nspden_updn,n3pi,shift 109 real(dp) :: cutoff,fact,fact0 110 real(dp) :: rloc,rr2,rx,ry,rz 111 real(dp) :: rshp,r2shp 112 real(dp) :: ucvol,xx,yy,zz 113 type(pawrad_type)::vale_mesh 114 !arrays 115 logical :: perx,pery,perz,gox,goy,goz 116 real(dp) :: hh(3) !fine grid spacing for wavelets 117 real(dp) :: gmet(3,3),gprimd(3,3),rcart(3),rmet(3,3),xcart(3,natom) 118 character(len=500) :: message ! to be uncommented, if needed 119 !allocatable arrays 120 integer,allocatable :: ifftsph_tmp(:),iindex(:) 121 real(dp),allocatable:: raux(:),raux2(:) 122 real(dp),allocatable:: rr(:)!,rred(:,:) 123 #endif 124 125 ! ************************************************************************* 126 127 DBG_ENTER("COLL") 128 129 #if defined HAVE_BIGDFT 130 131 !PENDING: PARALLELIZATION OVER ATOMS 132 133 write(message,'(a,a)') ch10,& 134 & ' wvl_initro: Initialize valence density from atomic data by splines' 135 call wrtout(std_out,message,'COLL') 136 137 !initialize 138 rhor(:,:)=zero 139 140 if(nspden==4)then 141 write(std_out,*)' initro : might work yet for nspden=4 (not checked)' 142 write(std_out,*)'spinat',spinat(1:3,1:natom) 143 ! stop 144 end if 145 146 !Check whether the values of spinat are acceptable 147 if(nspden==2)then 148 do itypat=1,ntypat 149 do iat=1,nattyp(itypat) 150 iatm=iatm+1;iatom=atindx1(iatm) 151 iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom) 152 153 if( sqrt(spinat(1,iatom)**2+spinat(2,iatom)**2+spinat(3,iatom)**2) & 154 & > abs(zion(itypat))*(1.0_dp + epsilon(0.0_dp)) ) then 155 write(message, '(a,a,a,a,i4,a,a,3es11.4,a,a,a,es11.4)' ) ch10,& 156 & ' initro : WARNING - ',ch10,& 157 & ' For atom number ',iatom,ch10,& 158 & ' input spinat=',spinat(:,iatom),' is larger, in magnitude,',ch10,& 159 & ' than zion(ia)=',zion(itypat) 160 call wrtout(std_out,message,'COLL') 161 call wrtout(ab_out,message,'COLL') 162 end if 163 end do 164 ia1=ia2+1 165 end do 166 end if 167 168 !Fine grid 169 hh(:)=0.5d0*h(:) 170 171 !mpi: 172 !Obtain n3pi, BigDFT quantity: 173 n3pi=wvl_den%denspot%dpbox%n3pi 174 i3s=wvl_den%denspot%dpbox%nscatterarr(me,3)+1-wvl_den%denspot%dpbox%nscatterarr(me,4) 175 shift=n1i*n2i*wvl_den%denspot%dpbox%nscatterarr(me,4) 176 177 !Compute xcart from xred 178 call xred2xcart(natom,rprimd,xcart,xred) 179 180 !Compute metric tensors and ucvol from rprimd 181 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 182 183 !Conditions for periodicity in the three directions 184 perx=(geocode /= 'F') 185 pery=(geocode == 'P') 186 perz=(geocode /= 'F') 187 188 !Compute values of external buffers 189 call ext_buffers(perx,nbl1,nbr1) 190 call ext_buffers(pery,nbl2,nbr2) 191 call ext_buffers(perz,nbl3,nbr3) 192 193 iatm=0 194 !Big loop on atom types 195 do itypat=1,ntypat 196 ! 197 rloc=psppar(0,0,itypat) 198 cutoff=10.d0*rloc 199 200 ! Create mesh_core object 201 ! since tnvale_mesh_size can be bigger than pawrad%mesh_size, 202 msz=pawtab(itypat)%tnvale_mesh_size 203 call pawrad_init(vale_mesh,mesh_size=msz,mesh_type=pawrad(itypat)%mesh_type,& 204 & rstep=pawrad(itypat)%rstep,lstep=pawrad(itypat)%lstep) 205 ! 206 ! Set radius size: 207 rshp=vale_mesh%rmax 208 r2shp=1.0000001_dp*rshp**2 209 210 ! allocate arrays 211 if (n3pi > 0) then 212 ! sphere: cycle i1,i2,i3 213 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 214 ! ncmax=1+int(1.1_dp*nfft*four_pi/(three*ucvol)*rshp**3) 215 ! 1+int(1.1* factors are included just for cautioness 216 ! circle: cycle only i1 and i2 217 ! ncmax=1+int(1.1d0*((rshp/hh(1))*(rshp/hh(2))*pi)) 218 ! line: 219 ncmax=1+int(1.1_dp*rshp/hh(1)*2.d0) 220 else 221 ncmax=1 222 end if 223 ! 224 ABI_MALLOC(ifftsph_tmp,(ncmax)) 225 ABI_MALLOC(iindex,(ncmax)) 226 ABI_MALLOC(rr,(ncmax)) 227 ABI_MALLOC(raux,(ncmax)) 228 if(nspden==2) then 229 ABI_MALLOC(raux2,(ncmax)) 230 end if 231 232 ! Big loop on atoms 233 do iat=1,nattyp(itypat) 234 iatm=iatm+1;iatom=atindx1(iatm) 235 iatom_tot=iatom; !if (mpi_enreg%nproc_atom>1) iatom_tot=mpi_enreg%atom_indx(iatom) 236 237 ! Spin 238 if(nspden==2) then 239 fact0=half/zion(itypat) 240 fact=fact0*(zion(itypat)+spinat(3,iatom)) 241 end if 242 243 ! 244 ! Define a "box" around each atom 245 rx=xcart(1,iatom_tot) 246 ry=xcart(2,iatom_tot) 247 rz=xcart(3,iatom_tot) 248 ! 249 isx=floor((rx-cutoff)/hh(1)) 250 isy=floor((ry-cutoff)/hh(2)) 251 isz=floor((rz-cutoff)/hh(3)) 252 253 iex=ceiling((rx+cutoff)/hh(1)) 254 iey=ceiling((ry+cutoff)/hh(2)) 255 iez=ceiling((rz+cutoff)/hh(3)) 256 ! 257 do i3=isz,iez 258 zz=real(i3,kind=8)*hh(3)-rz 259 call ind_positions(perz,i3,n3,j3,goz) 260 j3=j3+nbl3+1 261 ! 262 do i2=isy,iey 263 yy=real(i2,kind=8)*hh(2)-ry 264 call ind_positions(pery,i2,n2,j2,goy) 265 ! 266 ! Initialize counters 267 nfgd=0 268 ! nfgd_r0=0 269 ! 270 do i1=isx,iex 271 xx=real(i1,kind=8)*hh(1)-rx 272 call ind_positions(perx,i1,n1,j1,gox) 273 rr2=xx**2+yy**2+zz**2 274 if (j3 >= i3s .and. j3 <= i3s+n3pi-1 .and. goy .and. gox ) then 275 ! 276 if(rr2<=r2shp) then 277 if(rr2>tol5) then 278 ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i 279 nfgd=nfgd+1 280 rcart=[xx,yy,zz] 281 rr(nfgd)=(rr2)**0.5 282 ifftsph_tmp(nfgd)=shift+ind 283 ! DEBUG 284 ! write(itmp,'(i10,3(f13.7,x))')ind,xx+rx,yy+ry,zz+rz 285 ! write(itmp,'(6(f13.7,x))')rcart,rred(:,nfgd) 286 ! ENDDEBUG 287 ! else 288 ! ! We save r=0 vectors 289 ! ind=j1+1+nbl1+(j2+nbl2)*n1i+(j3-i3s)*n1i*n2i 290 ! ! We reuse the same variable "ifftshp_tmp", 291 ! ! but we start from the higher index 292 ! nfgd_r0=nfgd_r0+1 293 ! ifftsph_tmp(ncmax-nfgd_r0+1)=shift+ind 294 end if !rr2>tol5 295 end if !rr2<r2shp 296 end if !j3.. 297 end do !i1 298 299 ! All of the following could be done inside or outside the loops (i2,i1,i3) 300 ! Outside the loops: the memory consuption increases. 301 ! Inside the inner loop: the time of calculation increases. 302 303 if(nfgd==0) cycle 304 305 ! Evaluate spline fit of 1st der of core charge density 306 ! from tcoredens(:,2) and tcoredens(:,4) 307 do ii=1,nfgd 308 iindex(ii)=ii 309 end do 310 ! write(600,'(i4,x,9999f14.7)')nfgd, rr(1:nfgd) 311 call sort_dp(nfgd,rr(1:nfgd),iindex(1:nfgd),tol16) 312 call splint(msz,vale_mesh%rad,& 313 & pawtab(itypat)%tvalespl(:,1),pawtab(itypat)%tvalespl(:,2),& 314 & nfgd,rr(1:nfgd),raux(1:nfgd)) 315 316 317 ! Accumulate contributions to valence density on the entire cell 318 rhor(ifftsph_tmp(1:nfgd),1)=rhor(ifftsph_tmp(1:nfgd),1)+raux(iindex(1:nfgd)) 319 320 if(nspden==2) then 321 raux2(1:nfgd)=raux(iindex(1:nfgd))*fact 322 rhor(ifftsph_tmp(1:nfgd),2)=rhor(ifftsph_tmp(1:nfgd),1)+raux2(1:nfgd) 323 end if 324 ! DEBUG 325 ! do ii=1,msz 326 ! write(itmp,'(2(f15.7,1x))')vale_mesh%rad(ii),pawtab(itypat)%tvalespl(ii,1) 327 ! end do 328 ! do ii=1,nfgd 329 ! write(itmp,'(2i10)')ii,iindex(ii) 330 ! write(itmp,'(2(f15.7,1x))')rr(iindex(ii)),rhor(ifftsph_tmp(ii),1)!,raux(iindex(ii)) 331 ! end do 332 ! END DEBUG 333 334 end do !i2 335 end do !i1 336 end do !iat 337 338 ! Deallocate 339 call pawrad_free(vale_mesh) 340 ABI_FREE(ifftsph_tmp) 341 ABI_FREE(iindex) 342 ABI_FREE(rr) 343 ABI_FREE(raux) 344 if(nspden==2) then 345 ABI_FREE(raux2) 346 end if 347 348 end do !itypat 349 350 !nspden_updn: 1 for non-polarized, 2 for polarized 351 nspden_updn=min(nspden,2) 352 353 !Make the density positive everywhere 354 call mkdenpos(iwarn,nfft,nspden_updn,1,rhor(:,1:nspden_updn),xc_denpos) 355 356 !There seems to be a bug in the intel11 compiler 357 !rhor = reshape(wvl_den%denspot%rhov, shape(rhor)) 358 do ispden=1,nspden 359 do ifft=1,nfft 360 ii=ifft+nfft*(ispden-1) 361 ! rhor(ifft,ispden)=wvl_den%denspot%rhov(ii) 362 wvl_den%denspot%rhov(ii)=rhor(ifft,ispden) 363 end do 364 end do 365 wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY 366 write(message, '(a,a,a,a)' ) ch10, ' wvl_initro : but why are you copying me :..o(' 367 call wrtout(std_out,message,'COLL') 368 369 370 #else 371 BIGDFT_NOTENABLED_ERROR() 372 if (.false.) write(std_out,*) me,natom,ntypat,nfft,nspden,n1,n2,n1i,n2i,n3,h(1),& 373 & pawrad(1)%mesh_size,pawtab(1)%mesh_size,spinat(1,1),zion(1),rhor(1,1),xc_denpos,& 374 & geocode,wvl_den%symObj,atindx1(1),nattyp(1),psppar(1,1,1),rprimd(1,1),xred(1,1) 375 #endif 376 377 DBG_EXIT("COLL") 378 379 end subroutine wvl_initro
ABINIT/wvl_mkrho [ Functions ]
NAME
wvl_mkrho
FUNCTION
This method is just a wrapper around the BigDFT routine to compute the density from the wavefunctions.
INPUTS
dtset <type(dataset_type)>=input variables. mpi_enreg=information about MPI parallelization occ(dtset%mband)=occupation numbers. psps <type(pseudopotential_type)>=variables related to pseudopotentials wvl_wfs <type(wvl_projector_type)>=wavefunctions information for wavelets.
OUTPUT
rhor(dtset%nfft)=electron density in r space
SIDE EFFECTS
proj <type(wvl_projector_type)>=projectors information for wavelets. | proj(OUT)=computed projectors.
SOURCE
406 subroutine wvl_mkrho(dtset, irrzon, mpi_enreg, phnons, rhor, wvl_wfs, wvl_den) 407 408 #if defined HAVE_BIGDFT 409 use BigDFT_API, only : sumrho, symmetry_data, ELECTRONIC_DENSITY, communicate_density 410 #endif 411 412 !Arguments ------------------------------- 413 !scalars 414 type(MPI_type),intent(in) :: mpi_enreg 415 type(dataset_type),intent(in) :: dtset 416 type(wvl_wf_type),intent(inout) :: wvl_wfs 417 type(wvl_denspot_type), intent(inout) :: wvl_den 418 !arrays 419 real(dp),intent(inout) :: rhor(dtset%nfft,dtset%nspden) 420 integer, target, intent(in) :: irrzon(dtset%nfft**(1-1/dtset%nsym),2, & 421 & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 422 real(dp), target, intent(in) :: phnons(2,(dtset%ngfft(1)*dtset%ngfft(2)*dtset%ngfft(3))**(1-1/dtset%nsym), & 423 & (dtset%nspden/dtset%nsppol)-3*(dtset%nspden/4)) 424 425 !Local variables------------------------------- 426 #if defined HAVE_BIGDFT 427 !scalars 428 character(len=500) :: message 429 integer :: comm,me,nproc 430 type(symmetry_data) :: sym 431 !for debugging: 432 !integer::ifile,ierr 433 #endif 434 435 ! ************************************************************************* 436 437 DBG_ENTER("COLL") 438 439 #if defined HAVE_BIGDFT 440 comm=mpi_enreg%comm_wvl 441 me=xmpi_comm_rank(comm) 442 nproc=xmpi_comm_size(comm) 443 444 sym%symObj = wvl_den%symObj 445 sym%irrzon => irrzon 446 sym%phnons => phnons 447 448 call sumrho(wvl_den%denspot%dpbox,wvl_wfs%ks%orbs,wvl_wfs%ks%Lzd,& 449 & wvl_wfs%GPU,sym,wvl_den%denspot%rhod,wvl_den%denspot%xc,& 450 & wvl_wfs%ks%psi,wvl_den%denspot%rho_psi) 451 452 call communicate_density(wvl_den%denspot%dpbox,wvl_wfs%ks%orbs%nspin,& 453 & wvl_den%denspot%rhod,wvl_den%denspot%rho_psi,wvl_den%denspot%rhov,.false.) 454 455 wvl_den%denspot%rhov_is = ELECTRONIC_DENSITY 456 write(message, '(a,a,a,a)' ) ch10, ' wvl_mkrho : but why are you copying me :..o(' 457 call wrtout(std_out,message,'COLL') 458 459 call wvl_rho_abi2big(2,rhor,wvl_den) 460 461 #else 462 BIGDFT_NOTENABLED_ERROR() 463 if (.false.) write(std_out,*) mpi_enreg%me,dtset%nstep,wvl_wfs%ks,wvl_den%symObj,& 464 & rhor(1,1),irrzon(1,1,1),phnons(1,1,1) 465 #endif 466 467 DBG_EXIT("COLL") 468 469 end subroutine wvl_mkrho
ABINIT/wvl_prcref [ Functions ]
NAME
wvl_prcref
FUNCTION
FIXME: add description.
INPUTS
argin(sizein)=description
OUTPUT
argout(sizeout)=description
SIDE EFFECTS
NOTES
SOURCE
491 subroutine wvl_prcref(dielar,iprcel,my_natom,nfftprc,npawmix,nspden,pawrhoij,& 492 & rhoijrespc,usepaw,vresid,vrespc) 493 494 !Arguments ------------------------------------ 495 integer , intent(in) :: iprcel,nfftprc,my_natom,npawmix,nspden,usepaw 496 real(dp), intent(in) :: dielar(7) 497 real(dp), intent(in) :: vresid(nfftprc,nspden) 498 real(dp),intent(out) :: rhoijrespc(npawmix) 499 real(dp),intent(out) :: vrespc(nfftprc,nspden) 500 type(pawrhoij_type),intent(inout) :: pawrhoij(my_natom*usepaw) 501 502 !Local variables------------------------------- 503 integer :: iatom,index,ispden,klmn,kmix 504 real(dp):: diemix,diemixmag,mixfac_eff 505 character(len=500) :: message ! to be uncommented, if needed 506 507 ! ************************************************************************* 508 509 DBG_ENTER("COLL") 510 511 !PENDING: 512 !optres==1, or 0, for density or potential mixing. 513 !for potential mixing, we have to average over spins. 514 !check prcref.F90 and moddiel.F90 515 516 517 #if defined HAVE_BIGDFT 518 #endif 519 520 if(iprcel .ne. 0) then 521 write(message, '(a,i3,a,a,a,a)' )& 522 & ' From the calling routine, iprcel=',iprcel,ch10,& 523 & ' For wavelets, the only allowed value is 0.',ch10,& 524 & ' Action : correct your input file.' 525 ABI_ERROR(message) 526 end if 527 528 !call wvl_moddiel !PENDING 529 diemix=dielar(4) 530 !dielng=dielar(2) ; diemac=dielar(3) ; diemix=dielar(4) ; 531 diemixmag=abs(dielar(7)) 532 vrespc(:,1)=diemix*vresid(:,1) 533 if (nspden/=1) vrespc(:,2:nspden)=diemixmag*vresid(:,2:nspden) 534 535 !3) PAW only : precondition the rhoij quantities (augmentation 536 !occupancies) residuals. Use a simple preconditionning 537 !with the same mixing factor as the model dielectric function. 538 539 if (usepaw==1.and.my_natom>0) then 540 ABI_CHECK(pawrhoij(1)%qphase==1,'wvl_prcref: not available with qphase=1!') 541 ! mixfac=dielar(4);mixfacmag=abs(dielar(7)) 542 if (pawrhoij(1)%cplex_rhoij==1) then 543 index=0 544 do iatom=1,my_natom 545 do ispden=1,pawrhoij(iatom)%nspden 546 mixfac_eff=diemix;if (ispden>1) mixfac_eff=diemixmag 547 do kmix=1,pawrhoij(iatom)%lmnmix_sz 548 index=index+1;klmn=pawrhoij(iatom)%kpawmix(kmix) 549 rhoijrespc(index)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn,ispden) 550 end do 551 end do 552 end do 553 else 554 index=-1 555 do iatom=1,my_natom 556 do ispden=1,pawrhoij(iatom)%nspden 557 mixfac_eff=diemix;if (ispden>1) mixfac_eff=diemixmag 558 do kmix=1,pawrhoij(iatom)%lmnmix_sz 559 index=index+2;klmn=2*pawrhoij(iatom)%kpawmix(kmix)-1 560 rhoijrespc(index:index+1)=mixfac_eff*pawrhoij(iatom)%rhoijres(klmn:klmn+1,ispden) 561 end do 562 end do 563 end do 564 end if 565 end if 566 567 568 569 DBG_EXIT("COLL") 570 571 end subroutine wvl_prcref