TABLE OF CONTENTS


ABINIT/m_paw_atomorb [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_atomorb

FUNCTION

  This module provides the definition of the atomorb_type used
  to store atomic orbitals on a radial mesh as well
  as methods to operate on it.

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG)
 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

19 #if defined HAVE_CONFIG_H
20 #include "config.h"
21 #endif
22 
23 #include "abi_common.h"
24 
25 MODULE m_paw_atomorb
26 
27  use defs_basis
28  use m_abicore
29  use m_errors
30  use m_splines
31 
32  use m_io_tools,      only : open_file
33  use m_fstrings,      only : tolower
34  use m_paw_lmn,       only : make_indlmn, make_indklmn, make_kln2ln
35  use m_pawrad,        only : pawrad_type, pawrad_init, &
36 &                            pawrad_free, pawrad_print, pawrad_isame, pawrad_ifromr, simp_gen
37 
38  implicit none
39 
40  private

m_paw_atomorb/atomorb_type [ Types ]

[ Top ] [ m_paw_atomorb ] [ Types ]

NAME

FUNCTION

  Defines the atomorb_type datastructure type.
  It contains the atomic orbitals on a radial mesh for a given type of atom.

NOTES

  * Should the radial mesh included in the object or not?
  * We might have different meshes, useful for deep states in heavy atoms!
  * Methods to be added: corekin

SOURCE

 57  type, public :: atomorb_type
 58 
 59 !scalars
 60   integer :: ixc
 61    ! Exchange and correlation functional used to generate the orbitals
 62 
 63   integer :: method
 64   ! 1 for restricted, compatible only with nsppol=1.
 65   ! 2 for spin unrestricted, compatible only with nsppol=2.
 66 
 67   integer :: nspden
 68    ! Number of spin-density components.
 69 
 70   integer :: nsppol
 71    ! Number of independent spin-components.
 72    ! FIXME: here a lot of quantities might depend on nsppol in the
 73    ! case of magnetic atoms!
 74 
 75   integer :: nspinor
 76    ! Number of spinorial components
 77    ! TODO this is a quite delicate issue, for S.O. one should use J = L+S instead of L!
 78    ! If we use scalar relativistic then...
 79 
 80   integer :: l_max
 81    ! Maximum value of angular momentum l+1
 82 
 83   integer :: l_size
 84   ! Maximum value of l+1 leading to non zero Gaunt coeffs
 85   ! l_size=2*l_max-1
 86 
 87   integer :: ln_size
 88   ! Number of (l,n) components.
 89 
 90   integer :: ln2_size
 91   ! ln2_size=ln_size*(ln_size+1)/2
 92   ! where ln_size is the number of (l,n) elements for core orbitals.
 93 
 94   integer :: lmn_size
 95   ! Number of (l,m,n) elements.
 96 
 97   integer :: lmn2_size
 98    ! lmn2_size=lmn_size*(lmn_size+1)/2
 99    ! where lmn_size is the number of (l,m,n) elements for core orbitals.
