TABLE OF CONTENTS


ABINIT/m_paw_atom [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_atom

FUNCTION

  atompaw related operations

COPYRIGHT

  Copyright (C) 2012-2022 ABINIT group (T. Rangel, MT, JWZ, GJ)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt

SOURCE

20 #include "libpaw.h"
21 
22 module m_paw_atom
23 
24  USE_DEFS
25  USE_MSG_HANDLING
26  USE_MEMORY_PROFILING
27 
28  use m_paw_numeric, only : paw_jbessel, paw_solvbes, paw_spline, paw_splint
29  use m_pawrad,      only : pawrad_type, simp_gen, poisson, pawrad_deducer0, bound_deriv, pawrad_ifromr
30  use m_pawtab,      only : pawtab_type
31 
32  implicit none
33 
34  private
35 
36  public:: atompaw_shpfun
37  public:: atompaw_shapebes
38  public:: atompaw_vhnzc
39  public:: atompaw_dij0
40  public:: atompaw_kij

m_paw_atom/atompaw_atompaw_shapebes [ Functions ]

[ Top ] [ m_paw_atom ] [ Functions ]

NAME

 atompaw_shapebes

FUNCTION

    Find al and ql parameters for a "Bessel" shape function:
    Shape(r)=al1.jl(ql1.r)+al2.jl(ql2.r)
      such as Shape(r) and 2 derivatives are zero at r=rc
              Intg_0_rc[Shape(r).r^(l+2).dr]=1

INPUTS

  ll= l quantum number
  rc= cut-off radius

OUTPUT

  al(2)= al coefficients
  ql(2)= ql factors

SOURCE

189  subroutine atompaw_shapebes(al,ql,ll,rc)
190 
191 !Arguments ------------------------------------
192 !scalars
193  integer :: ll
194  real(dp) :: rc
195 !arrays
196  real(dp) :: al(2),ql(2)
197 
198 !Local variables-------------------------------
199 !scalars
200  integer :: ii
201  real(dp) :: alpha,beta,det,jbes,jbesp,jbespp,qr
202 !arrays
203  real(dp) :: amat(2,2),bb(2)
204 
205 ! *************************************************************************
206 
207  alpha=1._dp;beta=0._dp
208  call paw_solvbes(ql,alpha,beta,ll,2)
209  ql(1:2)=ql(1:2)/rc
210 
211  do ii=1,2
212    qr=ql(ii)*rc
213    call paw_jbessel(jbes,jbesp,jbespp,ll,1,qr)
214    amat(1,ii)=jbesp*ql(ii)
215    call paw_jbessel(jbes,jbesp,jbespp,ll+1,0,qr)
216    amat(2,ii)=jbes*rc**(ll+2)/ql(ii)  !  Intg_0_rc[jl(qr).r^(l+2).dr]
217  end do
218 
219  bb(1)=zero;bb(2)=one
220 
221  det=amat(1,1)*amat(2,2)-amat(1,2)*amat(2,1)
222  al(1)=(amat(2,2)*bb(1)-amat(1,2)*bb(2))/det
223  al(2)=(amat(1,1)*bb(2)-amat(2,1)*bb(1))/det
224 
225 end subroutine atompaw_shapebes

m_paw_atom/atompaw_dij0 [ Functions ]

[ Top ] [ m_paw_atom ] [ Functions ]

NAME

 atompaw_dij0

FUNCTION

  PAW: Compute "frozen" values of pseudopotential strengths Dij = Dij0

INPUTS

  indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn
  kij(pawtab%lmn2_size)= kinetic part of Dij
  lmnmax=max number of (l,m,n) components over all type of psps
  ncore(:)=atomic core density
  opt_init=flag defining the storage of PAW atomic data
           0: PAW atomic data have not been initialized (in pawtab)
           1: PAW atomic data have been initialized (in pawtab)
  pawtab <type(pawtab_type)>=paw tabulated starting data
  radmesh <type(pawrad_type)>=paw radial mesh (and related data)
  radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities
  radmesh_vloc <type(pawrad_type)>=radial mesh (and related data) for the local potential (VH(tnZc))
  vhtnzc(:)= local potential VH(tnZc)
  znucl= valence and total charge of the atomic species

OUTPUT

  pawtab%dij0(pawtab%lmn2_size)= Frozen part of the Dij term

SOURCE

310  subroutine atompaw_dij0(indlmn,kij,lmnmax,ncore,opt_init,pawtab,&
311 &                        radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl)
312 
313 !Arguments ---------------------------------------------
314 !scalars
315  integer,intent(in) :: lmnmax,opt_init
316  real(dp),intent(in) :: znucl
317  type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc
318  type(pawtab_type),intent(inout) :: pawtab
319 !arrays
320  integer,intent(in) :: indlmn(6,lmnmax)
321  real(dp),intent(in) :: kij(pawtab%lmn2_size)
322  real(dp),intent(in) :: ncore(:),vhtnzc(:)
323 !real(dp),optional,intent(in) :: vminushalf(:)
324 
325 !Local variables ---------------------------------------
326  integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size,meshsz,meshsz_core
327  integer :: meshsz_vhtnzc,meshsz_vmh
328  real(dp) :: intg,intvh,yp1,ypn
329  real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:)
330 
331 ! *********************************************************************
332 
333  lmn2_size=pawtab%lmn2_size
334  meshsz_vhtnzc=size(vhtnzc)
335  meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc)
336  LIBPAW_ALLOCATE(ff,(meshsz))
337 
338 !Retrieve VH(tnZc) on the correct radial mesh
339  LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz))
340  if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.&
341 &    (radmesh%rstep    /=radmesh_vloc%rstep)    .or.&
342 &    (radmesh%lstep    /=radmesh_vloc%lstep)) then
343    call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn)
344    LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc))
345    LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc))
346    call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1)
347    call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph)
348    LIBPAW_DEALLOCATE(work1)
349    LIBPAW_DEALLOCATE(work2)
350  else
351    vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz)
352  end if
353 
354 !Kinetic part of Dij0
355 !====================
356  pawtab%dij0(1:lmn2_size)=kij(1:lmn2_size)
357 
358 !Computation of <phi_i|vh(nZc)|phi_j> on the PAW sphere
359 !======================================================
360  meshsz_core=size(ncore)
361  LIBPAW_ALLOCATE(vhnzc,(meshsz_core))
362  call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
363  do jlmn=1,pawtab%lmn_size
364    j0lmn=jlmn*(jlmn-1)/2
365    jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
366    do ilmn=1,jlmn
367      klmn=j0lmn+ilmn
368      ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
369      if (jlm==ilm) then
370        ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz)
371        call simp_gen(intg,ff,radmesh)
372        pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
373      end if
374    end do
375  end do
376  LIBPAW_DEALLOCATE(vhnzc)
377 
378 !Computation of -<tphi_i|vh(tnZc)|tphi_j> on the PAW sphere
379 !==========================================================
380  do jlmn=1,pawtab%lmn_size
381    j0lmn=jlmn*(jlmn-1)/2
382    jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
383    do ilmn=1,jlmn
384      klmn=j0lmn+ilmn
385      ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
386      if (jlm==ilm) then
387        ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz)
388        call simp_gen(intg,ff,radmesh)
389        pawtab%dij0(klmn)=pawtab%dij0(klmn)-intg
390      end if
391    end do
392  end do
393 
394 !Computation of <phi_i|vminushalf|phi_j>  (if any)
395 !=================================================
396  if(pawtab%has_vminushalf==1) then
397    if(size(pawtab%vminushalf)>=1) then
398      meshsz_vmh=min(meshsz,size(pawtab%vminushalf))
399      do jlmn=1,pawtab%lmn_size
400        j0lmn=jlmn*(jlmn-1)/2
401        jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
402        do ilmn=1,jlmn
403          klmn=j0lmn+ilmn
404          ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
405          if (jlm==ilm) then
406            ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh)
407            call simp_gen(intg,ff(1:meshsz_vmh),radmesh)
408            pawtab%dij0(klmn)=pawtab%dij0(klmn)+intg
409          end if
410        end do
411      end do
412    end if
413  end if
414 
415 !Computation of -int[vh(tnzc)*Qijhat(r)dr]
416 !=========================================
417  if (opt_init==0) then
418    LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size))
419    call atompaw_shpfun(0,radmesh,intg,pawtab,shpf)
420    if (pawtab%shape_type==3) then
421      LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz))
422      r2k=zero
423      r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2
424      if(radmesh%mesh_type==5) then
425        call simp_gen(intg,r2k,radmesh)
426      else
427        call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp)
428      end if
429      shpf(1:meshsz)=shpf(1:meshsz)/intg
430      LIBPAW_DEALLOCATE(r2k)
431    end if
432    ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2
433    LIBPAW_DEALLOCATE(shpf)
434    call simp_gen(intvh,ff,radmesh)
435    do jlmn=1,pawtab%lmn_size
436      j0lmn=jlmn*(jlmn-1)/2
437      jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
438      do ilmn=1,jlmn
439        klmn=j0lmn+ilmn
440        il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
441        if (ilm==jlm) then
442          ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)&
443 &                     -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln))
444          call simp_gen(intg,ff,radmesh)
445          pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg
446        end if
447      end do
448    end do
449  else
450    ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2
451    call simp_gen(intvh,ff,radmesh)
452    do jlmn=1,pawtab%lmn_size
453      j0lmn=jlmn*(jlmn-1)/2
454      jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
455      do ilmn=1,jlmn
456        klmn=j0lmn+ilmn
457        il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
458        if (ilm==jlm) then
459          intg=pawtab%qijl(1,klmn)*sqrt(four_pi)
460          pawtab%dij0(klmn)=pawtab%dij0(klmn)-intvh*intg
461        end if
462      end do
463    end do
464  end if
465 
466  LIBPAW_DEALLOCATE(ff)
467  LIBPAW_DEALLOCATE(vhtnzc_sph)
468 
469  end subroutine atompaw_dij0

