TABLE OF CONTENTS
ABINIT/psolver_rhohxc [ Functions ]
NAME
psolver_rhohxc
FUNCTION
Given rho(r), compute Hartree potential considering the system as an isolated one. This potential is obtained from the convolution of 1/r and rho(r), treated in Fourier space. This method is a wrapper around Psolver() developped for BigDFT. It can compute the xc energy and potential if required. This computation is built on the drivexc() routine of ABINIT but access it directly from real space. The present routine is a real space counter part to rhotoxc().
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR,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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset mpi_enreg=MPI-parallelisation information. rhor(nfft,nspden)=electron density in real space in electrons/bohr**3
OUTPUT
enhartr=returned Hartree energy (hartree). enxc=returned exchange and correlation energy (hartree). envxc=returned energy of the Vxc potential (hartree). vhartr(nfft)=Hartree potential. vxc(nfft,nspden)=xc potential vxcavg=<Vxc>=unit cell average of Vxc = (1/ucvol) Int [Vxc(r) d^3 r]. NOTE In psolver, with nspden == 2, rhor(:,1) = density up and rhor(:,2) = density down. But in ABINIT (dtset%usewvl != 1) rhor(:,1) = total density and rhor(:,2) = density up . In ABINIT (dtset%usewvl != 1), the same convention is used as in psolver.
PARENTS
energy,rhotov,scfcv,setvtr
CHILDREN
h_potential,mean_fftr,metric,mkdenpos,psolver_kernel,wrtout wvl_rhov_abi2big,xc_potential
SOURCE
51 #if defined HAVE_CONFIG_H 52 #include "config.h" 53 #endif 54 55 #include "abi_common.h" 56 57 58 subroutine psolver_rhohxc(enhartr, enxc, envxc, icoulomb, ixc, & 59 & mpi_enreg, nfft, ngfft, nhat,nhatdim,& 60 & nscforder, nspden, n3xccc, rhor, rprimd,& 61 & usexcnhat,usepaw,usewvl,vhartr, vxc, vxcavg, wvl,wvl_den,wvl_e,& 62 & xccc3d,xclevel,xc_denpos) 63 64 ! use defs_basis,only: std_out,std_out_default 65 use defs_basis 66 use defs_abitypes 67 use defs_wvltypes 68 use m_profiling_abi 69 use m_errors 70 use m_abi2big 71 use m_cgtools 72 73 use m_xmpi, only: xmpi_comm_rank,xmpi_comm_size,xmpi_sum 74 use m_geometry, only : metric 75 76 #if defined HAVE_BIGDFT 77 use BigDFT_API, only : XC_potential,ELECTRONIC_DENSITY,coulomb_operator 78 use poisson_solver, only : H_potential 79 #endif 80 81 !This section has been created automatically by the script Abilint (TD). 82 !Do not modify the following lines by hand. 83 #undef ABI_FUNC 84 #define ABI_FUNC 'psolver_rhohxc' 85 use interfaces_14_hidewrite 86 use interfaces_41_xc_lowlevel 87 use interfaces_62_poisson, except_this_one => psolver_rhohxc 88 !End of the abilint section 89 90 implicit none 91 92 !Arguments ------------------------------------ 93 !scalars 94 integer, intent(in) :: nhatdim,nspden,n3xccc 95 integer, intent(in) :: nfft, icoulomb, ixc, nscforder, usewvl 96 integer,intent(in) :: usexcnhat,usepaw,xclevel 97 real(dp),intent(in) :: rprimd(3,3) 98 real(dp), intent(in) :: xc_denpos 99 real(dp), intent(out) :: enxc, envxc, enhartr, vxcavg 100 type(mpi_type), intent(in) :: mpi_enreg 101 type(wvl_internal_type), intent(in) :: wvl 102 type(wvl_denspot_type), intent(inout) :: wvl_den 103 type(wvl_energy_terms), intent(inout) :: wvl_e 104 !arrays 105 integer, intent(in) :: ngfft(18) 106 real(dp),intent(in) :: xccc3d(n3xccc) 107 real(dp),intent(in) :: nhat(nfft,nspden*nhatdim) 108 real(dp),intent(inout) :: rhor(nfft, nspden) 109 real(dp),intent(out) :: vhartr(nfft) 110 real(dp),intent(out) :: vxc(nfft, nspden) 111 112 !Local variables------------------------------- 113 #if defined HAVE_BIGDFT 114 ! n_c and \hat{n} can be added/rested inside bigdft by passing 115 ! them as pointers (rhocore and rhohat): 116 logical, parameter :: add_n_c_here=.true. !Add n_c here or inside bigdft 117 logical, parameter :: rest_hat_n_here=.true. !Rest \hat{n} here or inside bigdft 118 !scalars 119 integer :: me,nproc,comm 120 integer :: ifft,ispin 121 integer :: iwarn, opt_mkdenpos 122 integer :: nfftot,ngrad 123 integer :: n1i,n2i,n3d,n3i 124 real(dp) :: tmpDown, tmpUp, tmpPot,ucvol,ucvol_local 125 logical :: sumpion,test_nhat,use_psolver=.false. 126 character(len=500) :: message 127 character(len = 1) :: datacode, bndcode 128 !arrays 129 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 130 real(dp) :: hgrid(3) 131 real(dp) :: vxcmean(1) 132 real(dp), pointer :: rhocore(:,:,:,:),rhohat(:,:,:,:) 133 real(dp), pointer :: pot_ion(:,:,:,:),rhonow(:,:) 134 real(dp), dimension(6) :: xcstr 135 type(coulomb_operator) :: kernel 136 #endif 137 138 ! ********************************************************************* 139 140 DBG_ENTER("COLL") 141 142 #if defined HAVE_BIGDFT 143 144 nfftot=PRODUCT(ngfft(1:3)) 145 comm=mpi_enreg%comm_fft 146 if(usewvl==1) comm=mpi_enreg%comm_wvl 147 me=xmpi_comm_rank(comm) 148 nproc=xmpi_comm_size(comm) 149 150 if(n3xccc>0) then 151 if(nfft .ne. n3xccc)then 152 write(message,'(a,a,a,2(i0,1x))')& 153 & 'nfft and n3xccc should be equal,',ch10,& 154 & 'however, nfft and n3xccc=',nfft,n3xccc 155 MSG_BUG(message) 156 end if 157 end if 158 if(nspden==4) then 159 MSG_ERROR('nspden==4 not coded yet') 160 end if 161 162 if (ixc==0) then 163 vxcavg=zero 164 test_nhat=.false. 165 166 ! No xc at all is applied (usually for testing) 167 MSG_WARNING('Note that no xc is applied (ixc=0).') 168 169 else if (ixc/=20) then 170 171 ! ngrad=1 is for LDAs or LSDs, ngrad=2 is for GGAs 172 ngrad=1;if(xclevel==2)ngrad=2 173 ! ixc 31 to 34 are for mgga test purpose only (fake functionals based on LDA but need the gradients too) 174 if(ixc>=31 .and. ixc<=34)ngrad=2 175 ! Test: has a compensation density to be added/substracted (PAW) ? 176 ! test_nhat=((nhatdim==1).and.(usexcnhat==0.or.(ngrad==2.and.nhatgrdim==1))) 177 test_nhat=((nhatdim==1).and.(usexcnhat==0)) 178 end if 179 180 181 !Compute different geometric tensor, as well as ucvol, from rprimd 182 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 183 184 if (icoulomb == 0) then 185 ! The kernel is built with 'P'eriodic boundary counditions. 186 bndcode = 'P' 187 else if (icoulomb == 1) then 188 ! The kernel is built with 'F'ree boundary counditions. 189 bndcode = 'F' 190 else if (icoulomb == 2) then 191 ! The kernel is built with 'S'urface boundary counditions. 192 bndcode = 'S' 193 end if 194 195 !This makes the tests fail? 196 !For NC and n_c=0, call psolver, which uses less memory: 197 !if(usepaw==0 .and. n3xccc==0) use_psolver=.true. 198 199 if(nspden > 2)then 200 write(message, '(a,a,a,i0)' )& 201 & 'Only non-spin-polarised or collinear spin is allowed,',ch10,& 202 & 'while the argument nspden = ', nspden 203 MSG_ERROR(message) 204 end if 205 206 !We do the computation. 207 write(message, "(A,A,A,3I6)") "psolver_rhohxc(): compute potentials (Vhartree and Vxc)...", ch10, & 208 & " | dimension:", ngfft(1:3) 209 call wrtout(std_out, message,'COLL') 210 211 if(usewvl==1) then 212 hgrid=(/wvl_den%denspot%dpbox%hgrids(1),wvl_den%denspot%dpbox%hgrids(2),wvl_den%denspot%dpbox%hgrids(3)/) 213 else 214 hgrid=(/ rprimd(1,1) / ngfft(1), rprimd(2,2) / ngfft(2), rprimd(3,3) / ngfft(3) /) 215 end if 216 217 if (usewvl == 0) then 218 ! We get the kernel. 219 call psolver_kernel( hgrid, 2, icoulomb, me, kernel, comm, ngfft, nproc, nscforder) 220 elseif(usewvl==1) then 221 ! In this case, the kernel is already computed. 222 ! We just shallow copy it. 223 kernel = wvl_den%denspot%pkernel 224 end if 225 226 if(usewvl==1) then 227 if(wvl_den%denspot%rhov_is .ne. ELECTRONIC_DENSITY) then 228 message= "psolver_rhohxc: rhov should contain the electronic density" 229 MSG_ERROR(message) 230 end if 231 end if 232 233 if(usewvl==1) then 234 n1i=wvl%Glr%d%n1i; n2i=wvl%Glr%d%n2i; n3i=wvl%Glr%d%n3i 235 n3d=wvl_den%denspot%dpbox%n3d 236 else 237 n1i=ngfft(1); n2i=ngfft(2) ; n3i=ngfft(3) 238 n3d=ngfft(13) 239 end if 240 241 if (usewvl == 0) then 242 ! ucvol_local=product(hgrid)*half**3*real(n1i*n2i*n3i,dp) 243 ! write(*,*)'hgrid, n1i,n2i,n3i',hgrid,ngfft(1:3) 244 ! write(*,*)'ucvol_local',ucvol_local 245 ucvol_local = ucvol 246 ! write(*,*)'ucvol_local',ucvol_local 247 else 248 ! We need to tune the volume when wavelets are used because, not 249 ! all FFT points are used. 250 ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 251 ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(product(wvl_den%denspot%dpbox%ndims), dp) 252 end if 253 254 !Core density array 255 if(n3xccc==0 .or. add_n_c_here) then 256 nullify(rhocore) 257 ! Pending, next line should follow the same logic that the rest 258 if(usewvl==1 .and. usepaw==0) rhocore=> wvl_den%denspot%rho_C 259 else 260 if(usepaw==1) then 261 ABI_ALLOCATE(rhocore,(n1i,n2i,n3d,1)) !not spin dependent 262 call wvl_rhov_abi2big(1,xccc3d,rhocore) 263 264 ! Make rhocore positive to avoid numerical instabilities in V_xc 265 iwarn=0 ; opt_mkdenpos=0 266 call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhocore, tol20 ) 267 end if 268 end if 269 270 !write(*,*)'psolver_rhohxc, erase me, set rhocore=0' 271 !if( associated(wvl_den%denspot%rho_C))wvl_den%denspot%rho_C=zero 272 !if(associated(rhocore))rhocore=zero 273 274 !Rhohat array: 275 if(test_nhat .and. .not. rest_hat_n_here) then 276 ! rhohat => nhat !do not know how to point 4 index to 2 index 277 ! here we have to copy since convention for spin changes. 278 ABI_ALLOCATE(rhohat,(n1i,n2i,n3d,nspden)) 279 call wvl_rhov_abi2big(1,nhat,rhohat) 280 else 281 nullify(rhohat) 282 end if 283 284 !Data are always distributed when using the wavelets, even if nproc = 1. 285 !The size given is the complete size of the box, not the distributed size 286 !stored in ngfft. 287 if (nproc > 1 .or. usewvl > 0) then 288 datacode = 'D' 289 else 290 datacode = 'G' 291 end if 292 293 !If usewvl=1, vpsp(or v_ext) will be summed to vhartree 294 if(usewvl==1) then 295 pot_ion=>wvl_den%denspot%V_ext 296 sumpion=.false. 297 ! Note: 298 ! if sumpion==.true. 299 ! call wvl_newvtr in setvtr and rhotov 300 ! if sumpion==.false. 301 ! modify setvtr and rhotov to not use wvl_newvtr and follow the normal ABINIT flow. 302 else 303 ! This is not allowed 304 ! pot_ion=>vxc !this is a dummy variable here 305 sumpion=.false. 306 end if 307 308 309 !To make this work, make sure that xc_init has been called 310 !in gstate. 311 if(.not. use_psolver) then 312 ! T.Rangel: 313 ! Use this approach for PAW and sometimes for NC since 314 ! in psolver() the core density is not added. 315 ! 316 ! PAW case: 317 ! It is important to call H_potential before XC_potential: 318 ! In XC_potential, if test_nhat, we do: 319 ! 1) rhor=rhor-rhohat, 320 ! 2) makepositive(rhor,tol20) 321 ! 3) after Psolver, we do rhor=rhor+rhohat, 322 ! I found that rhor at input and output are slightly different, 323 ! These differences lead to a difference of ~0.01 hartree in V_hartree. 324 ! If PAW, substract compensation density from effective density: 325 ! - if GGA, because nhat gradients are computed separately 326 ! - if nhat does not have to be included in XC 327 328 ! save rhor in rhonow to avoid modifying it. 329 ABI_ALLOCATE(rhonow,(nfft,nspden)) 330 ! copy rhor into rhonow: 331 ! ABINIT convention is followed: (ispin=1: for spin up + spin down) 332 do ispin=1,nspden 333 do ifft=1,nfft 334 rhonow(ifft,ispin)=abs(rhor(ifft,ispin))+1.0d-20 335 end do 336 end do 337 338 if(usewvl==1) then 339 call H_potential(datacode,& 340 & kernel,rhonow,pot_ion,enhartr,& 341 & zero,sumpion) 342 else 343 ! Vxc is passed as a dummy argument 344 call H_potential(datacode,& 345 & kernel,rhonow,vxc,enhartr,& 346 & zero,sumpion) 347 end if 348 ! 349 do ifft=1,nfft 350 vhartr(ifft)=rhonow(ifft,1) 351 end do 352 ! write(*,*)'erase me psolver_rhohxc l350, set vhartr=0' 353 ! vhartr=zero ; enhartr=zero 354 ! 355 ! Since rhonow was modified inside H_potential: 356 ! copy rhor again into rhonow following the BigDFT convention: 357 call wvl_rhov_abi2big(1,rhor,rhonow) 358 359 ! Add n_c here: 360 if(n3xccc>0 .and. add_n_c_here) then 361 do ispin=1,nspden 362 do ifft=1,nfft 363 rhonow(ifft,ispin)=rhonow(ifft,ispin)+xccc3d(ifft) 364 end do 365 end do 366 end if 367 ! Remove \hat{n} here: 368 if(test_nhat .and. rest_hat_n_here) then 369 do ispin=1,nspden 370 do ifft=1,nfft 371 rhonow(ifft,ispin)=rhonow(ifft,ispin)-nhat(ifft,ispin) 372 end do 373 end do 374 end if 375 376 ! Make the density positive everywhere (but do not care about gradients) 377 iwarn=0 ; opt_mkdenpos=0 378 call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhonow, xc_denpos) 379 ! do ispin=1,nspden 380 ! do ifft=1,nfft 381 ! rhonow(ifft,ispin)=abs(rhonow(ifft,ispin))+1.0d-20 382 ! end do 383 ! end do 384 385 ! If PAW, substract compensation density from effective density: 386 ! - if GGA, because nhat gradients are computed separately 387 ! - if nhat does not have to be included in XC 388 if (test_nhat .and. .not. rest_hat_n_here) then 389 390 call XC_potential(bndcode,datacode,me,nproc,comm,& 391 & n1i,n2i,n3i,& 392 & wvl_den%denspot%xc,hgrid(1),hgrid(2),hgrid(3),& 393 & rhonow,enxc,envxc,nspden,rhocore,& 394 & vxc,xcstr,rhohat=rhohat) 395 396 else 397 398 call XC_potential(bndcode,datacode,me,nproc,comm,& 399 & n1i,n2i,n3i,& 400 & wvl_den%denspot%xc,hgrid(1),hgrid(2),hgrid(3),& 401 & rhonow,enxc,envxc,nspden,rhocore,& 402 & vxc,xcstr) 403 404 end if 405 406 ! write(*,*)'psolver_rhohxc: erase me, set vxc=0' 407 ! vxc=zero 408 ! enxc=zero 409 ! envxc=zero 410 411 ! deallocate temporary array 412 ABI_DEALLOCATE(rhonow) 413 414 else 415 ! NC case: here we optimize memory, and we reuse vhartree to store rhor: 416 417 ! We save total rhor in vhartr 418 do ifft=1,nfft 419 vhartr(ifft) = rhor(ifft, 1) 420 end do 421 422 ! In non-wavelet case, we change the rhor values. 423 if (nspden == 2) then 424 do ifft = 1, nfft 425 ! We change rhor for psolver call. 426 tmpDown = rhor(ifft, 1) - rhonow(ifft, 2) 427 tmpUp = rhor(ifft, 2) 428 rhor(ifft, 1) = tmpUp 429 rhor(ifft, 2) = tmpDown 430 end do 431 end if 432 ! Make the density positive everywhere (but do not care about gradients) 433 iwarn=0 ; opt_mkdenpos=0 434 call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhor, xc_denpos) 435 ! do ispin=1,nspden 436 ! do ifft=1,nfft 437 ! rhor(ifft,ispin)=abs(rhor(ifft,ispin))+1.0d-20 438 ! end do 439 ! end do 440 441 ! Call Poisson solver, here rhor(:,1) will contain Vhartree at output 442 ! This does not compile, check mklocl_realspace where it do work. 443 ! call psolver(bndcode, datacode, me, nproc, n1i, & 444 !& n2i,n3i, ixc, hgrid(1), hgrid(2), hgrid(3), & 445 !& rhor, kernel, vxc, enhartr, enxc, envxc, 0.d0, .false., nspden) 446 447 ! PSolver work in place, we set back the rhor values. 448 do ifft = 1, nfft, 1 449 tmpPot = rhor(ifft, 1) 450 ! Rhor total was saved in vhartr and current rhor(:,2) is down spin 451 rhor(ifft, 1) = vhartr(ifft) 452 if (nspden == 2) rhor(ifft, 2) = rhor(ifft, 1) - rhor(ifft, 2) 453 vhartr(ifft) = tmpPot 454 end do 455 end if 456 457 !Pass vhartr and vxc to BigDFT objects (useless?) 458 !if(usewvl==1) then 459 ! write(message, '(a,a,a,a)' ) ch10, ' rhotoxc_wvlpaw : but why are you copying me :..o(' 460 ! call wvl_vhartr_abi2big(1,vhartr,wvl_den) 461 ! (this can be commented out, since we do not use denspot%v_xc 462 ! call wvl_vxc_abi2big(1,vxc,wvl_den) 463 !end if 464 465 !Compute vxcavg 466 call mean_fftr(vxc, vxcmean, nfft, nfftot, nspden,mpi_comm_sphgrid=comm) 467 vxcavg = vxcmean(1) 468 469 !Pass energies to wvl object: 470 if(usewvl==1) then 471 wvl_e%energs%eh = enhartr 472 wvl_e%energs%exc = enxc 473 wvl_e%energs%evxc= envxc 474 end if 475 476 !Nullify pointers and deallocate arrays 477 if(test_nhat .and. .not. rest_hat_n_here) then 478 ! if(nspden==2) ABI_DEALLOCATE(rhohat) 479 ABI_DEALLOCATE(rhohat) 480 if(associated(rhohat)) nullify(rhohat) 481 end if 482 if( n3xccc>0 .and. .not. add_n_c_here) then 483 if(usepaw==1) then 484 ABI_DEALLOCATE(rhocore) 485 end if 486 end if 487 if(associated(rhocore)) nullify(rhocore) 488 if(associated(pot_ion)) nullify(pot_ion) 489 490 #else 491 BIGDFT_NOTENABLED_ERROR() 492 if (.false.) write(std_out,*) nhatdim,nspden,n3xccc,nfft,icoulomb,ixc,nscforder,usewvl,& 493 & usexcnhat,usepaw,xclevel,rprimd(1,1),xc_denpos,enxc,envxc,enhartr,vxcavg,mpi_enreg%nproc,& 494 & wvl%h(1),wvl_den%symObj,wvl_e%energs,ngfft(1),xccc3d(1),nhat(1,1),rhor(1,1),vhartr(1),vxc(1,1) 495 #endif 496 497 DBG_EXIT("COLL") 498 499 end subroutine psolver_rhohxc