100 
101   integer :: mesh_size
102   ! Dimension of the radial mesh.
103 
104   real(dp) :: rcore
105   ! Radius of the sphere used to describe core electrons.
106   ! It should be <= rpaw
107 
108   real(dp) :: zion
109    ! zionpsp
110    ! The ionic pseudo-charge, (giving raise to a long-range coulomb potential)
111 
112   ! TODO alchemy?
113   !real(dp) :: ziontypat
114    ! ziontypat
115    !  For each type of atom (might be alchemy wrt psps), the ionic pseudo-charge
116    ! (giving raise to a long-range coulomb potential)
117 
118   real(dp) :: znucl
119    ! The atomic number of the atom.
120 
121   ! TODO alchemy?
122   !real(dp) :: znucltypat
123    ! znucltypat
124    ! The atomic number of each type of atom (might be alchemy wrt psps)
125 
126   character(len=fnlen) :: fname
127    ! The filename for temporary storage.
128 
129 !arrays
130   integer, allocatable :: indlmn(:,:)
131   ! indlmn(6,lmn_size)
132   ! Array giving l,m,n,lm,ln,spin for i=lmn.
133 
134   integer, allocatable :: indln(:,:)
135   ! indln(2,ln_size)
136   ! Array giving l and n for i=ln
137 
138   integer, allocatable :: indklmn(:,:)
139    ! indklmn(8,lmn2_size)
140    ! Array giving klm, kln, abs(il-jl), (il+jl), ilm and jlm, ilmn and jlmn for each klmn=(ilmn,jlmn)
141    ! Note: ilmn=(il,im,in) and ilmn<=jlmn
142 
143   !integer, allocatable :: klm2lm TODO add
144    !  klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm, jlm for each klm=(ilm,jlm)
145    !  where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection.
146 
147   integer, allocatable :: klm_diag(:)
148    ! klm_diag(lmn2_size)
149    ! 1 il==jl and im==jm, 0 otherwise.
150 
151   integer, allocatable :: klmntomn(:,:) 
152    ! klmntomn(4,lmn2_size)
153    ! Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn)
154    ! Note: ilmn=(il,im,in) and ilmn<=jlmn
155    ! NB: klmntomn is an application and not a bijection
156 
157   integer, allocatable :: kln2ln(:,:) 
158    ! kln2ln(6,ln2_size)
159    ! Table giving il, jl ,in, jn, iln, jln for each kln=(iln,jln)
160    ! where iln=(il,in) and iln<=jln. NB: kln2ln is an application and not a bijection
161 
162   integer, allocatable :: mode(:,:) 
163   ! mode(ln_size,nsppol)
164   ! Flag defining how the orbital is treated.
165   ! During the pseudopotential generation we can have: ORB_FROZEN or ORB_VALENCE
166   ! For calculations in extended systems we can have:  ORB_FROZEN or ORB_RELAXED_CORE
167   ! Namely different treatment depending of the degree of localization.
168   ! For example the 1s in plutonium might be treated as ORB_FROZEN during
169   ! a relaxed core calculation.
170   ! TODO define function to test the type, much safer!
171 
172   real(dp), allocatable :: eig(:,:) 
173   ! eig(ln_size,nsppol)
174   ! Eigenvalues for each ln channel and spin.
175 
176   real(dp), allocatable :: occ(:,:)
177   ! occ(ln_size,nsppol)
178   ! Occupation for each ln channel and spin.
179 
180   real(dp), allocatable :: phi(:,:,:)
181   ! phi(mesh_size,ln_size,nsppol)
182   ! Here we might have different meshes, useful for deep states in heavy atoms!
183 
184   ! this might be retrieved with  a method get_atomden
185   !real(dp), allocatable :: density(:)
186    ! density(mesh_size,nspden)
187    ! Gives the core density of the atom for each spin channel
188    ! Total charge in first dimension,up component in second one (if present)
189 
190  end type atomorb_type
191 
192 ! public procedures.
193  public :: init_atomorb
194  public :: destroy_atomorb
195  public :: print_atomorb
196  public :: get_overlap

m_paw_atomorb/destroy_atomorb [ Functions ]

[ Top ] [ m_paw_atomorb ] [ Functions ]

NAME

  destroy_atomorb

FUNCTION

  Free the dynamic memory allocated in a structure of type atomorb_type.

SIDE EFFECTS

  Atm <type(atomorb_type)>=datastructure containing atomic orbitals for a given type of atom.

PARENTS

      m_paw_slater

CHILDREN

      wrtout

SOURCE

229 subroutine destroy_atomorb(Atm)
230 
231 
232 !This section has been created automatically by the script Abilint (TD).
233 !Do not modify the following lines by hand.
234 #undef ABI_FUNC
235 #define ABI_FUNC 'destroy_atomorb'
236 !End of the abilint section
237 
238  implicit none
239 
240 !Arguments ------------------------------------
241 !scalars
242  type(atomorb_type),intent(inout) :: Atm
243 
244 !************************************************************************
245 
246  !@atomorb_type
247 
248  ! integers
249  if (allocated(Atm%indlmn)) then
250    ABI_FREE(Atm%indlmn)
251  end if
252  if (allocated(Atm%indln)) then
253    ABI_FREE(Atm%indln)
254  end if
255  if (allocated(Atm%indklmn)) then
256    ABI_FREE(Atm%indklmn)
257  end if
258  if (allocated(Atm%klm_diag)) then
259    ABI_FREE(Atm%klm_diag)
260  end if
261  if (allocated(Atm%klmntomn)) then
262    ABI_FREE(Atm%klmntomn)
263  end if
264  if (allocated(Atm%kln2ln)) then
265    ABI_FREE(Atm%kln2ln)
266  end if
267  if (allocated(Atm%mode)) then
268    ABI_FREE(Atm%mode)
269  end if
270 
271  !real
272  if (allocated(Atm%eig)) then
273    ABI_FREE(Atm%eig)
274  end if
275  if (allocated(Atm%occ)) then
276    ABI_FREE(Atm%occ)
277  end if
278  if (allocated(Atm%phi)) then
279    ABI_FREE(Atm%phi)
280  end if
281 
282 end subroutine destroy_atomorb

m_paw_atomorb/get_atomorb_charge [ Functions ]

[ Top ] [ m_paw_atomorb ] [ Functions ]

NAME

  get_atomorb_charge

FUNCTION

  Get core charge from a structure of type atomorb_type
  and optionally core density.

INPUTS

  Atm<atomorb_type>=Structure defining the set of core orbitals.
  Radmesh<pawrad_type>=Info oh the Radial mesh used for core electrons.

