TABLE OF CONTENTS


ABINIT/pawdensities [ Functions ]

[ Top ] [ Functions ]

NAME

 pawdensities

FUNCTION

 Compute PAW on-site densities (all-electron ,pseudo and compensation) for a given atom

COPYRIGHT

 Copyright (C) 1998-2018 ABINIT group (MT)
 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

  cplex: if 1, on-site densities are REAL, if 2, COMPLEX (response function only)
  iatom=index of current atom (note: this is the absolute index, not the index on current proc)
  lm_size=number of (l,m) moments
  lmselectin(lm_size)=flags selecting the non-zero LM-moments of on-site densities
                      (value of these flags at input; must be .TRUE. for nzlmopt/=1)
  nspden=number of spin-density components
  nzlmopt=if -1, compute all LM-moments of densities (lmselectin=.true. forced)
                 initialize "lmselectout" (index of non-zero LM-moments of densities)
          if  0, compute all LM-moments of densities (lmselectin=.true. forced)
                 force "lmselectout" to .true. (index of non-zero LM-moments of densities)
          if  1, compute only non-zero LM-moments of densities (stored before in "lmselectin")
  one_over_rad2(mesh_size)= contains 1/r**2 for each point of the radial grid -optional argument-
  opt_compch=flag controlling the accumulation of compensation charge density moments
             inside PAW spheres (compch_sph)
  opt_dens=flag controlling which on-site density(ies) is (are) computed
           0: all on-site densities (all-electron, pseudo and compensation)
           1: all-electron and pseudo densities (no compensation)
           2: only all-electron density
  opt_l=controls which l-moment(s) contribute to the density:
        <0 : all l contribute
        >=0: only l=opt_l contributes
        Note: opt_l>=0 is only compatible with opt_dens=2
  opt_print=1 if the densities moments have to be printed out (only if pawprtvol>=2)
  pawang <type(pawang_type)>=paw angular mesh and related data
  pawprtvol=control print volume and debugging output for PAW
  pawrad <type(pawrad_type)>=paw radial mesh and related data (for the current atom type)
  pawrhoij <type(pawrhoij_type)>= paw rhoij occupancies and related data (for the current atom)
  pawtab <type(pawtab_type)>=paw tabulated starting data (for the current atom type)

OUTPUT

  nhat1(cplex*mesh_size,lm_size,nspden)= compensation charge on-site density for current atom
  rho1(cplex*mesh_size,lm_size,nspden)= all electron on-site density for current atom
  trho1(cplex*mesh_size,lm_size,nspden)= pseudo on-site density for current atom
  ==== if nzlmopt/=1
    lmselectout(lm_size)=flags selecting the non-zero LM-moments of on-site densities
                         (value of these flags at output if updated, i.e. if nzlmopt<1)

SIDE EFFECTS

  ==== if opt_compch==1
    compch_sph=compensation charge integral inside spheres computed over spherical meshes
               updated with the contribution of current atom

PARENTS

      make_efg_onsite,pawdenpot,pawdfptenergy,poslifetime,posratecore

CHILDREN

      pawrad_deducer0,simp_gen,wrtout

SOURCE

 67 #if defined HAVE_CONFIG_H
 68 #include "config.h"
 69 #endif
 70 
 71 #include "abi_common.h"
 72 
 73 subroutine pawdensities(compch_sph,cplex,iatom,lmselectin,lmselectout,lm_size,nhat1,nspden,nzlmopt,&
 74 &          opt_compch,opt_dens,opt_l,opt_print,pawang,pawprtvol,pawrad,pawrhoij,pawtab,rho1,trho1,&
 75 &          one_over_rad2) ! optional
 76 
 77  use m_profiling_abi
 78  use defs_basis
 79  use m_errors
 80  use m_pawang, only : pawang_type
 81  use m_pawrad, only : pawrad_type, pawrad_deducer0, simp_gen
 82  use m_pawtab, only : pawtab_type
 83  use m_pawrhoij, only : pawrhoij_type
 84 
 85 !This section has been created automatically by the script Abilint (TD).
 86 !Do not modify the following lines by hand.
 87 #undef ABI_FUNC
 88 #define ABI_FUNC 'pawdensities'
 89  use interfaces_14_hidewrite
 90 !End of the abilint section
 91 
 92  implicit none
 93 
 94 !Arguments ---------------------------------------------
 95 !scalars
 96  integer,intent(in) :: cplex,iatom,lm_size,nspden,nzlmopt,opt_compch,opt_dens,opt_l,opt_print,pawprtvol
 97 ! jmb  real(dp),intent(out) :: compch_sph
 98  real(dp),intent(inout) :: compch_sph
 99  type(pawang_type),intent(in) :: pawang