m_paw_atom/atompaw_kij [ Functions ]

[ Top ] [ m_paw_atom ] [ Functions ]

NAME

 atompaw_kij

FUNCTION

 PAW: deduce kinetic part of psp strength (Dij) from the knowledge of frozen Dij (Dij0)

INPUTS

  indlmn(6,lmnmax)= array giving l,m,n,lm,ln,s for i=lmn
  lmnmax=max number of (l,m,n) components over all type of psps
  ncore(:)=atomic core density
  opt_init=flag defining the storage of PAW atomic data
           0: PAW atomic data have not been initialized (in pawtab)
           1: PAW atomic data have been initialized (in pawtab)
  opt_vhnzc=flag defining the inclusion of VH(nZc) in computation
            0: VH(nZc) is not taken into account
            1: VH(nZc) is taken into account
  pawtab <type(pawtab_type)>=paw tabulated starting data
  radmesh <type(pawrad_type)>=paw radial mesh (and related data)
  radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities
  radmesh_vloc <type(pawrad_type)>=radial mesh (and related data) for the local potential (VH(tnZc))
  vhtnzc(:)= local potential VH(tnZc)
  znucl= valence and total charge of the atomic species

OUTPUT

  kij(pawtab%lmn2_size)= kinetic part of Dij