OUTPUT

  nele=core charge
  raddens(mesh_size)=core density (optional)

PARENTS

      m_paw_atomorb

CHILDREN

      wrtout

SOURCE

682 subroutine get_atomorb_charge(Atm,Radmesh,nele,radens)
683 
684 
685 !This section has been created automatically by the script Abilint (TD).
686 !Do not modify the following lines by hand.
687 #undef ABI_FUNC
688 #define ABI_FUNC 'get_atomorb_charge'
689 !End of the abilint section
690 
691  implicit none
692 
693 !Arguments ------------------------------------
694 !scalars
695  real(dp),intent(out) :: nele
696  type(atomorb_type),intent(in) :: Atm
697  type(pawrad_type),intent(in) :: Radmesh
698 !arrays
699  real(dp),optional,intent(out) :: radens(Atm%mesh_size,Atm%nspden)
700 
701 !Local variables-------------------------------
702 !scalars
703  integer :: iln,isppol
704  real(dp) :: intg,focc
705  real(dp),allocatable :: phi2nl(:)
706 
707 !************************************************************************
708 
709  if (Atm%nsppol==2) then
710    MSG_ERROR("nsppol==2 is Working in progress")
711  end if
712 
713  ABI_MALLOC(phi2nl,(Atm%mesh_size))
714  if (PRESENT(radens)) radens = zero
715 
716  nele = zero
717  do isppol=1,Atm%nsppol
718    do iln=1,Atm%ln_size
719      !Atm%mode(iln,isppol) TODO add option to select particular states
720      focc   = Atm%occ(iln,isppol)
721      if (ABS(focc) > tol16) then
722        phi2nl = Atm%phi(:,iln,isppol)**2
723        call simp_gen(intg,phi2nl,Radmesh)
724        nele = nele + focc*intg
725 !       if (PRESENT(radens)) then  !FIXME maybe it is better to rr**2 radens
726 !        radens(2:Atm%mesh_size) = radens(2:Atm%mesh_size) &
727 !&         + focc * phi2nl(2:Atm%mesh_size)/(four_pi*Radmesh%rad(2:Atm%mesh_size)**2)
728 !       end if
729      end if
730    end do
731  end do
732 
733  ABI_FREE(phi2nl)
734 
735 end subroutine get_atomorb_charge

m_paw_atomorb/get_overlap [ Functions ]

[ Top ] [ m_paw_atomorb ] [ Functions ]

NAME

  get_overlap

FUNCTION

  Get overlap between core and valence states

INPUTS

  Atm<atomorb_type>=Structure defining the set of core states
  Atmesh<pawrad_type>=Info oh the Radial mesh used for core states
  isppol=index for spin component
  nphi=number of core states
  phi(Radmesh2%mesh_size,nphi)=valence states
  phi_indln(nphi)=Array giving l and and n for i=1,nphi
  Radmesh2<pawrad_type>=Info oh the Radial mesh used for valence states

OUTPUT

  overlap(ln_size,nphi)=core-valence overlap matrix

PARENTS

      m_paw_slater

CHILDREN

      wrtout

SOURCE

