TABLE OF CONTENTS

ABINIT/mkdenpos [ Functions ]

[ Top ] [ Functions ]

NAME

 mkdenpos

FUNCTION

 Make a ground-state density positive everywhere :
 when the density (or spin-density) is smaller than xc_denpos,
 set it to the value of xc_denpos

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR)
 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

  nfft=(effective) number of FFT grid points (for this processor)
  nspden=number of spin-density components (max. 2)
  option=0 if density rhonow is stored as (up,dn)
         1 if density rhonow is stored as (up+dn,up)
         Active only when nspden=2
  xc_denpos= lowest allowed density (usually for the computation of the XC functionals)

OUTPUT

  (see side effects)

SIDE EFFECTS

  Input/output
  iwarn=At input: iwarn=0 a warning will be printed when rho is negative
                  iwarn>0 no warning will be printed out
        At output: iwarn is increased by 1
  rhonow(nfft,nspden)=electron (spin)-density in real space,
     either on the unshifted grid (if ishift==0,
     then equal to rhor),or on the shifted grid

NOTES

  At this stage, rhonow(:,1:nspden) contains the density in real space,
  on the unshifted or shifted grid. Now test for negative densities
  Note that, ignoring model core charge, as long as boxcut>=2
  the shifted density is derivable from the square of a Fourier
  interpolated charge density => CANNOT go < 0.
  However, actually can go < 0 to within machine precision;
  do not print useless warnings in this case, just fix it.
  Fourier interpolated core charge can go < 0 due to Gibbs
  oscillations; could avoid this by recomputing the model core
  charge at the new real space grid points (future work).

PARENTS

      bethe_salpeter,m_pawxc,mkcore_wvl,posdoppler,poslifetime,posratecore
      psolver_rhohxc,rhohxcpositron,rhotoxc,wvl_initro

CHILDREN

SOURCE

 58 #if defined HAVE_CONFIG_H
 59 #include "config.h"
 60 #endif
 61 
 62 #include "abi_common.h"
 63 
 64 
 65 subroutine mkdenpos(iwarn,nfft,nspden,option,rhonow,xc_denpos)
 66 
 67  use defs_basis
 68  use m_profiling_abi
 69  use m_errors
 70 
 71 !This section has been created automatically by the script Abilint (TD).
 72 !Do not modify the following lines by hand.
 73 #undef ABI_FUNC
 74 #define ABI_FUNC 'mkdenpos'
 75 !End of the abilint section
 76 
 77  implicit none
 78 
 79 !Arguments ------------------------------------
 80 !scalars
 81  integer,intent(in) :: nfft,nspden,option
 82  integer,intent(inout) :: iwarn
 83  real(dp),intent(in) :: xc_denpos
 84 !arrays
 85  real(dp),intent(inout) :: rhonow(nfft,nspden)
 86 
 87 !Local variables-------------------------------
 88 !scalars
 89  integer :: ifft,ispden,numneg
 90  real(dp) :: rhotmp,worst
 91  character(len=500) :: message
 92 !arrays
 93  real(dp) :: rho(2)
 94 
 95 ! *************************************************************************
 96 
 97  numneg=0
 98  worst=zero
 99 
100  if(nspden==1)then
101 
102 !  Non spin-polarized
103 !$OMP PARALLEL DO PRIVATE(ifft,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) SHARED(nfft,rhonow)
104    do ifft=1,nfft
105      rhotmp=rhonow(ifft,1)
106      if(rhotmp<xc_denpos)then
107        if(rhotmp<-xc_denpos)then
108 !        This case is probably beyond machine precision considerations
109          worst=min(worst,rhotmp)
110          numneg=numneg+1
111        end if
112        rhonow(ifft,1)=xc_denpos
113      end if
114    end do
115  else if (nspden==2) then
116 
117 !  Spin-polarized
118 
119 !  rhonow is stored as (up,dn)
120    if (option==0) then
121 
122 !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) REDUCTION(MIN:worst) REDUCTION(+:numneg) &
123 !$OMP&SHARED(nfft,nspden,rhonow)
124      do ifft=1,nfft
125 !      For polarized case, rho(1) is spin-up density, rho(2) is spin-down density
126        rho(1)=rhonow(ifft,1)
127        rho(2)=rhonow(ifft,2)
128        do ispden=1,nspden
129          if (rho(ispden)<xc_denpos) then
130            if (rho(ispden)<-xc_denpos) then
131 !            This case is probably beyond machine precision considerations
132              worst=min(worst,rho(ispden))
133              numneg=numneg+1
134            end if
135            rhonow(ifft,ispden)=xc_denpos
136          end if
137        end do
138      end do
139 
140 !    rhonow is stored as (up+dn,up)
141    else if (option==1) then
142 
143 !$OMP PARALLEL DO PRIVATE(ifft,ispden,rho,rhotmp) &
144 !$OMP&REDUCTION(MIN:worst) REDUCTION(+:numneg) &
145 !$OMP&SHARED(nfft,nspden,rhonow)
146      do ifft=1,nfft
147 !      For polarized case, rho(1) is spin-up density, rho(2) is spin-down density
148        rho(1)=rhonow(ifft,2)
149        rho(2)=rhonow(ifft,1)-rho(1)
150        do ispden=1,nspden
151          if (rho(ispden)<xc_denpos) then
152            if (rho(ispden)<-xc_denpos) then
153 !            This case is probably beyond machine precision considerations
154              worst=min(worst,rho(ispden))
155              numneg=numneg+1
156            end if
157            rho(ispden)=xc_denpos
158            rhonow(ifft,1)=rho(1)+rho(2)
159            rhonow(ifft,2)=rho(1)
160          end if
161        end do
162      end do
163 
164    end if  ! option
165 
166  else
167    MSG_BUG('nspden>2 not allowed !')
168  end if ! End choice between non-spin polarized and spin-polarized.
169 
170  if (numneg>0) then
171    if (iwarn==0) then
172      write(message,'(a,i0,a,a,a,es10.2,a,e10.2,a,a,a,a)')&
173 &     'Density went too small (lower than xc_denpos) at ',numneg,' points',ch10,&
174 &     'and was set to xc_denpos = ',xc_denpos,'. Lowest was ',worst,'.',ch10,&
175 &     'Likely due to too low boxcut or too low ecut for',' pseudopotential core charge.'
176      MSG_WARNING(message)
177    end if
178    iwarn=iwarn+1
179  end if
180 
181 end subroutine mkdenpos