TABLE OF CONTENTS
ABINIT/m_psolver [ Modules ]
NAME
m_psolver
FUNCTION
Poisson solver
COPYRIGHT
Copyright (C) 1998-2024 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 .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_psolver 23 24 use defs_basis 25 use defs_wvltypes 26 use m_abicore 27 use m_errors 28 use m_abi2big 29 use m_cgtools 30 use m_xmpi 31 32 use defs_abitypes, only : mpi_type 33 use m_geometry, only : metric 34 use m_drivexc, only : mkdenpos 35 36 implicit none 37 38 private
ABINIT/Psolver_hartree [ Functions ]
NAME
Psolver_hartree
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 does not compute the xc energy nor potential. See psolver_rhohxc() to do it. WARNING : the XC energy and potential computation capability has been for spin-polarized case, as everything is done as if nspden=1
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). vhartr(nfft)=Hartree potential. 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.
SOURCE
537 subroutine psolver_hartree(enhartr, hgrid, icoulomb, me, mpi_comm, nfft, ngfft, nproc, & 538 & nscforder, nspden, rhor, vhartr, usewvl) 539 540 #if defined HAVE_BIGDFT 541 use BigDFT_API, only : coulomb_operator 542 use poisson_solver, only : H_potential 543 #endif 544 implicit none 545 546 !Arguments ------------------------------------ 547 !scalars 548 integer, intent(in) :: nfft, nspden, icoulomb, usewvl, mpi_comm, me, nproc, nscforder 549 real(dp), intent(out) :: enhartr 550 !arrays 551 integer, intent(in) :: ngfft(3) 552 real(dp),intent(in) :: hgrid(3) 553 real(dp),intent(in) :: rhor(nfft,nspden) 554 real(dp),intent(out) :: vhartr(nfft) 555 556 !Local variables------------------------------- 557 #if defined HAVE_BIGDFT 558 !scalars 559 character(len=500) :: message 560 character(len = 1) :: datacode, bndcode 561 !arrays 562 real(dp), dimension(1) :: pot_ion_dummy 563 type(coulomb_operator):: kernel 564 #endif 565 566 ! ********************************************************************* 567 568 #if defined HAVE_BIGDFT 569 570 if (icoulomb == 0) then 571 ! The kernel is built with 'P'eriodic boundary counditions. 572 bndcode = 'P' 573 else if (icoulomb == 1) then 574 ! The kernel is built with 'F'ree boundary counditions. 575 bndcode = 'F' 576 else if (icoulomb == 2) then 577 ! The kernel is built with 'S'urface boundary counditions. 578 bndcode = 'S' 579 end if 580 581 if(nspden > 2 .and. usewvl/=0 )then 582 write(message, '(a,a,a,i0)' )& 583 & 'Only non-spin-polarised or collinear spin is allowed for wavelets,',ch10,& 584 & 'while the argument nspden = ', nspden 585 ABI_BUG(message) 586 end if 587 588 !We do the computation. 589 write(message, "(A,A,A,3I6)") "Psolver_hartree(): compute potential (Vhartree)...", ch10, & 590 & " | dimension:", ngfft(1:3) 591 call wrtout(std_out, message,'COLL') 592 593 if (usewvl == 0) then 594 vhartr(:) = rhor(:, 1) 595 596 datacode = 'G' 597 ! This may not work with MPI in the planewave code... 598 else 599 if(nspden==1)vhartr(:) = rhor(:, 1) 600 if(nspden==2)vhartr(:) = rhor(:, 1) + rhor(:, 2) 601 ! The data are 'D'istributed in the wavelet case or 'G'lobal otherwise. 602 if (nproc > 1) then 603 datacode = 'D' 604 else 605 datacode = 'G' 606 end if 607 end if 608 609 !We get the kernel. 610 call psolver_kernel( hgrid, 2, icoulomb, me, kernel, mpi_comm, ngfft, nproc, nscforder) 611 612 613 !We attack PSolver with the total density contained in vhartr. 614 !This is also valid for spin-polarized (collinear and non-collinear) 615 !systems. Thus we enter nspden (last arg of PSolver) as being 1. 616 !Warning : enxc and evxc are meaningless. 617 ! call psolver(bndcode, datacode, me, nproc, ngfft(1), ngfft(2), ngfft(3),& 618 !& 0, hgrid(1), hgrid(2), hgrid(3), vhartr, kernel%co%kernel, pot_ion_dummy, & 619 !& enhartr, enxc, evxc, 0.d0, .false., 1) 620 621 call H_potential(datacode,kernel,vhartr,pot_ion_dummy,& 622 & enhartr,zero,.false.) 623 624 625 #else 626 BIGDFT_NOTENABLED_ERROR() 627 if (.false.) write(std_out,*) nfft,nspden,icoulomb,usewvl,mpi_comm,me,nproc,nscforder,enhartr,& 628 & ngfft(1),hgrid(1),rhor(1,1),vhartr(1) 629 #endif 630 631 end subroutine psolver_hartree
ABINIT/psolver_kernel [ Functions ]
NAME
psolver_kernel
FUNCTION
Build, get or free the kernel matrix used by the Poisson solver to compute the the convolution between 1/r and rho. The kernel is a saved variable. If this routine is called for building while a kernel already exists, it is not recomputed if all parameters (grid step and data size) are unchanged. Otherwise the kernel is freed and recompute again. The build action has a returned variable which is a pointer on the kernel. The get action also returns the kernel, or NULL if none has been associated.
INPUTS
iaction=0 to free all kernel allocated array,
1 to compute the kernel (parallel case),
2 to get it (parallel case),
3 to compute the kernel (sequential),
4 to get the sequential kernel.
rprimd(3,3)=dimensional primitive translations in real space (bohr)
OUTPUT
kernel= associated kernel on build (iaction = 1) and get action (iaction = 2).
SOURCE
660 subroutine psolver_kernel(hgrid, iaction, icoulomb, & 661 & iproc, kernel, mpi_comm, ngfft, nproc, nscforder) 662 663 #if defined HAVE_BIGDFT 664 use BigDFT_API, only : coulomb_operator,nullify_coulomb_operator, & 665 & deallocate_coulomb_operator,mpi_environment 666 use poisson_solver, only : pkernel_init,pkernel_set 667 #else 668 use defs_wvltypes, only : coulomb_operator 669 #endif 670 implicit none 671 672 !Arguments ------------------------------------ 673 !scalars 674 integer,intent(in) :: iaction, icoulomb, mpi_comm, nscforder, iproc, nproc 675 !arrays 676 integer, intent(in) :: ngfft(3) 677 type(coulomb_operator),intent(inout)::kernel 678 real(dp),intent(in) :: hgrid(3) 679 680 !Local variables------------------------- 681 #if defined HAVE_BIGDFT 682 !scalars 683 integer,parameter :: igpu=0 !no GPUs 684 !arrays 685 integer, save :: kernel_scfOrder 686 integer, save :: kernel_icoulomb 687 integer, save :: data_size(3) = (/ -2, -2, -2 /) 688 real(dp), save :: kernel_hgrid(3) ! Grid step used when creating the kernel. 689 character(len = 1) :: geocode 690 character(len=500) :: message 691 integer :: current_size(3) 692 type(coulomb_operator),save :: pkernel, kernelseq 693 type(mpi_environment) :: mpi_env 694 #endif 695 696 ! ************************************************************************* 697 698 #if defined HAVE_BIGDFT 699 700 if (icoulomb == 0) then 701 ! The kernel is built with 'P'eriodic boundary counditions. 702 geocode = 'P' 703 else if (icoulomb == 1) then 704 ! The kernel is built with 'F'ree boundary counditions. 705 geocode = 'F' 706 else if (icoulomb == 2) then 707 ! The kernel is built with 'S'urface boundary counditions. 708 geocode = 'S' 709 end if 710 current_size(:) = ngfft(1:3) 711 712 !Initialise kernel_array pointer. 713 if (maxval(data_size) == -2) then 714 call nullify_coulomb_operator(pkernel) 715 call nullify_coulomb_operator(kernelseq) 716 end if 717 718 !If iaction == 0, we free the kernel. 719 if (iaction == 0) then 720 if (associated(pkernel%kernel)) then 721 write(message, "(A)") "Psolver_kernel() : deallocating pkernel..." 722 call wrtout(std_out, message,'COLL') 723 724 call deallocate_coulomb_operator(pkernel) 725 end if 726 if (associated(kernelseq%kernel)) then 727 write(message, "(A)") "Psolver_kernel() : deallocating kernelseq..." 728 call wrtout(std_out, message,'COLL') 729 730 call deallocate_coulomb_operator(kernelseq) 731 end if 732 data_size = (/ -1, -1, -1 /) 733 return 734 end if 735 736 737 !Action is build or get. We check the sizes before doing anything else. 738 739 !!$!Get the size depending on wavelets calculations or not 740 !!$ if (dtset%usewvl == 0) then 741 !!$ hgrid(1) = rprimd(1, 1) / ngfft(1) 742 !!$ hgrid(2) = rprimd(2, 2) / ngfft(2) 743 !!$ hgrid(3) = rprimd(3, 3) / ngfft(3) 744 !!$ 745 !!$ else 746 !!$ hgrid(:) = 0.5d0 * wvl%h(:) 747 !!$ current_size(1:3) = (/ wvl%Glr%d%n1i, wvl%Glr%d%n2i, wvl%Glr%d%n3i /) 748 !!$ end if 749 750 !Compute a new kernel if grid size has changed or if the kernel 751 !has never been computed. 752 if ((iaction == 1 .and. .not. associated(pkernel%kernel)) .or. & 753 & (iaction == 3 .and. .not. associated(kernelseq%kernel)) .or. & 754 & kernel_icoulomb /= icoulomb .or. & 755 & data_size(1) /= current_size(1) .or. & 756 & data_size(2) /= current_size(2) .or. & 757 & data_size(3) /= current_size(3) .or. & 758 & kernel_hgrid(1) /= hgrid(1) .or. & 759 & kernel_hgrid(2) /= hgrid(2) .or. & 760 & kernel_hgrid(3) /= hgrid(3) .or. & 761 & kernel_scfOrder /= nscforder) then 762 write(message, "(A,A,A,3I6)") "Psolver_kernel() : building kernel...", ch10, & 763 & " | data dimensions:", current_size 764 call wrtout(std_out, message, 'COLL') 765 766 if (iaction == 1 .or. iaction == 2) then 767 if (associated(pkernel%kernel)) then 768 call deallocate_coulomb_operator(pkernel) 769 end if 770 mpi_env%mpi_comm = mpi_comm 771 mpi_env%iproc = iproc 772 mpi_env%nproc = nproc 773 mpi_env%igroup = 0 ! no task groups 774 mpi_env%ngroup = 1 ! no task groups 775 pkernel= pkernel_init(.True.,iproc,nproc,igpu,geocode,& 776 & current_size,hgrid,nscforder, mpi_env = mpi_env) 777 call pkernel_set(pkernel,.True.) 778 end if 779 780 if (iaction == 3 .or. iaction == 4) then 781 if (associated(kernelseq%kernel)) then 782 call deallocate_coulomb_operator(kernelseq) 783 end if 784 mpi_env%mpi_comm = mpi_comm 785 mpi_env%iproc = 0 786 mpi_env%nproc = 1 787 mpi_env%igroup = 0 ! no task groups 788 mpi_env%ngroup = 1 ! no task groups 789 kernelseq= pkernel_init(.True.,iproc,nproc,igpu,geocode,& 790 & current_size,hgrid,nscforder, mpi_env = mpi_env) 791 call pkernel_set(kernelseq,.True.) 792 end if 793 794 ! Storing variables which were used to make the kernel 795 kernel_icoulomb = icoulomb 796 data_size(:) = current_size(:) 797 kernel_hgrid(:) = hgrid(:) 798 kernel_scfOrder = nscforder 799 end if 800 801 ! Shallow copy if kernel has been associated. 802 if (iaction == 1 .or. iaction == 2) then 803 kernel = pkernel 804 end if 805 if (iaction == 3 .or. iaction == 4) then 806 kernel = kernelseq 807 end if 808 809 #else 810 BIGDFT_NOTENABLED_ERROR() 811 if (.false.) write(std_out,*) iaction,icoulomb,mpi_comm,nscforder,iproc,nproc,ngfft(1),kernel%co,hgrid(1) 812 #endif 813 814 end subroutine psolver_kernel
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().
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.
SOURCE
84 subroutine psolver_rhohxc(enhartr, enxc, envxc, icoulomb, ixc, & 85 & mpi_enreg, nfft, ngfft, nhat,nhatdim,& 86 & nscforder, nspden, n3xccc, rhor, rprimd,& 87 & usexcnhat,usepaw,usewvl,vhartr, vxc, vxcavg, wvl,wvl_den,wvl_e,& 88 & xccc3d,xclevel,xc_denpos) 89 90 #if defined HAVE_BIGDFT 91 use BigDFT_API, only : XC_potential,ELECTRONIC_DENSITY,coulomb_operator 92 use poisson_solver, only : H_potential 93 #endif 94 implicit none 95 96 !Arguments ------------------------------------ 97 !scalars 98 integer, intent(in) :: nhatdim,nspden,n3xccc 99 integer, intent(in) :: nfft, icoulomb, ixc, nscforder, usewvl 100 integer,intent(in) :: usexcnhat,usepaw,xclevel 101 real(dp),intent(in) :: rprimd(3,3) 102 real(dp), intent(in) :: xc_denpos 103 real(dp), intent(out) :: enxc, envxc, enhartr, vxcavg 104 type(mpi_type), intent(in) :: mpi_enreg 105 type(wvl_internal_type), intent(in) :: wvl 106 type(wvl_denspot_type), intent(inout) :: wvl_den 107 type(wvl_energy_terms), intent(inout) :: wvl_e 108 !arrays 109 integer, intent(in) :: ngfft(18) 110 real(dp),intent(in) :: xccc3d(n3xccc) 111 real(dp),intent(in) :: nhat(nfft,nspden*nhatdim) 112 real(dp),intent(inout) :: rhor(nfft, nspden) 113 real(dp),intent(out) :: vhartr(nfft) 114 real(dp),intent(out) :: vxc(nfft, nspden) 115 116 !Local variables------------------------------- 117 #if defined HAVE_BIGDFT 118 ! n_c and \hat{n} can be added/rested inside bigdft by passing 119 ! them as pointers (rhocore and rhohat): 120 logical, parameter :: add_n_c_here=.true. !Add n_c here or inside bigdft 121 logical, parameter :: rest_hat_n_here=.true. !Rest \hat{n} here or inside bigdft 122 !scalars 123 integer :: me,nproc,comm 124 integer :: ifft,ispin 125 integer :: iwarn, opt_mkdenpos 126 integer :: nfftot,ngrad 127 integer :: n1i,n2i,n3d,n3i 128 real(dp) :: tmpDown, tmpUp, tmpPot,ucvol,ucvol_local 129 logical :: sumpion,test_nhat,use_psolver=.false. 130 character(len=500) :: message 131 character(len = 1) :: datacode, bndcode 132 !arrays 133 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 134 real(dp) :: hgrid(3) 135 real(dp) :: vxcmean(1) 136 real(dp), pointer :: rhocore(:,:,:,:),rhohat(:,:,:,:) 137 real(dp), pointer :: pot_ion(:,:,:,:),rhonow(:,:) 138 real(dp), dimension(6) :: xcstr 139 type(coulomb_operator) :: kernel 140 #endif 141 142 ! ********************************************************************* 143 144 DBG_ENTER("COLL") 145 146 #if defined HAVE_BIGDFT 147 148 nfftot=PRODUCT(ngfft(1:3)) 149 comm=mpi_enreg%comm_fft 150 if(usewvl==1) comm=mpi_enreg%comm_wvl 151 me=xmpi_comm_rank(comm) 152 nproc=xmpi_comm_size(comm) 153 154 if(n3xccc>0) then 155 if(nfft .ne. n3xccc)then 156 write(message,'(a,a,a,2(i0,1x))')& 157 & 'nfft and n3xccc should be equal,',ch10,& 158 & 'however, nfft and n3xccc=',nfft,n3xccc 159 ABI_BUG(message) 160 end if 161 end if 162 if(nspden==4) then 163 ABI_ERROR('nspden==4 not coded yet') 164 end if 165 166 if (ixc==0) then 167 vxcavg=zero 168 test_nhat=.false. 169 170 ! No xc at all is applied (usually for testing) 171 ABI_WARNING('Note that no xc is applied (ixc=0).') 172 173 else if (ixc/=20) then 174 175 ! ngrad=1 is for LDAs or LSDs, ngrad=2 is for GGAs 176 ngrad=1;if(xclevel==2)ngrad=2 177 ! ixc 31 to 35 are for mgga test purpose only (fake functionals based on LDA but need the gradients too) 178 if(ixc>=31 .and. ixc<=35)ngrad=2 179 ! Test: has a compensation density to be added/substracted (PAW) ? 180 ! test_nhat=((nhatdim==1).and.(usexcnhat==0.or.(ngrad==2.and.nhatgrdim==1))) 181 test_nhat=((nhatdim==1).and.(usexcnhat==0)) 182 end if 183 184 185 !Compute different geometric tensor, as well as ucvol, from rprimd 186 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 187 188 if (icoulomb == 0) then 189 ! The kernel is built with 'P'eriodic boundary counditions. 190 bndcode = 'P' 191 else if (icoulomb == 1) then 192 ! The kernel is built with 'F'ree boundary counditions. 193 bndcode = 'F' 194 else if (icoulomb == 2) then 195 ! The kernel is built with 'S'urface boundary counditions. 196 bndcode = 'S' 197 end if 198 199 !This makes the tests fail? 200 !For NC and n_c=0, call psolver, which uses less memory: 201 !if(usepaw==0 .and. n3xccc==0) use_psolver=.true. 202 203 if(nspden > 2)then 204 write(message, '(a,a,a,i0)' )& 205 & 'Only non-spin-polarised or collinear spin is allowed,',ch10,& 206 & 'while the argument nspden = ', nspden 207 ABI_ERROR(message) 208 end if 209 210 !We do the computation. 211 write(message, "(A,A,A,3I6)") "psolver_rhohxc(): compute potentials (Vhartree and Vxc)...", ch10, & 212 & " | dimension:", ngfft(1:3) 213 call wrtout(std_out, message,'COLL') 214 215 if(usewvl==1) then 216 hgrid=(/wvl_den%denspot%dpbox%hgrids(1),wvl_den%denspot%dpbox%hgrids(2),wvl_den%denspot%dpbox%hgrids(3)/) 217 else 218 hgrid=(/ rprimd(1,1) / ngfft(1), rprimd(2,2) / ngfft(2), rprimd(3,3) / ngfft(3) /) 219 end if 220 221 if (usewvl == 0) then 222 ! We get the kernel. 223 call psolver_kernel( hgrid, 2, icoulomb, me, kernel, comm, ngfft, nproc, nscforder) 224 elseif(usewvl==1) then 225 ! In this case, the kernel is already computed. 226 ! We just shallow copy it. 227 kernel = wvl_den%denspot%pkernel 228 end if 229 230 if(usewvl==1) then 231 if(wvl_den%denspot%rhov_is .ne. ELECTRONIC_DENSITY) then 232 message= "psolver_rhohxc: rhov should contain the electronic density" 233 ABI_ERROR(message) 234 end if 235 end if 236 237 if(usewvl==1) then 238 n1i=wvl%Glr%d%n1i; n2i=wvl%Glr%d%n2i; n3i=wvl%Glr%d%n3i 239 n3d=wvl_den%denspot%dpbox%n3d 240 else 241 n1i=ngfft(1); n2i=ngfft(2) ; n3i=ngfft(3) 242 n3d=ngfft(13) 243 end if 244 245 if (usewvl == 0) then 246 ! ucvol_local=product(hgrid)*half**3*real(n1i*n2i*n3i,dp) 247 ! write(*,*)'hgrid, n1i,n2i,n3i',hgrid,ngfft(1:3) 248 ! write(*,*)'ucvol_local',ucvol_local 249 ucvol_local = ucvol 250 ! write(*,*)'ucvol_local',ucvol_local 251 else 252 ! We need to tune the volume when wavelets are used because, not 253 ! all FFT points are used. 254 ! ucvol_local = (half * dtset%wvl_hgrid) ** 3 * ngfft(1)*ngfft(2)*ngfft(3) 255 ucvol_local = product(wvl_den%denspot%dpbox%hgrids) * real(product(wvl_den%denspot%dpbox%ndims), dp) 256 end if 257 258 !Core density array 259 if(n3xccc==0 .or. add_n_c_here) then 260 nullify(rhocore) 261 ! Pending, next line should follow the same logic that the rest 262 if(usewvl==1 .and. usepaw==0) rhocore=> wvl_den%denspot%rho_C 263 else 264 if(usepaw==1) then 265 ABI_MALLOC(rhocore,(n1i,n2i,n3d,1)) !not spin dependent 266 call wvl_rhov_abi2big(1,xccc3d,rhocore) 267 268 ! Make rhocore positive to avoid numerical instabilities in V_xc 269 iwarn=0 ; opt_mkdenpos=0 270 call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhocore, tol20 ) 271 end if 272 end if 273 274 !write(*,*)'psolver_rhohxc, erase me, set rhocore=0' 275 !if( associated(wvl_den%denspot%rho_C))wvl_den%denspot%rho_C=zero 276 !if(associated(rhocore))rhocore=zero 277 278 !Rhohat array: 279 if(test_nhat .and. .not. rest_hat_n_here) then 280 ! rhohat => nhat !do not know how to point 4 index to 2 index 281 ! here we have to copy since convention for spin changes. 282 ABI_MALLOC(rhohat,(n1i,n2i,n3d,nspden)) 283 call wvl_rhov_abi2big(1,nhat,rhohat) 284 else 285 nullify(rhohat) 286 end if 287 288 !Data are always distributed when using the wavelets, even if nproc = 1. 289 !The size given is the complete size of the box, not the distributed size 290 !stored in ngfft. 291 if (nproc > 1 .or. usewvl > 0) then 292 datacode = 'D' 293 else 294 datacode = 'G' 295 end if 296 297 !If usewvl=1, vpsp(or v_ext) will be summed to vhartree 298 if(usewvl==1) then 299 pot_ion=>wvl_den%denspot%V_ext 300 sumpion=.false. 301 ! Note: 302 ! if sumpion==.true. 303 ! call wvl_newvtr in setvtr and rhotov 304 ! if sumpion==.false. 305 ! modify setvtr and rhotov to not use wvl_newvtr and follow the normal ABINIT flow. 306 else 307 ! This is not allowed 308 ! pot_ion=>vxc !this is a dummy variable here 309 sumpion=.false. 310 end if 311 312 313 !To make this work, make sure that xc_init has been called 314 !in gstate. 315 if(.not. use_psolver) then 316 ! T.Rangel: 317 ! Use this approach for PAW and sometimes for NC since 318 ! in psolver() the core density is not added. 319 ! 320 ! PAW case: 321 ! It is important to call H_potential before XC_potential: 322 ! In XC_potential, if test_nhat, we do: 323 ! 1) rhor=rhor-rhohat, 324 ! 2) makepositive(rhor,tol20) 325 ! 3) after Psolver, we do rhor=rhor+rhohat, 326 ! I found that rhor at input and output are slightly different, 327 ! These differences lead to a difference of ~0.01 hartree in V_hartree. 328 ! If PAW, substract compensation density from effective density: 329 ! - if GGA, because nhat gradients are computed separately 330 ! - if nhat does not have to be included in XC 331 332 ! save rhor in rhonow to avoid modifying it. 333 ABI_MALLOC(rhonow,(nfft,nspden)) 334 ! copy rhor into rhonow: 335 ! ABINIT convention is followed: (ispin=1: for spin up + spin down) 336 do ispin=1,nspden 337 do ifft=1,nfft 338 rhonow(ifft,ispin)=abs(rhor(ifft,ispin))+1.0d-20 339 end do 340 end do 341 342 if(usewvl==1) then 343 call H_potential(datacode,& 344 & kernel,rhonow,pot_ion,enhartr,& 345 & zero,sumpion) 346 else 347 ! Vxc is passed as a dummy argument 348 call H_potential(datacode,& 349 & kernel,rhonow,vxc,enhartr,& 350 & zero,sumpion) 351 end if 352 ! 353 do ifft=1,nfft 354 vhartr(ifft)=rhonow(ifft,1) 355 end do 356 ! write(*,*)'erase me psolver_rhohxc l350, set vhartr=0' 357 ! vhartr=zero ; enhartr=zero 358 ! 359 ! Since rhonow was modified inside H_potential: 360 ! copy rhor again into rhonow following the BigDFT convention: 361 call wvl_rhov_abi2big(1,rhor,rhonow) 362 363 ! Add n_c here: 364 if(n3xccc>0 .and. add_n_c_here) then 365 do ispin=1,nspden 366 do ifft=1,nfft 367 rhonow(ifft,ispin)=rhonow(ifft,ispin)+xccc3d(ifft) 368 end do 369 end do 370 end if 371 ! Remove \hat{n} here: 372 if(test_nhat .and. rest_hat_n_here) then 373 do ispin=1,nspden 374 do ifft=1,nfft 375 rhonow(ifft,ispin)=rhonow(ifft,ispin)-nhat(ifft,ispin) 376 end do 377 end do 378 end if 379 380 ! Make the density positive everywhere (but do not care about gradients) 381 iwarn=0 ; opt_mkdenpos=0 382 call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhonow, xc_denpos) 383 ! do ispin=1,nspden 384 ! do ifft=1,nfft 385 ! rhonow(ifft,ispin)=abs(rhonow(ifft,ispin))+1.0d-20 386 ! end do 387 ! end do 388 389 ! If PAW, substract compensation density from effective density: 390 ! - if GGA, because nhat gradients are computed separately 391 ! - if nhat does not have to be included in XC 392 if (test_nhat .and. .not. rest_hat_n_here) then 393 394 call XC_potential(bndcode,datacode,me,nproc,comm,& 395 & n1i,n2i,n3i,& 396 & wvl_den%denspot%xc,hgrid(1),hgrid(2),hgrid(3),& 397 & rhonow,enxc,envxc,nspden,rhocore,& 398 & vxc,xcstr,rhohat=rhohat) 399 400 else 401 402 call XC_potential(bndcode,datacode,me,nproc,comm,& 403 & n1i,n2i,n3i,& 404 & wvl_den%denspot%xc,hgrid(1),hgrid(2),hgrid(3),& 405 & rhonow,enxc,envxc,nspden,rhocore,& 406 & vxc,xcstr) 407 408 end if 409 410 ! write(*,*)'psolver_rhohxc: erase me, set vxc=0' 411 ! vxc=zero 412 ! enxc=zero 413 ! envxc=zero 414 415 ! deallocate temporary array 416 ABI_FREE(rhonow) 417 418 else 419 ! NC case: here we optimize memory, and we reuse vhartree to store rhor: 420 421 ! We save total rhor in vhartr 422 do ifft=1,nfft 423 vhartr(ifft) = rhor(ifft, 1) 424 end do 425 426 ! In non-wavelet case, we change the rhor values. 427 if (nspden == 2) then 428 do ifft = 1, nfft 429 ! We change rhor for psolver call. 430 tmpDown = rhor(ifft, 1) - rhonow(ifft, 2) 431 tmpUp = rhor(ifft, 2) 432 rhor(ifft, 1) = tmpUp 433 rhor(ifft, 2) = tmpDown 434 end do 435 end if 436 ! Make the density positive everywhere (but do not care about gradients) 437 iwarn=0 ; opt_mkdenpos=0 438 call mkdenpos(iwarn, nfft, nspden, opt_mkdenpos, rhor, xc_denpos) 439 ! do ispin=1,nspden 440 ! do ifft=1,nfft 441 ! rhor(ifft,ispin)=abs(rhor(ifft,ispin))+1.0d-20 442 ! end do 443 ! end do 444 445 ! Call Poisson solver, here rhor(:,1) will contain Vhartree at output 446 ! This does not compile, check mklocl_realspace where it do work. 447 ! call psolver(bndcode, datacode, me, nproc, n1i, & 448 !& n2i,n3i, ixc, hgrid(1), hgrid(2), hgrid(3), & 449 !& rhor, kernel, vxc, enhartr, enxc, envxc, 0.d0, .false., nspden) 450 451 ! PSolver work in place, we set back the rhor values. 452 do ifft = 1, nfft, 1 453 tmpPot = rhor(ifft, 1) 454 ! Rhor total was saved in vhartr and current rhor(:,2) is down spin 455 rhor(ifft, 1) = vhartr(ifft) 456 if (nspden == 2) rhor(ifft, 2) = rhor(ifft, 1) - rhor(ifft, 2) 457 vhartr(ifft) = tmpPot 458 end do 459 end if 460 461 !Pass vhartr and vxc to BigDFT objects (useless?) 462 !if(usewvl==1) then 463 ! write(message, '(a,a,a,a)' ) ch10, ' rhotoxc_wvlpaw : but why are you copying me :..o(' 464 ! call wvl_vhartr_abi2big(1,vhartr,wvl_den) 465 ! (this can be commented out, since we do not use denspot%v_xc 466 ! call wvl_vxc_abi2big(1,vxc,wvl_den) 467 !end if 468 469 !Compute vxcavg 470 call mean_fftr(vxc, vxcmean, nfft, nfftot, nspden,mpi_comm_sphgrid=comm) 471 vxcavg = vxcmean(1) 472 473 !Pass energies to wvl object: 474 if(usewvl==1) then 475 wvl_e%energs%eh = enhartr 476 wvl_e%energs%exc = enxc 477 wvl_e%energs%evxc= envxc 478 end if 479 480 !Nullify pointers and deallocate arrays 481 if(test_nhat .and. .not. rest_hat_n_here) then 482 ! if(nspden==2) ABI_FREE(rhohat) 483 ABI_FREE(rhohat) 484 if(associated(rhohat)) nullify(rhohat) 485 end if 486 if( n3xccc>0 .and. .not. add_n_c_here) then 487 if(usepaw==1) then 488 ABI_FREE(rhocore) 489 end if 490 end if 491 if(associated(rhocore)) nullify(rhocore) 492 if(associated(pot_ion)) nullify(pot_ion) 493 494 #else 495 BIGDFT_NOTENABLED_ERROR() 496 if (.false.) write(std_out,*) nhatdim,nspden,n3xccc,nfft,icoulomb,ixc,nscforder,usewvl,& 497 & usexcnhat,usepaw,xclevel,rprimd(1,1),xc_denpos,enxc,envxc,enhartr,vxcavg,mpi_enreg%nproc,& 498 & wvl%h(1),wvl_den%symObj,wvl_e%energs,ngfft(1),xccc3d(1),nhat(1,1),rhor(1,1),vhartr(1),vxc(1,1) 499 #endif 500 501 DBG_EXIT("COLL") 502 503 end subroutine psolver_rhohxc