767 subroutine get_overlap(Atm,Atmesh,Radmesh2,isppol,nphi,phi,phi_indln,overlap)
768 
769 
770 !This section has been created automatically by the script Abilint (TD).
771 !Do not modify the following lines by hand.
772 #undef ABI_FUNC
773 #define ABI_FUNC 'get_overlap'
774 !End of the abilint section
775 
776  implicit none
777 
778 !Arguments ------------------------------------
779 !scalars
780  integer,intent(in) :: nphi,isppol
781  type(atomorb_type),intent(in) :: Atm
782  type(pawrad_type),target,intent(in) :: Atmesh,Radmesh2
783 !arrays
784  integer,intent(in) :: phi_indln(2,nphi)
785  real(dp),target,intent(in) :: phi(Radmesh2%mesh_size,nphi)
786  real(dp),intent(out) :: overlap(Atm%ln_size,nphi)
787 
788 !Local variables-------------------------------
789 !scalars
790  integer :: iln_atm,iphi,ll_phi,ll_atm,do_spline,iln
791  integer :: whichdenser,size4spl,my_mesh_size
792  real(dp) :: ybcbeg,ybcend,intg
793  logical :: hasameq
794 !arrays
795  real(dp),pointer :: ff_spl(:,:)
796  real(dp),allocatable :: der(:),ypp(:),func(:)
797  real(dp),pointer :: rad4spl(:),my_pts(:)
798 
799 !************************************************************************
800 
801  ABI_CHECK(isppol>0.and.isppol<=Atm%nsppol,"Wrong isppol")
802 
803  call pawrad_isame(Atmesh,Radmesh2,hasameq,whichdenser)
804 
805  do_spline= 0; if (.not.hasameq) do_spline=1
806 
807  my_mesh_size = MIN(Atmesh%mesh_size,Radmesh2%mesh_size)
808  ff_spl => phi
809 
810  ! === Spline valence onto Atom mesh (natural spline) ===
811  if (do_spline==1) then
812    MSG_COMMENT("Splining in overlap")
813 
814    my_mesh_size  =  Atmesh%mesh_size
815    my_pts        => Atmesh%rad(1:my_mesh_size)
816    ABI_MALLOC(ff_spl,(my_mesh_size,nphi))
817 
818    size4spl =  Radmesh2%mesh_size
819    rad4spl  => Radmesh2%rad
820    ABI_MALLOC(der,(size4spl))
821    ABI_MALLOC(ypp,(size4spl))
822 
823    do iln=1,nphi
824      ypp(:) = zero; ybcbeg = zero; ybcend = zero
825      call spline(rad4spl,phi(:,iln),size4spl,ybcbeg,ybcend,ypp)
826 
827      call splint(size4spl,rad4spl,phi(:,iln),ypp,my_mesh_size,my_pts,ff_spl(:,iln))
828    end do
829 
830    ABI_FREE(der)
831    ABI_FREE(ypp)
832  end if
833 
834  ABI_MALLOC(func,(my_mesh_size))
835  overlap = zero
836 
837  do iphi=1,nphi
838    ll_phi = phi_indln(1,iphi)
839    do iln_atm=1,Atm%ln_size
840      ll_atm = Atm%indln(1,iln_atm)
841 
842      if (ll_atm == ll_phi) then ! selection rule on l
843        func(:) = Atm%phi(1:my_mesh_size,iln_atm,isppol) * ff_spl(1:my_mesh_size,iphi)
844        call simp_gen(intg,func,Atmesh)
845        overlap(iln_atm,iphi)=intg
846        write(std_out,*)"overlap <phic_i|phi_j> for ll_phi",ll_phi,"ll_phic",ll_atm,"=",intg
847      end if
848 
849    end do
850  end do
851  ABI_FREE(func)
852 
853  if (do_spline==1)  then
854    ABI_FREE(ff_spl)
855  end if
856 
857 end subroutine get_overlap

m_paw_atomorb/init_atomorb [ Functions ]

[ Top ] [ m_paw_atomorb ] [ Functions ]

NAME

  init_atomorb

FUNCTION

  Initialize a structure of type atomorb_type from file.

INPUTS

  filename=Name of the file containing core electrons

OUTPUT

  Atmrad<pawrad_type>=Info oh the Radial mesh used for core electrons.
  Atm<atomorb_type>=Structure defining the set of core orbitals.
  ierr=Status error.
   * 1 if error during the opening of the file.
   * 2 for generic error during the reading.

PARENTS

      m_paw_slater

CHILDREN

      wrtout

SOURCE