SOURCE

503  subroutine atompaw_kij(indlmn,kij,lmnmax,ncore,opt_init,opt_vhnzc,pawtab, &
504 &                       radmesh,radmesh_core,radmesh_vloc,vhtnzc,znucl)
505 
506 !Arguments ---------------------------------------------
507 !scalars
508  integer,intent(in) :: lmnmax,opt_init,opt_vhnzc
509  real(dp),intent(in) :: znucl
510  type(pawrad_type),intent(in) :: radmesh,radmesh_core,radmesh_vloc
511  type(pawtab_type),intent(in) :: pawtab
512 !arrays
513  integer,intent(in) :: indlmn(6,lmnmax)
514  real(dp),intent(out) :: kij(pawtab%lmn2_size)
515  real(dp),intent(in) :: ncore(:)
516  real(dp),intent(in) :: vhtnzc(:)
517 
518 !Local variables ---------------------------------------
519  integer :: il,ilm,iln,ilmn,j0lmn,jl,jlm,jln,jlmn,klmn,lmn2_size
520  integer :: meshsz,meshsz_core,meshsz_vhtnzc,meshsz_vmh
521  real(dp) :: intg,intvh,yp1,ypn
522  real(dp),allocatable :: ff(:),r2k(:),shpf(:),vhnzc(:),vhtnzc_sph(:),work1(:),work2(:)
523 
524 ! *********************************************************************
525 
526  lmn2_size=pawtab%lmn2_size
527  meshsz_vhtnzc=size(vhtnzc)
528  meshsz=min(radmesh%mesh_size,radmesh_core%mesh_size,radmesh_vloc%mesh_size,meshsz_vhtnzc)
529  LIBPAW_ALLOCATE(ff,(meshsz))
530 
531 !Retrieve VH(tnZc) on the correct radial mesh
532  LIBPAW_ALLOCATE(vhtnzc_sph,(meshsz))
533  if ((radmesh%mesh_type/=radmesh_vloc%mesh_type).or.&
534 &    (radmesh%rstep    /=radmesh_vloc%rstep)    .or.&
535 &    (radmesh%lstep    /=radmesh_vloc%lstep)) then
536    call bound_deriv(vhtnzc,radmesh_vloc,meshsz_vhtnzc,yp1,ypn)
537    LIBPAW_ALLOCATE(work1,(meshsz_vhtnzc))
538    LIBPAW_ALLOCATE(work2,(meshsz_vhtnzc))
539    call paw_spline(radmesh_vloc%rad,vhtnzc,meshsz_vhtnzc,yp1,ypn,work1)
540    call paw_splint(meshsz_vhtnzc,radmesh_vloc%rad,vhtnzc,work1,meshsz,radmesh%rad(1:meshsz),vhtnzc_sph)
541    LIBPAW_DEALLOCATE(work1)
542    LIBPAW_DEALLOCATE(work2)
543  else
544    vhtnzc_sph(1:meshsz)=vhtnzc(1:meshsz)
545  end if
546 
547 !Initialize Kij with Dij0
548 !=========================
549  kij(1:lmn2_size)=pawtab%dij0(1:lmn2_size)
550 
551 !Substraction of -<phi_i|vh(nZc)|phi_j> on the PAW sphere
552 !========================================================
553  if (opt_vhnzc/=0) then
554    meshsz_core=size(ncore)
555    LIBPAW_ALLOCATE(vhnzc,(meshsz_core))
556    call atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
557    do jlmn=1,pawtab%lmn_size
558      j0lmn=jlmn*(jlmn-1)/2
559      jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
560      do ilmn=1,jlmn
561        klmn=j0lmn+ilmn
562        ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
563        if (jlm==ilm) then
564          ff(1:meshsz)=pawtab%phi(1:meshsz,iln)*pawtab%phi(1:meshsz,jln)*vhnzc(1:meshsz)
565          call simp_gen(intg,ff,radmesh)
566          kij(klmn)=kij(klmn)-intg
567        end if
568      end do
569    end do
570    LIBPAW_DEALLOCATE(vhnzc)
571  end if
572 
573 !Substraction of <tphi_i|vh(tnZc)|tphi_j> on the PAW sphere
574 !==========================================================
575  do jlmn=1,pawtab%lmn_size
576    j0lmn=jlmn*(jlmn-1)/2
577    jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
578    do ilmn=1,jlmn
579      klmn=j0lmn+ilmn
580      ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
581      if (jlm==ilm) then
582        ff(1:meshsz)=pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln)*vhtnzc_sph(1:meshsz)
583        call simp_gen(intg,ff,radmesh)
584        kij(klmn)=kij(klmn)+intg
585      end if
586    end do
587  end do
588 
589 !Computation of <phi_i|vminushalf|phi_j>  (if any)
590 !=================================================
591  if(pawtab%has_vminushalf==1) then
592    if(size(pawtab%vminushalf)>=1) then
593      meshsz_vmh=min(meshsz,size(pawtab%vminushalf))
594      do jlmn=1,pawtab%lmn_size
595        j0lmn=jlmn*(jlmn-1)/2
596        jlm=indlmn(4,jlmn);jln=indlmn(5,jlmn)
597        do ilmn=1,jlmn
598          klmn=j0lmn+ilmn
599          ilm=indlmn(4,ilmn);iln=indlmn(5,ilmn)
600          if (jlm==ilm) then
601            ff(1:meshsz_vmh)=pawtab%phi(1:meshsz_vmh,iln)*pawtab%phi(1:meshsz_vmh,jln)*pawtab%vminushalf(1:meshsz_vmh)
602            call simp_gen(intg,ff(1:meshsz_vmh),radmesh)
603            kij(klmn)=kij(klmn)-intg
604          end if
605        end do
606      end do
607    end if
608  end if
609 
610 !Computation of int[vh(tnzc)*Qijhat(r)dr]
611 !==========================================
612  if (opt_init==0) then
613    LIBPAW_ALLOCATE(shpf,(radmesh%mesh_size))
614    call atompaw_shpfun(0,radmesh,intg,pawtab,shpf)
615    if (pawtab%shape_type==3) then
616      LIBPAW_ALLOCATE(r2k,(radmesh%int_meshsz))
617      r2k=zero
618      r2k(2:radmesh%int_meshsz)=shpf(2:radmesh%int_meshsz)*radmesh%rad(2:radmesh%int_meshsz)**2
619      if(radmesh%mesh_type==5) then
620        call simp_gen(intg,r2k,radmesh)
621      else
622        call simp_gen(intg,r2k,radmesh,r_for_intg=pawtab%rshp)
623      end if
624      shpf(1:meshsz)=shpf(1:meshsz)/intg
625      LIBPAW_DEALLOCATE(r2k)
626    end if
627    ff(1:meshsz)=vhtnzc_sph(1:meshsz)*shpf(1:meshsz)*radmesh%rad(1:meshsz)**2
628    LIBPAW_DEALLOCATE(shpf)
629    call simp_gen(intvh,ff,radmesh)
630    do jlmn=1,pawtab%lmn_size
631      j0lmn=jlmn*(jlmn-1)/2
632      jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
633      do ilmn=1,jlmn
634        klmn=j0lmn+ilmn
635        il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
636        if (ilm==jlm) then
637          ff(1:meshsz)=(pawtab%phi (1:meshsz,iln)*pawtab%phi (1:meshsz,jln)&
638 &                     -pawtab%tphi(1:meshsz,iln)*pawtab%tphi(1:meshsz,jln))
639          call simp_gen(intg,ff,radmesh)
640          kij(klmn)=kij(klmn)+intvh*intg
641        end if
642      end do
643    end do
644  else
645    ff(1:meshsz)=vhtnzc_sph(1:meshsz)*pawtab%shapefunc(1:meshsz,1)*radmesh%rad(1:meshsz)**2
646    call simp_gen(intvh,ff,radmesh)
647    do jlmn=1,pawtab%lmn_size
648      j0lmn=jlmn*(jlmn-1)/2
649      jl=indlmn(1,jlmn);jln=indlmn(5,jlmn);jlm=indlmn(4,jlmn)
650      do ilmn=1,jlmn
651        klmn=j0lmn+ilmn
652        il=indlmn(1,ilmn);iln=indlmn(5,ilmn);ilm=indlmn(4,ilmn)
653        if (ilm==jlm) then
654          intg=pawtab%qijl(1,klmn)*sqrt(four_pi)
655          kij(klmn)=kij(klmn)+intvh*intg
656        end if
657      end do
658    end do
659  end if
660 
661  LIBPAW_DEALLOCATE(ff)
662  LIBPAW_DEALLOCATE(vhtnzc_sph)
663 
664  end subroutine atompaw_kij

