Branch data 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 : 1692 : 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 [ + + ]: 1692 : if(present(nremit))then
120 [ + + ]: 1525 : 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 : 167 : 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 [ + + + + ]: 8460 : iamax(:)=-1;ibmax(:)=-1
133 [ + + + + ]: 8460 : norm2_min(:)=-1.d0;ratio_percent_max(:)=-1.d0
134 : 1692 : iovl(:)=0
135 : :
136 : : !Loop on "overlapping" atoms with the maximum overlap
137 [ + + ]: 6242 : do ia=1,natom
138 : :
139 : 4550 : rcuta(1)=pawtab(typat(ia))%rpaw
140 : 4550 : rcuta(2)=pawtab(typat(ia))%rshp
141 : :
142 [ + + ]: 21849 : do ib=ia,natom
143 : :
144 : 15607 : rcutb(1)=pawtab(typat(ib))%rpaw
145 : 15607 : rcutb(2)=pawtab(typat(ib))%rshp
146 : 15607 : r2cut(1)=(rcuta(1)+rcutb(1))**2
147 : 15607 : r2cut(2)=(rcuta(2)+rcutb(2))**2
148 : :
149 : : ! Visit the box and its first images:
150 [ + + ]: 66978 : do t3=-1,1
151 [ + + ]: 202891 : do t2=-1,1
152 [ + + ]: 608673 : do t1=-1,1
153 : :
154 : 421389 : dif1=xred(1,ia)-(xred(1,ib)+dble(t1))
155 : 421389 : dif2=xred(2,ia)-(xred(2,ib)+dble(t2))
156 : 421389 : dif3=xred(3,ia)-(xred(3,ib)+dble(t3))
157 : 421389 : norm2=sqnrm_pawovlp(dif1,dif2,dif3)
158 : :
159 [ + + ]: 1404630 : do ii=1,2
160 : :
161 [ + + + + ]: 1264167 : 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 [ - + + - : 6768 : stop_on_error=(abs(pawovlp)<=tol6.or.(pawovlp>tol6.and.(maxval(ratio_percent_max(1:2))>pawovlp)))
+ + + - +
+ - + ]
192 : 1692 : decrease_nremit=0
193 : :
194 : : !Print adapted message with overlap value
195 [ + + ]: 1692 : 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 : 421389 : 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 : 421389 : & +rmet(1,3)*u1*u3+rmet(2,3)*u2*u3+rmet(3,3)*u3*u3
289 : :
290 : 421389 : 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 : 3080 : subroutine pawprt(dtset,my_natom,paw_ij,pawrhoij,pawtab,&
336 : : & electronpositron,& ! optional argument
337 [ + + + + : 1540 : & 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 : 1540 : 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 : 1540 : integer,allocatable :: jatom(:)
367 : 1540 : integer,pointer :: my_atmtab(:)
368 : : real(dp) :: rdum2(0,0),rdum4(0,0,0,0)
369 : 1540 : complex(dp),allocatable :: noccmmp_ylm(:,:,:),noccmmp_jmj(:,:),noccmmp_slm(:,:,:)
370 : 1540 : type(paw_ij_type), contiguous, pointer :: paw_ij_all(:)
371 : 1540 : type(pawrhoij_type),contiguous, pointer :: pawrhoij_all(:)
372 : : ! *********************************************************************
373 : :
374 : : DBG_ENTER("COLL")
375 : :
376 : : !Set up parallelism over atoms
377 [ + + ]: 1525 : paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom))
378 [ + + ]: 1540 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
379 [ + + ]: 1540 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
380 : 1540 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom)
381 : 1540 : me_atom=xmpi_comm_rank(my_comm_atom)
382 : :
383 : : !Continue only if comm_atom contains the master of the output comm
384 [ + + ]: 1540 : 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 [ + + ]: 1424 : natprt=natmax;if (dtset%natom==1) natprt=1
399 [ + + ]: 1424 : if (dtset%pawprtvol<0) natprt=dtset%natom
400 [ - + + - : 4272 : ABI_MALLOC(jatom,(natprt))
+ - + - -
+ - + ]
401 [ + + ]: 1424 : if (natprt==1) then
402 : 196 : jatom(1)=1
403 [ + + ]: 1228 : else if (natprt==2) then
404 : 1216 : 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 [ + + + + ]: 3473 : usepawu=(count(pawtab(:)%usepawu/=0)>0)
413 [ + + + + ]: 3473 : useexexch=(count(pawtab(:)%useexexch/=0)>0)
414 : 1424 : ipositron=0
415 [ + + ]: 1424 : if (present(electronpositron)) then
416 [ + + ]: 1409 : if (associated(electronpositron)) ipositron=electronpositron%calctype
417 : : end if
418 : :
419 : : !Main title
420 : 1424 : write(msg, '(2a)' ) ch10,&
421 : 2848 : & ' ==== Results concerning PAW augmentation regions ===='
422 : 1424 : call wrtout(ab_out,msg,'COLL')
423 : 1424 : call wrtout(std_out,msg,'COLL')
424 : 1424 : msg=' '
425 : 1424 : call wrtout(ab_out,msg,'COLL')
426 : 1424 : call wrtout(std_out,msg,'COLL')
427 : :
428 : : !If atomic data are distributed, retrieve all Dij on master proc
429 [ + + ]: 1424 : 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 : 982 : paw_ij_all => paw_ij
439 : : end if
440 : :
441 : : !Print out pseudopotential strength
442 : : !----------------------------------
443 [ + + ]: 1424 : if (me_atom==0) then
444 [ + + ]: 3363 : do i_unitfi=1,2
445 [ + + ]: 2242 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
446 [ + + ]: 7847 : do unt=1,2
447 [ + + + + : 4484 : 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 [ + + + - ]: 2470 : if ((unt==1).and.(dtset%enunit==0.or.dtset%enunit==2)) then
450 : 2228 : 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 : 2470 : call wrtout(unitfi,msg,'COLL')
455 [ + + ]: 2470 : 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 : 2470 : valmx=100._dp;if (ipositron>0) valmx=-1._dp
464 [ + + ]: 6982 : do iat=1,natprt
465 : 4512 : 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 : 8996 : & test_value=valmx,unit=unitfi,Ha_or_eV=unt,opt_prtvol=dtset%pawprtvol)
469 : : end do
470 : : end if
471 : 4484 : msg=' '
472 : 6726 : call wrtout(unitfi,msg,'COLL')
473 : : end do
474 : : end do
475 : : end if
476 [ + + + + ]: 1424 : 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 [ + + ]: 1424 : 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 : 982 : pawrhoij_all => pawrhoij
494 : : end if
495 : :
496 : : !Print out SYMMETRIZED occupancies of the partial waves
497 : : !------------------------------------------------------
498 [ + + ]: 1424 : if (me_atom==0) then
499 [ + + ]: 3363 : do i_unitfi=1,2
500 [ + + ]: 2242 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
501 : 2242 : write(msg,'(a)') ' Augmentation waves occupancies Rhoij:'
502 : 2242 : call wrtout(unitfi,msg,'COLL')
503 [ + + ]: 2242 : 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 [ + + + + : 2242 : 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 [ + + ]: 2242 : valmx=25._dp;if (ipositron>0) valmx=-1._dp
518 [ + + ]: 6370 : do iat=1,natprt
519 : 4128 : 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 : 6370 : & test_value=valmx,unit=unitfi,opt_prtvol=dtset%pawprtvol)
524 : : end do
525 : 2242 : msg=' '
526 : 3363 : 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 [ + + + + ]: 1424 : if ((usepawu.or.useexexch).and.ipositron/=1.and.me_atom==0) then
533 [ + + ]: 465 : do i_unitfi=1,2
534 [ + + ]: 310 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
535 [ + + ]: 310 : if(useexexch) write(msg,'(a)') &
536 : 2 : & ' "Local exact-exchange" part of augmentation waves occupancies Rhoij:'
537 [ + + ]: 310 : if(usepawu) write(msg,'(a)') &
538 : 308 : & ' "PAW+U" part of augmentation waves occupancies Rhoij:'
539 : 310 : call wrtout(unitfi,msg,'COLL')
540 [ + + ]: 1002 : do iatom=1,dtset%natom
541 : 692 : itypat=pawrhoij_all(iatom)%itypat
542 : 692 : nspden=pawrhoij_all(iatom)%nspden
543 [ + + ]: 692 : ll=-1;if (pawtab(itypat)%usepawu/=0) ll=pawtab(itypat)%lpawu
544 [ + + ]: 692 : llp=-1;if (pawtab(itypat)%useexexch/=0) llp=pawtab(itypat)%lexexch
545 [ + + - + ]: 692 : if (ll/=llp.and.ll/=-1.and.llp/=-1) then
546 : 0 : ABI_BUG("lpawu/=lexexch forbidden!")
547 : : end if
548 : 692 : ll=max(ll,llp)
549 [ + + ]: 1002 : 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 : 434 : & unit=unitfi,opt_prtvol=dtset%pawprtvol)
555 : : end if
556 : : end do ! iatom
557 : 310 : msg=' '
558 : 465 : call wrtout(unitfi,msg,'COLL')
559 : : end do
560 : : end if
561 : :
562 : : !PAW+U: print out occupations for correlated orbitals
563 : : !----------------------------------------------------
564 [ + + + + ]: 1424 : if (usepawu.and.ipositron/=1.and.me_atom==0) then
565 [ + + ]: 462 : do i_unitfi=1,2
566 [ + + ]: 308 : unitfi=ab_out;if (i_unitfi==2) unitfi=std_out
567 : : write(msg,'(3a)') &
568 : 308 : & ' ---------- DFT+U DATA --------------------------------------------------- ',ch10
569 : 308 : call wrtout(unitfi,msg,'COLL')
570 [ + + ]: 1146 : do iatom=1,dtset%natom
571 : 684 : itypat=dtset%typat(iatom);ll=pawtab(itypat)%lpawu
572 : 684 : nspden=paw_ij_all(iatom)%nspden;ndij=paw_ij_all(iatom)%ndij
573 : 684 : cplex_dij=paw_ij_all(iatom)%cplex_dij
574 [ + + + - ]: 992 : if ((ll>=0).and.(pawtab(itypat)%usepawu/=0)) then
575 : 430 : write(msg,fmt='(a,i5,a,i4,a)') " ====== For Atom ", iatom,&
576 : 860 : & ", occupations for correlated orbitals. lpawu =",ll,ch10
577 : 430 : call wrtout(unitfi,msg,'COLL')
578 [ + + ]: 430 : 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 [ + + ]: 430 : if(nspden==2) then
585 [ + + ]: 900 : do ispden=1,nspden
586 : 600 : write(msg,fmt='(a,i4,a,i3,a,f10.5)') " Atom", iatom,&
587 : 1200 : & ". Occ. for lpawu and for spin",ispden," =",paw_ij_all(iatom)%nocctot(ispden)
588 : 900 : call wrtout(unitfi,msg,'COLL')
589 : : end do
590 : 300 : localm=paw_ij_all(iatom)%nocctot(2)-paw_ij_all(iatom)%nocctot(1)
591 : 300 : write(msg,fmt='(a,i4,a,2x,f12.6)') " => On atom",iatom,&
592 : 600 : & ", local Mag. for lpawu is ",localm
593 : 300 : call wrtout(unitfi,msg,'COLL')
594 : : end if
595 [ + + ]: 430 : 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 : 430 : write(msg,'(3a)') ch10," == Occupation matrix for correlated orbitals:",ch10
617 : 430 : call wrtout(unitfi,msg,'COLL')
618 [ + + ]: 1370 : do ispden=1,ndij
619 [ + + + - ]: 940 : if (nspden==1.and.ndij/=4.and.(cplex_dij==1)) write(msg,fmt='(a)') " Up component only..."
620 [ + + ]: 940 : if (nspden==2) write(msg,fmt='(a,i3)')" Occupation matrix for spin",ispden
621 [ + + ]: 940 : if (ndij==4.or.(cplex_dij==2)) &
622 [ + - ]: 280 : & write(msg,fmt='(2a)') " Occupation matrix for component ",trim(dspin1(ispden+2*(ndij/4)))
623 : 940 : call wrtout(unitfi,msg,'COLL')
624 [ + + ]: 5684 : do im1=1,ll*2+1
625 [ + + ]: 4744 : if(cplex_dij==1)&
626 [ + - + + ]: 20340 : & write(msg,'(12(1x,9(1x,f10.5)))') (paw_ij_all(iatom)%noccmmp(1,im1,im2,ispden),im2=1,ll*2+1)
627 [ + + ]: 4744 : 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 : 5684 : call wrtout(unitfi,msg,'COLL')
631 : : end do
632 : 940 : write(msg,'(2a)') ch10,' '
633 : 1370 : call wrtout(unitfi,msg,'COLL')
634 : : end do
635 : : ! Transformation matrices: real->complex spherical harmonics
636 [ + + ]: 430 : 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 [ + + + - ]: 1424 : 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 : 1424 : msg=' '
741 : 1424 : call wrtout(ab_out,msg,'COLL')
742 : 1424 : call wrtout(std_out,msg,'COLL')
743 : :
744 : : !Destroy temporary stored atomic data
745 [ - + ]: 1424 : ABI_FREE(jatom)
746 : 1424 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
747 [ + + ]: 1424 : 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 [ + + - + : 3080 : end subroutine pawprt
- + - + -
+ - + - +
- + - + -
+ - + ]
759 : : !!***
760 : :
761 : : !----------------------------------------------------------------------
762 : :
763 : : END MODULE m_paw_tools
764 : : !!***
|