312 subroutine init_atomorb(Atm,Atmrad,rcut,filename,prtvol,ierr)
313 
314 
315 !This section has been created automatically by the script Abilint (TD).
316 !Do not modify the following lines by hand.
317 #undef ABI_FUNC
318 #define ABI_FUNC 'init_atomorb'
319 !End of the abilint section
320 
321  implicit none
322 
323 !Arguments ------------------------------------
324 !scalars
325  integer,intent(in) :: prtvol
326  integer,intent(out) :: ierr
327  real(dp),intent(in) :: rcut
328  character(len=*),intent(in) :: filename
329  type(atomorb_type),intent(out) :: Atm
330  type(pawrad_type),intent(out) :: Atmrad
331 
332 !Local variables-------------------------------
333 !scalars
334  integer :: creatorid
335  integer :: imainmesh,jlmn,jl,k0lmn,ilmn,il
336  integer :: lcutdens,version,klmn,nn,ll
337  integer :: lmax,pspdat
338  integer :: ii,unt,iln,ir,nmesh,imsh
339  integer :: iread1,pspcod,msz,isppol
340  integer :: mesh_type,mesh_size
341  real(dp) :: charge
342  real(dp) :: my_rcore
343  real(dp) :: rstep,lstep
344  real(dp):: occ,ene
345  character(len=80) :: line,title
346  character(len=500) :: msg
347 !arrays
348  integer,allocatable :: orbitals(:)
349  real(dp),allocatable :: phitmp(:,:,:),radens(:,:)
350  type(pawrad_type),allocatable :: Radmesh(:)
351 
352 ! ************************************************************************
353 
354  ! File format of formatted core orbitals.
355 
356  !(1) title (character) line
357  !(2) method, nspinor,nsppol
358  !(3) znucl, zion, pspdat
359  !(4) ixc, lmax
360  !(5) version, creatorID
361  !(6) ln_size, lmn_size
362  !(7) orbitals (for i=1 to ln_size)
363  !(8) number_of_meshes
364  !
365  !For imsh=1 to number_of_meshes
366  !(9)  mesh_index, mesh_type ,mesh_size, rad_step[, log_step]
367  !(10) rcore (SPH)
368  !
369  !For iln=1 to basis_size
370  !  (11) comment(character)
371  !  (12) radial mesh index for phi
372  !  (13) nn, ll,
373  !  (14) phi(r) (for ir=1 to phi_meshsz)
374  !
375  !Comments:
376  ! * allowed values for method are:
377  !    1 for restricted, compatible only with nsppol=1.
378  !    2 for spin unrestricted, compatible only with nsppol=2.
379  !* psp_version= ID of PAW_psp version
380  !    4 characters string of the form 'pawn' (with n varying)
381  !* creatorID= ID of psp generator
382  !  creatorid=1xyz : psp generated from Holzwarth AtomPAW generator version x.yz
383  !  creatorid=2xyz : psp generated from Vanderbilt ultra-soft generator version x.yz
384  !  creatorid=-1: psp for tests (for developpers only)
385  !* mesh_type= type of radial mesh
386  !    mesh_type=1 (regular grid): rad(i)=(i-1)*AA
387  !    mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
388  !    mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
389  !    mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
390  ! --------------------------------------------------------------------------------
391 
392  !@atomorb_type
393  ierr=0
394  if (open_file(filename,msg,newunit=unt,form="formatted",status="old",action="read") /=0) then
395    MSG_WARNING(msg)
396    ierr=1; RETURN
397  end if
398 
399  Atm%fname = filename
400  imainmesh = 1
401 
402  !1) title
403  read(unt,'(a80)',ERR=10)title
404  write(msg,'(a,1x,a)')ch10,TRIM(title)
405  call wrtout(std_out,msg,'COLL')
406 
407  read(unt,*,ERR=10)Atm%method, Atm%nspinor, Atm%nsppol
408  write(msg,'(3(i2,2x),22x,a)' )Atm%method, Atm%nspinor, Atm%nsppol,' method, nspinor, nsppol.'
409  call wrtout(std_out,msg,'COLL')
410 
411  Atm%nspden   = 1 !FIXME pass it as input
412 
413  !2)
414  read(unt,*,ERR=10) Atm%znucl, Atm%zion, pspdat
415  write(msg,'(2f10.5,2x,i8,2x,a)' )Atm%znucl, Atm%zion, pspdat,' znucl, zion, pspdat'
416  call wrtout(std_out,msg,'COLL')
417 
418  !3)
419  read(unt,*,ERR=10)pspcod,Atm%ixc,lmax
420  write(msg,'(2i5,2x,2x,a)')Atm%ixc,lmax,'ixc,lmax'
421  call wrtout(std_out,msg,'COLL')
422 
423  Atm%l_max  =  lmax+1
424  Atm%l_size =2*lmax+1
425 
426  !4)
427 !Read psp version in line 4 of the header
428  version=1
429  read(unt,'(a80)',ERR=10) line
430  line=ADJUSTL(line)
431 
432  if (tolower(line(1:4))=="core") read(unit=line(5:80),fmt=*) version
433  if (version/=1) then
434    write(msg,'(a,i2,3a)')&
435 &   'This version of core psp file (',version,') is not compatible with',ch10,&
436 &   'the current version of Abinit.'
437    MSG_ERROR(msg)
438  end if
439 
440  read(unit=line(6:80),fmt=*,ERR=10) creatorid
441 
442  !==========================================================
443 
444  ! 5)
445  read(unt,*,ERR=10)Atm%ln_size, Atm%lmn_size
446 
447  Atm%ln2_size  = Atm%ln_size *(Atm%ln_size +1)/2
448  Atm%lmn2_size = Atm%lmn_size*(Atm%lmn_size+1)/2
449 
450  ! 6)
451  ABI_MALLOC(orbitals,(Atm%ln_size))
452  read(unt,*,ERR=10) (orbitals(iln), iln=1,Atm%ln_size)
453 
454  lmax = MAXVAL(orbitals)
455  if (lmax+1/=Atm%l_max) then
456    write(msg,'(a)')" lmax read from file does not agree with orbitals. "
457    MSG_ERROR(msg)
458  end if
459 
460  ! 7)
461  read(unt,*)nmesh
462  ABI_DT_MALLOC(Radmesh,(nmesh))
463  do imsh=1,nmesh
464    lstep=zero
465    read(unt,'(a80)') line
466    read(unit=line,fmt=*,err=20,end=20) ii,mesh_type,mesh_size,rstep,lstep
467    20 continue
468    ABI_CHECK(ii<=nmesh,"Index of mesh out of range!")
469    call pawrad_init(Radmesh(ii),mesh_size,mesh_type,rstep,lstep,-one)
470  end do
471 
472  ! 8)
473  read(unt,*,ERR=10) my_rcore
474 
475  Atm%rcore = my_rcore
476  if (rcut > tol6) then
477    Atm%rcore = rcut
478    write(msg,'(a,f8.5,a,f8.5,a)')&
479 &    " Truncating radial mesh for core orbitals using new rcore = ",Atm%rcore," (",my_rcore,")"
480    call wrtout(std_out,msg,'COLL')
481    if (rcut > my_rcore) then
482      MSG_ERROR("rcut should not exceed my_rcore")
483    end if
484  end if
485 
486  !==========================================================
487  ! Mirror pseudopotential parameters to the output and log files
488 
489  write(msg,'(a,i1)')' Pseudopotential format is: core ',version
490  call wrtout(std_out,msg,'COLL')
491  write(msg,'(2(a,i3),a,64i4)')' (ln_size)= ',Atm%ln_size,' (lmn_size= ',Atm%lmn_size,'), orbitals= ',orbitals(1:Atm%ln_size)
492  call wrtout(std_out,msg,'COLL')
493  write(msg,'(a,f11.8)')' Radius used for core orbitals: rc_sph= ',Atm%rcore
494  call wrtout(std_out,msg,'COLL')
495  write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:'
496  call wrtout(std_out,msg,'COLL')
497 
498  do imsh=1,nmesh
499    call pawrad_print(Radmesh(imsh),prtvol=prtvol)
500  end do
501 
502  !==========================================================
503 
504  !---------------------------------
505  ! (11-12-13) Read wave-functions (phi)
506  ! * Initialize also the mesh.
507  ABI_MALLOC(Atm%indln,(2,Atm%ln_size))
508  ABI_MALLOC(Atm%eig,(Atm%ln_size,Atm%nsppol))
509  ABI_MALLOC(Atm%occ,(Atm%ln_size,Atm%nsppol))
510  Atm%occ = zero
511 
512  do isppol=1,Atm%nsppol
513    do iln=1,Atm%ln_size
514      read(unt,*,ERR=10) line
515      read(unt,*,ERR=10) iread1
516      if (iln==1.and.isppol==1) then
517        imainmesh=iread1
518        Atm%mesh_size = Radmesh(iread1)%mesh_size
519        ABI_MALLOC(Atm%phi,(Atm%mesh_size,Atm%ln_size,Atm%nsppol))
520      else if (iread1/=imainmesh) then
521        write(msg,'(3a)')&
522 &       ' All Phi core must be given on the same radial mesh !',ch10,&
523 &       ' Action: check your pseudopotential file.'
524        MSG_ERROR(msg)
525      end if
526      read(unt,*,ERR=10)nn,ll,ii
527      ABI_CHECK(ii==isppol,"Wrong spin index")
528      read(unt,*,ERR=10)ene,occ
529      Atm%indln(1,iln) = ll
530      Atm%indln(2,iln) = nn
531      Atm%occ(iln,isppol) = occ
532      Atm%eig(iln,isppol) = ene
533      read(unt,*,ERR=10) (Atm%phi(ir,iln,isppol),ir=1,Radmesh(imainmesh)%mesh_size)
534      !do ir=1,Radmesh(imainmesh)%mesh_size
535      !  write(100+iln,*)Radmesh(imainmesh)%rad(ir),Atm%phi(ir,iln,isppol)
536      !end do
537    end do
538  end do !isppol
539 
540  close(unt)
541 ! -------------------- END OF READING -------------------------------
542 
543  if ( rcut < tol16) then
544    ! * Use full mesh reported on file.
545    call pawrad_init(Atmrad,mesh_size=Radmesh(imainmesh)%mesh_size,mesh_type=Radmesh(imainmesh)%mesh_type,&
546 &                   rstep=Radmesh(imainmesh)%rstep,lstep=Radmesh(imainmesh)%lstep, &
547 !&                  r_for_intg=Atm%rcore) ! check this
548 &                   r_for_intg=-one)
549  else
550    ! * Truncate orbitals and radial mesh within rcut.
551    msz = MIN(pawrad_ifromr(Radmesh(imainmesh),rcut)+6, Radmesh(imainmesh)%mesh_size) ! add six more points
552 
553    write(msg,'(a,f8.5,a,f8.5,a)')&
554 &   " Truncating radial mesh for core orbitals ",Radmesh(imainmesh)%rad(msz),"(",my_rcore,")"
555    call wrtout(std_out,msg,'COLL')
556    Atm%rcore = Radmesh(imainmesh)%rad(msz)
557 
558    ABI_MALLOC(phitmp,(msz,Atm%ln_size,Atm%nsppol))
559    phitmp = Atm%phi(1:msz,:,:)
560 
561    ABI_FREE(Atm%phi)
562    ABI_MALLOC(Atm%phi,(msz,Atm%ln_size,Atm%nsppol))
563    Atm%phi = phitmp
564    ABI_FREE(phitmp)
565 
566    ! * Compute new mesh for core orbitals.
567    mesh_type = Radmesh(imainmesh)%mesh_type
568    mesh_size = msz !radmesh(imainmesh)%mesh_size
569    rstep     = Radmesh(imainmesh)%rstep
570    lstep     = Radmesh(imainmesh)%lstep
571 
572    call pawrad_init(Atmrad,mesh_size=mesh_size,mesh_type=mesh_type,rstep=rstep,lstep=lstep,&
573 !&                  r_for_intg=Atm%rcore)
574 &                   r_for_intg=-one)
575 
576    !do isppol=1,Atm%nsppol
577    !  do iln=1,Atm%ln_size
578    !    do ir=1,Atmrad%mesh_size
579    !      write(200+iln,*)Atmrad%rad(ir),Atm%phi(ir,iln,isppol)
580    !    end do
581    !  end do
582    !end do
583 
584  end if
585 
586  Atm%mesh_size = Atmrad%mesh_size
587  call pawrad_print(Atmrad,header="Final mesh",prtvol=prtvol)
588 
589  ABI_MALLOC(radens,(Atm%mesh_size,Atm%nspden))
590  call get_atomorb_charge(Atm,Atmrad,charge,radens=radens)
591 
592  write(std_out,*)"core charge  = ",charge
593  !do ii=1,Atmrad%mesh_size
594  ! write(77,*)Atmrad%rad(ii),(radens(ii,isppol),isppol=1,Atm%nspden)
595  !end do
596  ABI_FREE(radens)
597 
598  !==========================================================
599  ! Free temporary allocated space
600 
601  call pawrad_free(Radmesh)
602  ABI_DT_FREE(Radmesh)
603 
604  ! * Setup of indlmn
605  ABI_MALLOC(Atm%indlmn,(6,Atm%lmn_size))
606  call make_indlmn(Atm%ln_size, Atm%lmn_size, orbitals, Atm%indlmn)
607 
608  ABI_FREE(orbitals)
609 
610  ! * Setup of indklmn and klm_diag.
611  lcutdens=HUGE(1)
612  ABI_MALLOC(Atm%indklmn,(8,Atm%lmn2_size))
613  ABI_MALLOC(Atm%klm_diag,(Atm%lmn2_size))
614 
615  call make_indklmn(lcutdens, Atm%lmn_size, Atm%lmn2_size, Atm%indlmn, Atm%indklmn, Atm%klm_diag)
616 
617  ! * Setup of klmntomn.
618  ABI_MALLOC(Atm%klmntomn,(4,Atm%lmn2_size))
619 
620  do jlmn=1,Atm%lmn_size
621    jl= Atm%indlmn(1,jlmn)
622    k0lmn=jlmn*(jlmn-1)/2
623    do ilmn=1,jlmn
624      il= Atm%indlmn(1,ilmn)
625      klmn=k0lmn+ilmn
626      Atm%klmntomn(1,klmn) = Atm%indlmn(2,ilmn)+il+1 ! im
627      Atm%klmntomn(2,klmn) = Atm%indlmn(2,jlmn)+jl+1 ! jm
628      Atm%klmntomn(3,klmn) = Atm%indlmn(3,ilmn)      ! in
629      Atm%klmntomn(4,klmn) = Atm%indlmn(3,jlmn)      ! jn
630      !write(std_out,*) jlmn,ilmn,Atm%klmntomn(:,klmn)
631    end do
632  end do
633 
634  ! * Setup of kln2ln.
635  !TODO this has to be tested
636  ABI_MALLOC(Atm%kln2ln,(6,Atm%ln2_size))
637  call make_kln2ln(Atm%lmn_size,Atm%lmn2_size,Atm%ln2_size,Atm%indlmn,Atm%indklmn,Atm%kln2ln)
638 
639  ! * Treat all states as core.
640  ABI_MALLOC(Atm%mode,(Atm%ln_size,Atm%nsppol))
641  Atm%mode = ORB_FROZEN
642 
643  !call print_atomorb(Atm,prtvol=prtvol)
644 
645  return
646  !
647  ! === Propagate the error ===
648 10 continue
649  ierr=2
650  MSG_WARNING("Wrong file format")
651  return
652 
653 end subroutine init_atomorb