m_paw_atom/atompaw_shpfun [ Functions ]

[ Top ] [ m_paw_atom ] [ Functions ]

NAME

 atompaw_shpfun

FUNCTION

 Compute shape function used in the definition
 of compensation density (PAW)

INPUTS

  ll= l quantum number
  mesh <type(pawrad_type)>=data containing radial grid information
  pawtab <type(pawtab_type)>=paw tabulated starting data

OUTPUT

  norm= factor for shape function normalization

 SIDE effects
  shapefunc(:)=shape function g(r)
    In case of numerical shape function (shape_type=-1), shapefunc
    array contains the shape function read in psp file at input.

NOTES

  Types of shape functions:
   type -1: numerical shape function, given in psp file
   type  1: g(r)=k(r).r^l; k(r)=exp(-(r/sigma)^lambda)
   type  2: g(r)=k(r).r^l; k(r)=[sin(Pi.r/rshp)/(Pi.r/rshp)]^2
   type  3: g(r)=alpha1.jl(q1.r)+alpha2.jl(q2.r)

SOURCE

 76 subroutine atompaw_shpfun(ll,mesh,norm,pawtab,shapefunc)
 77 
 78 !Arguments ---------------------------------------------
 79 !scalars
 80  integer,intent(in) :: ll
 81  real(dp),intent(out) :: norm
 82  type(pawrad_type),intent(in) :: mesh
 83  type(pawtab_type),intent(in) :: pawtab
 84 !arrays
 85  real(dp),intent(inout) :: shapefunc(:)
 86 
 87 !Local variables ------------------------------
 88 !scalars
 89  integer :: ir,ishp,mesh_size
 90  real(dp) :: arg,besp,bespp,jbes1,jbes2
 91 !arrays
 92  real(dp) :: alpha(2),qq(2)
 93  real(dp),allocatable :: r2k(:)
 94 !no_abirules
 95 
 96 !***************************************************************************
 97 
 98  mesh_size=size(shapefunc)
 99  if (mesh_size>mesh%mesh_size) then
