TABLE OF CONTENTS


ABINIT/m_psolver [ Modules ]

[ Top ] [ Modules ]

NAME

  m_psolver

FUNCTION

  Poisson solver

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 .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_psolver
28 
29  use defs_basis
30  use defs_abitypes
31  use defs_wvltypes
32  use m_abicore
33  use m_errors
34  use m_abi2big
35  use m_cgtools
36  use m_xmpi
37 
38  use m_geometry, only : metric
39  use m_drivexc,  only : mkdenpos
40 
41  implicit none
42 
43  private

ABINIT/Psolver_hartree [ Functions ]

[ Top ] [ 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.

PARENTS

      mklocl_realspace,nres2vres

CHILDREN

      h_potential,psolver_kernel,wrtout

SOURCE

562 subroutine psolver_hartree(enhartr, hgrid, icoulomb, me, mpi_comm, nfft, ngfft, nproc, &
563      & nscforder, nspden, rhor, vhartr, usewvl)
564 
565 #if defined HAVE_BIGDFT
566  use BigDFT_API,     only : coulomb_operator
567  use poisson_solver, only : H_potential
568 #endif
569 
570 !This section has been created automatically by the script Abilint (TD).
571 !Do not modify the following lines by hand.
572 #undef ABI_FUNC
573 #define ABI_FUNC 'psolver_hartree'
574 !End of the abilint section
575 
576   implicit none
577 
578   !Arguments ------------------------------------
579   !scalars
580   integer, intent(in)           :: nfft, nspden, icoulomb, usewvl, mpi_comm, me, nproc, nscforder
581   real(dp), intent(out)         :: enhartr
582   !arrays
583   integer, intent(in)    :: ngfft(3)
584   real(dp),intent(in)    :: hgrid(3)
585   real(dp),intent(in)    :: rhor(nfft,nspden)
586   real(dp),intent(out)   :: vhartr(nfft)
587 
588   !Local variables-------------------------------
589 #if defined HAVE_BIGDFT
590   !scalars
591   character(len=500) :: message
592   character(len = 1) :: datacode, bndcode
593   !arrays
594   real(dp), dimension(1) :: pot_ion_dummy
595   type(coulomb_operator):: kernel
596 #endif
597 
598 ! *********************************************************************
599 
600 #if defined HAVE_BIGDFT
601 
602  if (icoulomb == 0) then
603 !  The kernel is built with 'P'eriodic boundary counditions.
604    bndcode = 'P'
605  else if (icoulomb == 1) then
606 !  The kernel is built with 'F'ree boundary counditions.
607    bndcode = 'F'
608  else if (icoulomb == 2) then
609 !  The kernel is built with 'S'urface boundary counditions.
610    bndcode = 'S'
611  end if
612 
613  if(nspden > 2 .and. usewvl/=0 )then
614    write(message, '(a,a,a,i0)' )&
615 &   'Only non-spin-polarised or collinear spin is allowed for wavelets,',ch10,&
616 &   'while the argument nspden = ', nspden
617    MSG_BUG(message)
618  end if
619 
620 !We do the computation.
621  write(message, "(A,A,A,3I6)") "Psolver_hartree(): compute potential (Vhartree)...", ch10, &
622 & " | dimension:", ngfft(1:3)
623  call wrtout(std_out, message,'COLL')
624 
625  if (usewvl == 0) then
626    vhartr(:)  = rhor(:, 1)
627 
628    datacode = 'G'
629 !  This may not work with MPI in the planewave code...
630  else
631    if(nspden==1)vhartr(:)  = rhor(:, 1)
632    if(nspden==2)vhartr(:)  = rhor(:, 1) + rhor(:, 2)
633 !  The data are 'D'istributed in the wavelet case or 'G'lobal otherwise.
634    if (nproc > 1) then
635      datacode = 'D'
636    else
637      datacode = 'G'
638    end if
639  end if
640 
641 !We get the kernel.
642  call psolver_kernel( hgrid, 2, icoulomb, me, kernel, mpi_comm, ngfft, nproc, nscforder)
643 
644 
645 !We attack PSolver with the total density contained in vhartr.
646 !This is also valid for spin-polarized (collinear and non-collinear)
647 !systems. Thus we enter nspden (last arg of PSolver) as being 1.
648 !Warning : enxc and evxc are meaningless.
649 ! call psolver(bndcode, datacode, me, nproc, ngfft(1), ngfft(2), ngfft(3),&
650 !& 0, hgrid(1), hgrid(2), hgrid(3), vhartr, kernel%co%kernel, pot_ion_dummy, &
651 !& enhartr, enxc, evxc, 0.d0, .false., 1)
652 
653  call H_potential(datacode,kernel,vhartr,pot_ion_dummy,&
654 & enhartr,zero,.false.)
655 
656 
657 #else
658  BIGDFT_NOTENABLED_ERROR()
659  if (.false.) write(std_out,*)  nfft,nspden,icoulomb,usewvl,mpi_comm,me,nproc,nscforder,enhartr,&
660 & ngfft(1),hgrid(1),rhor(1,1),vhartr(1)
661 #endif
662 
663 end subroutine psolver_hartree

ABINIT/psolver_kernel [ Functions ]

[ Top ] [ 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).

PARENTS

      gstate,mklocl_realspace,psolver_hartree,psolver_rhohxc
      wvl_wfsinp_reformat

CHILDREN

      deallocate_coulomb_operator,nullify_coulomb_operator,pkernel_set,wrtout

SOURCE

699 subroutine psolver_kernel(hgrid, iaction,  icoulomb, &
700 & iproc, kernel, mpi_comm, ngfft, nproc, nscforder)
701 
702 #if defined HAVE_BIGDFT
703  use BigDFT_API,     only  : coulomb_operator,nullify_coulomb_operator, &
704 &                            deallocate_coulomb_operator,mpi_environment
705  use poisson_solver, only : pkernel_init,pkernel_set
706 #else
707  use defs_wvltypes,  only : coulomb_operator
708 #endif
709 
710 !This section has been created automatically by the script Abilint (TD).
711 !Do not modify the following lines by hand.
712 #undef ABI_FUNC
713 #define ABI_FUNC 'psolver_kernel'
714 !End of the abilint section
715 
716   implicit none
717 
718 !Arguments ------------------------------------
719   !scalars
720   integer,intent(in) :: iaction, icoulomb, mpi_comm, nscforder, iproc, nproc
721   !arrays
722   integer, intent(in) :: ngfft(3)
723   type(coulomb_operator),intent(inout)::kernel
724   real(dp),intent(in) :: hgrid(3)
725 
726 !Local variables-------------------------
727 #if defined HAVE_BIGDFT
728   !scalars
729   integer,parameter :: igpu=0 !no GPUs
730   !arrays
731   integer, save :: kernel_scfOrder
732   integer, save :: kernel_icoulomb
733   integer, save :: data_size(3) = (/ -2, -2, -2 /)
734   real(dp), save :: kernel_hgrid(3)   ! Grid step used when creating the kernel.
735   character(len = 1) :: geocode
736   character(len=500) :: message
737   integer :: current_size(3)
738   type(coulomb_operator),save :: pkernel, kernelseq
739   type(mpi_environment) :: mpi_env
740 #endif
741 
742 ! *************************************************************************
743 
744 #if defined HAVE_BIGDFT
745 
746  if (icoulomb == 0) then
747 !  The kernel is built with 'P'eriodic boundary counditions.
748    geocode = 'P'
749  else if (icoulomb == 1) then
750 !  The kernel is built with 'F'ree boundary counditions.
751    geocode = 'F'
752  else if (icoulomb == 2) then
753 !  The kernel is built with 'S'urface boundary counditions.
754    geocode = 'S'
755  end if
756  current_size(:) = ngfft(1:3)
757 
758 !Initialise kernel_array pointer.
759  if (maxval(data_size) == -2) then
760    call nullify_coulomb_operator(pkernel)
761    call nullify_coulomb_operator(kernelseq)
762  end if
763 
764 !If iaction == 0, we free the kernel.
765  if (iaction == 0) then
766    if (associated(pkernel%kernel)) then
767      write(message, "(A)") "Psolver_kernel() : deallocating pkernel..."
768      call wrtout(std_out, message,'COLL')
769 
770      call deallocate_coulomb_operator(pkernel)
771    end if
772    if (associated(kernelseq%kernel)) then
773      write(message, "(A)") "Psolver_kernel() : deallocating kernelseq..."
774      call wrtout(std_out, message,'COLL')
775 
776      call deallocate_coulomb_operator(kernelseq)
777    end if
778    data_size = (/ -1, -1, -1 /)
779    return
780  end if
781 
782 
783 !Action is build or get. We check the sizes before doing anything else.
784 
785 !!$!Get the size depending on wavelets calculations or not
786 !!$ if (dtset%usewvl == 0) then
787 !!$   hgrid(1) = rprimd(1, 1) / ngfft(1)
788 !!$   hgrid(2) = rprimd(2, 2) / ngfft(2)
789 !!$   hgrid(3) = rprimd(3, 3) / ngfft(3)
790 !!$
791 !!$ else
792 !!$   hgrid(:) = 0.5d0 * wvl%h(:)
793 !!$   current_size(1:3) = (/ wvl%Glr%d%n1i, wvl%Glr%d%n2i, wvl%Glr%d%n3i /)
794 !!$ end if
795 
796 !Compute a new kernel if grid size has changed or if the kernel
797 !has never been computed.
798  if ((iaction == 1 .and. .not. associated(pkernel%kernel))    .or. &
799 & (iaction == 3 .and. .not. associated(kernelseq%kernel)) .or. &
800 & kernel_icoulomb /= icoulomb  .or. &
801 & data_size(1)    /= current_size(1) .or. &
802 & data_size(2)    /= current_size(2) .or. &
803 & data_size(3)    /= current_size(3) .or. &
804 & kernel_hgrid(1) /= hgrid(1)        .or. &
805 & kernel_hgrid(2) /= hgrid(2)        .or. &
806 & kernel_hgrid(3) /= hgrid(3)        .or. &
807 & kernel_scfOrder /= nscforder) then
808    write(message, "(A,A,A,3I6)") "Psolver_kernel() : building kernel...", ch10, &
809 &   " | data dimensions:", current_size
810    call wrtout(std_out, message, 'COLL')
811 
812    if (iaction == 1 .or. iaction == 2) then
813      if (associated(pkernel%kernel)) then
814        call deallocate_coulomb_operator(pkernel)
815      end if
816      mpi_env%mpi_comm = mpi_comm
817      mpi_env%iproc    = iproc
818      mpi_env%nproc    = nproc
819      mpi_env%igroup   = 0 ! no task groups
820      mpi_env%ngroup   = 1 ! no task groups
821      pkernel= pkernel_init(.True.,iproc,nproc,igpu,geocode,&
822 &     current_size,hgrid,nscforder, mpi_env = mpi_env)
823      call pkernel_set(pkernel,.True.)
824    end if
825 
826    if (iaction == 3 .or. iaction == 4) then
827      if (associated(kernelseq%kernel)) then
828        call deallocate_coulomb_operator(kernelseq)
829      end if
830      mpi_env%mpi_comm = mpi_comm
831      mpi_env%iproc    = 0
832      mpi_env%nproc    = 1
833      mpi_env%igroup   = 0 ! no task groups
834      mpi_env%ngroup   = 1 ! no task groups
835      kernelseq= pkernel_init(.True.,iproc,nproc,igpu,geocode,&
836 &     current_size,hgrid,nscforder, mpi_env = mpi_env)
837      call pkernel_set(kernelseq,.True.)
838    end if
839 
840 !  Storing variables which were used to make the kernel
841    kernel_icoulomb = icoulomb
842    data_size(:)    = current_size(:)
843    kernel_hgrid(:) = hgrid(:)
844    kernel_scfOrder = nscforder
845  end if
846 
847  ! Shallow copy if kernel has been associated.
848  if (iaction == 1 .or. iaction == 2) then
849    kernel = pkernel
850  end if
851  if (iaction == 3 .or. iaction == 4) then
852    kernel = kernelseq
853  end if
854 
855 #else
856  BIGDFT_NOTENABLED_ERROR()
857  if (.false.) write(std_out,*)  iaction,icoulomb,mpi_comm,nscforder,iproc,nproc,ngfft(1),kernel%co,hgrid(1)
858 #endif
859 
860 end subroutine psolver_kernel

ABINIT/psolver_rhohxc [ Functions ]

[ Top ] [ 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.

PARENTS

      energy,rhotov,scfcv,setvtr

CHILDREN

      h_potential,mean_fftr,metric,mkdenpos,psolver_kernel,wrtout
      wvl_rhov_abi2big,xc_potential

SOURCE

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