m_paw_atomorb/my_mode2str [ Functions ]

[ Top ] [ m_paw_atomorb ] [ Functions ]

NAME

  my_mode2str

FUNCTION

  Converts an integer flags defining the way an orbital is treated to a string.

INPUTS

  mode=Integer

OUTPUT

  str=mode. Either "Frozen", "Relazed Core", "Valence"

PARENTS

CHILDREN

SOURCE

 976 function my_mode2str(mode) result(str)
 977 
 978 
 979 !This section has been created automatically by the script Abilint (TD).
 980 !Do not modify the following lines by hand.
 981 #undef ABI_FUNC
 982 #define ABI_FUNC 'my_mode2str'
 983 !End of the abilint section
 984 
 985  implicit none
 986 
 987 !Arguments ------------------------------------
 988 !scalars
 989  integer,intent(in) :: mode
 990  character(len=50) :: str
 991 
 992 !Local variables
 993  character(len=500) :: msg
 994 
 995 !************************************************************************
 996 
 997  select case (mode)
 998  case (ORB_FROZEN)
 999    str="Frozen Orbital"
1000  case (ORB_RELAXED_CORE)
1001    str="Relaxed Core Orbital"
1002  case (ORB_VALENCE)
1003    str="Valence Orbital"
1004  case default
1005    write(msg,'(a,i3)')" Wrong mode= ",mode
1006    MSG_BUG(msg)
1007  end select
1008 
1009 end function my_mode2str