100    LIBPAW_BUG('wrong size!')
101  end if
102 
103 !Index for shape function cut-off radius
104  ishp=pawrad_ifromr(mesh,pawtab%rshp)-1
105 
106 !Computation of non-normalized shape function
107  if (pawtab%shape_type==-1) then
108    shapefunc(1:ishp)=pawtab%shapefunc(1:ishp,1+ll)
109  else if (pawtab%shape_type==1) then
110    if (ll==0) then
111      shapefunc(1)=one
112      do ir=2,ishp
113        shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda)
114      end do
115    else
116      shapefunc(1)=zero
117      do ir=2,ishp
118        shapefunc(ir)=exp(-(mesh%rad(ir)/pawtab%shape_sigma)**pawtab%shape_lambda)*mesh%rad(ir)**ll
119      end do
120    end if
121  else if (pawtab%shape_type==2) then
122    if (ll==0) then
123      shapefunc(1)=one
124      do ir=2,ishp
125        arg=pi*mesh%rad(ir)/pawtab%rshp
126        shapefunc(ir)=(sin(arg)/arg)**2
127      end do
128    else
129      shapefunc(1)=zero
130      do ir=2,ishp
131        arg=pi*mesh%rad(ir)/pawtab%rshp
132        shapefunc(ir)=(sin(arg)/(arg))**2 *mesh%rad(ir)**ll
133      end do
134    end if
135  else if (pawtab%shape_type==3) then
136    alpha(1:2)=pawtab%shape_alpha(1:2,1+ll)
137    qq(1:2)=pawtab%shape_q(1:2,1+ll)
138    do ir=1,ishp
139      call paw_jbessel(jbes1,besp,bespp,ll,0,qq(1)*mesh%rad(ir))
140      call paw_jbessel(jbes2,besp,bespp,ll,0,qq(2)*mesh%rad(ir))
141      shapefunc(ir)=alpha(1)*jbes1+alpha(2)*jbes2
142    end do
143  end if
144 
145  if (ishp<mesh_size) shapefunc(ishp+1:mesh_size)=zero
146 
147 !Shape function normalization
148  if (pawtab%shape_type==-1.or.pawtab%shape_type==1.or.pawtab%shape_type==2) then
149    LIBPAW_ALLOCATE(r2k,(mesh_size))
150    r2k=zero
151    r2k(2:ishp)=shapefunc(2:ishp)*mesh%rad(2:ishp)**(2+ll)
152    if (mesh%mesh_type==5) then
153      call simp_gen(norm,r2k,mesh);norm=one/norm
154    else
155      call simp_gen(norm,r2k,mesh,r_for_intg=pawtab%rshp);norm=one/norm
156    end if
157    shapefunc(1:ishp)=shapefunc(1:ishp)*norm
158    if (pawtab%shape_type==-1) norm=one
159    LIBPAW_DEALLOCATE(r2k)
160  else if (pawtab%shape_type==3) then
161    norm=one
162  end if
163 
164 end subroutine atompaw_shpfun

