TABLE OF CONTENTS
- ABINIT/m_atom
- m_atom/atom_type
- m_atom/destroy_atom
- m_atom/get_atomcharge
- m_atom/get_overlap
- m_atom/init_atom
- m_atom/my_mode2str
- m_atom/print_atom
ABINIT/m_atom [ Modules ]
NAME
m_atom
FUNCTION
This module provides the definition of the atom_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_atom 26 27 use defs_basis 28 use m_profiling_abi 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_lmn_indices, 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_atom/atom_type [ Types ]
NAME
FUNCTION
Defines the atom_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 :: atom_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 atom_type 191 192 ! public procedures. 193 public :: init_atom 194 public :: destroy_atom 195 public :: print_atom 196 !public :: get_atomcharge 197 public :: get_overlap
m_atom/destroy_atom [ Functions ]
[ Top ] [ m_atom ] [ Functions ]
NAME
destroy_atom
FUNCTION
Free the dynamic memory allocated in a structure of type atom_type.
SIDE EFFECTS
Atm <type(atom_type)>=datastructure containing atomic orbitals for a given type of atom.
PARENTS
m_paw_slater
CHILDREN
wrtout
SOURCE
230 subroutine destroy_atom(Atm) 231 232 233 !This section has been created automatically by the script Abilint (TD). 234 !Do not modify the following lines by hand. 235 #undef ABI_FUNC 236 #define ABI_FUNC 'destroy_atom' 237 !End of the abilint section 238 239 implicit none 240 241 !Arguments ------------------------------------ 242 !scalars 243 type(atom_type),intent(inout) :: Atm 244 245 !************************************************************************ 246 247 !@atom_type 248 249 ! integers 250 if (allocated(Atm%indlmn)) then 251 ABI_FREE(Atm%indlmn) 252 end if 253 if (allocated(Atm%indln)) then 254 ABI_FREE(Atm%indln) 255 end if 256 if (allocated(Atm%indklmn)) then 257 ABI_FREE(Atm%indklmn) 258 end if 259 if (allocated(Atm%klm_diag)) then 260 ABI_FREE(Atm%klm_diag) 261 end if 262 if (allocated(Atm%klmntomn)) then 263 ABI_FREE(Atm%klmntomn) 264 end if 265 if (allocated(Atm%kln2ln)) then 266 ABI_FREE(Atm%kln2ln) 267 end if 268 if (allocated(Atm%mode)) then 269 ABI_FREE(Atm%mode) 270 end if 271 272 !real 273 if (allocated(Atm%eig)) then 274 ABI_FREE(Atm%eig) 275 end if 276 if (allocated(Atm%occ)) then 277 ABI_FREE(Atm%occ) 278 end if 279 if (allocated(Atm%phi)) then 280 ABI_FREE(Atm%phi) 281 end if 282 283 end subroutine destroy_atom
m_atom/get_atomcharge [ Functions ]
[ Top ] [ m_atom ] [ Functions ]
NAME
get_atomcharge
FUNCTION
Get core charge from a structure of type atom_type and optionally core density.
INPUTS
Atm<atom_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_atom
CHILDREN
wrtout
SOURCE
684 subroutine get_atomcharge(Atm,Radmesh,nele,radens) 685 686 687 !This section has been created automatically by the script Abilint (TD). 688 !Do not modify the following lines by hand. 689 #undef ABI_FUNC 690 #define ABI_FUNC 'get_atomcharge' 691 !End of the abilint section 692 693 implicit none 694 695 !Arguments ------------------------------------ 696 !scalars 697 real(dp),intent(out) :: nele 698 type(atom_type),intent(in) :: Atm 699 type(pawrad_type),intent(in) :: Radmesh 700 !arrays 701 real(dp),optional,intent(out) :: radens(Atm%mesh_size,Atm%nspden) 702 703 !Local variables------------------------------- 704 !scalars 705 integer :: iln,isppol 706 real(dp) :: intg,focc 707 real(dp),allocatable :: phi2nl(:) 708 709 !************************************************************************ 710 711 if (Atm%nsppol==2) then 712 MSG_ERROR("nsppol==2 is Working in progress") 713 end if 714 715 ABI_MALLOC(phi2nl,(Atm%mesh_size)) 716 if (PRESENT(radens)) radens = zero 717 718 nele = zero 719 do isppol=1,Atm%nsppol 720 do iln=1,Atm%ln_size 721 !Atm%mode(iln,isppol) TODO add option to select particular states 722 focc = Atm%occ(iln,isppol) 723 if (ABS(focc) > tol16) then 724 phi2nl = Atm%phi(:,iln,isppol)**2 725 call simp_gen(intg,phi2nl,Radmesh) 726 nele = nele + focc*intg 727 ! if (PRESENT(radens)) then !FIXME maybe it is better to rr**2 radens 728 ! radens(2:Atm%mesh_size) = radens(2:Atm%mesh_size) & 729 !& + focc * phi2nl(2:Atm%mesh_size)/(four_pi*Radmesh%rad(2:Atm%mesh_size)**2) 730 ! end if 731 end if 732 end do 733 end do 734 735 ABI_FREE(phi2nl) 736 737 end subroutine get_atomcharge
m_atom/get_overlap [ Functions ]
[ Top ] [ m_atom ] [ Functions ]
NAME
get_overlap
FUNCTION
Get overlap between core and valence states
INPUTS
Atm<atom_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
769 subroutine get_overlap(Atm,Atmesh,Radmesh2,isppol,nphi,phi,phi_indln,overlap) 770 771 772 !This section has been created automatically by the script Abilint (TD). 773 !Do not modify the following lines by hand. 774 #undef ABI_FUNC 775 #define ABI_FUNC 'get_overlap' 776 !End of the abilint section 777 778 implicit none 779 780 !Arguments ------------------------------------ 781 !scalars 782 integer,intent(in) :: nphi,isppol 783 type(atom_type),intent(in) :: Atm 784 type(pawrad_type),target,intent(in) :: Atmesh,Radmesh2 785 !arrays 786 integer,intent(in) :: phi_indln(2,nphi) 787 real(dp),target,intent(in) :: phi(Radmesh2%mesh_size,nphi) 788 real(dp),intent(out) :: overlap(Atm%ln_size,nphi) 789 790 !Local variables------------------------------- 791 !scalars 792 integer :: iln_atm,iphi,ll_phi,ll_atm,do_spline,iln 793 integer :: whichdenser,size4spl,my_mesh_size 794 real(dp) :: ybcbeg,ybcend,intg 795 logical :: hasameq 796 !arrays 797 real(dp),pointer :: ff_spl(:,:) 798 real(dp),allocatable :: der(:),ypp(:),func(:) 799 real(dp),pointer :: rad4spl(:),my_pts(:) 800 801 !************************************************************************ 802 803 ABI_CHECK(isppol>0.and.isppol<=Atm%nsppol,"Wrong isppol") 804 805 call pawrad_isame(Atmesh,Radmesh2,hasameq,whichdenser) 806 807 do_spline= 0; if (.not.hasameq) do_spline=1 808 809 my_mesh_size = MIN(Atmesh%mesh_size,Radmesh2%mesh_size) 810 ff_spl => phi 811 812 ! === Spline valence onto Atom mesh (natural spline) === 813 if (do_spline==1) then 814 MSG_COMMENT("Splining in overlap") 815 816 my_mesh_size = Atmesh%mesh_size 817 my_pts => Atmesh%rad(1:my_mesh_size) 818 ABI_MALLOC(ff_spl,(my_mesh_size,nphi)) 819 820 size4spl = Radmesh2%mesh_size 821 rad4spl => Radmesh2%rad 822 ABI_MALLOC(der,(size4spl)) 823 ABI_MALLOC(ypp,(size4spl)) 824 825 do iln=1,nphi 826 ypp(:) = zero; ybcbeg = zero; ybcend = zero 827 call spline(rad4spl,phi(:,iln),size4spl,ybcbeg,ybcend,ypp) 828 829 call splint(size4spl,rad4spl,phi(:,iln),ypp,my_mesh_size,my_pts,ff_spl(:,iln)) 830 end do 831 832 ABI_FREE(der) 833 ABI_FREE(ypp) 834 end if 835 836 ABI_MALLOC(func,(my_mesh_size)) 837 overlap = zero 838 839 do iphi=1,nphi 840 ll_phi = phi_indln(1,iphi) 841 do iln_atm=1,Atm%ln_size 842 ll_atm = Atm%indln(1,iln_atm) 843 844 if (ll_atm == ll_phi) then ! selection rule on l 845 func(:) = Atm%phi(1:my_mesh_size,iln_atm,isppol) * ff_spl(1:my_mesh_size,iphi) 846 call simp_gen(intg,func,Atmesh) 847 overlap(iln_atm,iphi)=intg 848 write(std_out,*)"overlap <phic_i|phi_j> for ll_phi",ll_phi,"ll_phic",ll_atm,"=",intg 849 end if 850 851 end do 852 end do 853 ABI_FREE(func) 854 855 if (do_spline==1) then 856 ABI_FREE(ff_spl) 857 end if 858 859 end subroutine get_overlap
m_atom/init_atom [ Functions ]
[ Top ] [ m_atom ] [ Functions ]
NAME
init_atom
FUNCTION
Initialize a structure of type atom_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<atom_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
313 subroutine init_atom(Atm,Atmrad,rcut,filename,prtvol,ierr) 314 315 316 !This section has been created automatically by the script Abilint (TD). 317 !Do not modify the following lines by hand. 318 #undef ABI_FUNC 319 #define ABI_FUNC 'init_atom' 320 use interfaces_14_hidewrite 321 !End of the abilint section 322 323 implicit none 324 325 !Arguments ------------------------------------ 326 !scalars 327 integer,intent(in) :: prtvol 328 integer,intent(out) :: ierr 329 real(dp),intent(in) :: rcut 330 character(len=*),intent(in) :: filename 331 type(atom_type),intent(out) :: Atm 332 type(pawrad_type),intent(out) :: Atmrad 333 334 !Local variables------------------------------- 335 !scalars 336 integer :: creatorid 337 integer :: imainmesh,jlmn,jl,k0lmn,ilmn,il 338 integer :: lcutdens,version,klmn,nn,ll 339 integer :: lmax,pspdat 340 integer :: ii,unt,iln,ir,nmesh,imsh 341 integer :: iread1,pspcod,msz,isppol 342 integer :: mesh_type,mesh_size 343 real(dp) :: charge 344 real(dp) :: my_rcore 345 real(dp) :: rstep,lstep 346 real(dp):: occ,ene 347 character(len=80) :: line,title 348 character(len=500) :: msg 349 !arrays 350 integer,allocatable :: orbitals(:) 351 real(dp),allocatable :: phitmp(:,:,:),radens(:,:) 352 type(pawrad_type),allocatable :: Radmesh(:) 353 354 ! ************************************************************************ 355 356 ! File format of formatted core orbitals. 357 358 !(1) title (character) line 359 !(2) method, nspinor,nsppol 360 !(3) znucl, zion, pspdat 361 !(4) ixc, lmax 362 !(5) version, creatorID 363 !(6) ln_size, lmn_size 364 !(7) orbitals (for i=1 to ln_size) 365 !(8) number_of_meshes 366 ! 367 !For imsh=1 to number_of_meshes 368 !(9) mesh_index, mesh_type ,mesh_size, rad_step[, log_step] 369 !(10) rcore (SPH) 370 ! 371 !For iln=1 to basis_size 372 ! (11) comment(character) 373 ! (12) radial mesh index for phi 374 ! (13) nn, ll, 375 ! (14) phi(r) (for ir=1 to phi_meshsz) 376 ! 377 !Comments: 378 ! * allowed values for method are: 379 ! 1 for restricted, compatible only with nsppol=1. 380 ! 2 for spin unrestricted, compatible only with nsppol=2. 381 !* psp_version= ID of PAW_psp version 382 ! 4 characters string of the form 'pawn' (with n varying) 383 !* creatorID= ID of psp generator 384 ! creatorid=1xyz : psp generated from Holzwarth AtomPAW generator version x.yz 385 ! creatorid=2xyz : psp generated from Vanderbilt ultra-soft generator version x.yz 386 ! creatorid=-1: psp for tests (for developpers only) 387 !* mesh_type= type of radial mesh 388 ! mesh_type=1 (regular grid): rad(i)=(i-1)*AA 389 ! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1) 390 ! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0 391 ! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n 392 ! -------------------------------------------------------------------------------- 393 394 !@atom_type 395 ierr=0 396 if (open_file(filename,msg,newunit=unt,form="formatted",status="old",action="read") /=0) then 397 MSG_WARNING(msg) 398 ierr=1; RETURN 399 end if 400 401 Atm%fname = filename 402 imainmesh = 1 403 404 !1) title 405 read(unt,'(a80)',ERR=10)title 406 write(msg,'(a,1x,a)')ch10,TRIM(title) 407 call wrtout(std_out,msg,'COLL') 408 409 read(unt,*,ERR=10)Atm%method, Atm%nspinor, Atm%nsppol 410 write(msg,'(3(i2,2x),22x,a)' )Atm%method, Atm%nspinor, Atm%nsppol,' method, nspinor, nsppol.' 411 call wrtout(std_out,msg,'COLL') 412 413 Atm%nspden = 1 !FIXME pass it as input 414 415 !2) 416 read(unt,*,ERR=10) Atm%znucl, Atm%zion, pspdat 417 write(msg,'(2f10.5,2x,i8,2x,a)' )Atm%znucl, Atm%zion, pspdat,' znucl, zion, pspdat' 418 call wrtout(std_out,msg,'COLL') 419 420 !3) 421 read(unt,*,ERR=10)pspcod,Atm%ixc,lmax 422 write(msg,'(2i5,2x,2x,a)')Atm%ixc,lmax,'ixc,lmax' 423 call wrtout(std_out,msg,'COLL') 424 425 Atm%l_max = lmax+1 426 Atm%l_size =2*lmax+1 427 428 !4) 429 !Read psp version in line 4 of the header 430 version=1 431 read(unt,'(a80)',ERR=10) line 432 line=ADJUSTL(line) 433 434 if (tolower(line(1:4))=="core") read(unit=line(5:80),fmt=*) version 435 if (version/=1) then 436 write(msg,'(a,i2,3a)')& 437 & 'This version of core psp file (',version,') is not compatible with',ch10,& 438 & 'the current version of Abinit.' 439 MSG_ERROR(msg) 440 end if 441 442 read(unit=line(6:80),fmt=*,ERR=10) creatorid 443 444 !========================================================== 445 446 ! 5) 447 read(unt,*,ERR=10)Atm%ln_size, Atm%lmn_size 448 449 Atm%ln2_size = Atm%ln_size *(Atm%ln_size +1)/2 450 Atm%lmn2_size = Atm%lmn_size*(Atm%lmn_size+1)/2 451 452 ! 6) 453 ABI_MALLOC(orbitals,(Atm%ln_size)) 454 read(unt,*,ERR=10) (orbitals(iln), iln=1,Atm%ln_size) 455 456 lmax = MAXVAL(orbitals) 457 if (lmax+1/=Atm%l_max) then 458 write(msg,'(a)')" lmax read from file does not agree with orbitals. " 459 MSG_ERROR(msg) 460 end if 461 462 ! 7) 463 read(unt,*)nmesh 464 ABI_DT_MALLOC(Radmesh,(nmesh)) 465 do imsh=1,nmesh 466 lstep=zero 467 read(unt,'(a80)') line 468 read(unit=line,fmt=*,err=20,end=20) ii,mesh_type,mesh_size,rstep,lstep 469 20 continue 470 ABI_CHECK(ii<=nmesh,"Index of mesh out of range!") 471 call pawrad_init(Radmesh(ii),mesh_size,mesh_type,rstep,lstep,-one) 472 end do 473 474 ! 8) 475 read(unt,*,ERR=10) my_rcore 476 477 Atm%rcore = my_rcore 478 if (rcut > tol6) then 479 Atm%rcore = rcut 480 write(msg,'(a,f8.5,a,f8.5,a)')& 481 & " Truncating radial mesh for core orbitals using new rcore = ",Atm%rcore," (",my_rcore,")" 482 call wrtout(std_out,msg,'COLL') 483 if (rcut > my_rcore) then 484 MSG_ERROR("rcut should not exceed my_rcore") 485 end if 486 end if 487 488 !========================================================== 489 ! Mirror pseudopotential parameters to the output and log files 490 491 write(msg,'(a,i1)')' Pseudopotential format is: core ',version 492 call wrtout(std_out,msg,'COLL') 493 write(msg,'(2(a,i3),a,64i4)')' (ln_size)= ',Atm%ln_size,' (lmn_size= ',Atm%lmn_size,'), orbitals= ',orbitals(1:Atm%ln_size) 494 call wrtout(std_out,msg,'COLL') 495 write(msg,'(a,f11.8)')' Radius used for core orbitals: rc_sph= ',Atm%rcore 496 call wrtout(std_out,msg,'COLL') 497 write(msg,'(a,i1,a)')' ',nmesh,' radial meshes are used:' 498 call wrtout(std_out,msg,'COLL') 499 500 do imsh=1,nmesh 501 call pawrad_print(Radmesh(imsh),prtvol=prtvol) 502 end do 503 504 !========================================================== 505 506 !--------------------------------- 507 ! (11-12-13) Read wave-functions (phi) 508 ! * Initialize also the mesh. 509 ABI_MALLOC(Atm%indln,(2,Atm%ln_size)) 510 ABI_MALLOC(Atm%eig,(Atm%ln_size,Atm%nsppol)) 511 ABI_MALLOC(Atm%occ,(Atm%ln_size,Atm%nsppol)) 512 Atm%occ = zero 513 514 do isppol=1,Atm%nsppol 515 do iln=1,Atm%ln_size 516 read(unt,*,ERR=10) line 517 read(unt,*,ERR=10) iread1 518 if (iln==1.and.isppol==1) then 519 imainmesh=iread1 520 Atm%mesh_size = Radmesh(iread1)%mesh_size 521 ABI_MALLOC(Atm%phi,(Atm%mesh_size,Atm%ln_size,Atm%nsppol)) 522 else if (iread1/=imainmesh) then 523 write(msg,'(3a)')& 524 & ' All Phi core must be given on the same radial mesh !',ch10,& 525 & ' Action: check your pseudopotential file.' 526 MSG_ERROR(msg) 527 end if 528 read(unt,*,ERR=10)nn,ll,ii 529 ABI_CHECK(ii==isppol,"Wrong spin index") 530 read(unt,*,ERR=10)ene,occ 531 Atm%indln(1,iln) = ll 532 Atm%indln(2,iln) = nn 533 Atm%occ(iln,isppol) = occ 534 Atm%eig(iln,isppol) = ene 535 read(unt,*,ERR=10) (Atm%phi(ir,iln,isppol),ir=1,Radmesh(imainmesh)%mesh_size) 536 !do ir=1,Radmesh(imainmesh)%mesh_size 537 ! write(100+iln,*)Radmesh(imainmesh)%rad(ir),Atm%phi(ir,iln,isppol) 538 !end do 539 end do 540 end do !isppol 541 542 close(unt) 543 ! -------------------- END OF READING ------------------------------- 544 545 if ( rcut < tol16) then 546 ! * Use full mesh reported on file. 547 call pawrad_init(Atmrad,mesh_size=Radmesh(imainmesh)%mesh_size,mesh_type=Radmesh(imainmesh)%mesh_type,& 548 & rstep=Radmesh(imainmesh)%rstep,lstep=Radmesh(imainmesh)%lstep, & 549 !& r_for_intg=Atm%rcore) ! check this 550 & r_for_intg=-one) 551 else 552 ! * Truncate orbitals and radial mesh within rcut. 553 msz = MIN(pawrad_ifromr(Radmesh(imainmesh),rcut)+6, Radmesh(imainmesh)%mesh_size) ! add six more points 554 555 write(msg,'(a,f8.5,a,f8.5,a)')& 556 & " Truncating radial mesh for core orbitals ",Radmesh(imainmesh)%rad(msz),"(",my_rcore,")" 557 call wrtout(std_out,msg,'COLL') 558 Atm%rcore = Radmesh(imainmesh)%rad(msz) 559 560 ABI_MALLOC(phitmp,(msz,Atm%ln_size,Atm%nsppol)) 561 phitmp = Atm%phi(1:msz,:,:) 562 563 ABI_FREE(Atm%phi) 564 ABI_MALLOC(Atm%phi,(msz,Atm%ln_size,Atm%nsppol)) 565 Atm%phi = phitmp 566 ABI_FREE(phitmp) 567 568 ! * Compute new mesh for core orbitals. 569 mesh_type = Radmesh(imainmesh)%mesh_type 570 mesh_size = msz !radmesh(imainmesh)%mesh_size 571 rstep = Radmesh(imainmesh)%rstep 572 lstep = Radmesh(imainmesh)%lstep 573 574 call pawrad_init(Atmrad,mesh_size=mesh_size,mesh_type=mesh_type,rstep=rstep,lstep=lstep,& 575 !& r_for_intg=Atm%rcore) 576 & r_for_intg=-one) 577 578 !do isppol=1,Atm%nsppol 579 ! do iln=1,Atm%ln_size 580 ! do ir=1,Atmrad%mesh_size 581 ! write(200+iln,*)Atmrad%rad(ir),Atm%phi(ir,iln,isppol) 582 ! end do 583 ! end do 584 !end do 585 586 end if 587 588 Atm%mesh_size = Atmrad%mesh_size 589 call pawrad_print(Atmrad,header="Final mesh",prtvol=prtvol) 590 591 ABI_MALLOC(radens,(Atm%mesh_size,Atm%nspden)) 592 call get_atomcharge(Atm,Atmrad,charge,radens=radens) 593 594 write(std_out,*)"core charge = ",charge 595 !do ii=1,Atmrad%mesh_size 596 ! write(77,*)Atmrad%rad(ii),(radens(ii,isppol),isppol=1,Atm%nspden) 597 !end do 598 ABI_FREE(radens) 599 600 !========================================================== 601 ! Free temporary allocated space 602 603 call pawrad_free(Radmesh) 604 ABI_DT_FREE(Radmesh) 605 606 ! * Setup of indlmn 607 ABI_MALLOC(Atm%indlmn,(6,Atm%lmn_size)) 608 call make_indlmn(Atm%ln_size, Atm%lmn_size, orbitals, Atm%indlmn) 609 610 ABI_FREE(orbitals) 611 612 ! * Setup of indklmn and klm_diag. 613 lcutdens=HUGE(1) 614 ABI_MALLOC(Atm%indklmn,(8,Atm%lmn2_size)) 615 ABI_MALLOC(Atm%klm_diag,(Atm%lmn2_size)) 616 617 call make_indklmn(lcutdens, Atm%lmn_size, Atm%lmn2_size, Atm%indlmn, Atm%indklmn, Atm%klm_diag) 618 619 ! * Setup of klmntomn. 620 ABI_MALLOC(Atm%klmntomn,(4,Atm%lmn2_size)) 621 622 do jlmn=1,Atm%lmn_size 623 jl= Atm%indlmn(1,jlmn) 624 k0lmn=jlmn*(jlmn-1)/2 625 do ilmn=1,jlmn 626 il= Atm%indlmn(1,ilmn) 627 klmn=k0lmn+ilmn 628 Atm%klmntomn(1,klmn) = Atm%indlmn(2,ilmn)+il+1 ! im 629 Atm%klmntomn(2,klmn) = Atm%indlmn(2,jlmn)+jl+1 ! jm 630 Atm%klmntomn(3,klmn) = Atm%indlmn(3,ilmn) ! in 631 Atm%klmntomn(4,klmn) = Atm%indlmn(3,jlmn) ! jn 632 !write(std_out,*) jlmn,ilmn,Atm%klmntomn(:,klmn) 633 end do 634 end do 635 636 ! * Setup of kln2ln. 637 !TODO this has to be tested 638 ABI_MALLOC(Atm%kln2ln,(6,Atm%ln2_size)) 639 call make_kln2ln(Atm%lmn_size,Atm%lmn2_size,Atm%ln2_size,Atm%indlmn,Atm%indklmn,Atm%kln2ln) 640 641 ! * Treat all states as core. 642 ABI_MALLOC(Atm%mode,(Atm%ln_size,Atm%nsppol)) 643 Atm%mode = ORB_FROZEN 644 645 !call print_atom(Atm,prtvol=prtvol) 646 647 return 648 ! 649 ! === Propagate the error === 650 10 continue 651 ierr=2 652 MSG_WARNING("Wrong file format") 653 return 654 655 end subroutine init_atom
m_atom/my_mode2str [ Functions ]
[ Top ] [ m_atom ] [ 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
979 function my_mode2str(mode) result(str) 980 981 982 !This section has been created automatically by the script Abilint (TD). 983 !Do not modify the following lines by hand. 984 #undef ABI_FUNC 985 #define ABI_FUNC 'my_mode2str' 986 !End of the abilint section 987 988 implicit none 989 990 !Arguments ------------------------------------ 991 !scalars 992 integer,intent(in) :: mode 993 character(len=50) :: str 994 995 !Local variables 996 character(len=500) :: msg 997 998 !************************************************************************ 999 1000 select case (mode) 1001 case (ORB_FROZEN) 1002 str="Frozen Orbital" 1003 case (ORB_RELAXED_CORE) 1004 str="Relaxed Core Orbital" 1005 case (ORB_VALENCE) 1006 str="Valence Orbital" 1007 case default 1008 write(msg,'(a,i3)')" Wrong mode= ",mode 1009 MSG_BUG(msg) 1010 end select 1011 1012 end function my_mode2str
m_atom/print_atom [ Functions ]
[ Top ] [ m_atom ] [ Functions ]
NAME
print_atom
FUNCTION
Reports info on a structure of type atom_type.
INPUTS
Atm <type(atom_type)>=datastructure containing atomic orbitals for a given type of atom.
OUTPUT
PARENTS
m_paw_slater
CHILDREN
wrtout
SOURCE
884 subroutine print_atom(Atm,header,unit,prtvol,mode_paral) 885 886 887 !This section has been created automatically by the script Abilint (TD). 888 !Do not modify the following lines by hand. 889 #undef ABI_FUNC 890 #define ABI_FUNC 'print_atom' 891 use interfaces_14_hidewrite 892 !End of the abilint section 893 894 implicit none 895 896 !Arguments ------------------------------------ 897 !scalars 898 type(atom_type),intent(in) :: Atm 899 integer,optional,intent(in) :: prtvol,unit 900 character(len=*),optional,intent(in) :: header 901 character(len=4),optional,intent(in) :: mode_paral 902 903 !Local variables------------------------------- 904 integer :: my_unt,my_prtvol,iln,ll,nn,isppol 905 character(len=4) :: my_mode 906 character(len=500) :: msg 907 ! ************************************************************************ 908 909 !@atom_type 910 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 911 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 912 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 913 914 msg=' ==== Info on the atom_type ==== ' 915 if (PRESENT(header)) msg=header 916 call wrtout(my_unt,msg,my_mode) 917 918 select case (Atm%method) 919 case (1) 920 msg = " Spin restricted" 921 case(2) 922 msg = " Spin unrestricted" 923 case default 924 write(msg,'(a,i3)')" Wrong method= ",Atm%method 925 MSG_BUG(msg) 926 end select 927 call wrtout(my_unt,msg,my_mode) 928 929 write(msg,'(7(a,i5,a),(a,f8.5,a))')& 930 & ' Number of spinorial components ...... ',Atm%nspinor,ch10,& 931 & ' Number of ind. spin polarizations ... ',Atm%nsppol,ch10,& 932 & ' Number of spin-density components ... ',Atm%nspden,ch10,& 933 & ' Maximum angular momentum + 1 ........ ',Atm%l_max,ch10,& 934 & ' Number of (l,n) orbitals ........... ',Atm%ln_size,ch10,& 935 & ' Number of (l,m,n) orbitals ......... ',Atm%lmn_size,ch10,& 936 & ' Dimensions of radial mesh ........... ',Atm%mesh_size,ch10,& 937 & ' Core Radius ........................ ',Atm%rcore,ch10 938 call wrtout(my_unt,msg,my_mode) 939 940 write(msg,'(2(a,f8.5,a))')& 941 & ' Ionic charge ........................ ',Atm%zion,ch10,& 942 & ' Atomic number ....................... ',Atm%znucl,ch10 943 call wrtout(my_unt,msg,my_mode) 944 945 do isppol=1,Atm%nsppol 946 do iln=1,Atm%ln_size 947 ll = Atm%indln(1,iln) 948 nn = Atm%indln(2,iln) 949 write(msg,'(" n=",i2,", l=",i2,", spin=",i2,", nocc=",f15.7,", energy=",f15.7,2x,"(",a,")")')& 950 & nn,ll,isppol,Atm%occ(iln,isppol),Atm%eig(iln,isppol),TRIM(my_mode2str(Atm%mode(iln,isppol))) 951 call wrtout(my_unt,msg,my_mode) 952 end do 953 end do 954 955 end subroutine print_atom