m_paw_atomorb/print_atomorb [ Functions ]

[ Top ] [ m_paw_atomorb ] [ Functions ]

NAME

  print_atomorb

FUNCTION

  Reports info on a structure of type atomorb_type.

INPUTS

  Atm <type(atomorb_type)>=datastructure containing atomic orbitals for a given type of atom.

OUTPUT

PARENTS

      m_paw_slater

CHILDREN

      wrtout

SOURCE

882 subroutine print_atomorb(Atm,header,unit,prtvol,mode_paral)
883 
884 
885 !This section has been created automatically by the script Abilint (TD).
886 !Do not modify the following lines by hand.
887 #undef ABI_FUNC
888 #define ABI_FUNC 'print_atomorb'
889 !End of the abilint section
890 
891  implicit none
892 
893 !Arguments ------------------------------------
894 !scalars
895  type(atomorb_type),intent(in) :: Atm
896  integer,optional,intent(in) :: prtvol,unit
897  character(len=*),optional,intent(in) :: header
898  character(len=4),optional,intent(in) :: mode_paral
899 
900 !Local variables-------------------------------
901  integer :: my_unt,my_prtvol,iln,ll,nn,isppol
902  character(len=4) :: my_mode
903  character(len=500) :: msg
904 ! ************************************************************************
905 
906  !@atomorb_type
907  my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
908  my_prtvol=0      ; if (PRESENT(prtvol    )) my_prtvol=prtvol
909  my_mode  ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
910 
911  msg=' ==== Info on the atomorb_type ==== '
912  if (PRESENT(header)) msg=header
913  call wrtout(my_unt,msg,my_mode)
914 
915  select case (Atm%method)
916  case (1)
917    msg = "  Spin restricted"
918  case(2)
919    msg = "  Spin unrestricted"
920  case default
921    write(msg,'(a,i3)')" Wrong method= ",Atm%method
922    MSG_BUG(msg)
923  end select
924  call wrtout(my_unt,msg,my_mode)
925 
926  write(msg,'(7(a,i5,a),(a,f8.5,a))')&
927 & '  Number of spinorial components ...... ',Atm%nspinor,ch10,&
928 & '  Number of ind. spin polarizations ... ',Atm%nsppol,ch10,&
929 & '  Number of spin-density components ... ',Atm%nspden,ch10,&
930 & '  Maximum angular momentum + 1 ........ ',Atm%l_max,ch10,&
931 & '  Number of (l,n) orbitals  ........... ',Atm%ln_size,ch10,&
932 & '  Number of (l,m,n) orbitals  ......... ',Atm%lmn_size,ch10,&
933 & '  Dimensions of radial mesh ........... ',Atm%mesh_size,ch10,&
934 & '  Core Radius  ........................ ',Atm%rcore,ch10
935  call wrtout(my_unt,msg,my_mode)
936 
937  write(msg,'(2(a,f8.5,a))')&
938 & '  Ionic charge ........................ ',Atm%zion,ch10,&
939 & '  Atomic number ....................... ',Atm%znucl,ch10
940  call wrtout(my_unt,msg,my_mode)
941 
942  do isppol=1,Atm%nsppol
943    do iln=1,Atm%ln_size
944      ll = Atm%indln(1,iln)
945      nn = Atm%indln(2,iln)
946      write(msg,'(" n=",i2,", l=",i2,", spin=",i2,", nocc=",f15.7,", energy=",f15.7,2x,"(",a,")")')&
947 &      nn,ll,isppol,Atm%occ(iln,isppol),Atm%eig(iln,isppol),TRIM(my_mode2str(Atm%mode(iln,isppol)))
948      call wrtout(my_unt,msg,my_mode)
949    end do
950  end do
951 
952 end subroutine print_atomorb