TABLE OF CONTENTS
ABINIT/m_paw_tools [ Modules ]
NAME
m_paw_tools
FUNCTION
This module contains miscelaneous routines used in the PAW context.
COPYRIGHT
Copyright (C) 2018-2024 ABINIT group (FJ,MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 MODULE m_paw_tools 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_xmpi 28 use m_dtset 29 30 use m_paral_atom, only : get_my_atmtab, free_my_atmtab 31 use m_electronpositron, only : electronpositron_type,electronpositron_calctype,EP_POSITRON 32 use m_pawang, only : pawang_type 33 use m_pawtab, only : pawtab_type 34 use m_paw_ij, only : paw_ij_type, paw_ij_free, paw_ij_nullify, paw_ij_gather 35 use m_pawdij, only : pawdij_print_dij 36 use m_pawrhoij, only : pawrhoij_type, pawrhoij_free, pawrhoij_gather, pawrhoij_nullify, & 37 & pawrhoij_print_rhoij 38 use m_paw_io, only : pawio_print_ij 39 use m_paw_sphharm, only : mat_mlms2jmj, mat_slm2ylm 40 use m_paw_correlations, only : setnoccmmp 41 42 implicit none 43 44 private 45 46 !public procedures. 47 public :: chkpawovlp 48 public :: pawprt 49 50 CONTAINS !========================================================================================
m_paw_tools/chkpawovlp [ Functions ]
[ Top ] [ m_paw_tools ] [ Functions ]
NAME
chkpawovlp
FUNCTION
Verify that the PAW spheres are not overlapping
INPUTS
natom=number of atoms in cell. nremit [optional] = if non-zero initialize the number of possible remits before stop ntypat=number of types of atoms in unit cell. pawovlp=percentage of voluminal overlap ratio allowed to continue execution (if negative value, execution always continues) pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data: rmet(3,3)=real space metric ($\textrm{bohr}^{2}$). typat(natom)=type (integer) for each atom xred(3,natom)=reduced dimensionless atomic coordinates
OUTPUT
(only checking)
NOTES
SOURCE
80 subroutine chkpawovlp(natom,ntypat,pawovlp,pawtab,rmet,typat,xred,nremit) 81 82 !Arguments --------------------------------------------- 83 !scalars 84 integer,intent(in) :: natom,ntypat 85 integer,intent(in),optional :: nremit 86 real(dp) :: pawovlp 87 !arrays 88 integer,intent(in) :: typat(natom) 89 real(dp),intent(in) :: rmet(3,3),xred(3,natom) 90 type(pawtab_type),intent(in) :: pawtab(ntypat) 91 92 !Local variables --------------------------------------- 93 !scalars 94 integer :: decrease_nremit,ia,ib,ii,t1,t2,t3 95 integer,save :: nremit_counter=0 96 logical :: stop_on_error 97 real(dp) :: dd,dif1,dif2,dif3,ha,hb,norm2 98 real(dp) :: ratio_percent,va,vb,vv 99 character(len=750) :: message 100 !arrays 101 integer :: iamax(2),ibmax(2),iovl(2) 102 real(dp) :: norm2_min(2),r2cut(2),ratio_percent_max(2),rcuta(2),rcutb(2) 103 104 105 ! ************************************************************************* 106 107 DBG_ENTER("COLL") 108 109 ! if(present(nremit))then 110 ! if(nremit/=0)nremit_counter=abs(nremit) 111 ! else 112 ! nremit_counter=0 113 ! endif 114 115 !DEBUG 116 ! write(std_out,'(a,a,i4)')ch10,' m_paw_tools, chkpawovlp : enter, saved nremit_counter=',nremit_counter 117 !ENDDEBUG 118 119 if(present(nremit))then 120 if(nremit/=0)nremit_counter=abs(nremit) 121 !DEBUG 122 ! write(std_out,'(a,i4)')' m_paw_tools, chkpawovlp : optional arg nremit present, nremit=',nremit 123 !ENDDEBUG 124 else 125 nremit_counter=0 126 endif 127 !DEBUG 128 ! write(std_out,'(a,i4)')' m_paw_tools, chkpawovlp : after init, nremit_counter=',nremit_counter 129 !ENDDEBUG 130 131 132 iamax(:)=-1;ibmax(:)=-1 133 norm2_min(:)=-1.d0;ratio_percent_max(:)=-1.d0 134 iovl(:)=0 135 136 !Loop on "overlapping" atoms with the maximum overlap 137 do ia=1,natom 138 139 rcuta(1)=pawtab(typat(ia))%rpaw 140 rcuta(2)=pawtab(typat(ia))%rshp 141 142 do ib=ia,natom 143 144 rcutb(1)=pawtab(typat(ib))%rpaw 145 rcutb(2)=pawtab(typat(ib))%rshp 146 r2cut(1)=(rcuta(1)+rcutb(1))**2 147 r2cut(2)=(rcuta(2)+rcutb(2))**2 148 149 ! Visit the box and its first images: 150 do t3=-1,1 151 do t2=-1,1 152 do t1=-1,1 153 154 dif1=xred(1,ia)-(xred(1,ib)+dble(t1)) 155 dif2=xred(2,ia)-(xred(2,ib)+dble(t2)) 156 dif3=xred(3,ia)-(xred(3,ib)+dble(t3)) 157 norm2=sqnrm_pawovlp(dif1,dif2,dif3) 158 159 do ii=1,2 160 161 if(norm2>tol10.and.norm2<r2cut(ii)) then 162 163 iovl(ii)=iovl(ii)+1 164 165 ! Compute the overlap ratio: 166 dd=sqrt(norm2) 167 va=4._dp/3._dp*pi*rcuta(ii)**3 168 vb=4._dp/3._dp*pi*rcutb(ii)**3 169 ha=(rcutb(ii)**2-(dd-rcuta(ii))**2)/(two*dd) 170 hb=(rcuta(ii)**2-(dd-rcutb(ii))**2)/(two*dd) 171 vv=pi/3.d0*(ha**2*(three*rcuta(ii)-ha)+hb**2*(three*rcutb(ii)-hb)) 172 ratio_percent=100._dp*min(vv/min(va,vb),one) 173 if (ratio_percent>ratio_percent_max(ii)) then 174 ratio_percent_max(ii)=ratio_percent 175 norm2_min(ii)=norm2 176 iamax(ii)=ia;ibmax(ii)=ib 177 end if 178 179 end if 180 end do 181 end do 182 end do 183 end do 184 end do 185 end do 186 187 !DEBUG 188 !write(std_out,'(a,f8.4)')' chkpawovlp : maxval(ratio_percent_max(1:2))=',maxval(ratio_percent_max(1:2)) 189 !ENDDEBUG 190 191 stop_on_error=(abs(pawovlp)<=tol6.or.(pawovlp>tol6.and.(maxval(ratio_percent_max(1:2))>pawovlp))) 192 decrease_nremit=0 193 194 !Print adapted message with overlap value 195 if (iovl(1)+iovl(2)>0) then 196 197 do ii=1,2 198 199 if(ratio_percent_max(ii)>zero)then 200 if (ii==1) write(message,' (a)' ) 'PAW SPHERES ARE OVERLAPPING!' 201 if (ii==2) write(message, '(a)' ) 'PAW COMPENSATION DENSITIES ARE OVERLAPPING !' 202 203 if (iovl(ii)==1) then 204 write(message, '(3a)' ) trim(message),ch10,& 205 & ' There is one pair of overlapping atoms.' 206 else 207 write(message, '(3a,i5,a)' ) trim(message),ch10,& 208 & ' There are ', iovl(ii),' pairs of overlapping atoms.' 209 end if 210 write(message, '(3a,i3,a,i3,a)' ) trim(message),ch10,& 211 ' The maximum overlap percentage is obtained for the atoms ',iamax(ii),' and ',ibmax(ii),'.' 212 write(message, '(2a,2(a,i3),a,f9.5)' ) trim(message),ch10,& 213 & ' | Distance between atoms ',iamax(ii),' and ',ibmax(ii),' is : ',sqrt(norm2_min(ii)) 214 if(ii==1)then 215 write(message, '(2a,2(a,i3,a,f9.5,a))' ) trim(message),ch10,& 216 & ' | PAW radius of the sphere around atom ',iamax(ii),' is: ',pawtab(typat(iamax(ii)))%rpaw,ch10,& 217 & ' | PAW radius of the sphere around atom ',ibmax(ii),' is: ',pawtab(typat(ibmax(ii)))%rpaw,ch10 218 else if(ii==2)then 219 write(message, '(2a,2(a,i3,a,f9.5,a))' ) trim(message),ch10,& 220 & ' | Radius of the compensation sphere around atom ',iamax(ii),' is: ',pawtab(typat(iamax(ii)))%rshp,ch10,& 221 & ' | Radius of the compensation sphere around atom ',ibmax(ii),' is: ',pawtab(typat(ibmax(ii)))%rshp,ch10 222 endif 223 write(message, '(2a,f7.4,a)' ) trim(message),& 224 & ' | This leads to a (voluminal) overlap ratio of ',ratio_percent_max(ii),' %' 225 if (ii==1) then 226 write(message, '(3a)' ) trim(message),ch10,& 227 & 'THIS IS DANGEROUS, as PAW formalism assumes non-overlapping PAW spheres.' 228 else if (ii==2) then 229 write(message, '(3a)' ) trim(message),ch10,& 230 & 'THIS IS DANGEROUS, as PAW formalism assumes non-overlapping compensation densities.' 231 end if 232 if (stop_on_error .and. nremit_counter==0) then 233 ABI_ERROR_NOSTOP(message,ia) !ia is dummy 234 else 235 ABI_WARNING(message) 236 if(stop_on_error .and. nremit_counter/=0)decrease_nremit=1 237 end if 238 endif ! ratio_percent_max(ii)>zero 239 240 enddo ! ii 241 242 ! Print advice 243 if (stop_on_error) then 244 write(message, '(3a)' )& 245 & ' Action: 1- decrease cutoff radius of PAW dataset',ch10,& 246 & ' OR 2- ajust "pawovlp" input variable to allow overlap (risky)' 247 if(nremit_counter==0)then 248 ABI_ERROR(message) 249 endif 250 end if 251 252 nremit_counter=nremit_counter-decrease_nremit 253 254 ! Print last message if execution continues: 255 if (pawovlp<=tol6) then 256 write(message, '(6a)' ) & 257 & ' Results might be approximate,',ch10,& 258 & ' and even inaccurate (if overlap is too big) !',ch10,& 259 & ' Assume experienced user. Execution will continue.',ch10 260 call wrtout(std_out,message,'COLL') 261 else if (ratio_percent_max(1)<=pawovlp .and. ratio_percent_max(2)<=pawovlp) then 262 write(message, '(8a)' ) & 263 & ' Overlap ratio seems to be acceptable (less than value',ch10,& 264 & ' of "pawovlp" input parameter): execution will continue.',ch10,& 265 & ' But be aware that results might be approximate,',ch10,& 266 & ' and even inaccurate (depending on your physical system) !',ch10 267 call wrtout(std_out,message,'COLL') 268 else if(decrease_nremit==1)then 269 write(message, '(3a)' ) & 270 & ' First time that overlap is bigger than "pawovlp" input parameter.',ch10,& 271 & ' Execution will continue, but such overlap will not be tolerated twice.' 272 call wrtout(std_out,message,'COLL') 273 end if 274 275 end if !iovl>0 276 277 DBG_EXIT("COLL") 278 279 contains 280 281 function sqnrm_pawovlp(u1,u2,u3) 282 !squared norm of a vector 283 real(dp) :: sqnrm_pawovlp 284 real(dp),intent(in) :: u1,u2,u3 285 286 sqnrm_pawovlp=rmet(1,1)*u1*u1+rmet(2,1)*u2*u1+rmet(3,1)*u3*u1& 287 & +rmet(1,2)*u1*u2+rmet(2,2)*u2*u2+rmet(3,2)*u3*u2& 288 & +rmet(1,3)*u1*u3+rmet(2,3)*u2*u3+rmet(3,3)*u3*u3 289 290 end function sqnrm_pawovlp 291 292 end subroutine chkpawovlp
m_paw_tools/pawprt [ Functions ]
[ Top ] [ m_paw_tools ] [ Functions ]
NAME
pawprt
FUNCTION
Print out data concerning PAW formalism (pseudopotential strength, augmentation occupancies...) To be called at the end of the SCF cycle
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (FJ,MT,BA) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
INPUTS
dtset <type(dataset_type)>=all input variables in this dataset | enunit=parameter determining units of output energies | kptopt=option for the generation of k points | natom=number of atoms in cell | ntypat = number of atom types | pawprtvol= printing volume | pawspnorb=flag: 1 if spin-orbit coupling is activated | typat(natom)=type of each atom electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument) mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc comm_atom=--optional-- MPI communicator over atoms my_natom=number of atoms treated by current processor paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
(only printing)
SOURCE
335 subroutine pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,& 336 & electronpositron,& ! optional argument 337 & mpi_atmtab,comm_atom) ! optional arguments (parallelism) 338 339 !Arguments ------------------------------------ 340 !scalars 341 integer,intent(in) :: my_natom 342 integer,optional,intent(in) :: comm_atom 343 type(dataset_type),intent(in) :: dtset 344 type(electronpositron_type),pointer,optional :: electronpositron 345 !arrays 346 integer,optional,target,intent(in) :: mpi_atmtab(:) 347 type(paw_ij_type),target,intent(inout) :: paw_ij(my_natom) 348 type(pawrhoij_type),target,intent(inout) :: pawrhoij(my_natom) 349 type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat) 350 351 !Local variables------------------------------- 352 !scalars 353 integer,parameter :: natmax=2 354 integer :: cplex_dij,group1,group2,iat,iatom,ierr,ii,im1,im2,ipositron,ispden 355 integer :: i_unitfi,itypat,ll,llp,me_atom,my_comm_atom,natprt,ndij,nspden,nsppol 356 integer :: unitfi,unt 357 real(dp) :: mnorm,mx,my,mz,ntot,valmx,localm 358 logical :: my_atmtab_allocated,paral_atom,useexexch,usepawu 359 type(pawang_type):: pawang_dum 360 character(len=7),parameter :: dspin1(6)=(/"up ","down ","up-up ","dwn-dwn","up-dwn ","dwn-up "/) 361 character(len=8),parameter :: dspin2(6)=(/"up ","down ","dens (n)","magn (x)","magn (y)","magn (z)"/) 362 character(len=500) :: msg 363 !arrays 364 integer :: idum(1) 365 integer :: idum1(0),idum3(0,0,0) 366 integer,allocatable :: jatom(:) 367 integer,pointer :: my_atmtab(:) 368 real(dp) :: rdum2(0,0),rdum4(0,0,0,0) 369 complex(dpc),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:) 370 type(paw_ij_type), ABI_CONTIGUOUS pointer :: paw_ij_all(:) 371 type(pawrhoij_type),ABI_CONTIGUOUS pointer :: pawrhoij_all(:) 372 373 ! ********************************************************************* 374 375 DBG_ENTER("COLL") 376 377 !Set up parallelism over atoms 378 paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom)) 379 nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab 380 my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom 381 call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom) 382 me_atom=xmpi_comm_rank(my_comm_atom) 383 384 !Continue only if comm_atom contains the master of the output comm 385 if (paral_atom) then 386 call xmpi_comm_group(abinit_comm_output,group1,ierr) 387 call xmpi_comm_group(my_comm_atom,group2,ierr) 388 call xmpi_group_translate_ranks(group1,1,(/0/),group2,idum,ierr) 389 call xmpi_group_free(group1) 390 call xmpi_group_free(group2) 391 if (idum(1)==xmpi_undefined) then 392 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 393 return 394 end if 395 end if 396 397 !Initializations 398 natprt=natmax;if (dtset%natom==1) natprt=1 399 if (dtset%pawprtvol<0) natprt=dtset%natom 400 ABI_MALLOC(jatom,(natprt)) 401 if (natprt==1) then 402 jatom(1)=1 403 else if (natprt==2) then 404 jatom(1)=1;jatom(2)=dtset%natom 405 else if (natprt==dtset%natom) then 406 do iat=1,dtset%natom 407 jatom(iat)=iat 408 end do 409 else 410 ABI_BUG("invalid value of natprt!") 411 end if 412 usepawu=(count(pawtab(:)%usepawu/=0)>0) 413 useexexch=(count(pawtab(:)%useexexch/=0)>0) 414 ipositron=0 415 if (present(electronpositron)) then 416 if (associated(electronpositron)) ipositron=electronpositron%calctype 417 end if 418 419 !Main title 420 write(msg, '(2a)' ) ch10,& 421 & ' ==== Results concerning PAW augmentation regions ====' 422 call wrtout(ab_out,msg,'COLL') 423 call wrtout(std_out,msg,'COLL') 424 msg=' ' 425 call wrtout(ab_out,msg,'COLL') 426 call wrtout(std_out,msg,'COLL') 427 428 !If atomic data are distributed, retrieve all Dij on master proc 429 if (paral_atom) then 430 if (me_atom==0) then 431 ABI_MALLOC(paw_ij_all,(dtset%natom)) 432 call paw_ij_nullify(paw_ij_all) 433 else 434 ABI_MALLOC(paw_ij_all,(0)) 435 end if 436 call paw_ij_gather(paw_ij,paw_ij_all,0,my_comm_atom) 437 else 438 paw_ij_all => paw_ij 439 end if 440 441 !Print out pseudopotential strength 442 !---------------------------------- 443 if (me_atom==0) then 444 do i_unitfi=1,2 445 unitfi=ab_out;if (i_unitfi==2) unitfi=std_out 446 do unt=1,2 447 if (((unt==1).and.(dtset%enunit==0.or.dtset%enunit==2)).or.& 448 & ((unt==2).and.(dtset%enunit==1.or.dtset%enunit==2))) then 449 if ((unt==1).and.(dtset%enunit==0.or.dtset%enunit==2)) then 450 write(msg,'(a)') ' Total pseudopotential strength Dij (hartree):' 451 else if ((unt==2).and.(dtset%enunit==1.or.dtset%enunit==2)) then 452 write(msg,'(a)') ' Total pseudopotential strength Dij (eV):' 453 end if 454 call wrtout(unitfi,msg,'COLL') 455 if (ipositron>0) then 456 if (electronpositron%has_pos_ham==0) then 457 write(msg,'(a)') ' -Note: these are the electronic Dij' 458 else 459 write(msg,'(a)') ' -Note: these are the positronic Dij' 460 end if 461 call wrtout(unitfi,msg,'COLL') 462 end if 463 valmx=100._dp;if (ipositron>0) valmx=-1._dp 464 do iat=1,natprt 465 iatom=jatom(iat) 466 call pawdij_print_dij(paw_ij_all(iatom)%dij,paw_ij_all(iatom)%cplex_dij,& 467 & paw_ij_all(iatom)%qphase,iatom,dtset%natom,paw_ij_all(iatom)%nspden,& 468 & test_value=valmx,unit=unitfi,Ha_or_eV=unt,opt_prtvol=dtset%pawprtvol) 469 end do 470 end if 471 msg=' ' 472 call wrtout(unitfi,msg,'COLL') 473 end do 474 end do 475 end if 476 if (paral_atom.and.(.not.usepawu).and.(.not.useexexch)) then 477 call paw_ij_free(paw_ij_all) 478 ABI_FREE(paw_ij_all) 479 end if 480 481 !If atomic data are distributed, retrieve all Rhoij on master proc 482 if (paral_atom) then 483 if (me_atom==0) then 484 ABI_MALLOC(pawrhoij_all,(dtset%natom)) 485 else 486 ABI_MALLOC(pawrhoij_all,(0)) 487 end if 488 call pawrhoij_nullify(pawrhoij_all) 489 call pawrhoij_gather(pawrhoij,pawrhoij_all,0,my_comm_atom,& 490 & with_grhoij=.false.,with_lmnmix=.false.,& 491 & with_rhoij_=.false.,with_rhoijres=.false.) 492 else 493 pawrhoij_all => pawrhoij 494 end if 495 496 !Print out SYMMETRIZED occupancies of the partial waves 497 !------------------------------------------------------ 498 if (me_atom==0) then 499 do i_unitfi=1,2 500 unitfi=ab_out;if (i_unitfi==2) unitfi=std_out 501 write(msg,'(a)') ' Augmentation waves occupancies Rhoij:' 502 call wrtout(unitfi,msg,'COLL') 503 if (ipositron>0) then 504 if (electronpositron%particle==EP_POSITRON) then 505 write(msg,'(a)') ' -Note: these are the electronic Rhoij' 506 else 507 write(msg,'(a)') ' -Note: these are the positronic Rhoij' 508 end if 509 call wrtout(unitfi,msg,'COLL') 510 end if 511 if (dtset%pawspnorb>0.and.pawrhoij_all(1)%cplex_rhoij==1.and.dtset%kptopt/=1.and.dtset%kptopt/=2) then 512 write(msg,'(6a)') ' pawprt: - WARNING:',ch10,& 513 & ' Spin-orbit coupling is activated but only real part of Rhoij occupancies',ch10,& 514 & ' has been computed; they could have an imaginary part (not printed here).' 515 call wrtout(unitfi,msg,'COLL') 516 end if 517 valmx=25._dp;if (ipositron>0) valmx=-1._dp 518 do iat=1,natprt 519 iatom=jatom(iat);nspden=pawrhoij_all(iatom)%nspden 520 call pawrhoij_print_rhoij(pawrhoij_all(iatom)%rhoijp,pawrhoij_all(iatom)%cplex_rhoij,& 521 & pawrhoij_all(iatom)%qphase,iatom,dtset%natom,& 522 & rhoijselect=pawrhoij_all(iatom)%rhoijselect,& 523 & test_value=valmx,unit=unitfi,opt_prtvol=dtset%pawprtvol) 524 end do 525 msg=' ' 526 call wrtout(unitfi,msg,'COLL') 527 end do 528 end if 529 530 !PAW+U or local exact-exchange: print out +U components of occupancies 531 !--------------------------------------------------------------------- 532 if ((usepawu.or.useexexch).and.ipositron/=1.and.me_atom==0) then 533 do i_unitfi=1,2 534 unitfi=ab_out;if (i_unitfi==2) unitfi=std_out 535 if(useexexch) write(msg,'(a)') & 536 & ' "Local exact-exchange" part of augmentation waves occupancies Rhoij:' 537 if(usepawu) write(msg,'(a)') & 538 & ' "PAW+U" part of augmentation waves occupancies Rhoij:' 539 call wrtout(unitfi,msg,'COLL') 540 do iatom=1,dtset%natom 541 itypat=pawrhoij_all(iatom)%itypat 542 nspden=pawrhoij_all(iatom)%nspden 543 ll=-1;if (pawtab(itypat)%usepawu/=0) ll=pawtab(itypat)%lpawu 544 llp=-1;if (pawtab(itypat)%useexexch/=0) llp=pawtab(itypat)%lexexch 545 if (ll/=llp.and.ll/=-1.and.llp/=-1) then 546 ABI_BUG("lpawu/=lexexch forbidden!") 547 end if 548 ll=max(ll,llp) 549 if (ll>=0) then 550 call pawrhoij_print_rhoij(pawrhoij_all(iatom)%rhoijp,pawrhoij_all(iatom)%cplex_rhoij,& 551 & pawrhoij_all(iatom)%qphase,iatom,dtset%natom,& 552 & rhoijselect=pawrhoij_all(iatom)%rhoijselect,& 553 & l_only=ll,indlmn=pawtab(itypat)%indlmn,& 554 & unit=unitfi,opt_prtvol=dtset%pawprtvol) 555 end if 556 end do ! iatom 557 msg=' ' 558 call wrtout(unitfi,msg,'COLL') 559 end do 560 end if 561 562 !PAW+U: print out occupations for correlated orbitals 563 !---------------------------------------------------- 564 if (usepawu.and.ipositron/=1.and.me_atom==0) then 565 do i_unitfi=1,2 566 unitfi=ab_out;if (i_unitfi==2) unitfi=std_out 567 write(msg,'(3a)') & 568 & ' ---------- DFT+U DATA --------------------------------------------------- ',ch10 569 call wrtout(unitfi,msg,'COLL') 570 do iatom=1,dtset%natom 571 itypat=dtset%typat(iatom);ll=pawtab(itypat)%lpawu 572 nspden=paw_ij_all(iatom)%nspden;ndij=paw_ij_all(iatom)%ndij 573 cplex_dij=paw_ij_all(iatom)%cplex_dij 574 if ((ll>=0).and.(pawtab(itypat)%usepawu/=0)) then 575 write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom ", iatom,& 576 & ", occupations for correlated orbitals. lpawu =",ll,ch10 577 call wrtout(unitfi,msg,'COLL') 578 if(pawtab(itypat)%usepawu>=10) then 579 write(msg,fmt='(a)') " (This is PAW atomic orbital occupations)" 580 call wrtout(unitfi,msg,'COLL') 581 write(msg,fmt='(a)') " (For Wannier orbital occupations, refer to DFT+DMFT occupations above)" 582 call wrtout(unitfi,msg,'COLL') 583 end if 584 if(nspden==2) then 585 do ispden=1,nspden 586 write(msg,fmt='(a,i4,a,i3,a,f10.5)') " Atom", iatom,& 587 & ". Occ. for lpawu and for spin",ispden," =",paw_ij_all(iatom)%nocctot(ispden) 588 call wrtout(unitfi,msg,'COLL') 589 end do 590 localm=paw_ij_all(iatom)%nocctot(2)-paw_ij_all(iatom)%nocctot(1) 591 write(msg,fmt='(a,i4,a,2x,f12.6)') " => On atom",iatom,& 592 & ", local Mag. for lpawu is ",localm 593 call wrtout(unitfi,msg,'COLL') 594 end if 595 if(ndij==4) then 596 ntot=paw_ij_all(iatom)%nocctot(1) 597 mx=paw_ij_all(iatom)%nocctot(2) 598 my=paw_ij_all(iatom)%nocctot(3) 599 mz=paw_ij_all(iatom)%nocctot(4) 600 mnorm=sqrt(mx*mx+my*my+mz*mz) 601 write(msg,'(a,i4,a,2x,e15.8)') " => On atom",iatom,", for lpawu, local Mag. x is ",mx 602 call wrtout(unitfi,msg,'COLL') 603 write(msg,'(14x,a,2x,e15.8)') " local Mag. y is ",my 604 call wrtout(unitfi,msg,'COLL') 605 write(msg,'(14x,a,2x,e15.8)') " local Mag. z is ",mz 606 call wrtout(unitfi,msg,'COLL') 607 write(msg,'(14x,a,2x,e15.8)') " norm of Mag. is ",mnorm 608 call wrtout(unitfi,msg,'COLL') 609 write(msg,fmt='(8x,a,2x,f10.5)') " (along mag axis) occ. for majority spin is = ",& 610 & half*(ntot+mnorm) 611 call wrtout(unitfi,msg,'COLL') 612 write(msg,fmt='(8x,a,2x,f10.5)') " (along mag axis) occ. for minority spin is = ",& 613 & half*(ntot-mnorm) 614 call wrtout(unitfi,msg,'COLL') 615 end if 616 write(msg,'(3a)') ch10," == Occupation matrix for correlated orbitals:",ch10 617 call wrtout(unitfi,msg,'COLL') 618 do ispden=1,ndij 619 if (nspden==1.and.ndij/=4.and.(cplex_dij==1)) write(msg,fmt='(a)') " Up component only..." 620 if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden 621 if (ndij==4.or.(cplex_dij==2)) & 622 & write(msg,fmt='(2a)') " Occupation matrix for component ",trim(dspin1(ispden+2*(ndij/4))) 623 call wrtout(unitfi,msg,'COLL') 624 do im1=1,ll*2+1 625 if(cplex_dij==1)& 626 & write(msg,'(12(1x,9(1x,f10.5)))') (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1) 627 if(cplex_dij==2)& 628 & write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')& 629 & (paw_ij_all(iatom)%noccmmp(:,im1,im2,ispden),im2=1,ll*2+1) 630 call wrtout(unitfi,msg,'COLL') 631 end do 632 write(msg,'(2a)') ch10,' ' 633 call wrtout(unitfi,msg,'COLL') 634 end do 635 ! Transformation matrices: real->complex spherical harmonics 636 if(paw_ij_all(iatom)%ndij==4) then 637 ABI_MALLOC(noccmmp_ylm,(2*ll+1,2*ll+1,paw_ij_all(iatom)%ndij)) 638 noccmmp_ylm=czero 639 ABI_MALLOC(noccmmp_slm,(2*ll+1,2*ll+1,paw_ij_all(iatom)%ndij)) 640 noccmmp_slm=czero 641 ! Go from real notation for complex noccmmp to complex notation in noccmmp_slm 642 noccmmp_slm(:,:,:)=cmplx(paw_ij_all(iatom)%noccmmp(1,:,:,:),paw_ij_all(iatom)%noccmmp(2,:,:,:)) 643 ii=std_out;if (unitfi==ab_out) ii=-1 644 call mat_slm2ylm(ll,noccmmp_slm,noccmmp_ylm,paw_ij_all(iatom)%ndij,& 645 & 1,1,dtset%pawprtvol,ii,'COLL') ! optspin=1 because up spin are first 646 do ispden=1,paw_ij_all(iatom)%ndij 647 write(msg,'(3a)') ch10,& 648 & "== Occupation matrix in the complex harmonics basis for component ",& 649 & trim(dspin1(ispden+2*(paw_ij_all(iatom)%ndij/4))) 650 call wrtout(unitfi,msg,'COLL') 651 do im1=1,ll*2+1 652 write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') & 653 & (noccmmp_ylm(im1,im2,ispden),im2=1,ll*2+1) 654 call wrtout(unitfi,msg,'COLL') 655 end do 656 end do 657 write(msg,'(a)') ch10 658 call wrtout(unitfi,msg,'COLL') 659 if (dtset%pawspnorb>0) then 660 ABI_MALLOC(noccmmp_jmj,(2*(2*ll+1),2*(2*ll+1))) 661 noccmmp_jmj=czero 662 ii=std_out;if (unitfi==ab_out) ii=-1 663 call mat_mlms2jmj(ll,noccmmp_ylm,noccmmp_jmj,paw_ij_all(iatom)%ndij,& 664 & 1,1,dtset%pawprtvol,-1,'COLL') ! optspin=1: up spin are first 665 write(msg,'(3a)') ch10,"== Occupation matrix in the J (= L-1/2, L+1/2) and M_J basis" 666 call wrtout(unitfi,msg,'COLL') 667 do im1=1,2*(ll*2+1) 668 write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') & 669 & (noccmmp_jmj(im1,im2),im2=1,2*(ll*2+1)) 670 call wrtout(unitfi,msg,'COLL') 671 end do 672 write(msg,'(a)') ch10 673 call wrtout(unitfi,msg,'COLL') 674 ABI_FREE(noccmmp_jmj) 675 end if ! pawspnorb 676 ABI_FREE(noccmmp_ylm) 677 ABI_FREE(noccmmp_slm) 678 end if ! ndij==4 679 end if ! ((ll>=0).and.(pawtab(itypat)%usepawu/=0)) 680 end do 681 end do 682 end if 683 684 !Exact exchange: print out occupations for correlated orbitals 685 !------------------------------------------------------------- 686 if (useexexch.and.ipositron/=1.and.me_atom==0) then 687 nspden=paw_ij_all(1)%nspden;nsppol=paw_ij_all(1)%nsppol;ndij=paw_ij_all(1)%ndij 688 do iatom=1,dtset%natom 689 itypat=dtset%typat(iatom);ll=pawtab(itypat)%lexexch 690 cplex_dij=paw_ij_all(iatom)%cplex_dij 691 if (ll>=0.and.pawtab(itypat)%useexexch/=0) then 692 ABI_MALLOC(paw_ij_all(iatom)%noccmmp,(cplex_dij,2*ll+1,2*ll+1,ndij)) 693 ABI_MALLOC(paw_ij_all(iatom)%nocctot,(nspden)) 694 end if 695 end do 696 call setnoccmmp(1,0,rdum4,0,0,idum3,dtset%natom,dtset%natom,0,1,nsppol,0,dtset%ntypat,& 697 & paw_ij_all,pawang_dum,dtset%pawprtvol,pawrhoij_all,pawtab,rdum2,idum1,dtset%typat,1,0) 698 do i_unitfi=1,2 699 unitfi=ab_out;if (i_unitfi==2) unitfi=std_out 700 write(msg, '(3a)' ) & 701 & ' ---------- Exact Exchange --------------------------------------------------- ',ch10 702 call wrtout(unitfi,msg,'COLL') 703 do iatom=1,dtset%natom 704 itypat=dtset%typat(iatom);ll=pawtab(itypat)%lexexch 705 cplex_dij=paw_ij_all(iatom)%cplex_dij 706 if ((ll>=0).and.(pawtab(itypat)%useexexch/=0)) then 707 write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom",iatom,& 708 & ", occupations for correlated orbitals. l =",ll,ch10 709 call wrtout(unitfi,msg,'COLL') 710 do ispden=1,ndij 711 if (nspden==1.and.ndij/=4) write(msg,fmt='(a)') " Up component only..." 712 if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden 713 if (ndij==4) write(msg,fmt='(2a)') " Occupation matrix for component ",& 714 & trim(dspin2(ispden+2*(ndij/4))) 715 call wrtout(unitfi,msg,'COLL') 716 do im1=1,ll*2+1 717 if(cplex_dij==1)& 718 & write(msg,'(12(1x,9(1x,f10.5)))')& 719 & (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1) 720 if(cplex_dij==2)& 721 & write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') & 722 & (paw_ij_all(iatom)%noccmmp(:,im1,im2,ispden),im2=1,ll*2+1) 723 call wrtout(unitfi,msg,'COLL') 724 end do 725 call wrtout(unitfi,' ','COLL') 726 end do 727 end if 728 end do 729 end do 730 do iatom=1,dtset%natom 731 if (allocated(paw_ij_all(iatom)%noccmmp)) then 732 ABI_FREE(paw_ij_all(iatom)%noccmmp) 733 end if 734 if (allocated(paw_ij_all(iatom)%nocctot)) then 735 ABI_FREE(paw_ij_all(iatom)%nocctot) 736 end if 737 end do 738 end if 739 740 msg=' ' 741 call wrtout(ab_out,msg,'COLL') 742 call wrtout(std_out,msg,'COLL') 743 744 !Destroy temporary stored atomic data 745 ABI_FREE(jatom) 746 call free_my_atmtab(my_atmtab,my_atmtab_allocated) 747 if (paral_atom) then 748 call pawrhoij_free(pawrhoij_all) 749 ABI_FREE(pawrhoij_all) 750 if (usepawu.or.useexexch) then 751 call paw_ij_free(paw_ij_all) 752 ABI_FREE(paw_ij_all) 753 end if 754 end if 755 756 DBG_EXIT("COLL") 757 758 end subroutine pawprt