m_paw_atom/atompaw_vhnzc [ Functions ]

[ Top ] [ m_paw_atom ] [ Functions ]

NAME

 atompaw_vhnzc

FUNCTION

 PAW: compute Hartree potential for n_{Zc}

INPUTS

  ncore(:)=atomic core density
  radmesh_core <type(pawrad_type)>=radial mesh (and related data) for the core densities
  znucl= valence and total charge of the atomic species

OUTPUT

  vhnzc(:)=Hartree potential due to Z_nc

SOURCE

247  subroutine atompaw_vhnzc(ncore,radmesh_core,vhnzc,znucl)
248 
249 !Arguments ---------------------------------------------
250 !scalars
251  real(dp),intent(in) :: znucl
252  type(pawrad_type),intent(in) :: radmesh_core
253 !arrays
254  real(dp),intent(in) :: ncore(:)
255  real(dp), intent(out) :: vhnzc(:)
256 
257 !Local variables ---------------------------------------
258  integer :: mesh_size
259  real(dp),allocatable :: nwk(:)
260 
261 ! *********************************************************************
262 
263  mesh_size=size(ncore)
264  if (mesh_size/=size(vhnzc).or.mesh_size>radmesh_core%mesh_size) then
265    LIBPAW_BUG('wrong sizes!')
266  end if
267 
268  LIBPAW_ALLOCATE(nwk,(mesh_size))
269 
270  nwk(:)=ncore(:)*four_pi*radmesh_core%rad(:)**2
271  call poisson(nwk,0,radmesh_core,vhnzc)
272  vhnzc(2:mesh_size)=(vhnzc(2:mesh_size)-znucl)/radmesh_core%rad(2:mesh_size)
273  call pawrad_deducer0(vhnzc,mesh_size,radmesh_core)
274 
275  LIBPAW_DEALLOCATE(nwk)
276 
277  end subroutine atompaw_vhnzc