TABLE OF CONTENTS
- ABINIT/m_paw_atomorb
- m_paw_atomorb/atomorb_type
- m_paw_atomorb/destroy_atomorb
- m_paw_atomorb/get_atomorb_charge
- m_paw_atomorb/get_overlap
- m_paw_atomorb/init_atomorb
- m_paw_atomorb/my_mode2str
- m_paw_atomorb/print_atomorb
ABINIT/m_paw_atomorb [ 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-2024 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
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 MODULE m_paw_atomorb 25 26 use defs_basis 27 use m_abicore 28 use m_errors 29 use m_splines 30 31 use m_io_tools, only : open_file 32 use m_fstrings, only : tolower 33 use m_paw_lmn, only : make_indlmn, make_indklmn, make_kln2ln 34 use m_pawrad, only : pawrad_type, pawrad_init, & 35 & pawrad_free, pawrad_print, pawrad_isame, pawrad_ifromr, simp_gen 36 37 implicit none 38 39 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
56 type, public :: atomorb_type 57 58 !scalars 59 integer :: ixc 60 ! Exchange and correlation functional used to generate the orbitals 61 62 integer :: method 63 ! 1 for restricted, compatible only with nsppol=1. 64 ! 2 for spin unrestricted, compatible only with nsppol=2. 65 66 integer :: nspden 67 ! Number of spin-density components. 68 69 integer :: nsppol 70 ! Number of independent spin-components. 71 ! FIXME: here a lot of quantities might depend on nsppol in the 72 ! case of magnetic atoms! 73 74 integer :: nspinor 75 ! Number of spinorial components 76 ! TODO this is a quite delicate issue, for S.O. one should use J = L+S instead of L! 77 ! If we use scalar relativistic then... 78 79 integer :: l_max 80 ! Maximum value of angular momentum l+1 81 82 integer :: l_size 83 ! Maximum value of l+1 leading to non zero Gaunt coeffs 84 ! l_size=2*l_max-1 85 86 integer :: ln_size 87 ! Number of (l,n) components. 88 89 integer :: ln2_size 90 ! ln2_size=ln_size*(ln_size+1)/2 91 ! where ln_size is the number of (l,n) elements for core orbitals. 92 93 integer :: lmn_size 94 ! Number of (l,m,n) elements. 95 96 integer :: lmn2_size 97 ! lmn2_size=lmn_size*(lmn_size+1)/2 98 ! where lmn_size is the number of (l,m,n) elements for core orbitals. 99 100 integer :: mesh_size 101 ! Dimension of the radial mesh. 102 103 real(dp) :: rcore 104 ! Radius of the sphere used to describe core electrons. 105 ! It should be <= rpaw 106 107 real(dp) :: zion 108 ! zionpsp 109 ! The ionic pseudo-charge, (giving raise to a long-range coulomb potential) 110 111 ! TODO alchemy? 112 !real(dp) :: ziontypat 113 ! ziontypat 114 ! For each type of atom (might be alchemy wrt psps), the ionic pseudo-charge 115 ! (giving raise to a long-range coulomb potential) 116 117 real(dp) :: znucl 118 ! The atomic number of the atom. 119 120 ! TODO alchemy? 121 !real(dp) :: znucltypat 122 ! znucltypat 123 ! The atomic number of each type of atom (might be alchemy wrt psps) 124 125 character(len=fnlen) :: fname 126 ! The filename for temporary storage. 127 128 !arrays 129 integer, allocatable :: indlmn(:,:) 130 ! indlmn(6,lmn_size) 131 ! Array giving l,m,n,lm,ln,spin for i=lmn. 132 133 integer, allocatable :: indln(:,:) 134 ! indln(2,ln_size) 135 ! Array giving l and n for i=ln 136 137 integer, allocatable :: indklmn(:,:) 138 ! indklmn(8,lmn2_size) 139 ! Array giving klm, kln, abs(il-jl), (il+jl), ilm and jlm, ilmn and jlmn for each klmn=(ilmn,jlmn) 140 ! Note: ilmn=(il,im,in) and ilmn<=jlmn 141 142 !integer, allocatable :: klm2lm TODO add 143 ! klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm, jlm for each klm=(ilm,jlm) 144 ! where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection. 145 146 integer, allocatable :: klm_diag(:) 147 ! klm_diag(lmn2_size) 148 ! 1 il==jl and im==jm, 0 otherwise. 149 150 integer, allocatable :: klmntomn(:,:) 151 ! klmntomn(4,lmn2_size) 152 ! Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn) 153 ! Note: ilmn=(il,im,in) and ilmn<=jlmn 154 ! NB: klmntomn is an application and not a bijection 155 156 integer, allocatable :: kln2ln(:,:) 157 ! kln2ln(6,ln2_size) 158 ! Table giving il, jl ,in, jn, iln, jln for each kln=(iln,jln) 159 ! where iln=(il,in) and iln<=jln. NB: kln2ln is an application and not a bijection 160 161 integer, allocatable :: mode(:,:) 162 ! mode(ln_size,nsppol) 163 ! Flag defining how the orbital is treated. 164 ! During the pseudopotential generation we can have: ORB_FROZEN or ORB_VALENCE 165 ! For calculations in extended systems we can have: ORB_FROZEN or ORB_RELAXED_CORE 166 ! Namely different treatment depending of the degree of localization. 167 ! For example the 1s in plutonium might be treated as ORB_FROZEN during 168 ! a relaxed core calculation. 169 ! TODO define function to test the type, much safer! 170 171 real(dp), allocatable :: eig(:,:) 172 ! eig(ln_size,nsppol) 173 ! Eigenvalues for each ln channel and spin. 174 175 real(dp), allocatable :: occ(:,:) 176 ! occ(ln_size,nsppol) 177 ! Occupation for each ln channel and spin. 178 179 real(dp), allocatable :: phi(:,:,:) 180 ! phi(mesh_size,ln_size,nsppol) 181 ! Here we might have different meshes, useful for deep states in heavy atoms! 182 183 ! this might be retrieved with a method get_atomden 184 !real(dp), allocatable :: density(:) 185 ! density(mesh_size,nspden) 186 ! Gives the core density of the atom for each spin channel 187 ! Total charge in first dimension,up component in second one (if present) 188 189 end type atomorb_type 190 191 ! public procedures. 192 public :: init_atomorb 193 public :: destroy_atomorb 194 public :: print_atomorb 195 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.
SOURCE
222 subroutine destroy_atomorb(Atm) 223 224 implicit none 225 226 !Arguments ------------------------------------ 227 !scalars 228 type(atomorb_type),intent(inout) :: Atm 229 230 !************************************************************************ 231 232 !@atomorb_type 233 234 ! integers 235 if (allocated(Atm%indlmn)) then 236 ABI_FREE(Atm%indlmn) 237 end if 238 if (allocated(Atm%indln)) then 239 ABI_FREE(Atm%indln) 240 end if 241 if (allocated(Atm%indklmn)) then 242 ABI_FREE(Atm%indklmn) 243 end if 244 if (allocated(Atm%klm_diag)) then 245 ABI_FREE(Atm%klm_diag) 246 end if 247 if (allocated(Atm%klmntomn)) then 248 ABI_FREE(Atm%klmntomn) 249 end if 250 if (allocated(Atm%kln2ln)) then 251 ABI_FREE(Atm%kln2ln) 252 end if 253 if (allocated(Atm%mode)) then 254 ABI_FREE(Atm%mode) 255 end if 256 257 !real 258 if (allocated(Atm%eig)) then 259 ABI_FREE(Atm%eig) 260 end if 261 if (allocated(Atm%occ)) then 262 ABI_FREE(Atm%occ) 263 end if 264 if (allocated(Atm%phi)) then 265 ABI_FREE(Atm%phi) 266 end if 267 268 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)
SOURCE
649 subroutine get_atomorb_charge(Atm,Radmesh,nele,radens) 650 651 implicit none 652 653 !Arguments ------------------------------------ 654 !scalars 655 real(dp),intent(out) :: nele 656 type(atomorb_type),intent(in) :: Atm 657 type(pawrad_type),intent(in) :: Radmesh 658 !arrays 659 real(dp),optional,intent(out) :: radens(Atm%mesh_size,Atm%nspden) 660 661 !Local variables------------------------------- 662 !scalars 663 integer :: iln,isppol 664 real(dp) :: intg,focc 665 real(dp),allocatable :: phi2nl(:) 666 667 !************************************************************************ 668 669 if (Atm%nsppol==2) then 670 ABI_ERROR("nsppol==2 is Working in progress") 671 end if 672 673 ABI_MALLOC(phi2nl,(Atm%mesh_size)) 674 if (PRESENT(radens)) radens = zero 675 676 nele = zero 677 do isppol=1,Atm%nsppol 678 do iln=1,Atm%ln_size 679 !Atm%mode(iln,isppol) TODO add option to select particular states 680 focc = Atm%occ(iln,isppol) 681 if (ABS(focc) > tol16) then 682 phi2nl = Atm%phi(:,iln,isppol)**2 683 call simp_gen(intg,phi2nl,Radmesh) 684 nele = nele + focc*intg 685 ! if (PRESENT(radens)) then !FIXME maybe it is better to rr**2 radens 686 ! radens(2:Atm%mesh_size) = radens(2:Atm%mesh_size) & 687 !& + focc * phi2nl(2:Atm%mesh_size)/(four_pi*Radmesh%rad(2:Atm%mesh_size)**2) 688 ! end if 689 end if 690 end do 691 end do 692 693 ABI_FREE(phi2nl) 694 695 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
SOURCE
721 subroutine get_overlap(Atm,Atmesh,Radmesh2,isppol,nphi,phi,phi_indln,overlap) 722 723 implicit none 724 725 !Arguments ------------------------------------ 726 !scalars 727 integer,intent(in) :: nphi,isppol 728 type(atomorb_type),intent(in) :: Atm 729 type(pawrad_type),target,intent(in) :: Atmesh,Radmesh2 730 !arrays 731 integer,intent(in) :: phi_indln(2,nphi) 732 real(dp),target,intent(in) :: phi(Radmesh2%mesh_size,nphi) 733 real(dp),intent(out) :: overlap(Atm%ln_size,nphi) 734 735 !Local variables------------------------------- 736 !scalars 737 integer :: iln_atm,iphi,ll_phi,ll_atm,do_spline,iln 738 integer :: whichdenser,size4spl,my_mesh_size 739 real(dp) :: ybcbeg,ybcend,intg 740 logical :: hasameq 741 !arrays 742 real(dp),pointer :: ff_spl(:,:) 743 real(dp),allocatable :: der(:),ypp(:),func(:) 744 real(dp),pointer :: rad4spl(:),my_pts(:) 745 746 !************************************************************************ 747 748 ABI_CHECK(isppol>0.and.isppol<=Atm%nsppol,"Wrong isppol") 749 750 call pawrad_isame(Atmesh,Radmesh2,hasameq,whichdenser) 751 752 do_spline= 0; if (.not.hasameq) do_spline=1 753 754 my_mesh_size = MIN(Atmesh%mesh_size,Radmesh2%mesh_size) 755 ff_spl => phi 756 757 ! === Spline valence onto Atom mesh (natural spline) === 758 if (do_spline==1) then 759 ABI_COMMENT("Splining in overlap") 760 761 my_mesh_size = Atmesh%mesh_size 762 my_pts => Atmesh%rad(1:my_mesh_size) 763 ABI_MALLOC(ff_spl,(my_mesh_size,nphi)) 764 765 size4spl = Radmesh2%mesh_size 766 rad4spl => Radmesh2%rad 767 ABI_MALLOC(der,(size4spl)) 768 ABI_MALLOC(ypp,(size4spl)) 769 770 do iln=1,nphi 771 ypp(:) = zero; ybcbeg = zero; ybcend = zero 772 call spline(rad4spl,phi(:,iln),size4spl,ybcbeg,ybcend,ypp) 773 774 call splint(size4spl,rad4spl,phi(:,iln),ypp,my_mesh_size,my_pts,ff_spl(:,iln)) 775 end do 776 777 ABI_FREE(der) 778 ABI_FREE(ypp) 779 end if 780 781 ABI_MALLOC(func,(my_mesh_size)) 782 overlap = zero 783 784 do iphi=1,nphi 785 ll_phi = phi_indln(1,iphi) 786 do iln_atm=1,Atm%ln_size 787 ll_atm = Atm%indln(1,iln_atm) 788 789 if (ll_atm == ll_phi) then ! selection rule on l 790 func(:) = Atm%phi(1:my_mesh_size,iln_atm,isppol) * ff_spl(1:my_mesh_size,iphi) 791 call simp_gen(intg,func,Atmesh) 792 overlap(iln_atm,iphi)=intg 793 write(std_out,*)"overlap <phic_i|phi_j> for ll_phi",ll_phi,"ll_phic",ll_atm,"=",intg 794 end if 795 796 end do 797 end do 798 ABI_FREE(func) 799 800 if (do_spline==1) then 801 ABI_FREE(ff_spl) 802 end if 803 804 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.
SOURCE
292 subroutine init_atomorb(Atm,Atmrad,rcut,filename,prtvol,ierr) 293 294 implicit none 295 296 !Arguments ------------------------------------ 297 !scalars 298 integer,intent(in) :: prtvol 299 integer,intent(out) :: ierr 300 real(dp),intent(in) :: rcut 301 character(len=*),intent(in) :: filename 302 type(atomorb_type),intent(out) :: Atm 303 type(pawrad_type),intent(out) :: Atmrad 304 305 !Local variables------------------------------- 306 !scalars 307 integer :: creatorid 308 integer :: imainmesh,jlmn,jl,k0lmn,ilmn,il 309 integer :: lcutdens,version,klmn,nn,ll 310 integer :: lmax,pspdat 311 integer :: ii,unt,iln,ir,nmesh,imsh 312 integer :: iread1,pspcod,msz,isppol 313 integer :: mesh_type,mesh_size 314 real(dp) :: corecharge 315 real(dp) :: my_rcore 316 real(dp) :: rstep,lstep 317 real(dp):: occ,ene 318 character(len=80) :: line,title 319 character(len=500) :: msg 320 !arrays 321 integer,allocatable :: orbitals(:) 322 real(dp),allocatable :: phitmp(:,:,:),radens(:,:) 323 type(pawrad_type),allocatable :: Radmesh(:) 324 325 ! ************************************************************************ 326 327 ! File format of formatted core orbitals. 328 329 !(1) title (character) line 330 !(2) method, nspinor,nsppol 331 !(3) znucl, zion, pspdat 332 !(4) ixc, lmax 333 !(5) version, creatorID 334 !(6) ln_size, lmn_size 335 !(7) orbitals (for i=1 to ln_size) 336 !(8) number_of_meshes 337 ! 338 !For imsh=1 to number_of_meshes 339 !(9) mesh_index, mesh_type ,mesh_size, rad_step[, log_step] 340 !(10) rcore (SPH) 341 ! 342 !For iln=1 to basis_size 343 ! (11) comment(character) 344 ! (12) radial mesh index for phi 345 ! (13) nn, ll, 346 ! (14) phi(r) (for ir=1 to phi_meshsz) 347 ! 348 !Comments: 349 ! * allowed values for method are: 350 ! 1 for restricted, compatible only with nsppol=1. 351 ! 2 for spin unrestricted, compatible only with nsppol=2. 352 !* psp_version= ID of PAW_psp version 353 ! 4 characters string of the form 'pawn' (with n varying) 354 !* creatorID= ID of psp generator 355 ! creatorid=1xyz : psp generated from Holzwarth AtomPAW generator version x.yz 356 ! creatorid=2xyz : psp generated from Vanderbilt ultra-soft generator version x.yz 357 ! creatorid=-1: psp for tests (for developpers only) 358 !* mesh_type= type of radial mesh 359 ! mesh_type=1 (regular grid): rad(i)=(i-1)*AA 360 ! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) 361 ! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 362 ! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n 363 ! -------------------------------------------------------------------------------- 364 365 !@atomorb_type 366 ierr=0 367 if (open_file(filename,msg,newunit=unt,form="formatted",status="old",action="read") /=0) then 368 ABI_WARNING(msg) 369 ierr=1; RETURN 370 end if 371 372 Atm%fname = filename 373 imainmesh = 1 374 375 !1) title 376 read(unt,'(a80)',ERR=10)title 377 write(msg,'(a,1x,a)')ch10,TRIM(title) 378 call wrtout(std_out,msg,'COLL') 379 380 read(unt,*,ERR=10)Atm%method, Atm%nspinor, Atm%nsppol 381 write(msg,'(3(i2,2x),22x,a)' )Atm%method, Atm%nspinor, Atm%nsppol,' method, nspinor, nsppol.' 382 call wrtout(std_out,msg,'COLL') 383 384 Atm%nspden = 1 !FIXME pass it as input 385 386 !2) 387 read(unt,*,ERR=10) Atm%znucl, Atm%zion, pspdat 388 write(msg,'(2f10.5,2x,i8,2x,a)' )Atm%znucl, Atm%zion, pspdat,' znucl, zion, pspdat' 389 call wrtout(std_out,msg,'COLL') 390 391 !3) 392 read(unt,*,ERR=10)pspcod,Atm%ixc,lmax 393 write(msg,'(2i5,2x,2x,a)')Atm%ixc,lmax,'ixc,lmax' 394 call wrtout(std_out,msg,'COLL') 395 396 Atm%l_max = lmax+1 397 Atm%l_size =2*lmax+1 398 399 !4) 400 !Read psp version in line 4 of the header 401 version=1 402 read(unt,'(a80)',ERR=10) line 403 line=ADJUSTL(line) 404 405 if (tolower(line(1:4))=="core") read(unit=line(5:80),fmt=*) version 406 if (version/=1) then 407 write(msg,'(a,i2,3a)')& 408 & 'This version of core psp file (',version,') is not compatible with',ch10,& 409 & 'the current version of Abinit.' 410 ABI_ERROR(msg) 411 end if 412 413 read(unit=line(6:80),fmt=*,ERR=10) creatorid 414 415 !========================================================== 416 417 ! 5) 418 read(unt,*,ERR=10)Atm%ln_size, Atm%lmn_size 419 420 Atm%ln2_size = Atm%ln_size *(Atm%ln_size +1)/2 421 Atm%lmn2_size = Atm%lmn_size*(Atm%lmn_size+1)/2 422 423 ! 6) 424 ABI_MALLOC(orbitals,(Atm%ln_size)) 425 read(unt,*,ERR=10) (orbitals(iln), iln=1,Atm%ln_size) 426 427 lmax = MAXVAL(orbitals) 428 if (lmax+1/=Atm%l_max) then 429 write(msg,'(a)')" lmax read from file does not agree with orbitals. " 430 ABI_ERROR(msg) 431 end if 432 433 ! 7) 434 read(unt,*)nmesh 435 ABI_MALLOC(Radmesh,(nmesh)) 436 do imsh=1,nmesh 437 lstep=zero 438 read(unt,'(a80)') line 439 read(unit=line,fmt=*,err=20,end=20) ii,mesh_type,mesh_size,rstep,lstep 440 20 continue 441 ABI_CHECK(ii<=nmesh,"Index of mesh out of range!") 442 call pawrad_init(Radmesh(ii),mesh_size,mesh_type,rstep,lstep,-one) 443 end do 444 445 ! 8) 446 read(unt,*,ERR=10) my_rcore 447 448 Atm%rcore = my_rcore 449 if (rcut > tol6) then 450 Atm%rcore = rcut 451 write(msg,'(a,f8.5,a,f8.5,a)')& 452 & " Truncating radial mesh for core orbitals using new rcore = ",Atm%rcore," (",my_rcore,")" 453 call wrtout(std_out,msg,'COLL') 454 if (rcut > my_rcore) then 455 ABI_ERROR("rcut should not exceed my_rcore") 456 end if 457 end if 458 459 !========================================================== 460 ! Mirror pseudopotential parameters to the output and log files 461 462 write(msg,'(a,i1)')' Pseudopotential format is: core ',version 463 call wrtout(std_out,msg,'COLL') 464 write(msg,'(2(a,i3),a,64i4)')' (ln_size)= ',Atm%ln_size,' (lmn_size= ',Atm%lmn_size,'), orbitals= ',orbitals(1:Atm%ln_size) 465 call wrtout(std_out,msg,'COLL') 466 write(msg,'(a,f11.8)')' Radius used for core orbitals: rc_sph= ',Atm%rcore 467 call wrtout(std_out,msg,'COLL') 468 write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:' 469 call wrtout(std_out,msg,'COLL') 470 471 do imsh=1,nmesh 472 call pawrad_print(Radmesh(imsh),prtvol=prtvol) 473 end do 474 475 !========================================================== 476 477 !--------------------------------- 478 ! (11-12-13) Read wave-functions (phi) 479 ! * Initialize also the mesh. 480 ABI_MALLOC(Atm%indln,(2,Atm%ln_size)) 481 ABI_MALLOC(Atm%eig,(Atm%ln_size,Atm%nsppol)) 482 ABI_MALLOC(Atm%occ,(Atm%ln_size,Atm%nsppol)) 483 Atm%occ = zero 484 485 do isppol=1,Atm%nsppol 486 do iln=1,Atm%ln_size 487 read(unt,*,ERR=10) line 488 read(unt,*,ERR=10) iread1 489 if (iln==1.and.isppol==1) then 490 imainmesh=iread1 491 Atm%mesh_size = Radmesh(iread1)%mesh_size 492 ABI_MALLOC(Atm%phi,(Atm%mesh_size,Atm%ln_size,Atm%nsppol)) 493 else if (iread1/=imainmesh) then 494 write(msg,'(3a)')& 495 & ' All Phi core must be given on the same radial mesh !',ch10,& 496 & ' Action: check your pseudopotential file.' 497 ABI_ERROR(msg) 498 end if 499 read(unt,*,ERR=10)nn,ll,ii 500 ABI_CHECK(ii==isppol,"Wrong spin index") 501 read(unt,*,ERR=10)ene,occ 502 Atm%indln(1,iln) = ll 503 Atm%indln(2,iln) = nn 504 Atm%occ(iln,isppol) = occ 505 Atm%eig(iln,isppol) = ene 506 read(unt,*,ERR=10) (Atm%phi(ir,iln,isppol),ir=1,Radmesh(imainmesh)%mesh_size) 507 !do ir=1,Radmesh(imainmesh)%mesh_size 508 ! write(100+iln,*)Radmesh(imainmesh)%rad(ir),Atm%phi(ir,iln,isppol) 509 !end do 510 end do 511 end do !isppol 512 513 close(unt) 514 ! -------------------- END OF READING ------------------------------- 515 516 if ( rcut < tol16) then 517 ! * Use full mesh reported on file. 518 call pawrad_init(Atmrad,mesh_size=Radmesh(imainmesh)%mesh_size,mesh_type=Radmesh(imainmesh)%mesh_type,& 519 & rstep=Radmesh(imainmesh)%rstep,lstep=Radmesh(imainmesh)%lstep, & 520 !& r_for_intg=Atm%rcore) ! check this 521 & r_for_intg=-one) 522 else 523 ! * Truncate orbitals and radial mesh within rcut. 524 msz = MIN(pawrad_ifromr(Radmesh(imainmesh),rcut)+6, Radmesh(imainmesh)%mesh_size) ! add six more points 525 526 write(msg,'(a,f8.5,a,f8.5,a)')& 527 & " Truncating radial mesh for core orbitals ",Radmesh(imainmesh)%rad(msz),"(",my_rcore,")" 528 call wrtout(std_out,msg,'COLL') 529 Atm%rcore = Radmesh(imainmesh)%rad(msz) 530 531 ABI_MALLOC(phitmp,(msz,Atm%ln_size,Atm%nsppol)) 532 phitmp = Atm%phi(1:msz,:,:) 533 534 ABI_FREE(Atm%phi) 535 ABI_MALLOC(Atm%phi,(msz,Atm%ln_size,Atm%nsppol)) 536 Atm%phi = phitmp 537 ABI_FREE(phitmp) 538 539 ! * Compute new mesh for core orbitals. 540 mesh_type = Radmesh(imainmesh)%mesh_type 541 mesh_size = msz !radmesh(imainmesh)%mesh_size 542 rstep = Radmesh(imainmesh)%rstep 543 lstep = Radmesh(imainmesh)%lstep 544 545 call pawrad_init(Atmrad,mesh_size=mesh_size,mesh_type=mesh_type,rstep=rstep,lstep=lstep,& 546 !& r_for_intg=Atm%rcore) 547 & r_for_intg=-one) 548 549 !do isppol=1,Atm%nsppol 550 ! do iln=1,Atm%ln_size 551 ! do ir=1,Atmrad%mesh_size 552 ! write(200+iln,*)Atmrad%rad(ir),Atm%phi(ir,iln,isppol) 553 ! end do 554 ! end do 555 !end do 556 557 end if 558 559 Atm%mesh_size = Atmrad%mesh_size 560 call pawrad_print(Atmrad,header="Final mesh",prtvol=prtvol) 561 562 ABI_MALLOC(radens,(Atm%mesh_size,Atm%nspden)) 563 call get_atomorb_charge(Atm,Atmrad,corecharge,radens=radens) 564 565 write(std_out,*)"core charge = ",corecharge 566 !do ii=1,Atmrad%mesh_size 567 ! write(77,*)Atmrad%rad(ii),(radens(ii,isppol),isppol=1,Atm%nspden) 568 !end do 569 ABI_FREE(radens) 570 571 !========================================================== 572 ! Free temporary allocated space 573 574 call pawrad_free(Radmesh) 575 ABI_FREE(Radmesh) 576 577 ! * Setup of indlmn 578 ABI_MALLOC(Atm%indlmn,(6,Atm%lmn_size)) 579 call make_indlmn(Atm%ln_size, Atm%lmn_size, orbitals, Atm%indlmn) 580 581 ABI_FREE(orbitals) 582 583 ! * Setup of indklmn and klm_diag. 584 lcutdens=HUGE(1) 585 ABI_MALLOC(Atm%indklmn,(8,Atm%lmn2_size)) 586 ABI_MALLOC(Atm%klm_diag,(Atm%lmn2_size)) 587 588 call make_indklmn(lcutdens, Atm%lmn_size, Atm%lmn2_size, Atm%indlmn, Atm%indklmn, Atm%klm_diag) 589 590 ! * Setup of klmntomn. 591 ABI_MALLOC(Atm%klmntomn,(4,Atm%lmn2_size)) 592 593 do jlmn=1,Atm%lmn_size 594 jl= Atm%indlmn(1,jlmn) 595 k0lmn=jlmn*(jlmn-1)/2 596 do ilmn=1,jlmn 597 il= Atm%indlmn(1,ilmn) 598 klmn=k0lmn+ilmn 599 Atm%klmntomn(1,klmn) = Atm%indlmn(2,ilmn)+il+1 ! im 600 Atm%klmntomn(2,klmn) = Atm%indlmn(2,jlmn)+jl+1 ! jm 601 Atm%klmntomn(3,klmn) = Atm%indlmn(3,ilmn) ! in 602 Atm%klmntomn(4,klmn) = Atm%indlmn(3,jlmn) ! jn 603 !write(std_out,*) jlmn,ilmn,Atm%klmntomn(:,klmn) 604 end do 605 end do 606 607 ! * Setup of kln2ln. 608 !TODO this has to be tested 609 ABI_MALLOC(Atm%kln2ln,(6,Atm%ln2_size)) 610 call make_kln2ln(Atm%lmn_size,Atm%lmn2_size,Atm%ln2_size,Atm%indlmn,Atm%indklmn,Atm%kln2ln) 611 612 ! * Treat all states as core. 613 ABI_MALLOC(Atm%mode,(Atm%ln_size,Atm%nsppol)) 614 Atm%mode = ORB_FROZEN 615 616 !call print_atomorb(Atm,prtvol=prtvol) 617 618 return 619 ! 620 ! === Propagate the error === 621 10 continue 622 ierr=2 623 ABI_WARNING("Wrong file format") 624 return 625 626 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"
SOURCE
906 function my_mode2str(mode) result(str) 907 908 implicit none 909 910 !Arguments ------------------------------------ 911 !scalars 912 integer,intent(in) :: mode 913 character(len=50) :: str 914 915 !Local variables 916 character(len=500) :: msg 917 918 !************************************************************************ 919 920 select case (mode) 921 case (ORB_FROZEN) 922 str="Frozen Orbital" 923 case (ORB_RELAXED_CORE) 924 str="Relaxed Core Orbital" 925 case (ORB_VALENCE) 926 str="Valence Orbital" 927 case default 928 write(msg,'(a,i3)')" Wrong mode= ",mode 929 ABI_BUG(msg) 930 end select 931 932 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
SOURCE
823 subroutine print_atomorb(Atm,header,unit,prtvol,mode_paral) 824 825 implicit none 826 827 !Arguments ------------------------------------ 828 !scalars 829 type(atomorb_type),intent(in) :: Atm 830 integer,optional,intent(in) :: prtvol,unit 831 character(len=*),optional,intent(in) :: header 832 character(len=4),optional,intent(in) :: mode_paral 833 834 !Local variables------------------------------- 835 integer :: my_unt,my_prtvol,iln,ll,nn,isppol 836 character(len=4) :: my_mode 837 character(len=500) :: msg 838 ! ************************************************************************ 839 840 !@atomorb_type 841 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 842 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 843 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 844 845 msg=' ==== Info on the atomorb_type ==== ' 846 if (PRESENT(header)) msg=header 847 call wrtout(my_unt,msg,my_mode) 848 849 select case (Atm%method) 850 case (1) 851 msg = " Spin restricted" 852 case(2) 853 msg = " Spin unrestricted" 854 case default 855 write(msg,'(a,i3)')" Wrong method= ",Atm%method 856 ABI_BUG(msg) 857 end select 858 call wrtout(my_unt,msg,my_mode) 859 860 write(msg,'(7(a,i5,a),(a,f8.5,a))')& 861 & ' Number of spinorial components ...... ',Atm%nspinor,ch10,& 862 & ' Number of ind. spin polarizations ... ',Atm%nsppol,ch10,& 863 & ' Number of spin-density components ... ',Atm%nspden,ch10,& 864 & ' Maximum angular momentum + 1 ........ ',Atm%l_max,ch10,& 865 & ' Number of (l,n) orbitals ........... ',Atm%ln_size,ch10,& 866 & ' Number of (l,m,n) orbitals ......... ',Atm%lmn_size,ch10,& 867 & ' Dimensions of radial mesh ........... ',Atm%mesh_size,ch10,& 868 & ' Core Radius ........................ ',Atm%rcore,ch10 869 call wrtout(my_unt,msg,my_mode) 870 871 write(msg,'(2(a,f8.5,a))')& 872 & ' Ionic charge ........................ ',Atm%zion,ch10,& 873 & ' Atomic number ....................... ',Atm%znucl,ch10 874 call wrtout(my_unt,msg,my_mode) 875 876 do isppol=1,Atm%nsppol 877 do iln=1,Atm%ln_size 878 ll = Atm%indln(1,iln) 879 nn = Atm%indln(2,iln) 880 write(msg,'(" n=",i2,", l=",i2,", spin=",i2,", nocc=",f15.7,", energy=",f15.7,2x,"(",a,")")')& 881 & nn,ll,isppol,Atm%occ(iln,isppol),Atm%eig(iln,isppol),TRIM(my_mode2str(Atm%mode(iln,isppol))) 882 call wrtout(my_unt,msg,my_mode) 883 end do 884 end do 885 886 end subroutine print_atomorb