Line data Source code
1 : !!****m* ABINIT/m_paw_tools
2 : !! NAME
3 : !! m_paw_tools
4 : !!
5 : !! FUNCTION
6 : !! This module contains miscelaneous routines used in the PAW context.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2018-2026 ABINIT group (FJ,MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
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 !========================================================================================
51 : !!***
52 :
53 : !----------------------------------------------------------------------
54 :
55 : !!****f* m_paw_tools/chkpawovlp
56 : !! NAME
57 : !! chkpawovlp
58 : !!
59 : !! FUNCTION
60 : !! Verify that the PAW spheres are not overlapping
61 : !!
62 : !! INPUTS
63 : !! natom=number of atoms in cell.
64 : !! nremit [optional] = if non-zero initialize the number of possible remits before stop
65 : !! ntypat=number of types of atoms in unit cell.
66 : !! pawovlp=percentage of voluminal overlap ratio allowed to continue execution
67 : !! (if negative value, execution always continues)
68 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data:
69 : !! rmet(3,3)=real space metric ($\textrm{bohr}^{2}$).
70 : !! typat(natom)=type (integer) for each atom
71 : !! xred(3,natom)=reduced dimensionless atomic coordinates
72 : !!
73 : !! OUTPUT
74 : !! (only checking)
75 : !!
76 : !! NOTES
77 : !!
78 : !! SOURCE
79 :
80 1659 : 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 1659 : if(present(nremit))then
120 1500 : 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 159 : 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 8295 : iamax(:)=-1;ibmax(:)=-1
133 8295 : norm2_min(:)=-1.d0;ratio_percent_max(:)=-1.d0
134 1659 : iovl(:)=0
135 :
136 : !Loop on "overlapping" atoms with the maximum overlap
137 6147 : do ia=1,natom
138 :
139 4488 : rcuta(1)=pawtab(typat(ia))%rpaw
140 4488 : rcuta(2)=pawtab(typat(ia))%rshp
141 :
142 21657 : do ib=ia,natom
143 :
144 15510 : rcutb(1)=pawtab(typat(ib))%rpaw
145 15510 : rcutb(2)=pawtab(typat(ib))%rshp
146 15510 : r2cut(1)=(rcuta(1)+rcutb(1))**2
147 15510 : r2cut(2)=(rcuta(2)+rcutb(2))**2
148 :
149 : ! Visit the box and its first images:
150 66528 : do t3=-1,1
151 201630 : do t2=-1,1
152 604890 : do t1=-1,1
153 :
154 418770 : dif1=xred(1,ia)-(xred(1,ib)+dble(t1))
155 418770 : dif2=xred(2,ia)-(xred(2,ib)+dble(t2))
156 418770 : dif3=xred(3,ia)-(xred(3,ib)+dble(t3))
157 418770 : norm2=sqnrm_pawovlp(dif1,dif2,dif3)
158 :
159 1395900 : do ii=1,2
160 :
161 1256310 : if(norm2>tol10.and.norm2<r2cut(ii)) then
162 :
163 2114 : iovl(ii)=iovl(ii)+1
164 :
165 : ! Compute the overlap ratio:
166 2114 : dd=sqrt(norm2)
167 2114 : va=4._dp/3._dp*pi*rcuta(ii)**3
168 2114 : vb=4._dp/3._dp*pi*rcutb(ii)**3
169 2114 : ha=(rcutb(ii)**2-(dd-rcuta(ii))**2)/(two*dd)
170 2114 : hb=(rcuta(ii)**2-(dd-rcutb(ii))**2)/(two*dd)
171 2114 : vv=pi/3.d0*(ha**2*(three*rcuta(ii)-ha)+hb**2*(three*rcutb(ii)-hb))
172 2114 : ratio_percent=100._dp*min(vv/min(va,vb),one)
173 2114 : if (ratio_percent>ratio_percent_max(ii)) then
174 988 : ratio_percent_max(ii)=ratio_percent
175 988 : norm2_min(ii)=norm2
176 988 : 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 6636 : stop_on_error=(abs(pawovlp)<=tol6.or.(pawovlp>tol6.and.(maxval(ratio_percent_max(1:2))>pawovlp)))
192 1659 : decrease_nremit=0
193 :
194 : !Print adapted message with overlap value
195 1659 : if (iovl(1)+iovl(2)>0) then
196 :
197 1443 : do ii=1,2
198 :
199 1443 : if(ratio_percent_max(ii)>zero)then
200 779 : if (ii==1) write(message,' (a)' ) 'PAW SPHERES ARE OVERLAPPING!'
201 779 : if (ii==2) write(message, '(a)' ) 'PAW COMPENSATION DENSITIES ARE OVERLAPPING !'
202 :
203 779 : if (iovl(ii)==1) then
204 418 : write(message, '(3a)' ) trim(message),ch10,&
205 836 : & ' There is one pair of overlapping atoms.'
206 : else
207 361 : write(message, '(3a,i5,a)' ) trim(message),ch10,&
208 722 : & ' There are ', iovl(ii),' pairs of overlapping atoms.'
209 : end if
210 779 : write(message, '(3a,i4,a,i4,a)' ) trim(message),ch10,&
211 1558 : ' The maximum overlap percentage is obtained for the atoms ',iamax(ii),' and ',ibmax(ii),'.'
212 779 : write(message, '(2a,2(a,i3),a,f9.5)' ) trim(message),ch10,&
213 1558 : & ' | Distance between atoms ',iamax(ii),' and ',ibmax(ii),' is : ',sqrt(norm2_min(ii))
214 779 : if(ii==1)then
215 481 : write(message, '(2a,2(a,i3,a,f9.5,a))' ) trim(message),ch10,&
216 481 : & ' | PAW radius of the sphere around atom ',iamax(ii),' is: ',pawtab(typat(iamax(ii)))%rpaw,ch10,&
217 962 : & ' | PAW radius of the sphere around atom ',ibmax(ii),' is: ',pawtab(typat(ibmax(ii)))%rpaw,ch10
218 : else if(ii==2)then
219 298 : write(message, '(2a,2(a,i3,a,f9.5,a))' ) trim(message),ch10,&
220 298 : & ' | Radius of the compensation sphere around atom ',iamax(ii),' is: ',pawtab(typat(iamax(ii)))%rshp,ch10,&
221 596 : & ' | Radius of the compensation sphere around atom ',ibmax(ii),' is: ',pawtab(typat(ibmax(ii)))%rshp,ch10
222 : endif
223 779 : write(message, '(2a,f7.4,a)' ) trim(message),&
224 1558 : & ' | This leads to a (voluminal) overlap ratio of ',ratio_percent_max(ii),' %'
225 779 : if (ii==1) then
226 481 : write(message, '(3a)' ) trim(message),ch10,&
227 962 : & 'THIS IS DANGEROUS, as PAW formalism assumes non-overlapping PAW spheres.'
228 : else if (ii==2) then
229 298 : write(message, '(3a)' ) trim(message),ch10,&
230 596 : & 'THIS IS DANGEROUS, as PAW formalism assumes non-overlapping compensation densities.'
231 : end if
232 779 : if (stop_on_error .and. nremit_counter==0) then
233 0 : ABI_ERROR_NOSTOP(message,ia) !ia is dummy
234 : else
235 779 : ABI_WARNING(message)
236 779 : 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 481 : if (stop_on_error) then
244 : write(message, '(3a)' )&
245 0 : & ' Action: 1- decrease cutoff radius of PAW dataset',ch10,&
246 0 : & ' OR 2- ajust "pawovlp" input variable to allow overlap (risky)'
247 0 : if(nremit_counter==0)then
248 0 : ABI_ERROR(message)
249 : endif
250 : end if
251 :
252 481 : nremit_counter=nremit_counter-decrease_nremit
253 :
254 : ! Print last message if execution continues:
255 481 : if (pawovlp<=tol6) then
256 : write(message, '(6a)' ) &
257 4 : & ' Results might be approximate,',ch10,&
258 4 : & ' and even inaccurate (if overlap is too big) !',ch10,&
259 8 : & ' Assume experienced user. Execution will continue.',ch10
260 4 : call wrtout(std_out,message,'COLL')
261 477 : else if (ratio_percent_max(1)<=pawovlp .and. ratio_percent_max(2)<=pawovlp) then
262 : write(message, '(8a)' ) &
263 477 : & ' Overlap ratio seems to be acceptable (less than value',ch10,&
264 477 : & ' of "pawovlp" input parameter): execution will continue.',ch10,&
265 477 : & ' But be aware that results might be approximate,',ch10,&
266 954 : & ' and even inaccurate (depending on your physical system) !',ch10
267 477 : call wrtout(std_out,message,'COLL')
268 0 : else if(decrease_nremit==1)then
269 : write(message, '(3a)' ) &
270 0 : & ' First time that overlap is bigger than "pawovlp" input parameter.',ch10,&
271 0 : & ' Execution will continue, but such overlap will not be tolerated twice.'
272 0 : 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 418770 : 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 418770 : & +rmet(1,3)*u1*u3+rmet(2,3)*u2*u3+rmet(3,3)*u3*u3
289 :
290 418770 : end function sqnrm_pawovlp
291 :
292 : end subroutine chkpawovlp
293 : !!***
294 :
295 : !----------------------------------------------------------------------
296 :
297 : !!****f* m_paw_tools/pawprt
298 : !! NAME
299 : !! pawprt
300 : !!
301 : !! FUNCTION
302 : !! Print out data concerning PAW formalism
303 : !! (pseudopotential strength, augmentation occupancies...)
304 : !! To be called at the end of the SCF cycle
305 : !!
306 : !! COPYRIGHT
307 : !! Copyright (C) 1998-2026 ABINIT group (FJ,MT,BA)
308 : !! This file is distributed under the terms of the
309 : !! GNU General Public License, see ~abinit/COPYING
310 : !! or http://www.gnu.org/copyleft/gpl.txt .
311 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt.
312 : !!
313 : !! INPUTS
314 : !! dtset <type(dataset_type)>=all input variables in this dataset
315 : !! | enunit=parameter determining units of output energies
316 : !! | kptopt=option for the generation of k points
317 : !! | natom=number of atoms in cell
318 : !! | ntypat = number of atom types
319 : !! | pawprtvol= printing volume
320 : !! | pawspnorb=flag: 1 if spin-orbit coupling is activated
321 : !! | typat(natom)=type of each atom
322 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation (optional argument)
323 : !! mpi_atmtab(:)=--optional-- indexes of the atoms treated by current proc
324 : !! comm_atom=--optional-- MPI communicator over atoms
325 : !! my_natom=number of atoms treated by current processor
326 : !! paw_ij(my_natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels
327 : !! pawrhoij(my_natom) <type(pawrhoij_type)>= paw rhoij occupancies and related data
328 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
329 : !!
330 : !! OUTPUT
331 : !! (only printing)
332 : !!
333 : !! SOURCE
334 :
335 3024 : subroutine pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,&
336 : & electronpositron,& ! optional argument
337 1512 : & 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 1512 : 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),int_arr(1)
365 : integer :: idum1(0),idum3(0,0,0)
366 1512 : integer,allocatable :: jatom(:)
367 1512 : integer,pointer :: my_atmtab(:)
368 : real(dp) :: rdum2(0,0),rdum4(0,0,0,0)
369 1512 : complex(dp),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:)
370 1512 : type(paw_ij_type), ABI_CONTIGUOUS pointer :: paw_ij_all(:)
371 1512 : type(pawrhoij_type),ABI_CONTIGUOUS pointer :: pawrhoij_all(:)
372 : ! *********************************************************************
373 :
374 : DBG_ENTER("COLL")
375 :
376 : !Set up parallelism over atoms
377 1500 : paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom))
378 1512 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
379 1512 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
380 1512 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom)
381 1512 : me_atom=xmpi_comm_rank(my_comm_atom)
382 :
383 : !Continue only if comm_atom contains the master of the output comm
384 1512 : if (paral_atom) then
385 558 : call xmpi_comm_group(abinit_comm_output,group1,ierr)
386 558 : call xmpi_comm_group(my_comm_atom,group2,ierr)
387 558 : int_arr(1) = 0
388 558 : call xmpi_group_translate_ranks(group1,1,int_arr,group2,idum,ierr)
389 558 : call xmpi_group_free(group1)
390 558 : call xmpi_group_free(group2)
391 558 : if (idum(1)==xmpi_undefined) then
392 116 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
393 116 : return
394 : end if
395 : end if
396 :
397 : !Initializations
398 1396 : natprt=natmax;if (dtset%natom==1) natprt=1
399 1396 : if (dtset%pawprtvol<0) natprt=dtset%natom
400 4188 : ABI_MALLOC(jatom,(natprt))
401 1396 : if (natprt==1) then
402 192 : jatom(1)=1
403 1204 : else if (natprt==2) then
404 1192 : jatom(1)=1;jatom(2)=dtset%natom
405 12 : else if (natprt==dtset%natom) then
406 60 : do iat=1,dtset%natom
407 60 : jatom(iat)=iat
408 : end do
409 : else
410 0 : ABI_BUG("invalid value of natprt!")
411 : end if
412 3393 : usepawu=(count(pawtab(:)%usepawu/=0)>0)
413 3393 : useexexch=(count(pawtab(:)%useexexch/=0)>0)
414 1396 : ipositron=0
415 1396 : if (present(electronpositron)) then
416 1384 : if (associated(electronpositron)) ipositron=electronpositron%calctype
417 : end if
418 :
419 : !Main title
420 1396 : write(msg, '(2a)' ) ch10,&
421 2792 : & ' ==== Results concerning PAW augmentation regions ===='
422 1396 : call wrtout(ab_out,msg,'COLL')
423 1396 : call wrtout(std_out,msg,'COLL')
424 1396 : msg=' '
425 1396 : call wrtout(ab_out,msg,'COLL')
426 1396 : call wrtout(std_out,msg,'COLL')
427 :
428 : !If atomic data are distributed, retrieve all Dij on master proc
429 1396 : if (paral_atom) then
430 442 : if (me_atom==0) then
431 797 : ABI_MALLOC(paw_ij_all,(dtset%natom))
432 139 : call paw_ij_nullify(paw_ij_all)
433 : else
434 303 : ABI_MALLOC(paw_ij_all,(0))
435 : end if
436 442 : call paw_ij_gather(paw_ij,paw_ij_all,0,my_comm_atom)
437 : else
438 954 : paw_ij_all => paw_ij
439 : end if
440 :
441 : !Print out pseudopotential strength
442 : !----------------------------------
443 1396 : if (me_atom==0) then
444 3279 : do i_unitfi=1,2
445 2186 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
446 7651 : do unt=1,2
447 4372 : 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 2414 : if ((unt==1).and.(dtset%enunit==0.or.dtset%enunit==2)) then
450 2172 : write(msg,'(a)') ' Total pseudopotential strength Dij (hartree):'
451 242 : else if ((unt==2).and.(dtset%enunit==1.or.dtset%enunit==2)) then
452 242 : write(msg,'(a)') ' Total pseudopotential strength Dij (eV):'
453 : end if
454 2414 : call wrtout(unitfi,msg,'COLL')
455 2414 : if (ipositron>0) then
456 32 : if (electronpositron%has_pos_ham==0) then
457 12 : write(msg,'(a)') ' -Note: these are the electronic Dij'
458 : else
459 20 : write(msg,'(a)') ' -Note: these are the positronic Dij'
460 : end if
461 32 : call wrtout(unitfi,msg,'COLL')
462 : end if
463 2414 : valmx=100._dp;if (ipositron>0) valmx=-1._dp
464 6822 : do iat=1,natprt
465 4408 : 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 8780 : & test_value=valmx,unit=unitfi,Ha_or_eV=unt,opt_prtvol=dtset%pawprtvol)
469 : end do
470 : end if
471 4372 : msg=' '
472 6558 : call wrtout(unitfi,msg,'COLL')
473 : end do
474 : end do
475 : end if
476 1396 : if (paral_atom.and.(.not.usepawu).and.(.not.useexexch)) then
477 434 : call paw_ij_free(paw_ij_all)
478 434 : ABI_FREE(paw_ij_all)
479 : end if
480 :
481 : !If atomic data are distributed, retrieve all Rhoij on master proc
482 1396 : if (paral_atom) then
483 442 : if (me_atom==0) then
484 797 : ABI_MALLOC(pawrhoij_all,(dtset%natom))
485 : else
486 303 : ABI_MALLOC(pawrhoij_all,(0))
487 : end if
488 442 : call pawrhoij_nullify(pawrhoij_all)
489 : call pawrhoij_gather(pawrhoij,pawrhoij_all,0,my_comm_atom,&
490 : & with_grhoij=.false.,with_lmnmix=.false.,&
491 442 : & with_rhoij_=.false.,with_rhoijres=.false.)
492 : else
493 954 : pawrhoij_all => pawrhoij
494 : end if
495 :
496 : !Print out SYMMETRIZED occupancies of the partial waves
497 : !------------------------------------------------------
498 1396 : if (me_atom==0) then
499 3279 : do i_unitfi=1,2
500 2186 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
501 2186 : write(msg,'(a)') ' Augmentation waves occupancies Rhoij:'
502 2186 : call wrtout(unitfi,msg,'COLL')
503 2186 : if (ipositron>0) then
504 32 : if (electronpositron%particle==EP_POSITRON) then
505 24 : write(msg,'(a)') ' -Note: these are the electronic Rhoij'
506 : else
507 8 : write(msg,'(a)') ' -Note: these are the positronic Rhoij'
508 : end if
509 32 : call wrtout(unitfi,msg,'COLL')
510 : end if
511 2186 : if (dtset%pawspnorb>0.and.pawrhoij_all(1)%cplex_rhoij==1.and.dtset%kptopt/=1.and.dtset%kptopt/=2) then
512 2 : write(msg,'(6a)') ' pawprt: - WARNING:',ch10,&
513 2 : & ' Spin-orbit coupling is activated but only real part of Rhoij occupancies',ch10,&
514 4 : & ' has been computed; they could have an imaginary part (not printed here).'
515 2 : call wrtout(unitfi,msg,'COLL')
516 : end if
517 2186 : valmx=25._dp;if (ipositron>0) valmx=-1._dp
518 6210 : do iat=1,natprt
519 4024 : 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 6210 : & test_value=valmx,unit=unitfi,opt_prtvol=dtset%pawprtvol)
524 : end do
525 2186 : msg=' '
526 3279 : 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 1396 : if ((usepawu.or.useexexch).and.ipositron/=1.and.me_atom==0) then
533 462 : do i_unitfi=1,2
534 308 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
535 308 : if(useexexch) write(msg,'(a)') &
536 2 : & ' "Local exact-exchange" part of augmentation waves occupancies Rhoij:'
537 308 : if(usepawu) write(msg,'(a)') &
538 306 : & ' "PAW+U" part of augmentation waves occupancies Rhoij:'
539 308 : call wrtout(unitfi,msg,'COLL')
540 992 : do iatom=1,dtset%natom
541 684 : itypat=pawrhoij_all(iatom)%itypat
542 684 : nspden=pawrhoij_all(iatom)%nspden
543 684 : ll=-1;if (pawtab(itypat)%usepawu/=0) ll=pawtab(itypat)%lpawu
544 684 : llp=-1;if (pawtab(itypat)%useexexch/=0) llp=pawtab(itypat)%lexexch
545 684 : if (ll/=llp.and.ll/=-1.and.llp/=-1) then
546 0 : ABI_BUG("lpawu/=lexexch forbidden!")
547 : end if
548 684 : ll=max(ll,llp)
549 992 : 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 430 : & unit=unitfi,opt_prtvol=dtset%pawprtvol)
555 : end if
556 : end do ! iatom
557 308 : msg=' '
558 462 : call wrtout(unitfi,msg,'COLL')
559 : end do
560 : end if
561 :
562 : !PAW+U: print out occupations for correlated orbitals
563 : !----------------------------------------------------
564 1396 : if (usepawu.and.ipositron/=1.and.me_atom==0) then
565 459 : do i_unitfi=1,2
566 306 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
567 : write(msg,'(3a)') &
568 306 : & ' ---------- DFT+U DATA --------------------------------------------------- ',ch10
569 306 : call wrtout(unitfi,msg,'COLL')
570 1135 : do iatom=1,dtset%natom
571 676 : itypat=dtset%typat(iatom);ll=pawtab(itypat)%lpawu
572 676 : nspden=paw_ij_all(iatom)%nspden;ndij=paw_ij_all(iatom)%ndij
573 676 : cplex_dij=paw_ij_all(iatom)%cplex_dij
574 982 : if ((ll>=0).and.(pawtab(itypat)%usepawu/=0)) then
575 426 : write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom ", iatom,&
576 852 : & ", occupations for correlated orbitals. lpawu =",ll,ch10
577 426 : call wrtout(unitfi,msg,'COLL')
578 426 : if(pawtab(itypat)%usepawu>=10) then
579 112 : write(msg,fmt='(a)') " (This is PAW atomic orbital occupations)"
580 112 : call wrtout(unitfi,msg,'COLL')
581 112 : write(msg,fmt='(a)') " (For Wannier orbital occupations, refer to DFT+DMFT occupations above)"
582 112 : call wrtout(unitfi,msg,'COLL')
583 : end if
584 426 : if(nspden==2) then
585 888 : do ispden=1,nspden
586 592 : write(msg,fmt='(a,i4,a,i3,a,f10.5)') " Atom", iatom,&
587 1184 : & ". Occ. for lpawu and for spin",ispden," =",paw_ij_all(iatom)%nocctot(ispden)
588 888 : call wrtout(unitfi,msg,'COLL')
589 : end do
590 296 : localm=paw_ij_all(iatom)%nocctot(2)-paw_ij_all(iatom)%nocctot(1)
591 296 : write(msg,fmt='(a,i4,a,2x,f12.6)') " => On atom",iatom,&
592 592 : & ", local Mag. for lpawu is ",localm
593 296 : call wrtout(unitfi,msg,'COLL')
594 : end if
595 426 : if(ndij==4) then
596 70 : ntot=paw_ij_all(iatom)%nocctot(1)
597 70 : mx=paw_ij_all(iatom)%nocctot(2)
598 70 : my=paw_ij_all(iatom)%nocctot(3)
599 70 : mz=paw_ij_all(iatom)%nocctot(4)
600 70 : mnorm=sqrt(mx*mx+my*my+mz*mz)
601 70 : write(msg,'(a,i4,a,2x,e15.8)') " => On atom",iatom,", for lpawu, local Mag. x is ",mx
602 70 : call wrtout(unitfi,msg,'COLL')
603 70 : write(msg,'(14x,a,2x,e15.8)') " local Mag. y is ",my
604 70 : call wrtout(unitfi,msg,'COLL')
605 70 : write(msg,'(14x,a,2x,e15.8)') " local Mag. z is ",mz
606 70 : call wrtout(unitfi,msg,'COLL')
607 70 : write(msg,'(14x,a,2x,e15.8)') " norm of Mag. is ",mnorm
608 70 : call wrtout(unitfi,msg,'COLL')
609 70 : write(msg,fmt='(8x,a,2x,f10.5)') " (along mag axis) occ. for majority spin is = ",&
610 140 : & half*(ntot+mnorm)
611 70 : call wrtout(unitfi,msg,'COLL')
612 70 : write(msg,fmt='(8x,a,2x,f10.5)') " (along mag axis) occ. for minority spin is = ",&
613 140 : & half*(ntot-mnorm)
614 70 : call wrtout(unitfi,msg,'COLL')
615 : end if
616 426 : write(msg,'(3a)') ch10," == Occupation matrix for correlated orbitals:",ch10
617 426 : call wrtout(unitfi,msg,'COLL')
618 1358 : do ispden=1,ndij
619 932 : if (nspden==1.and.ndij/=4.and.(cplex_dij==1)) write(msg,fmt='(a)') " Up component only..."
620 932 : if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden
621 932 : if (ndij==4.or.(cplex_dij==2)) &
622 280 : & write(msg,fmt='(2a)') " Occupation matrix for component ",trim(dspin1(ispden+2*(ndij/4)))
623 932 : call wrtout(unitfi,msg,'COLL')
624 5636 : do im1=1,ll*2+1
625 4704 : if(cplex_dij==1)&
626 20100 : & write(msg,'(12(1x,9(1x,f10.5)))') (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1)
627 4704 : if(cplex_dij==2)&
628 : & write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))')&
629 23584 : & (paw_ij_all(iatom)%noccmmp(:,im1,im2,ispden),im2=1,ll*2+1)
630 5636 : call wrtout(unitfi,msg,'COLL')
631 : end do
632 932 : write(msg,'(2a)') ch10,' '
633 1358 : call wrtout(unitfi,msg,'COLL')
634 : end do
635 : ! Transformation matrices: real->complex spherical harmonics
636 426 : if(paw_ij_all(iatom)%ndij==4) then
637 350 : ABI_MALLOC(noccmmp_ylm,(2*ll+1,2*ll+1,paw_ij_all(iatom)%ndij))
638 9166 : noccmmp_ylm=czero
639 280 : ABI_MALLOC(noccmmp_slm,(2*ll+1,2*ll+1,paw_ij_all(iatom)%ndij))
640 9166 : noccmmp_slm=czero
641 : ! Go from real notation for complex noccmmp to complex notation in noccmmp_slm
642 9166 : noccmmp_slm(:,:,:)=cmplx(paw_ij_all(iatom)%noccmmp(1,:,:,:),paw_ij_all(iatom)%noccmmp(2,:,:,:))
643 70 : ii=std_out;if (unitfi==ab_out) ii=-1
644 : call mat_slm2ylm(ll,noccmmp_slm,noccmmp_ylm,paw_ij_all(iatom)%ndij,&
645 70 : & 1,1,dtset%pawprtvol,ii,'COLL') ! optspin=1 because up spin are first
646 350 : do ispden=1,paw_ij_all(iatom)%ndij
647 280 : write(msg,'(3a)') ch10,&
648 280 : & "== Occupation matrix in the complex harmonics basis for component ",&
649 560 : & trim(dspin1(ispden+2*(paw_ij_all(iatom)%ndij/4)))
650 280 : call wrtout(unitfi,msg,'COLL')
651 1782 : do im1=1,ll*2+1
652 : write(msg,'(12(1x,9(1x,"(",f9.5,",",f9.5,")")))') &
653 1432 : & (noccmmp_ylm(im1,im2,ispden),im2=1,ll*2+1)
654 1712 : call wrtout(unitfi,msg,'COLL')
655 : end do
656 : end do
657 70 : write(msg,'(a)') ch10
658 70 : call wrtout(unitfi,msg,'COLL')
659 70 : if (dtset%pawspnorb>0) then
660 248 : ABI_MALLOC(noccmmp_jmj,(2*(2*ll+1),2*(2*ll+1)))
661 7282 : noccmmp_jmj=czero
662 62 : ii=std_out;if (unitfi==ab_out) ii=-1
663 : call mat_mlms2jmj(ll,noccmmp_ylm,noccmmp_jmj,paw_ij_all(iatom)%ndij,&
664 62 : & 1,1,dtset%pawprtvol,-1,'COLL') ! optspin=1: up spin are first
665 62 : write(msg,'(3a)') ch10,"== Occupation matrix in the J (= L-1/2, L+1/2) and M_J basis"
666 62 : call wrtout(unitfi,msg,'COLL')
667 698 : do im1=1,2*(ll*2+1)
668 : write(msg,'(12(1x,18(1x,"(",f7.3,",",f7.3,")")))') &
669 636 : & (noccmmp_jmj(im1,im2),im2=1,2*(ll*2+1))
670 698 : call wrtout(unitfi,msg,'COLL')
671 : end do
672 62 : write(msg,'(a)') ch10
673 62 : call wrtout(unitfi,msg,'COLL')
674 62 : ABI_FREE(noccmmp_jmj)
675 : end if ! pawspnorb
676 70 : ABI_FREE(noccmmp_ylm)
677 70 : 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 1396 : if (useexexch.and.ipositron/=1.and.me_atom==0) then
687 1 : nspden=paw_ij_all(1)%nspden;nsppol=paw_ij_all(1)%nsppol;ndij=paw_ij_all(1)%ndij
688 5 : do iatom=1,dtset%natom
689 4 : itypat=dtset%typat(iatom);ll=pawtab(itypat)%lexexch
690 4 : cplex_dij=paw_ij_all(iatom)%cplex_dij
691 5 : if (ll>=0.and.pawtab(itypat)%useexexch/=0) then
692 12 : ABI_MALLOC(paw_ij_all(iatom)%noccmmp,(cplex_dij,2*ll+1,2*ll+1,ndij))
693 6 : 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 1 : & paw_ij_all,pawang_dum,dtset%pawprtvol,pawrhoij_all,pawtab,rdum2,idum1,dtset%typat,1,0)
698 3 : do i_unitfi=1,2
699 2 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
700 : write(msg, '(3a)' ) &
701 2 : & ' ---------- Exact Exchange --------------------------------------------------- ',ch10
702 2 : call wrtout(unitfi,msg,'COLL')
703 11 : do iatom=1,dtset%natom
704 8 : itypat=dtset%typat(iatom);ll=pawtab(itypat)%lexexch
705 8 : cplex_dij=paw_ij_all(iatom)%cplex_dij
706 10 : if ((ll>=0).and.(pawtab(itypat)%useexexch/=0)) then
707 4 : write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom",iatom,&
708 8 : & ", occupations for correlated orbitals. l =",ll,ch10
709 4 : call wrtout(unitfi,msg,'COLL')
710 12 : do ispden=1,ndij
711 8 : if (nspden==1.and.ndij/=4) write(msg,fmt='(a)') " Up component only..."
712 8 : if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden
713 8 : if (ndij==4) write(msg,fmt='(2a)') " Occupation matrix for component ",&
714 0 : & trim(dspin2(ispden+2*(ndij/4)))
715 8 : call wrtout(unitfi,msg,'COLL')
716 48 : do im1=1,ll*2+1
717 40 : if(cplex_dij==1)&
718 : & write(msg,'(12(1x,9(1x,f10.5)))')&
719 240 : & (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1)
720 40 : if(cplex_dij==2)&
721 : & write(msg,'(12(1x,9(1x,"(",f7.3,",",f7.3,")")))') &
722 0 : & (paw_ij_all(iatom)%noccmmp(:,im1,im2,ispden),im2=1,ll*2+1)
723 48 : call wrtout(unitfi,msg,'COLL')
724 : end do
725 12 : call wrtout(unitfi,' ','COLL')
726 : end do
727 : end if
728 : end do
729 : end do
730 5 : do iatom=1,dtset%natom
731 4 : if (allocated(paw_ij_all(iatom)%noccmmp)) then
732 2 : ABI_FREE(paw_ij_all(iatom)%noccmmp)
733 : end if
734 5 : if (allocated(paw_ij_all(iatom)%nocctot)) then
735 2 : ABI_FREE(paw_ij_all(iatom)%nocctot)
736 : end if
737 : end do
738 : end if
739 :
740 1396 : msg=' '
741 1396 : call wrtout(ab_out,msg,'COLL')
742 1396 : call wrtout(std_out,msg,'COLL')
743 :
744 : !Destroy temporary stored atomic data
745 1396 : ABI_FREE(jatom)
746 1396 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
747 1396 : if (paral_atom) then
748 442 : call pawrhoij_free(pawrhoij_all)
749 442 : ABI_FREE(pawrhoij_all)
750 442 : if (usepawu.or.useexexch) then
751 8 : call paw_ij_free(paw_ij_all)
752 8 : ABI_FREE(paw_ij_all)
753 : end if
754 : end if
755 :
756 : DBG_EXIT("COLL")
757 :
758 3024 : end subroutine pawprt
759 : !!***
760 :
761 : !----------------------------------------------------------------------
762 :
763 : END MODULE m_paw_tools
764 : !!***
|