100  type(pawrad_type),intent(in) :: pawrad
101  type(pawrhoij_type),intent(in) :: pawrhoij
102  type(pawtab_type),intent(in) :: pawtab
103 !arrays
104  logical,intent(in) :: lmselectin(lm_size)
105  logical,intent(inout) :: lmselectout(lm_size)
106  real(dp),intent(in),target,optional :: one_over_rad2(pawtab%mesh_size)
107  real(dp),intent(out) :: nhat1(cplex*pawtab%mesh_size,lm_size,nspden*(1-((opt_dens+1)/2)))
108  real(dp),intent(out) ::  rho1(cplex*pawtab%mesh_size,lm_size,nspden)
109  real(dp),intent(out) :: trho1(cplex*pawtab%mesh_size,lm_size,nspden*(1-(opt_dens/2)))
110 
111 !Local variables ---------------------------------------
112 !scalars
113  integer :: dplex,ii,ilm,iplex,ir,irhoij,isel,ispden,jrhoij
114  integer :: klm,klmn,kln,ll,lmax,lmin,mesh_size
115  real(dp) :: m1,mt1,rdum
116  character(len=500) :: msg
117 !arrays
118  real(dp) :: compchspha(cplex),compchsphb(cplex),ro(cplex),ro_ql(cplex),ro_rg(cplex)
119  real(dp),allocatable :: aa(:),bb(:)
120  real(dp),pointer :: one_over_rad2_(:)
121 
122 ! *************************************************************************
123 
124  DBG_ENTER("COLL")
125 
126 !Compatibility tests
127  if (opt_dens/=2.and.opt_l>=0) then
128    msg='  opt_dens/=2 incompatible with opt_l>=0 !'
129    MSG_BUG(msg)
130  end if
131  if(nzlmopt/=0.and.nzlmopt/=1.and.nzlmopt/=-1) then
132    msg='  invalid value for variable "nzlmopt".'
133    MSG_BUG(msg)
134  end if
135  if(nspden>pawrhoij%nspden) then
136    msg='  nspden must be <= pawrhoij%nspden !'
137    MSG_BUG(msg)
138  end if
139  if (cplex>pawrhoij%cplex) then
140    msg='  cplex must be <= pawrhoij%cplex !'
141    MSG_BUG(msg)
142  end if
143  if (nzlmopt/=1) then
144    if (any(.not.lmselectin(1:lm_size))) then
145      msg='  With nzlmopt/=1, lmselectin must be true !'
146      MSG_BUG(msg)
147    end if
148  end if
149  if (pawang%gnt_option==0) then
150    msg='  pawang%gnt_option=0 !'
151    MSG_BUG(msg)
152  end if
153 
154 !Various inits
155  rho1=zero
156  if (opt_dens <2) trho1=zero
157  if (opt_dens==0) nhat1=zero
158  mesh_size=pawtab%mesh_size;dplex=cplex-1
159  if (nzlmopt<1) lmselectout(1:lm_size)=.true.
160  if (present(one_over_rad2)) then
161    one_over_rad2_ => one_over_rad2
162  else
163    ABI_ALLOCATE(one_over_rad2_,(mesh_size))
164    one_over_rad2_(1)=zero
165    one_over_rad2_(2:mesh_size)=one/pawrad%rad(2:mesh_size)**2
166  end if
167 
168 !===== Compute "on-site" densities (n1, ntild1, nhat1) =====
169 !==========================================================
170 
171  do ispden=1,nspden
172 
173 !  -- Loop over ij channels (basis components)
174    jrhoij=1
175    do irhoij=1,pawrhoij%nrhoijsel
176      klmn=pawrhoij%rhoijselect(irhoij)
177      klm =pawtab%indklmn(1,klmn)
178      kln =pawtab%indklmn(2,klmn)
179      lmin=pawtab%indklmn(3,klmn)
180      lmax=pawtab%indklmn(4,klmn)
181 
182 
183 !    Retrieve rhoij
184      if (pawrhoij%nspden/=2) then
185        ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,ispden)
186      else
187        if (ispden==1) then
188          ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,1)&
189 &         +pawrhoij%rhoijp(jrhoij:jrhoij+dplex,2)
190        else if (ispden==2) then
191          ro(1:cplex)=pawrhoij%rhoijp(jrhoij:jrhoij+dplex,1)
192        end if
193      end if
194      ro(1:cplex)=pawtab%dltij(klmn)*ro(1:cplex)
195 
196 !    First option: all on-site densities are computed (opt_dens==0)
197 !    --------------------------------------------------------------
198      if (opt_dens==0) then
199        do ll=lmin,lmax,2
200          do ilm=ll**2+1,(ll+1)**2
201            if (lmselectin(ilm)) then
202              isel=pawang%gntselect(ilm,klm)
203              if (isel>0) then
204                ro_ql(1:cplex)=ro(1:cplex)*pawtab%qijl(ilm,klmn)
205                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
206 !              == nhat1(r=0)
207                nhat1(1:cplex,ilm,ispden)=nhat1(1:cplex,ilm,ispden) &
208 &               +ro_ql(1:cplex)*pawtab%shapefunc(1,ll+1) 
209 !              == rho1(r>0), trho1(r>0), nhat1(r>0)
210                do ir=2,mesh_size
211                  rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
212 &                 +ro_rg(1:cplex)*pawtab%phiphj  (ir,kln)*one_over_rad2_(ir)
213                  trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
214 &                 +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
215                  nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)=nhat1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
216 &                 +ro_ql(1:cplex)*pawtab%shapefunc(ir,ll+1)
217                end do
218              end if
219            end if
220          end do  ! End loops over ll,lm
221        end do
222 
223 !      2nd option: AE and pseudo densities are computed (opt_dens==1)
224 !      --------------------------------------------------------------
225      else if (opt_dens==1) then
226        do ll=lmin,lmax,2
227          do ilm=ll**2+1,(ll+1)**2
228            if (lmselectin(ilm)) then
229              isel=pawang%gntselect(ilm,klm)
230              if (isel>0) then
231                ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
232 !              == rho1(r>0), trho1(r>0)
233                do ir=2,mesh_size
234                  rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
235 &                 +ro_rg(1:cplex)*pawtab%phiphj  (ir,kln)*one_over_rad2_(ir)
236                  trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)=trho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
237 &                 +ro_rg(1:cplex)*pawtab%tphitphj(ir,kln)*one_over_rad2_(ir)
238                end do
239              end if
240            end if
241          end do  ! End loops over ll,lm
242        end do
243 
244 !      3rd option: only all-electron on-site density is computed (opt_dens==2)
245 !      -----------------------------------------------------------------------
246      else if (opt_dens==2) then
247        if (opt_l<0.or.(pawtab%indklmn(3,klmn)==0.and.pawtab%indklmn(4,klmn)==2*opt_l)) then
248          do ll=lmin,lmax,2
249            do ilm=ll**2+1,(ll+1)**2
250              if (lmselectin(ilm)) then
251                isel=pawang%gntselect(ilm,klm)
252                if (isel>0) then
253                  ro_rg(1:cplex)=ro(1:cplex)*pawang%realgnt(isel)
254 !                == rho1(r>0)
255                  do ir=2,mesh_size
256                    rho1(cplex*ir-dplex:ir*cplex,ilm,ispden) =rho1(cplex*ir-dplex:ir*cplex,ilm,ispden)&
257 &                   +ro_rg(1:cplex)*pawtab%phiphj(ir,kln)*one_over_rad2_(ir)
258                  end do
259                end if
260              end if
261            end do  ! End loops over ll, lm
262          end do
263        end if
264      end if
265 
266 !    -- End loop over ij channels
267      jrhoij=jrhoij+pawrhoij%cplex
268    end do
269 
270 !  Scale densities with 1/r**2 and compute rho1(r=0) and trho1(r=0)
271    if (cplex==2)  then
272      ABI_ALLOCATE(aa,(5))
273      ABI_ALLOCATE(bb,(5))
274    end if
275    if (opt_dens==0.or.opt_dens==1) then
276      do ll=0,pawtab%lcut_size-1
277        do ilm=ll**2+1,(ll+1)**2
278          if (lmselectin(ilm)) then
279            if (cplex==1) then
280              call pawrad_deducer0(rho1 (:,ilm,ispden),mesh_size,pawrad)
281              call pawrad_deducer0(trho1(:,ilm,ispden),mesh_size,pawrad)
282            else
283              do ii=0,1
284                do ir=2,5
285                  aa(ir)=rho1 (2*ir-ii,ilm,ispden)
286                  bb(ir)=trho1(2*ir-ii,ilm,ispden)
287                end do
288                call pawrad_deducer0(aa,5,pawrad)
289                call pawrad_deducer0(bb,5,pawrad)
290                rho1 (2-ii,ilm,ispden)=aa(1)
291                trho1(2-ii,ilm,ispden)=bb(1)
292              end do
293            end if
294          end if
295        end do
296      end do
297    else
298      do ll=0,pawtab%lcut_size-1
299        do ilm=ll**2+1,(ll+1)**2
300          if (lmselectin(ilm)) then
301            if (cplex==1) then
302              call pawrad_deducer0(rho1(:,ilm,ispden),mesh_size,pawrad)
303            else
304              do ii=0,1
305                do ir=2,5
306                  aa(ir)=rho1 (2*ir-ii,ilm,ispden)
307                end do
308                call pawrad_deducer0(aa,5,pawrad)
309                rho1(2-ii,ilm,ispden)=aa(1)
310              end do
311            end if
312          end if
313        end do
314      end do
315    end if
316    if (cplex==2)  then
317      ABI_DEALLOCATE(aa)
318      ABI_DEALLOCATE(bb)
319    end if
320 
321 !  -- Test moments of densities and store non-zero ones
322    if (nzlmopt==-1) then
323      do ll=0,pawtab%lcut_size-1
324        do ilm=ll**2+1,(ll+1)**2
325          m1=zero;mt1=zero
326          if (cplex==1) then
327            m1=maxval(abs(rho1 (1:mesh_size,ilm,ispden)))
328            if (opt_dens<2) mt1=maxval(abs(trho1(1:mesh_size,ilm,ispden)))
329          else
330            do ir=1,mesh_size
331              rdum=sqrt(rho1(2*ir-1,ilm,ispden)**2+rho1(2*ir,ilm,ispden)**2)
332              m1=max(m1,rdum)
333            end do
334            if (opt_dens<2) then
335              do ir=1,mesh_size
336                rdum=sqrt(trho1(2*ir-1,ilm,ispden)**2+trho1(2*ir,ilm,ispden)**2)
337                mt1=max(mt1,rdum)
338              end do
339            end if
340          end if
341          if (ispden==1) then
342            if ((ilm>1).and.(m1<tol16).and.(mt1<tol16)) then
343              lmselectout(ilm)=.false.
344            end if
345          else if (.not.(lmselectout(ilm))) then
346            lmselectout(ilm)=((m1>=tol16).or.(mt1>=tol16))
347          end if
348        end do
349      end do
350    end if
351 
352 !  -- Compute integral of (n1-tn1) inside spheres
353    if (opt_compch==1.and.ispden==1.and.opt_dens<2) then
354      ABI_ALLOCATE(aa,(mesh_size))
355      aa(1:mesh_size)=(rho1(1:mesh_size,1,1)-trho1(1:mesh_size,1,1)) &
356 &     *pawrad%rad(1:mesh_size)**2
357 ! jmb 
358 !    compchspha(1)=zero
359      call simp_gen(compchspha(1),aa,pawrad)
360      compch_sph=compch_sph+compchspha(1)*sqrt(four_pi)
361      ABI_DEALLOCATE(aa)
362    end if
363 
364 !  -- Print out moments of densities (if requested)
365    if (abs(pawprtvol)>=2.and.opt_print==1.and.opt_dens<2) then
366      ABI_ALLOCATE(aa,(cplex*mesh_size))
367      ABI_ALLOCATE(bb,(cplex*mesh_size))
368      if (opt_dens==0) then
369        write(msg,'(2a,i3,a,i1,3a)') ch10, &
370 &       ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
371 &       '  ******* Moment of (n1-tn1)  ** Moment of (n1-tn1-nhat1)'
372      else
373        write(msg,'(2a,i3,a,i1,3a)') ch10, &
374 &       ' Atom ',iatom,' (ispden=',ispden,'):',ch10,&
375 &       '  ******* Moment of (n1-tn1)'
376      end if
377      call wrtout(std_out,msg,'PERS')
378      do ll=0,pawtab%lcut_size-1
379        do ilm=ll**2+1,(ll+1)**2
380          if (lmselectin(ilm)) then
381            do iplex=1,cplex
382              if (opt_dens==0) then
383                do ir=1,mesh_size
384                  ii=cplex*(ir-1)+iplex
385                  ro(1)=pawrad%rad(ir)**(2+ll)
386                  aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
387                  bb(ir)=ro(1)*nhat1(ii,ilm,ispden)
388                end do
389                call simp_gen(compchspha(iplex),aa,pawrad)
390                call simp_gen(compchsphb(iplex),bb,pawrad)
391              else
392                do ir=1,mesh_size
393                  ii=cplex*(ir-1)+iplex
394                  ro(1)=pawrad%rad(ir)**(2+ll)
395                  aa(ir)=ro(1)*(rho1(ii,ilm,ispden)-trho1(ii,ilm,ispden))
396                end do
397                call simp_gen(compchspha(iplex),aa,pawrad)
398              end if
399            end do
400            if (opt_dens==0) then
401              if (cplex==1) then
402                write(msg,'(3x,a,2i2,2(a,es14.7))') &
403 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1),&
404 &               ' **    M=',compchspha(1)-compchsphb(1)
405              else
406                write(msg,'(3x,a,2i2,2(a,2es14.7))') &
407 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2),&
408 &               ' **    M=',compchspha(1:2)-compchsphb(1:2)
409              end if
410            else
411              if (cplex==1) then
412                write(msg,'(3x,a,2i2,a,es14.7)') &
413 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1)
414              else
415                write(msg,'(3x,a,2i2,a,2es14.7)') &
416 &               'l,m=',ll,ilm-(ll**2+ll+1),': M=',compchspha(1:2)
417              end if
418            end if
419            call wrtout(std_out,msg,'PERS')
420          end if
421        end do
422      end do
423      ABI_DEALLOCATE(aa)
424      ABI_DEALLOCATE(bb)
425    end if
426 
427 !  ----- End loop over spin components
428  end do
429 
430  if (.not.present(one_over_rad2))  then
431    ABI_DEALLOCATE(one_over_rad2_)
432  end if
433 
434  DBG_EXIT("COLL")
435 
436 end subroutine pawdensities