Line data Source code
1 : !!****m* ABINIT/m_epjdos
2 : !! NAME
3 : !! m_epjdos
4 : !!
5 : !! FUNCTION
6 : !! Tools for the computiation of electronic PJDOSes
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MVer, XG, SM, MT, BAmadon, MG, MB)
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_epjdos
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_xmpi
27 : use m_errors
28 : use m_htetra
29 : use m_splines
30 : use m_cgtools
31 : use m_atomdata
32 : use m_crystal
33 : use m_ebands
34 : use m_nctk
35 : use netcdf
36 : use m_hdr
37 : use m_mpinfo
38 : use m_sort
39 : use m_dtset
40 :
41 : use defs_abitypes, only : MPI_type
42 : use defs_datatypes, only : pseudopotential_type
43 : use m_occ, only : dos_hdr_write
44 : use m_time, only : cwtime, timab
45 : use m_io_tools, only : open_file
46 : use m_numeric_tools, only : simpson, simpson_int
47 : use m_fstrings, only : int2char4, strcat
48 : use m_special_funcs, only : jlspline_t
49 : use m_kpts, only : tetra_from_kptrlatt
50 : use m_kg, only : ph1d3d, getph
51 : use m_gsphere, only : getkpgnorm
52 : use m_fftcore, only : sphereboundary
53 : use m_fft, only : fftpac, fourwf, fourdp
54 : use m_pawrad, only : pawrad_type, simp_gen
55 : use m_pawtab, only : pawtab_type
56 : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_get, pawcprj_free
57 : use m_initylmg, only : initylmg
58 :
59 : implicit none
60 :
61 : private
62 : !!***
63 :
64 : public :: recip_ylm ! Project input wavefunctions (real space) on to Ylm
65 : public :: dens_in_sph ! Calculate integrated density in sphere around each atom
66 : !!***
67 :
68 : !----------------------------------------------------------------------
69 :
70 : !!****t* m_epjdos/epjdos_t
71 : !! NAME
72 : !! epjdos_t
73 : !!
74 : !! FUNCTION
75 : !! Stores different contributions to the electronic DOS.
76 : !!
77 : !! NOTES
78 : !! Please contact gmatteo if you plan to change the internal implementation
79 : !! or add new DOSes. These results are saved in a netcdf file (see fatbands_ncwrite)
80 : !! so that one can read it with python and plot fatbands and PJDOSEs.
81 : !! The python version is able to handle the different cases (L, LM, Spin ...) but
82 : !! any change in the internal Abinit implementation is likely to break the python interface.
83 : !!
84 : !! SOURCE
85 :
86 : type,public :: epjdos_t
87 :
88 : integer :: mbesslang
89 : ! Max L+1 used in LM-DOS (Bessel function expansion)
90 :
91 : integer :: ndosfraction
92 : ! Defines the last dimension of the dos arrays.
93 : ! Actual value depends on the other variables.
94 :
95 : integer :: prtdos
96 : ! 2 --> Standard DOS with tetra.
97 : ! 3 --> L-DOS with tetra (prtdosm>0 if LM is wanted in Ylm/Slm basis).
98 : ! 4 --> L-DOS with gaussian (prtdosm if LM is wanted in Ylm/Slm basis).
99 : ! 5 --> Spin-DOS
100 :
101 : integer :: prtdosm
102 : ! Option for the m-contributions to the partial DOS
103 : ! 1 if LM-projection is done onto complex Ylm
104 : ! 2 if LM-projection is done onto real Slm
105 :
106 : integer :: partial_dos_flag
107 :
108 : integer :: paw_dos_flag
109 : ! 1 if both PAW contributions are evaluated AND stored
110 :
111 : !integer :: pawfatbnd
112 : integer :: fatbands_flag
113 :
114 : integer :: nkpt, mband, nsppol
115 : ! Used to dimension arrays
116 :
117 : integer,allocatable :: mlang_type(:)
118 : ! mlang_type(ntypat + natsph_extra)
119 : ! Max L+1 used in LM-DOS for each atom type
120 :
121 : real(dp),allocatable :: fractions(:,:,:,:)
122 : ! fractions(nkpt,mband,nsppol,ndosfraction))
123 : ! TODO: replace nsppol with nspden = 1, 2 (nsppol==2) or 4 (nspinor==2)
124 :
125 : real(dp),allocatable :: fractions_m(:,:,:,:)
126 : ! fractions_m(nkpt,mband,nsppol,ndosfraction*mbesslang)
127 :
128 : real(dp),allocatable :: fractions_paw1(:,:,:,:)
129 : ! fractions_paw1(nkpt,mband,nsppol,ndosfraction)
130 :
131 : real(dp),allocatable :: fractions_pawt1(:,:,:,:)
132 : ! fractions_pawt1(nkpt,mband,nsppol,ndosfraction))
133 :
134 : contains
135 :
136 : procedure :: init => epjdos_init
137 : ! Create new object
138 :
139 : procedure :: calcnwrite => epjdos_calcnwrite
140 : ! Calculate DOS and write results to file(s).
141 :
142 : procedure :: partial_dos_fractions => partial_dos_fractions
143 : ! Calculate partial DOS fractions to feed to the tetrahedron method (PW part)
144 :
145 : procedure :: partial_dos_fractions_paw => partial_dos_fractions_paw
146 : ! Calculate PAW contributions to the partial DOS fractions.
147 :
148 : procedure :: prtfatbands => epjdos_prtfatbands
149 : ! Print PJDOS contributions in xmgrace format.
150 :
151 : procedure :: ncwrite => epjdos_ncwrite
152 : ! Write PJDOS contributions to netcdf file.
153 :
154 : procedure :: free => epjdos_free
155 : ! Free dynamic memory
156 :
157 : end type epjdos_t
158 :
159 : !----------------------------------------------------------------------
160 :
161 : contains !============================================================
162 : !!***
163 :
164 : !!****f* m_epjdos/epjdos_init
165 : !! NAME
166 : !! epjdos_init
167 : !!
168 : !! FUNCTION
169 : !! Create new object from dataset input variables.
170 : !!
171 : !! INPUTS
172 : !! dtset <type(dataset_type)>=all input variables for this dataset
173 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
174 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
175 : !!
176 : !! SOURCE
177 :
178 49 : subroutine epjdos_init(new, dtset, psps, pawtab)
179 :
180 : !Arguments ------------------------------------
181 : class(epjdos_t),intent(inout) :: new
182 : type(dataset_type),intent(in) :: dtset
183 : type(pseudopotential_type),intent(in) :: psps
184 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
185 :
186 : !Local variables-------------------------------
187 : !scalars
188 : integer :: ierr,itypat,iat
189 : ! *********************************************************************
190 :
191 : !write(std_out,*)' m_epjdos%epjdos_new, enter '
192 :
193 49 : new%nkpt = dtset%nkpt; new%mband = dtset%mband; new%nsppol = dtset%nsppol
194 :
195 49 : new%prtdos = dtset%prtdos
196 49 : new%partial_dos_flag = 0
197 : if (new%prtdos==2) new%partial_dos_flag = 0 ! Standard DOS with tetra.
198 49 : if (new%prtdos==3) new%partial_dos_flag = 1 ! L-DOS with tetra (prtdosm>0 if LM is wanted in Ylm/Slm basis).
199 49 : if (new%prtdos==4) new%partial_dos_flag = 1 ! L-DOS with gaussian (prtdosm if LM is wanted in Ylm/Slm basis).
200 49 : if (new%prtdos==5) new%partial_dos_flag = 2 ! Spin DOS
201 :
202 49 : new%prtdosm=0
203 49 : if (new%partial_dos_flag==1) new%prtdosm=dtset%prtdosm
204 : ! paw_dos_flag= 1 if both PAW contributions are evaluated AND stored
205 49 : new%paw_dos_flag=0
206 49 : if (dtset%usepaw==1 .and. new%partial_dos_flag==1 .and. dtset%pawprtdos==1) new%paw_dos_flag=1
207 :
208 49 : new%fatbands_flag=0
209 49 : if (dtset%pawfatbnd>0 .and. new%prtdosm==0) new%fatbands_flag=1
210 49 : if (new%prtdosm==1.and.dtset%pawfatbnd>0)then
211 : ! because they compute quantities in real and complex harmonics respectively
212 0 : ABI_ERROR('pawfatbnd>0 and prtdosm=1 are not compatible')
213 : end if
214 :
215 : ! mjv : initialization is needed as mbesslang is used for allocation below
216 : ! NOTE: 10/5/2010 the whole of this could be looped over ndosfraction,
217 : ! to store much less in memory. The DOS is accumulated in an array
218 : ! and then printed to file at the end.
219 49 : new%mbesslang = 1
220 49 : if (new%partial_dos_flag==1 .or. new%fatbands_flag==1) then
221 :
222 132 : ABI_MALLOC(new%mlang_type, (dtset%ntypat + dtset%natsph_extra))
223 110 : new%mlang_type = 0
224 :
225 : ! TODO: Could use mbesslang = 4 or compute it from psps/pawtab
226 : ! Increment by one (could underestimate if vloc = vlmax)
227 44 : if (dtset%usepaw == 0) then
228 52 : do iat=1,dtset%natsph
229 32 : itypat = dtset%typat(dtset%iatsph(iat))
230 173 : new%mlang_type(itypat) = 1 + maxval(psps%indlmn(1, :, itypat))
231 : end do
232 : else
233 62 : do iat=1,dtset%natsph
234 38 : itypat= dtset%typat(dtset%iatsph(iat))
235 62 : new%mlang_type(itypat) = 1 + (pawtab(itypat)%l_size - 1) / 2
236 : end do
237 : end if
238 :
239 : ! Up to l=g if we have natsph_extra.
240 56 : if (dtset%natsph_extra > 0) new%mlang_type(dtset%ntypat+1:) = 5
241 :
242 110 : new%mlang_type = 5 ! This is to preserve the old implementation
243 110 : new%mbesslang = maxval(new%mlang_type)
244 44 : new%ndosfraction = (dtset%natsph + dtset%natsph_extra) * new%mbesslang
245 :
246 5 : else if (new%partial_dos_flag == 2) then
247 1 : new%ndosfraction = 7
248 :
249 : else
250 4 : new%ndosfraction = 1
251 4 : new%mbesslang = 0
252 : end if
253 :
254 : ! Check allocations status as these arrays are not distributed and the wavefunctions are still in memory.
255 294 : ABI_MALLOC_OR_DIE(new%fractions, (dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr)
256 66018 : new%fractions = zero
257 :
258 49 : if (new%prtdosm>=1 .or. new%fatbands_flag==1) then
259 18 : ABI_MALLOC_OR_DIE(new%fractions_m,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction*new%mbesslang), ierr)
260 57874 : new%fractions_m = zero
261 : end if
262 :
263 49 : if (dtset%usepaw==1 .and. new%partial_dos_flag==1) then
264 138 : ABI_MALLOC_OR_DIE(new%fractions_paw1,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr)
265 138 : ABI_MALLOC_OR_DIE(new%fractions_pawt1,(dtset%nkpt,dtset%mband,dtset%nsppol,new%ndosfraction), ierr)
266 66663 : new%fractions_paw1 = zero; new%fractions_pawt1 = zero
267 : end if
268 :
269 : !write(std_out,*)' m_epjdos%epjdos_new, exit '
270 :
271 49 : end subroutine epjdos_init
272 : !!***
273 :
274 : !!****f* m_epjdos/epjdos_free
275 : !! NAME
276 : !! epjdos_free
277 : !!
278 : !! FUNCTION
279 : !! Deallocate memory
280 : !!
281 : !! SOURCE
282 :
283 1774 : subroutine epjdos_free(self)
284 :
285 : !Arguments ------------------------------------
286 : class(epjdos_t),intent(inout) :: self
287 : ! *********************************************************************
288 :
289 : ! integer
290 1774 : ABI_SFREE(self%mlang_type)
291 :
292 : ! real
293 1774 : ABI_SFREE(self%fractions)
294 1774 : ABI_SFREE(self%fractions_m)
295 1774 : ABI_SFREE(self%fractions_paw1)
296 1774 : ABI_SFREE(self%fractions_pawt1)
297 :
298 1774 : end subroutine epjdos_free
299 : !!***
300 :
301 : !!****f* m_epjdos/epjdos_calcnwrite
302 : !! NAME
303 : !! epjdos_calcnwrite
304 : !!
305 : !! FUNCTION
306 : !! calculate DOS and write results to file(s)
307 : !!
308 : !! INPUTS
309 : !! dos_fractions= projections of wavefunctions on each angular momentum Ylm
310 : !! which is the weight going into the DOS for an l-decomposed dos
311 : !! dos_fractions_m= same as dos_fractions, but m-decomposed not just l-
312 : !! dos_fractions_paw1= contribution to dos fractions from the PAW partial waves (phi)
313 : !! dos_fractions_pawt1= contribution to dos fractions from the PAW pseudo partial waves (phi_tild)
314 : !! dtset structured datatype, in particular one uses :
315 : !! kptrlatt(3,3)=lattice vectors for full kpoint grid
316 : !! nshiftk =number of kpoint grid shifts
317 : !! pawprtdos =option to output the individual contributions to the partial DOS (0, 1 or 2)
318 : !! shiftk(3,nshiftk)=kpoint shifts
319 : !! usepaw =option for PAW
320 : !! crystal<crystal_t>=Object defining the unit cell and its symmetries.
321 : !! ebands<ebands_t>=Band structure data.
322 : !! fermie=Fermi energy
323 : !! fildata=name of the DOS output file
324 : !! mbesslang=maximum angular momentum for Bessel function expansion
325 : !! prtdosm=option for the m-contributions to the partial DOS
326 : !! ndosfraction= number of types of DOS we are calculating, e.g. the number
327 : !! of l channels. Could be much more general, for other types of partial DOS
328 : !! paw_dos_flag= option for partial dos in PAW
329 : !! comm=MPI communicator.
330 : !!
331 : !! OUTPUT
332 : !! (no explicit output)
333 : !!
334 : !! SOURCE
335 :
336 46 : subroutine epjdos_calcnwrite(dos,dtset,crystal,ebands,fildata,comm)
337 :
338 : !Arguments ------------------------------------
339 : !scalars
340 : class(epjdos_t),intent(in) :: dos
341 : integer,intent(in) :: comm
342 : character(len=*),intent(in) :: fildata
343 : type(dataset_type),intent(in) :: dtset
344 : type(crystal_t),intent(in) :: crystal
345 : type(ebands_t),intent(in) :: ebands
346 :
347 : !Local variables-------------------------------
348 : !scalars
349 : integer,parameter :: bcorr0=0,master=0
350 : integer :: iat,iband,iene,ikpt,isppol,natsph,natsph_extra,nkpt,nsppol,i1,i2
351 : integer :: nene,prtdos,unitdos,ierr,prtdosm,paw_dos_flag,mbesslang,ndosfraction
352 : integer :: my_rank,nprocs,cnt,ifrac,ii
353 : real(dp),parameter :: dos_max=9999.9999_dp
354 : real(dp) :: buffer,deltaene,enemax,enemin,integral_DOS,max_occ
355 : real(dp) :: cpu,wall,gflops
356 : logical :: bigDOS,iam_master
357 : character(len=10) :: tag
358 : character(len=500) :: frmt,frmt_extra,msg
359 46 : type(htetra_t) :: tetra
360 : !arrays
361 46 : integer,allocatable :: unt_atsph(:)
362 : real(dp) :: list_dp(3)
363 46 : real(dp),allocatable :: tmp_eigen(:),total_dos(:,:,:),eig_dos(:,:)
364 46 : real(dp),allocatable :: dos_m(:,:,:),dos_paw1(:,:,:),dos_pawt1(:,:,:), wdt(:,:)
365 : ! *********************************************************************
366 :
367 : !write(std_out,*)' m_epjdos%dos_calcncwrite, enter '
368 :
369 46 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm); iam_master = (my_rank == master)
370 :
371 46 : prtdosm = dos%prtdosm; paw_dos_flag = dos%paw_dos_flag
372 46 : mbesslang = dos%mbesslang; ndosfraction = dos%ndosfraction
373 :
374 46 : nkpt = dtset%nkpt; nsppol = dtset%nsppol
375 :
376 : !m-decomposed DOS not compatible with PAW-decomposed DOS
377 46 : if (prtdosm>=1.and.paw_dos_flag==1) then
378 0 : msg = 'm-decomposed DOS (prtdosm>=1) not compatible with PAW-decomposed DOS (pawprtdos=1) !'
379 0 : ABI_ERROR(msg)
380 : end if
381 :
382 : !Refuse nband different for different kpoints
383 : !Note: This means we can pass ebands%eig(:,:,:) instead of eigen(mband*nkpt*nsppol) in packed form
384 98 : do isppol=1,nsppol
385 486 : do ikpt=1,nkpt
386 440 : if ( dtset%nband(nkpt*(isppol-1) + ikpt) /= dtset%nband(1) ) then
387 0 : write(std_out,*) 'tetrahedron: skip subroutine.'
388 0 : write(std_out,*) 'nband must be the same for all kpoints'
389 0 : write(std_out,*) 'nband=', dtset%nband
390 0 : ABI_WARNING('tetrahedron: skip subroutine. See message above')
391 0 : return
392 : end if
393 : end do
394 : end do
395 :
396 46 : call cwtime(cpu, wall, gflops, "start")
397 :
398 : tetra = tetra_from_kptrlatt(crystal, dtset%kptopt, dtset%kptrlatt, dtset%nshiftk, &
399 46 : dtset%shiftk, dtset%nkpt, dtset%kpt, comm, msg, ierr)
400 46 : if (ierr /= 0) then
401 0 : call tetra%free()
402 0 : ABI_WARNING(msg)
403 0 : return
404 : end if
405 :
406 46 : natsph=dtset%natsph; natsph_extra=dtset%natsph_extra
407 :
408 : ! Master opens the DOS files.
409 46 : if (iam_master) then
410 25 : if (any(dtset%prtdos == [2, 5])) then
411 5 : if (open_file(fildata, msg, newunit=unitdos, status='unknown', form='formatted', action="write") /= 0) then
412 0 : ABI_ERROR(msg)
413 : end if
414 :
415 20 : else if (dtset%prtdos == 3) then
416 : ! unt_atsph(0) is used for the total DOS.
417 60 : ABI_MALLOC(unt_atsph,(0:natsph+natsph_extra))
418 :
419 : ! Open file for total DOS as well.
420 20 : if (open_file(strcat(fildata, '_TOTAL'), msg, newunit=unt_atsph(0), &
421 : status='unknown', form='formatted', action="write") /= 0) then
422 20 : ABI_ERROR(msg)
423 : end if
424 :
425 50 : do iat=1,natsph
426 30 : call int2char4(dtset%iatsph(iat),tag)
427 30 : ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
428 30 : if (open_file(strcat(fildata, '_AT', tag), msg, newunit=unt_atsph(iat), &
429 20 : status='unknown', form='formatted', action="write") /= 0) then
430 30 : ABI_ERROR(msg)
431 : end if
432 : end do
433 : ! do extra spheres in vacuum too. Use _ATEXTRA[NUM] suffix
434 32 : do iat=1,natsph_extra
435 12 : call int2char4(iat,tag)
436 12 : ABI_CHECK((tag(1:1)/='#'),'Bug: string length too short!')
437 12 : if (open_file(strcat(fildata, '_ATEXTRA', tag), msg, newunit=unt_atsph(natsph+iat), &
438 20 : status='unknown', form='formatted', action="write") /= 0) then
439 12 : ABI_ERROR(msg)
440 : end if
441 : end do
442 : end if
443 : end if
444 :
445 : ! Write the header of the DOS file, and determine the energy range and spacing
446 46 : prtdos=dtset%prtdos
447 46 : buffer=0.01_dp ! Size of the buffer around the min and max ranges
448 :
449 : ! A Similar section is present is getnel. Should move all DOS stuff to m_ebands
450 : ! Choose the lower and upper energies
451 6323 : enemax = maxval(ebands%eig) + buffer
452 6323 : enemin = minval(ebands%eig) - buffer
453 :
454 : ! Extend the range to a nicer value
455 46 : enemax=0.1_dp*ceiling(enemax*10._dp)
456 46 : enemin=0.1_dp*floor(enemin*10._dp)
457 :
458 : ! Choose the energy increment
459 46 : if(abs(dtset%dosdeltae)<tol10)then
460 44 : deltaene=0.001_dp
461 44 : if(dtset%prtdos>=2)deltaene=0.0005_dp ! Higher resolution possible (and wanted) for tetrahedron
462 : else
463 2 : deltaene=dtset%dosdeltae
464 : end if
465 46 : nene=nint((enemax-enemin)/deltaene)+1
466 :
467 46 : call xmpi_bcast(nene, master, comm, ierr)
468 121 : if (iam_master) list_dp(1:3) = [deltaene, enemin, enemax]
469 46 : call xmpi_bcast(list_dp, master, comm, ierr)
470 46 : deltaene = list_dp(1); enemin = list_dp(2); enemax = list_dp(3)
471 :
472 46 : if (iam_master) then
473 25 : if (any(dtset%prtdos == [2, 5])) then
474 : call dos_hdr_write(deltaene,ebands%eig,enemax,enemin,ebands%fermie,ebands%fermih,&
475 : dtset%mband,dtset%nband,nene,nkpt,nsppol,dtset%occopt,prtdos,&
476 5 : dtset%tphysel,dtset%tsmear,unitdos)
477 20 : else if (dtset%prtdos == 3) then
478 82 : do iat=0,natsph+natsph_extra
479 : call dos_hdr_write(deltaene,ebands%eig,enemax,enemin,ebands%fermie,ebands%fermih,&
480 : dtset%mband,dtset%nband,nene,nkpt,nsppol,dtset%occopt,prtdos,&
481 82 : dtset%tphysel,dtset%tsmear,unt_atsph(iat))
482 : end do
483 : end if
484 : end if
485 :
486 : ! Tetra weights
487 138 : ABI_MALLOC(wdt, (nene, 2))
488 :
489 : ! Allocate arrays to store DOSes and fill with zeros.
490 : ! 1--> DOS , 2--> IDOS
491 230 : ABI_MALLOC(total_dos,(nene,ndosfraction,2))
492 92 : ABI_MALLOC(eig_dos, (nene, 2))
493 :
494 46 : if (paw_dos_flag==1) then
495 0 : ABI_MALLOC(dos_paw1,(nene,ndosfraction,2))
496 0 : ABI_MALLOC(dos_pawt1,(nene,ndosfraction,2))
497 : end if
498 46 : if (prtdosm>=1) then
499 10 : ABI_MALLOC(dos_m, (nene,ndosfraction*mbesslang,2))
500 : end if
501 :
502 : !Get maximum occupation value (2 or 1)
503 46 : max_occ = one; if (dtset%nspinor == 1 .and. nsppol == 1) max_occ = two
504 :
505 : !-------------------------------------------------------------------
506 : !For each spin polarisation and band, interpolate band over kpoints
507 : !calculate integration weights and DOS contib from
508 : !-------------------------------------------------------------------
509 :
510 : ! Workspace arrays.
511 138 : ABI_MALLOC(tmp_eigen,(nkpt))
512 :
513 46 : cnt = 0
514 98 : do isppol=1,nsppol
515 :
516 7702572 : total_dos = zero; eig_dos = zero
517 900458 : if (prtdosm>=1) dos_m = zero
518 52 : if (paw_dos_flag==1) then
519 0 : dos_paw1 = zero; dos_pawt1 = zero
520 : end if
521 :
522 440 : do ikpt=1,nkpt
523 6277 : do iband=1,ebands%nband(ikpt+(isppol-1)*ebands%nkpt)
524 5837 : cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! Mpi parallelism.
525 :
526 : ! Accumulate total DOS from eigenvalues (this is the **exact** total DOS)
527 24672 : tmp_eigen(:) = ebands%eig(iband, :, isppol)
528 2837 : call tetra%get_onewk(ikpt,bcorr0,nene,nkpt,tmp_eigen,enemin,enemax,max_occ,wdt)
529 48586585 : wdt = wdt*ebands%wtk(ikpt)
530 48589422 : eig_dos = eig_dos + wdt
531 :
532 : ! Accumulate L-DOS.
533 8511 : do ii=1,2
534 46713 : do ifrac=1,ndosfraction
535 383816878 : total_dos(:,ifrac,ii) = total_dos(:,ifrac,ii) + wdt(:,ii) * dos%fractions(ikpt,iband,isppol,ifrac)
536 : end do
537 : end do
538 :
539 2837 : if (paw_dos_flag==1) then
540 : ! Accumulate L-DOS (on-site terms).
541 0 : do ii=1,2
542 0 : do ifrac=1,ndosfraction
543 0 : dos_paw1(:,ifrac,ii) = dos_paw1(:,ifrac,ii) + wdt(:,ii) * dos%fractions_paw1(ikpt,iband,isppol,ifrac)
544 0 : dos_pawt1(:,ifrac,ii) = dos_pawt1(:,ifrac,ii) + wdt(:,ii) * dos%fractions_pawt1(ikpt,iband,isppol,ifrac)
545 : end do
546 : end do
547 : end if
548 :
549 3225 : if (prtdosm>=1) then
550 : ! Accumulate LM-DOS.
551 609 : do ii=1,2
552 10909 : do ifrac=1,ndosfraction*mbesslang
553 86141006 : dos_m(:,ifrac,ii) = dos_m(:,ifrac,ii) + wdt(:,ii) * dos%fractions_m(ikpt, iband, isppol, ifrac)
554 : end do
555 : end do
556 : end if
557 :
558 : end do ! ikpt
559 : end do ! iband
560 :
561 : ! Collect results on master
562 52 : call xmpi_sum_master(eig_dos, master, comm, ierr)
563 52 : call xmpi_sum_master(total_dos, master, comm, ierr)
564 3446676 : bigDOS=(maxval(total_dos(:,:,1))>999._dp)
565 :
566 52 : if (paw_dos_flag == 1) then
567 0 : call xmpi_sum_master(dos_paw1, master, comm, ierr)
568 0 : call xmpi_sum_master(dos_pawt1, master, comm, ierr)
569 : end if
570 52 : if (prtdosm >= 1) call xmpi_sum_master(dos_m, master, comm, ierr)
571 :
572 : ! Write the DOS value in the DOS file
573 : ! Print the data for this energy. Note the upper limit (dos_max), to be consistent with the format.
574 : ! The use of "E" format is not adequate, for portability of the self-testing procedure.
575 : ! header lines depend on the type of DOS (projected etc...) which is output
576 :
577 52 : if (.not. iam_master) goto 10
578 31 : call write_extra_headers()
579 :
580 31 : if (prtdos==2) then
581 : ! E, DOS, IDOS
582 23610 : do iene=1,nene
583 : write(unitdos, '(f11.5,1x,2(f10.4,1x))') &
584 47215 : enemin + (iene-1)*deltaene, min(total_dos(iene,:,1), dos_max), total_dos(iene,:,2)
585 : end do
586 :
587 26 : else if (prtdos==3) then
588 :
589 : ! Write E, DOS, IDOS
590 142330 : do iene=1,nene
591 : write(unt_atsph(0), '(f11.5,1x,2(f10.4,1x))') &
592 142330 : enemin + (iene-1)*deltaene, min(eig_dos(iene,1), dos_max), eig_dos(iene,2)
593 : end do
594 :
595 : ! E, DOS(L=1,LMAX), IDOS(L=1,LMAX)
596 : ! Here we assume mpsang = 5 in the format.
597 25 : if (paw_dos_flag/=1.or.dtset%pawprtdos==2) then
598 25 : frmt = '(f11.5,1x,5(f9.4,1x),10x,5(f8.2,1x),10x,25(f8.2,1x))'
599 25 : if (bigDOS) frmt = '(f11.5,1x,5(f10.4,1x),10x,5(f8.2,1x),10x,25(f8.2,1x))'
600 : ! for extra atoms in vacuum need more precision
601 25 : frmt_extra = '(f11.5,1x,5(f20.16,1x),10x,5(f20.16,1x),10x,25(f20.16,1x))'
602 :
603 60 : do iat=1,natsph
604 35 : i1 = (iat-1)*mbesslang+1; i2 = iat*mbesslang
605 60 : if (prtdosm==0) then
606 216542 : do iene=1,nene
607 216511 : write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, &
608 1515608 : min(total_dos(iene, i1:i2, 1), dos_max), total_dos(iene, i1:i2,2)
609 : end do
610 : else
611 18008 : do iene=1,nene
612 18004 : write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, &
613 108024 : min(total_dos(iene, i1:i2, 1), dos_max),&
614 18004 : total_dos(iene, i1:i2, 2),&
615 486112 : min(dos_m(iene,(iat-1)*mbesslang**2+1:iat*mbesslang**2,1), dos_max)
616 : end do
617 : end if
618 : end do
619 :
620 : ! Extra spheres.
621 37 : do iat=natsph+1,natsph+natsph_extra
622 12 : i1 = (iat-1)*mbesslang+1; i2 = iat*mbesslang
623 37 : if (prtdosm==0) then
624 14424 : do iene=1,nene
625 14412 : write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, &
626 14412 : total_dos(iene, i1:i2, 1), &
627 28836 : total_dos(iene, i1:i2, 2)
628 : end do
629 : else
630 0 : do iene=1,nene
631 0 : write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, &
632 0 : total_dos(iene, i1:i2, 1),&
633 0 : total_dos(iene, i1:i2, 2),&
634 0 : dos_m(iene,(iat-1)*mbesslang**2+1:iat*mbesslang**2, 1)
635 : end do
636 : end if
637 : end do
638 :
639 : else
640 0 : frmt = '(f11.5,1x,5(f9.4,1x),3(6x,5f9.4))'
641 0 : if (bigDOS) frmt = '(f11.5,1x,5(f10.4,1x),3(6x,5f10.4))'
642 : ! for extra atom spheres in vacuum need more precision
643 0 : frmt_extra = '(f11.5,1x,5(f20.16,1x),3(6x,5f20.16))'
644 :
645 0 : do iat=1,natsph
646 0 : i1 = iat*5-4; i2 = iat*5
647 0 : do iene=1,nene
648 0 : write(unt_atsph(iat), fmt=frmt) enemin + (iene-1)*deltaene, &
649 0 : min(total_dos(iene,i1:i2,1), dos_max),&
650 0 : min(total_dos(iene,i1:i2,1) - dos_paw1(iene,i1:i2,1) + dos_pawt1(iene,i1:i2,1), dos_max),&
651 0 : min(dos_paw1(iene,i1:i2,1), dos_max),&
652 0 : min(dos_pawt1(iene,i1:i2,1), dos_max)
653 : end do
654 : end do
655 :
656 : ! Extra spheres.
657 0 : do iat=natsph+1,natsph+natsph_extra
658 0 : i1 = iat*5-4; i2 = iat*5
659 0 : do iene=1,nene
660 0 : write(unt_atsph(iat), fmt=frmt_extra) enemin + (iene-1)*deltaene, &
661 0 : min(total_dos(iene,i1:i2,1), dos_max),&
662 0 : min(total_dos(iene,i1:i2,1) - dos_paw1(iene,i1:i2,1) + dos_pawt1(iene,i1:i2,1), dos_max),&
663 0 : min(dos_paw1(iene,i1:i2,1), dos_max),&
664 0 : min(dos_pawt1(iene,i1:i2,1), dos_max)
665 : end do
666 : end do
667 : end if
668 :
669 1 : else if (prtdos==5)then
670 : ! E, SPIN-DOS
671 1 : frmt = '(f11.5,1x,7(f9.4,1x),10x,7(f8.2,1x))'
672 1 : if (bigDOS) frmt = '(f11.5,1x,7(f10.4,1x),10x,7(f8.2,1x))'
673 9402 : do iene=1,nene
674 75209 : write(unitdos, fmt=frmt) enemin + (iene-1)*deltaene, min(total_dos(iene,1:7,1), dos_max), total_dos(iene,1:7,2)
675 : end do
676 : end if
677 :
678 : 10 continue
679 464 : integral_DOS=sum(total_dos(nene,:,2))
680 52 : write(msg, '(a,es16.8)' ) ' tetrahedron : integrate to',integral_DOS
681 98 : call wrtout(std_out,msg)
682 : end do ! isppol
683 :
684 : ! Close files.
685 46 : if (iam_master) then
686 25 : if (any(prtdos == [2, 5])) then
687 5 : close(unitdos)
688 20 : else if (prtdos == 3) then
689 82 : do iat=0,natsph+natsph_extra
690 82 : close(unt_atsph(iat))
691 : end do
692 20 : ABI_FREE(unt_atsph)
693 : end if
694 : end if
695 :
696 46 : ABI_FREE(tmp_eigen)
697 46 : ABI_FREE(total_dos)
698 46 : ABI_FREE(wdt)
699 46 : ABI_FREE(eig_dos)
700 :
701 46 : if (prtdosm>=1) then
702 2 : ABI_FREE(dos_m)
703 : end if
704 :
705 46 : if (paw_dos_flag==1) then
706 0 : ABI_FREE(dos_paw1)
707 0 : ABI_FREE(dos_pawt1)
708 : end if
709 :
710 46 : call tetra%free()
711 :
712 46 : call cwtime(cpu,wall,gflops,"stop")
713 46 : write(msg,'(2(a,f8.2),a)')" tetrahedron: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
714 138 : call wrtout(std_out,msg,"PERS")
715 :
716 : !write(std_out,*)' m_epjdos%dos_calcncwrite, exit '
717 :
718 : contains
719 :
720 31 : subroutine write_extra_headers()
721 :
722 31 : if (nsppol==2) then
723 12 : if(isppol==1) write(msg,'(a,16x,a)') '#','Spin-up DOS'
724 12 : if(isppol==2) write(msg,'(2a,16x,a)') ch10,'#','Spin-dn DOS'
725 : ! NB: dtset%prtdos == 5 should not happen for nsppol==2
726 :
727 12 : if (any(dtset%prtdos == [2, 5])) then
728 2 : write(unitdos, "(a)")trim(msg)
729 :
730 10 : else if (dtset%prtdos == 3) then
731 30 : do iat=0,natsph+natsph_extra
732 30 : write(unt_atsph(iat), "(a)")trim(msg)
733 : end do
734 : end if
735 :
736 : end if
737 :
738 31 : if (prtdos==2) then
739 5 : write(unitdos, '(a)' )'# energy(Ha) DOS integrated DOS'
740 :
741 26 : else if (prtdos==3) then
742 :
743 25 : write(unt_atsph(0), '(a)' )'# energy(Ha) DOS integrated DOS'
744 :
745 25 : if (paw_dos_flag/=1.or.dtset%pawprtdos==2) then
746 60 : do iat=1,natsph
747 : write(unt_atsph(iat), '(3a,i5,a,i5,a,a,es16.6,3a)' ) &
748 35 : '# Local DOS (columns 2-6) and integrated local DOS (columns 7-11),',ch10,&
749 35 : '# for atom number iat=',iat,' iatom=',dtset%iatsph(iat),ch10,&
750 70 : '# inside sphere of radius ratsph=',dtset%ratsph(dtset%typat(dtset%iatsph(iat))),' Bohr.',ch10,"#"
751 :
752 35 : if (dtset%usepaw==1.and.dtset%pawprtdos==2) then
753 : write(unt_atsph(iat), '(3a)' ) &
754 6 : '# PAW: note that only all-electron on-site part has been used to compute DOS !',ch10,"#"
755 : end if
756 35 : if (bigDOS) then
757 : write(msg, '(a,a)' ) &
758 14 : '# energy(Ha) l=0 l=1 l=2 l=3 l=4',&
759 28 : ' (integral=>) l=0 l=1 l=2 l=3 l=4'
760 : else
761 : write(msg, '(a,a)' ) &
762 21 : '# energy(Ha) l=0 l=1 l=2 l=3 l=4',&
763 42 : ' (integral=>) l=0 l=1 l=2 l=3 l=4'
764 : end if
765 35 : if (prtdosm>=1) then
766 4 : write(msg, '(7a)' ) trim(msg),' ',&
767 4 : ' lm=0 0',&
768 4 : ' lm=1-1 lm=1 0 lm=1 1',&
769 4 : ' lm=2-2 lm=2-1 lm=2 0 lm=2 1 lm=2 2',&
770 4 : ' lm=3-3 lm=3-2 lm=3-1 lm=3 0 lm=3 1 lm=3 2 lm=3 3',&
771 8 : ' lm=4-4 lm=4-3 lm=4-2 lm=4-1 lm=4 0 lm=4 1 lm=4 2 lm=4 3 lm=4 4'
772 : end if
773 60 : write(unt_atsph(iat), "(a)")trim(msg)
774 : end do
775 : else
776 0 : do iat=1,natsph
777 : write(unt_atsph(iat), '(9a,i5,a,i5,a,a,es16.6,3a)' ) &
778 0 : '# Local DOS (columns 2-6),',ch10,&
779 0 : '# plane-waves contrib. to DOS (columns 7-11),',ch10,&
780 0 : '# AE on-site contrib. to DOS (columns 12-16),',ch10,&
781 0 : '# -PS on-site contrib. to DOS (columns 17-21),',ch10,&
782 0 : '# for atom number iat=',iat,' iatom=',dtset%iatsph(iat),ch10,&
783 0 : '# inside sphere of radius ratsph=',dtset%ratsph(dtset%typat(dtset%iatsph(iat))),' Bohr.',ch10,"#"
784 0 : if (bigDOS) then
785 : write(msg, '(4a)' ) &
786 0 : '#energy(Ha) l=0 l=1 l=2 l=3 l=4',&
787 0 : ' (PW) l=0 l=1 l=2 l=3 l=4',&
788 0 : ' (Phi) l=0 l=1 l=2 l=3 l=4',&
789 0 : ' (tPhi) l=0 l=1 l=2 l=3 l=4'
790 : else
791 : write(msg, '(4a)' ) &
792 0 : '#energy(Ha) l=0 l=1 l=2 l=3 l=4',&
793 0 : ' (PW) l=0 l=1 l=2 l=3 l=4',&
794 0 : ' (Phi) l=0 l=1 l=2 l=3 l=4',&
795 0 : ' (tPhi) l=0 l=1 l=2 l=3 l=4'
796 : end if
797 0 : write(unt_atsph(iat), "(a)")trim(msg)
798 : end do
799 : end if
800 37 : do iat=1,natsph_extra
801 : write(unt_atsph(natsph+iat), '(3a,i5,2a,es16.6,3a)' ) &
802 12 : '# Local DOS (columns 2-6) and integrated local DOS (columns 7-11),',ch10,&
803 12 : '# for non-atomic sphere number iat=',iat,ch10,&
804 24 : '# of radius ratsph=',dtset%ratsph_extra,' Bohr.',ch10,"#"
805 12 : if (bigDOS) then
806 : write(msg, '(a,a)' ) &
807 0 : '# energy(Ha) l=0 l=1 l=2 l=3 l=4',&
808 0 : ' (integral=>) l=0 l=1 l=2 l=3 l=4'
809 : else
810 : write(msg, '(a,a)' ) &
811 12 : '# energy(Ha) l=0 l=1 l=2 l=3 l=4',&
812 24 : ' (integral=>) l=0 l=1 l=2 l=3 l=4'
813 : end if
814 12 : if (prtdosm>=1) then
815 0 : write(msg, '(7a)' ) trim(msg),' ',&
816 0 : ' lm=0 0',&
817 0 : ' lm=1-1 lm=1 0 lm=1 1',&
818 0 : ' lm=2-2 lm=2-1 lm=2 0 lm=2 1 lm=2 2',&
819 0 : ' lm=3-3 lm=3-2 lm=3-1 lm=3 0 lm=3 1 lm=3 2 lm=3 3',&
820 0 : ' lm=4-4 lm=4-3 lm=4-2 lm=4-1 lm=4 0 lm=4 1 lm=4 2 lm=4 3 lm=4 4'
821 : end if
822 37 : write(unt_atsph(natsph+iat), "(a)")trim(msg)
823 : end do
824 :
825 1 : else if (prtdos==5) then
826 : write(unitdos, '(a)' )&
827 1 : '# energy(Ha) DOS up,up up,dn dn,up dn,dn sigma_x sigma_y sigma_z and integrated DOS components'
828 : end if ! prtdos value
829 :
830 31 : end subroutine write_extra_headers
831 :
832 : end subroutine epjdos_calcnwrite
833 : !!***
834 :
835 : !!****f* m_epjdos/recip_ylm
836 : !! NAME
837 : !! recip_ylm
838 : !!
839 : !! FUNCTION
840 : !! Project input wavefunctions in reciprocal space on to Ylm
841 : !! (real or complex harmonics depending on rc_ylm).
842 : !!
843 : !! INPUTS
844 : !! bess_fit(mpw,nradintmax,ll) = Bessel functions for L, splined
845 : !! with arguments $2 \pi |k+G| \Delta r$, for all G vectors in sphere
846 : !! and all points on radial grid.
847 : !! cg_1band(2,npw_k)=wavefunction in recip space (note that nspinor is missing, see Notes).
848 : !! istwfk= storage mode of cg_1band
849 : !! nradint(natsph)=number of points on radial real-space grid for a given atom.
850 : !! nradintmax=dimension of rint array.
851 : !! me_g0=1 if this processor has G=0, 0 otherwise
852 : !! mlang=maximum angular momentum in Bessel functions.
853 : !! mpw=Maximum number of planewaves. Used to dimension bess_fit
854 : !! natsph=number of atoms around which ang mom projection has to be done
855 : !! typat_extra(natsph)=Type of each atom. ntypat + 1 if empty sphere
856 : !! mlang_type(ntypat + natsph_extra)=Max L+1 for each atom type
857 : !! npw_k=number of plane waves for this kpt
858 : !! nspinor=number of spinor components
859 : !! ph3d(2,npw_k,natsph)=3-dim structure factors, for each atom and plane wave.
860 : !! prtsphere= if 1, print a complete analysis of the angular momenta in atomic spheres
861 : !! rint(nradintmax) = points on radial real-space grid for integration
862 : !! rmax(natsph)=maximum radius for real space integration sphere
863 : !! rc_ylm= 1 for real spherical harmonics. 2 for complex spherical harmonics,
864 : !! ucvol=unit cell volume in bohr**3.
865 : !! ylm_k(npw_k,mlang**2)=real spherical harmonics for each G and LM.
866 : !! znucl_sph(natsph)=gives the nuclear number for each type of atom
867 : !!
868 : !! OUTPUT
869 : !! sum_1ll_1atom(mlang,natsph)= projected scalars for each atom and ang. mom.
870 : !! sum_1lm_1atom(mlang*mlang,natsph)= projected scalars for each atom and LM component.
871 : !! cplx_1lm_1atom(2,dtset%nspinor**2,dos%mbesslang**2,natsph_tot) = complex projection of wave function on atomic like orbital
872 : !!
873 : !! NOTES
874 : !! * ph3d atoms are ordered with natsph and must be provided by the caller in the correct order!
875 : !!
876 : !! * spinor components are not treated here. This facilitates the implementation of spinor parallelism
877 : !! because the caller can easily call the routine inside a loop over spinors and then sum the
878 : !! different contributions outside the loop thus reducing the number of MPI calls.
879 : !!
880 : !! SOURCE
881 :
882 12852 : subroutine recip_ylm (bess_fit, cg_1band, istwfk, mpi_enreg, nradint, nradintmax, mlang,&
883 4284 : & mpw, natsph, typat_extra, mlang_type, npw_k, nspinor, ph3d, prtsphere, rint, rmax,&
884 4284 : & rc_ylm, sum_1ll_1atom, sum_1lm_1atom, cplx_1lm_1atom, ucvol, ylm_k, znucl_sph)
885 :
886 : !Arguments ------------------------------------
887 : !scalars
888 : integer,intent(in) :: istwfk,mlang,mpw,natsph,npw_k,nradintmax
889 : integer,intent(in) :: nspinor
890 : integer,intent(in) :: prtsphere,rc_ylm
891 : real(dp),intent(in) :: ucvol
892 : !arrays
893 : integer,intent(in) :: nradint(natsph),typat_extra(natsph),mlang_type(:)
894 : real(dp),intent(in) :: bess_fit(mpw,nradintmax,mlang),cg_1band(:,:) !(2,my_nspinor*npw_k)
895 : real(dp),intent(in) :: ph3d(2,npw_k,natsph),rint(nradintmax)
896 : real(dp),intent(in) :: rmax(natsph),ylm_k(npw_k,mlang*mlang)
897 : real(dp),intent(in) :: znucl_sph(natsph)
898 : type(MPI_type),intent(in) :: mpi_enreg
899 : real(dp),intent(out) :: sum_1ll_1atom(nspinor**2, mlang, natsph)
900 : real(dp),intent(out) :: sum_1lm_1atom(nspinor**2, mlang*mlang,natsph)
901 : real(dp),intent(out) :: cplx_1lm_1atom(2,nspinor, mlang*mlang,natsph)
902 :
903 : !Local variables-------------------------------
904 : !scalars
905 : integer :: ilm,iat,ipw,ixint,ll,mm,il,jlm,ierr,lm_size,itypat, ispinor,ipauli, is, isp, my_nspinor
906 : real(dp),parameter :: invsqrt2=one/sqrt2
907 : real(dp) :: sum_all, dr, fact
908 : type(atomdata_t) :: atom
909 : character(len=500) :: msg
910 : !arrays
911 8568 : integer :: ilang(mlang**2)
912 4284 : integer, allocatable :: ispinors(:)
913 8568 : real(dp) :: c1(2),c2(2), sum_1atom(natsph),sum_1ll(mlang),sum_1lm(mlang**2)
914 8568 : real(dp) :: func(nradintmax), func_cplx(nradintmax,2)
915 8568 : complex(dp) :: vect(npw_k)
916 4284 : complex(dp),allocatable :: tmppsia(:,:),tmppsim(:,:),dotc(:), values(:,:,:,:)
917 : ! *************************************************************************
918 :
919 : ! Workspace array (used to reduce the number of MPI communications)
920 : ! One could reduce a bit the memory requirement by using non-blocking operations ...
921 25704 : ABI_MALLOC_OR_DIE(values, (nradintmax, nspinor, mlang**2, natsph), ierr)
922 81628070 : values = czero
923 :
924 4284 : my_nspinor = max(1,nspinor/mpi_enreg%nproc_spinor)
925 17136 : ABI_MALLOC(tmppsia, (npw_k,my_nspinor))
926 12852 : ABI_MALLOC(tmppsim, (npw_k,my_nspinor))
927 12852 : ABI_MALLOC(dotc, (my_nspinor))
928 12852 : ABI_MALLOC(ispinors, (my_nspinor))
929 4284 : if (my_nspinor == 2) then
930 496 : ispinors(1) = 1
931 496 : ispinors(2) = 2
932 : else
933 3788 : ispinors(1) = mpi_enreg%me_spinor+1
934 : end if
935 :
936 519720 : sum_1lm_1atom = zero
937 841520 : cplx_1lm_1atom = zero
938 :
939 25704 : do ll=0,mlang-1
940 132804 : do mm=-ll,ll
941 107100 : ilm = (ll+1)**2-ll+mm
942 128520 : ilang(ilm) = ll+1
943 : end do
944 : end do
945 :
946 : ! Big loop on all atoms
947 10720 : do iat=1,natsph
948 6436 : itypat = typat_extra(iat)
949 6436 : lm_size = mlang_type(itypat) ** 2
950 6436 : dr = rmax(iat) / (nradint(iat)-1)
951 :
952 : ! u(G) e^{i(k+G).Ra}
953 : ! tmppsia = Temporary array for part which depends only on iat
954 13768 : do ispinor=1,my_nspinor
955 1293880 : do ipw=1,npw_k
956 : tmppsia(ipw,ispinor) = dcmplx(cg_1band(1,ipw+(ispinor-1)*npw_k),cg_1band(2,ipw+(ispinor-1)*npw_k)) &
957 1287444 : & * dcmplx(ph3d(1,ipw,iat), ph3d(2,ipw,iat))
958 : end do
959 : end do
960 :
961 : ! tmppsim = temporary arrays for part of psi which does not depend on ixint = tmppsia * ylm.
962 : ! could remove this intermediate array to save memory...
963 : ! u(G) Y_LM^*(k+G) e^{i(k+G).Ra}
964 : ! Take into account the fact that ylm_k are REAL spherical harmonics, see initylmg.f
965 : ! For time-reversal states, detailed treatment show that only the real or imaginary
966 : ! part of tmppsia is needed here, depending on l being even or odd: only one of the coef is 1, the other 0
967 171620 : do ilm=1,lm_size
968 160900 : il = ilang(ilm)
969 160900 : ll = ilang(ilm) - 1
970 160900 : mm = ilm - (ll+1)**2 + ll
971 :
972 160900 : select case (rc_ylm)
973 : case (1)
974 : ! to get PDOS for real spherical harmonics, simply multiply here by ylm_k
975 2750 : do ispinor=1,my_nspinor
976 497025 : do ipw=1,npw_k
977 495650 : tmppsim(ipw,ispinor) = tmppsia(ipw,ispinor) * ylm_k(ipw,ilm)
978 : end do
979 : end do
980 :
981 : ! Handle time-reversal
982 : ! TODO: check if time reversal with spinors is special and may need some treatment.
983 : ! normally SOC will simply impose istwfk 1 in appropriate cases (kptopt 4).
984 1375 : if (istwfk /= 1) then
985 400 : if (mod(ll, 2) == 0) then
986 64485 : tmppsim(:,:) = dcmplx(real(tmppsim(:,:)),zero)
987 : else
988 42990 : tmppsim(:,:) = dcmplx(aimag(tmppsim(:,:)),zero)
989 : end if
990 : ! f2008 version:
991 : ! if (mod(ll, 2) == 0) then
992 : ! tmppsim(:,:)%im = zero
993 : ! else
994 : ! tmppsim(:,:)%re = tmppsim(:,:)%im
995 : ! tmppsim(:,:)%im = zero
996 : ! end if
997 : end if
998 :
999 : case (2)
1000 : ! to get PDOS for complex spherical harmonics, build linear combination of real ylm_k
1001 159525 : jlm = (ll+1)**2-ll-mm ! index of (l, -m)
1002 159525 : if (mm == 0) then
1003 5975930 : vect(:) = dcmplx(ylm_k(1:npw_k,ilm),zero)
1004 127620 : else if (mm > 0) then
1005 : !vect(1,:) = invsqrt2 * ylm_k(1:npw_k,ilm) * (-1)**mm
1006 : !vect(2,:) = +invsqrt2 * ylm_k(1:npw_k,jlm) * (-1)**mm
1007 63810 : c1 = sy(ll, mm, mm)
1008 63810 : c2 = sy(ll,-mm, mm)
1009 : vect(:) = dcmplx(c1(1) * ylm_k(1:npw_k,ilm) + c2(1) * ylm_k(1:npw_k,jlm), &
1010 11951860 : c1(2) * ylm_k(1:npw_k,ilm) + c2(2) * ylm_k(1:npw_k,jlm))
1011 :
1012 : else if (mm < 0) then
1013 : !vect(1,:) = invsqrt2 * ylm_k(1:npw_k,jlm) !* (-1)**mm
1014 : !vect(2,:) = -invsqrt2 * ylm_k(1:npw_k,ilm) !* (-1)**mm
1015 63810 : c1 = sy(ll, mm, mm)
1016 63810 : c2 = sy(ll,-mm, mm)
1017 : vect(:) = dcmplx(c1(1) * ylm_k(1:npw_k,ilm) + c2(1) * ylm_k(1:npw_k,jlm),&
1018 11951860 : c1(2) * ylm_k(1:npw_k,ilm) + c2(2) * ylm_k(1:npw_k,jlm))
1019 : end if
1020 29879650 : vect(:) = dcmplx(real(vect(:)), -aimag(vect(:)))
1021 : !vect(:)%im = -vect(:)%im
1022 :
1023 159525 : if (istwfk == 1) then
1024 334700 : do ispinor=1,my_nspinor
1025 31602950 : do ipw=1,npw_k
1026 31446800 : tmppsim(ipw,ispinor) = tmppsia(ipw,ispinor) * vect(ipw)
1027 : end do
1028 : end do
1029 : else
1030 : ! Handle time-reversal
1031 3375 : if (mod(ll, 2) == 0) then
1032 4050 : do ispinor=1,my_nspinor
1033 148215 : do ipw=1,npw_k
1034 146190 : tmppsim(ipw,ispinor) = real(tmppsia(ipw,ispinor)) * vect(ipw)
1035 : end do
1036 : end do
1037 : else
1038 2700 : do ispinor=1,my_nspinor
1039 98810 : do ipw=1,npw_k
1040 97460 : tmppsim(ipw,ispinor) = aimag(tmppsia(ipw,ispinor)) * vect(ipw)
1041 : end do
1042 : end do
1043 : end if
1044 : end if
1045 :
1046 : case default
1047 160900 : ABI_ERROR("Wrong value for rc_ylm")
1048 : end select
1049 :
1050 : ! Compute integral $ \int_0^{rc} dr r**2 ||\sum_G u(G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r)||**2 $
1051 : ! or more general spinor case integral
1052 : ! $ \int_0^{rc} dr r**2 dotc^*_s \sigma^x_{ss'} dotc_{s'}
1053 : ! where dotc_s = \sum_G u_s (G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r)
1054 64389286 : do ixint=1,nradint(iat)
1055 134103100 : dotc = czero
1056 134103100 : do ispinor=1, my_nspinor
1057 16359986925 : do ipw=1,npw_k
1058 16295764975 : dotc(ispinor) = dotc(ispinor) + bess_fit(ipw, ixint, il) * tmppsim(ipw, ispinor)
1059 : end do
1060 : end do
1061 64221950 : if (istwfk /= 1) then
1062 2813950 : dotc = two * dotc
1063 1406975 : if (istwfk == 2 .and. mpi_enreg%me_g0 == 1) then
1064 876700 : dotc(:) = dotc(:) - bess_fit(1, ixint, il) * tmppsim(1, :)
1065 : end if
1066 : end if
1067 :
1068 : ! Store results to reduce number of xmpi_sum calls if MPI
1069 134264000 : do ispinor=1, my_nspinor
1070 134103100 : values(ixint, ispinors(ispinor), ilm, iat) = dotc(ispinor)
1071 : end do
1072 : end do ! ixint
1073 :
1074 : end do ! ilm
1075 : end do ! iat
1076 :
1077 4284 : ABI_FREE(tmppsia)
1078 4284 : ABI_FREE(tmppsim)
1079 4284 : ABI_FREE(dotc)
1080 4284 : ABI_FREE(ispinors)
1081 :
1082 : ! Collect results in comm_pw (data are distributed over plane waves)
1083 4284 : call xmpi_sum(values, mpi_enreg%comm_bandfft, ierr)
1084 : ! ! Collect results in mpi_enreg%comm_spinor (data are distributed over spinor components)
1085 4284 : call xmpi_sum(values, mpi_enreg%comm_spinor, ierr)
1086 :
1087 : ! Multiply by r**2 and take norm, integrate
1088 10720 : do iat=1,natsph
1089 6436 : itypat = typat_extra(iat)
1090 6436 : lm_size = mlang_type(itypat) ** 2
1091 171620 : do ilm=1,lm_size
1092 :
1093 509000 : do ipauli=0,nspinor**2-1
1094 111067650 : do ixint=1,nradint(iat)
1095 110719550 : func(ixint) = zero
1096 283435900 : do is=1,nspinor
1097 580145850 : do isp=1,nspinor
1098 : func(ixint) = func(ixint) + real(conjg(values(ixint, is, ilm, iat)) * pauli_mat(is,isp,ipauli)*&
1099 469426300 : & values(ixint, isp, ilm, iat))
1100 : end do
1101 : end do
1102 111067650 : func(ixint) = rint(ixint)**2 * func(ixint)
1103 : end do
1104 : ! Here I should treat the case in which the last point /= rcut
1105 : ! NB: indexing is from 1 not 0 for spin matrix components
1106 509000 : sum_1lm_1atom (ipauli+1, ilm, iat) = simpson(dr, func(1:nradint(iat)))
1107 : end do ! ipauli
1108 :
1109 390636 : do is = 1, nspinor
1110 81456450 : func_cplx(:,1) = real(values(:, is, ilm, iat))
1111 81456450 : func_cplx(:,2) = aimag(values(:, is, ilm, iat))
1112 79944450 : do ixint=1,nradint(iat)
1113 239386750 : func_cplx(ixint,:) = rint(ixint)**2 * func_cplx(ixint,:)
1114 : end do
1115 223300 : cplx_1lm_1atom(1, is, ilm, iat) = simpson(dr, func_cplx(1:nradint(iat),1))
1116 384200 : cplx_1lm_1atom(2, is, ilm, iat) = simpson(dr, func_cplx(1:nradint(iat),2))
1117 : end do
1118 :
1119 : end do ! ilm
1120 : end do ! iat
1121 :
1122 : ! Normalize with unit cell volume and include 4pi term coming from Rayleigh expansion.
1123 4284 : fact = four_pi**2 / ucvol
1124 519720 : sum_1lm_1atom = fact * sum_1lm_1atom
1125 841520 : cplx_1lm_1atom = fact * cplx_1lm_1atom
1126 :
1127 : ! sum up the m-independent fractions
1128 112520 : sum_1ll_1atom = zero
1129 10720 : do iat=1,natsph
1130 6436 : itypat = typat_extra(iat)
1131 6436 : lm_size = mlang_type(itypat) ** 2
1132 171620 : do ilm=1,lm_size
1133 160900 : il = ilang(ilm)
1134 515436 : sum_1ll_1atom(:,il, iat) = sum_1ll_1atom(:,il, iat) + sum_1lm_1atom(:,ilm, iat)
1135 : end do
1136 : end do
1137 :
1138 4284 : ABI_FREE(values)
1139 :
1140 : ! Output
1141 8568 : if (prtsphere == 1) then
1142 240 : sum_1ll = zero
1143 1040 : sum_1lm = zero
1144 95 : sum_1atom = zero
1145 95 : do iat=1,natsph
1146 1430 : sum_1atom(iat) = sum(sum_1lm_1atom(1,:,iat))
1147 330 : sum_1ll(:)=sum_1ll(:)+sum_1ll_1atom(1,:,iat)
1148 1470 : sum_1lm(:)=sum_1lm(:)+sum_1lm_1atom(1,:,iat)
1149 : end do
1150 95 : sum_all = sum(sum_1atom)
1151 :
1152 40 : if (rc_ylm == 1) msg = " Angular analysis (real spherical harmonics)"
1153 40 : if (rc_ylm == 2) msg = " Angular analysis (complex spherical harmonics)"
1154 40 : call wrtout(std_out, msg)
1155 95 : do iat=1,natsph
1156 55 : call atomdata_from_znucl(atom, znucl_sph(iat))
1157 55 : call wrtout(std_out, " ")
1158 55 : write(msg,'(a,i3,a,a,a,f10.6)' )' Atom # ',iat, ' is ', atom%symbol,', in-sphere charge =',sum_1atom(iat)
1159 55 : call wrtout(std_out, msg)
1160 425 : do ll=0,mlang-1
1161 : write(msg,'(a,i1,a,f9.6,a,9f6.3)' )&
1162 275 : ' l=',ll,', charge=',sum_1ll_1atom(1,ll+1,iat),&
1163 550 : ', m=-l,l splitting:',sum_1lm_1atom(1,1+ll**2:(ll+1)**2,iat)
1164 330 : call wrtout(std_out, msg)
1165 : end do ! ll
1166 : end do ! iat
1167 40 : write(msg,'(a,a)') ch10,' Sum of angular contributions for all atomic spheres '
1168 40 : call wrtout(std_out, msg)
1169 240 : do ll=0,mlang-1
1170 : write(msg,'(a,i1,a,f9.6,a,f9.6)' )&
1171 200 : ' l=',ll,', charge =',sum_1ll(ll+1),' proportion =',sum_1ll(ll+1)/sum_all
1172 240 : call wrtout(std_out, msg)
1173 : end do
1174 40 : write(msg,'(a,a,f10.6)' ) ch10,' Total over all atoms and l=0 to 4 :',sum_all
1175 40 : call wrtout(std_out, msg)
1176 40 : call wrtout(std_out, " ")
1177 : end if
1178 :
1179 : contains
1180 :
1181 255240 : function sy(ll, mm, mp)
1182 : use m_paw_sphharm, only : ys
1183 : ! Computes the matrix element <Slm|Ylm'>
1184 : integer,intent(in) :: ll,mm, mp
1185 :
1186 : real(dp) :: sy(2)
1187 : complex(dp) :: ys_val
1188 :
1189 : ! Computes the matrix element <Yl'm'|Slm>
1190 255240 : call ys(ll,mp,ll,mm,ys_val)
1191 : !call ys(ll,mm,ll,mp,ys_val)
1192 255240 : sy(1) = real(ys_val)
1193 255240 : sy(2) = -aimag(ys_val)
1194 :
1195 : end function sy
1196 :
1197 : end subroutine recip_ylm
1198 : !!***
1199 :
1200 : !!****f* m_epjdos/dens_in_sph
1201 : !! NAME
1202 : !! dens_in_sph
1203 : !!
1204 : !! FUNCTION
1205 : !! Calculate integrated density in sphere around each atom
1206 : !!
1207 : !! INPUTS
1208 : !! cg = wavefunction coefficitents in recip space
1209 : !! gmet = metric in recip space
1210 : !! istwfk = storage mode for cg coefficients
1211 : !! kg_k = G vector indices
1212 : !! natom = number of atoms
1213 : !! mpi_enreg=information about MPI parallelization
1214 : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1215 : !! npw_k = number of plane waves for this kpoint
1216 : !! ph1d = phase factors for different atoms for all G vectors
1217 : !! rmax(natom) = max radius to integrate to (in bohr)
1218 : !!
1219 : !! OUTPUT
1220 : !! cmax = integrated density for each atom for a rmax-radius sphere
1221 : !!
1222 : !! WARNING
1223 : !! cg should not be modified by fourwf.
1224 : !!
1225 : !! SOURCE
1226 :
1227 40 : subroutine dens_in_sph(cmax,cg,gmet,istwfk,kg_k,natom,ngfft,mpi_enreg,npw_k,&
1228 40 : & ph1d,rmax,ucvol)
1229 :
1230 : !Arguments ------------------------------------
1231 : !scalars
1232 : integer,intent(in) :: istwfk,natom,npw_k
1233 : real(dp),intent(in) :: ucvol
1234 : type(MPI_type),intent(in) :: mpi_enreg
1235 : !arrays
1236 : integer,intent(in) :: kg_k(3,npw_k),ngfft(18)
1237 : real(dp),intent(in) :: gmet(3,3)
1238 : real(dp),intent(in) :: ph1d(2,(2*ngfft(1)+1+2*ngfft(2)+1+2*ngfft(3)+1)*natom)
1239 : real(dp),intent(in) :: rmax(natom)
1240 : real(dp),intent(inout) :: cg(2,npw_k)
1241 : real(dp),intent(out) :: cmax(natom)
1242 :
1243 : !Local variables -------------------------
1244 : !scalars
1245 : integer,parameter :: tim_fourwf=0
1246 : integer :: cplex,i1,i2,i3,iatom,id1,id2,id3,ifft,mgfft,n1,n2,n3,n4,n5,n6,nfft,nfftot
1247 : real(dp) :: cmaxr,g1,g2,g3,norm,weight
1248 : !arrays
1249 : integer :: ngfft_here(18)
1250 40 : integer,allocatable :: garr(:,:),gbound(:,:)
1251 40 : real(dp),allocatable :: denpot(:,:,:),fofgout(:,:),fofr(:,:,:,:),gnorm(:)
1252 40 : real(dp),allocatable :: ph3d(:,:,:),phkxred(:,:),rhog(:,:),rhor(:), sphrhog(:,:)
1253 : ! *********************************************************************
1254 :
1255 40 : n1=ngfft(1)
1256 40 : n2=ngfft(2)
1257 40 : n3=ngfft(3)
1258 40 : n4=ngfft(4)
1259 40 : n5=ngfft(5)
1260 40 : n6=ngfft(6)
1261 40 : nfftot = n1*n2*n3
1262 40 : nfft=n1*n2*n3
1263 40 : ngfft_here(:) = ngfft(:)
1264 : !fourwf doesnt work with other options for mode 0 (fft G -> r)
1265 40 : ngfft_here(7)=111
1266 40 : ngfft_here(8)=256
1267 160 : mgfft=maxval(ngfft_here(1:3))
1268 :
1269 40 : call sqnorm_g(norm,istwfk,npw_k,cg,mpi_enreg%me_g0,mpi_enreg%comm_fft)
1270 :
1271 40 : if (abs(one-norm) > tol6) then
1272 0 : write(std_out,'(a,f8.5)' ) ' dens_in_sph : this state is not normalized : norm=',norm
1273 : end if
1274 :
1275 : !-----------------------------------------------------------------
1276 : !inverse FFT of wavefunction to real space => density in real space
1277 : !-----------------------------------------------------------------
1278 160 : ABI_MALLOC(gbound,(2*mgfft+8,2))
1279 40 : call sphereboundary(gbound,istwfk,kg_k,mgfft,npw_k)
1280 :
1281 40 : weight = one
1282 40 : cplex=1
1283 200 : ABI_MALLOC(denpot,(cplex*n4,n5,n6))
1284 500360 : denpot(:,:,:)=zero
1285 120 : ABI_MALLOC(fofgout,(2,npw_k))
1286 200 : ABI_MALLOC(fofr,(2,n4,n5,n6))
1287 :
1288 : call fourwf(cplex,denpot,cg,fofgout,fofr,gbound,gbound, &
1289 : istwfk,kg_k,kg_k,mgfft,mpi_enreg,1,ngfft_here,npw_k,&
1290 40 : npw_k,n4,n5,n6,1,tim_fourwf,weight,weight)
1291 :
1292 40 : ABI_FREE(fofgout)
1293 40 : ABI_FREE(fofr)
1294 40 : ABI_FREE(gbound)
1295 :
1296 500360 : norm = sum(denpot(:,:,:))/nfftot
1297 40 : if (abs(one-norm) > tol6) then
1298 0 : write(std_out,'(a,f8.5)') ' dens_in_sph : this state is not normalized in real space : norm=',norm
1299 : end if
1300 :
1301 : !-----------------------------------------------------------------
1302 : !FFT of new density: we obtain n(G) in rhog(1,:)
1303 : !-----------------------------------------------------------------
1304 :
1305 : !Change the packing of the reciprocal space density
1306 120 : ABI_MALLOC(rhor,(nfft))
1307 40 : call fftpac(1,mpi_enreg,1,n1,n2,n3,n4,n5,n6,ngfft,rhor,denpot,1)
1308 :
1309 120 : ABI_MALLOC(rhog,(2,nfft))
1310 40 : call fourdp(1,rhog,rhor,-1,mpi_enreg,nfft,1,ngfft,0)
1311 :
1312 40 : ABI_FREE(rhor)
1313 40 : ABI_FREE(denpot)
1314 :
1315 439208 : do ifft=1,nfft
1316 1317544 : rhog(:,ifft) = rhog(:,ifft) / ucvol
1317 : end do
1318 :
1319 : !-----------------------------------------------------------------
1320 : !calculate norms of G vectors
1321 : !-----------------------------------------------------------------
1322 :
1323 120 : ABI_MALLOC(garr,(3,nfft))
1324 80 : ABI_MALLOC(gnorm,(nfft))
1325 40 : id3=ngfft(3)/2+2 ; id2=ngfft(2)/2+2 ; id1=ngfft(1)/2+2
1326 920 : do i3=1,n3
1327 880 : g3=i3-(i3/(id3))*ngfft(3)-one
1328 20472 : do i2=1,n2
1329 19552 : g2=i2-(i2/(id2))*ngfft(2)-one
1330 459600 : do i1=1,n1
1331 439168 : g1=i1-(i1/(id1))*ngfft(1)-one
1332 439168 : ifft=i1+(i2-1)*n1+(i3-1)*n1*n2
1333 439168 : garr(1,ifft)=nint(g1)
1334 439168 : garr(2,ifft)=nint(g2)
1335 439168 : garr(3,ifft)=nint(g3)
1336 : gnorm(ifft)=sqrt(gmet(1,1)*g1*g1 + &
1337 : & two*gmet(2,1)*g2*g1 + &
1338 : & two*gmet(3,1)*g3*g1 + &
1339 : & gmet(2,2)*g2*g2 + &
1340 : & gmet(3,2)*g3*g2 + &
1341 458720 : & gmet(3,3)*g3*g3)
1342 : end do
1343 : end do
1344 : end do
1345 :
1346 : !-----------------------------------------------------------------
1347 : !For each atom call sphericaldens to calculate
1348 : !n(G) * 1/|G|^3 * int_0^2*\pi*r_{max}*|G| 4 \pi y^2 j_0 (y) dy
1349 : !for all G vectors put into array sphrhog
1350 : !scalar product of phase factors with spherically convoluted density
1351 : !-----------------------------------------------------------------
1352 :
1353 : !largest mem occupation = nfft * (2(sphrog) +2*1(ph3d) +3(garr) +2(rhog) +1(gnorm)) = nfft * 10
1354 80 : ABI_MALLOC(sphrhog,(2,nfft))
1355 120 : ABI_MALLOC(phkxred,(2,natom))
1356 95 : phkxred(1,:)=one
1357 95 : phkxred(2,:)=zero
1358 80 : ABI_MALLOC(ph3d,(2,nfft,1))
1359 :
1360 95 : do iatom=1,natom
1361 :
1362 55 : call sphericaldens(rhog,gnorm,nfft,rmax(iatom),sphrhog)
1363 : ! -----------------------------------------------------------------
1364 : ! Compute the phases for the whole set of fft vectors
1365 : ! -----------------------------------------------------------------
1366 :
1367 55 : call ph1d3d(iatom,iatom,garr,natom,natom,nfft,ngfft(1),ngfft(2),ngfft(3),phkxred,ph1d,ph3d)
1368 :
1369 : ! For the phase factors, take the compex conjugate, before evaluating the scalar product
1370 542519 : do ifft=1,nfft
1371 542519 : ph3d(2,ifft,1)=-ph3d(2,ifft,1)
1372 : end do
1373 55 : cplex=2
1374 55 : call dotprod_v(cplex,cmaxr,nfft,1,0,ph3d,sphrhog,mpi_enreg%comm_fft)
1375 95 : cmax(iatom) = cmaxr
1376 : ! write(std_out,'(a,i4,a,es14.6,a,es12.6)' )' dens_in_sph : At ', iatom, ' has ',cmaxr, ' el.s in a sphere of rad ', rmax
1377 : end do
1378 :
1379 40 : ABI_FREE(rhog)
1380 40 : ABI_FREE(gnorm)
1381 40 : ABI_FREE(garr)
1382 40 : ABI_FREE(sphrhog)
1383 40 : ABI_FREE(ph3d)
1384 40 : ABI_FREE(phkxred)
1385 :
1386 40 : end subroutine dens_in_sph
1387 : !!***
1388 :
1389 : !!****f* m_epjdos/sphericaldens
1390 : !! NAME
1391 : !! sphericaldens
1392 : !!
1393 : !! FUNCTION
1394 : !! Compute the convolution of a function with
1395 : !! the unity constant function over a sphere of radius rmax .
1396 : !! The function is to be given in reciprocal space,
1397 : !! the resulting function is also given in reciprocal space.
1398 : !! The routine needs the norm of the reciprocal space vectors.
1399 : !!
1400 : !! The resulting function in reciprocal space can give the
1401 : !! integral of the density in any sphere of that radius, centered
1402 : !! on any point, by a simple scalar product.
1403 : !!
1404 : !! INPUTS
1405 : !! fofg(2,nfft)=initial function, in reciprocal space
1406 : !! gnorm(nfft)=norm of the reciprocal space vectors
1407 : !! nfft=(effective) number of FFT grid points (for this processor)
1408 : !! rmax=radius of the sphere
1409 : !!
1410 : !! OUTPUT
1411 : !! sphfofg(2,nfft)=convoluted function, in reciprocal space
1412 : !!
1413 : !! SOURCE
1414 :
1415 55 : subroutine sphericaldens(fofg,gnorm,nfft,rmax,sphfofg)
1416 :
1417 : !Arguments ------------------------------------
1418 : !scalars
1419 : integer,intent(in) :: nfft
1420 : real(dp),intent(in) :: rmax
1421 : !arrays
1422 : real(dp),intent(in) :: fofg(2,nfft),gnorm(nfft)
1423 : real(dp),intent(out) :: sphfofg(2,nfft)
1424 :
1425 : !Local variables-------------------------------
1426 : !scalars
1427 : integer :: ifft
1428 : real(dp) :: factor,int0yy,rmax_2pi,yy
1429 : ! *************************************************************************
1430 :
1431 55 : rmax_2pi=two_pi*rmax
1432 55 : factor=four_pi/(two_pi)**3
1433 :
1434 542519 : do ifft=1,nfft
1435 542519 : if(abs(gnorm(ifft)) < tol12)then
1436 55 : sphfofg(1,ifft)=fofg(1,ifft)*four_pi*third*rmax**3
1437 55 : sphfofg(2,ifft)=fofg(2,ifft)*four_pi*third*rmax**3
1438 : else
1439 542409 : yy=gnorm(ifft)*rmax_2pi
1440 542409 : int0yy=factor*(sin(yy)-yy*cos(yy))/(gnorm(ifft)**3)
1441 542409 : sphfofg(1,ifft)=fofg(1,ifft)*int0yy
1442 542409 : sphfofg(2,ifft)=fofg(2,ifft)*int0yy
1443 : end if
1444 : end do
1445 :
1446 55 : end subroutine sphericaldens
1447 : !!***
1448 :
1449 : !!****f* m_epjdos/epjdos_prtfatbands
1450 : !! NAME
1451 : !! epjdos_prtfatbands
1452 : !!
1453 : !! FUNCTION
1454 : !! Print dos_fractions_m in order to plot easily fatbands
1455 : !! if pawfatbnd=1 1 : fatbands are resolved in L.
1456 : !! if pawfatbnd=1 2 : fatbands are resolved in L and M.
1457 : !!
1458 : !! INPUTS
1459 : !! dos_fractions_m(nkpt,mband,nsppol,ndosfraction*mbesslang*m_dos_flag)
1460 : !! = m-resolved projected dos inside PAW sphere.
1461 : !! dtset = Input variables
1462 : !! ebands<ebands_t>=Band structure data.
1463 : !! pawfatbnd = keyword for fatbands
1464 : !! mbesslang =maximum angular momentum for Bessel function expansion
1465 : !! m_dos_flag =option for the m-contributions to the partial DOS
1466 : !! ndosfraction =natsph*mbesslang
1467 : !!
1468 : !! OUTPUT
1469 : !! (only writing)
1470 : !!
1471 : !! NOTES
1472 : !! This routine should be called by master only
1473 : !!
1474 : !! SOURCE
1475 :
1476 1 : subroutine epjdos_prtfatbands(dos,dtset,ebands,fildata,pawfatbnd,pawtab)
1477 :
1478 : !Arguments ------------------------------------
1479 : !scalars
1480 : class(epjdos_t),intent(in) :: dos
1481 : integer,intent(in) :: pawfatbnd
1482 : type(ebands_t),intent(in) :: ebands
1483 : type(dataset_type),intent(in) :: dtset
1484 : character(len=fnlen),intent(in) :: fildata
1485 : !arrays
1486 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat)
1487 :
1488 : !Local variables-------------------------------
1489 : !scalars
1490 : integer :: iall,il,iat,natsph,inbfatbands,iband,mband,ixfat,isppol,nkpt,lmax,ll,mm
1491 : integer :: ikpt,nband_k,ndosfraction,mbesslang
1492 : real(dp) :: xfatband,cpu,wall,gflops
1493 : character(len=1) :: tag_l,tag_1m,tag_is
1494 : character(len=2) :: tag_2m
1495 : character(len=10) :: tag_il,tag_at,tag_grace
1496 : character(len=1500) :: msg
1497 : character(len=fnlen) :: tmpfil
1498 : type(atomdata_t) :: atom
1499 : !arrays
1500 : integer :: units(2)
1501 1 : integer,allocatable :: unitfatbands_arr(:,:)
1502 1 : real(dp),allocatable :: eigenvalues(:,:,:)
1503 : character(len=2) :: symbol
1504 : !*************************************************************************
1505 :
1506 : DBG_ENTER("COLL")
1507 :
1508 3 : units = [std_out, ab_out]
1509 1 : ndosfraction = dos%ndosfraction; mbesslang = dos%mbesslang
1510 :
1511 1 : if(dos%prtdosm.ne.0) then
1512 : write(msg,'(3a)')&
1513 0 : 'm decomposed dos is activated',ch10, &
1514 0 : 'Action: deactivate it with prtdosm=0 !'
1515 0 : ABI_ERROR(msg)
1516 : end if
1517 :
1518 1 : if(dtset%nspinor==2) then
1519 0 : ABI_WARNING("Fatbands are not yet available in the case nspinor==2!")
1520 : end if
1521 :
1522 1 : ABI_CHECK(allocated(dos%fractions_m), "dos%fractions_m is not allocated!")
1523 :
1524 1 : natsph=dtset%natsph
1525 1 : nkpt=dtset%nkpt
1526 1 : mband=dtset%mband
1527 :
1528 1 : if(natsph>1000) then
1529 : write(msg,'(3a)')&
1530 0 : 'Too big number of fat bands!',ch10, &
1531 0 : 'Action: decrease natsph in input file !'
1532 0 : ABI_ERROR(msg)
1533 : end if
1534 :
1535 : !-------------- PRINTING IN LOG
1536 1 : call cwtime(cpu, wall, gflops, "start")
1537 1 : write(msg,'(a,a,a,a,i5,a,a,1000i5)') ch10," ***** Print of fatbands activated ****** ",ch10,&
1538 1 : " Number of atom: natsph = ",natsph,ch10, &
1539 6 : " atoms are = ",(dtset%iatsph(iat),iat=1,natsph)
1540 1 : call wrtout(units, msg)
1541 1 : iall=0;inbfatbands=0
1542 :
1543 1 : if(pawfatbnd==1) then
1544 0 : inbfatbands=mbesslang-1
1545 0 : write(msg,'(3a)')" (fatbands are in eV and are given for each value of L)",ch10
1546 1 : else if(pawfatbnd==2) then
1547 1 : write(msg,'(3a)')" (fatbands are in eV and are given for each value of L and M)",ch10
1548 1 : inbfatbands=(mbesslang-1)**2
1549 : end if
1550 1 : call wrtout(units, msg)
1551 :
1552 1 : write(msg,'(a,e12.5,a,e12.5,a)') " Fermi energy is ",ebands%fermie*Ha_eV," eV = ",ebands%fermie," Ha"
1553 1 : call wrtout(std_out,msg)
1554 :
1555 : !-------------- OPEN AND NAME FILES FOR FATBANDS
1556 4 : ABI_MALLOC(unitfatbands_arr,(natsph*inbfatbands,dtset%nsppol))
1557 66 : unitfatbands_arr = -3
1558 :
1559 5 : do iat=1,natsph
1560 4 : lmax=(pawtab(dtset%typat(dtset%iatsph(iat)))%l_size-1)/2
1561 4 : call int2char4(dtset%iatsph(iat),tag_at)
1562 4 : ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!')
1563 4 : call atomdata_from_znucl(atom,dtset%znucl(dtset%typat(dtset%iatsph(iat))))
1564 4 : symbol = atom%symbol
1565 69 : do il=1,inbfatbands
1566 64 : iall=iall+1
1567 64 : ll=int(sqrt(float(il-1))) ! compute l
1568 68 : if(ll.le.lmax) then ! print only angular momentum included in the PAW data
1569 72 : do isppol=1,dtset%nsppol
1570 36 : write(tag_is,'(i1)')isppol
1571 36 : if(pawfatbnd==1) then
1572 0 : call int2char4(il-1,tag_il)
1573 0 : ABI_CHECK((tag_il(1:1)/='#'),'Bug: string length too short!')
1574 0 : tmpfil = trim(fildata)//'_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//trim(tag_il)
1575 36 : else if (pawfatbnd==2) then
1576 36 : write(tag_l,'(i1)') ll
1577 36 : mm=il-(ll**2+ll+1) ! compute m
1578 36 : if(mm<0) write(tag_2m,'(i2)') mm
1579 36 : if(mm>=0) write(tag_1m,'(i1)') mm
1580 36 : if(mm<0) tmpfil = trim(fildata)// &
1581 12 : & '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//tag_l//'_m'//tag_2m
1582 36 : if(mm>=0) tmpfil = trim(fildata)// &
1583 24 : & '_at'//trim(tag_at)//'_'//trim(adjustl(symbol))//'_is'//tag_is//'_l'//tag_l//'_m+'//tag_1m
1584 : end if
1585 : !unitfatbands_arr(iall,isppol)=tmp_unit+100+iall-1+(natsph*inbfatbands)*(isppol-1)
1586 : !open (unit=unitfatbands_arr(iall,isppol),file=trim(tmpfil),status='unknown',form='formatted')
1587 36 : if (open_file(tmpfil, msg, newunit=unitfatbands_arr(iall,isppol), status='unknown',form='formatted') /= 0) then
1588 0 : ABI_ERROR(msg)
1589 : end if
1590 :
1591 36 : write(msg,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitfatbands_arr(iall,isppol)
1592 36 : call wrtout(std_out,msg)
1593 36 : write(msg,'(9a)') "# ",ch10,"# ABINIT package : FATBAND file ", ch10,&
1594 36 : "# It contains, for each band: the eigenvalues in eV (and the character of the band) as a function of the k-point",&
1595 72 : ch10,"# This file can be read with xmgrace (http://plasma-gate.weizmann.ac.il/Grace/) ",ch10,"# "
1596 36 : write(unitfatbands_arr(iall,isppol), "(a)")trim(msg)
1597 1116 : do iband=1,mband
1598 1080 : call int2char4(iband-1,tag_grace)
1599 1080 : ABI_CHECK((tag_grace(1:1)/='#'),'Bug: string length too short!')
1600 1080 : write(msg,'(16a)') ch10,"@ s",trim(tag_grace)," line color 1",&
1601 1080 : ch10,"@ s",trim(tag_grace)," errorbar color 2",&
1602 1080 : ch10,"@ s",trim(tag_grace)," errorbar riser linewidth 5.0", &
1603 2160 : ch10,"@ s",trim(tag_grace)," errorbar linestyle 0"
1604 1116 : write(unitfatbands_arr(iall,isppol), "(a)")trim(msg)
1605 : end do !iband
1606 72 : write(unitfatbands_arr(iall,isppol), '(a,a)') ch10,'@type xydy'
1607 : end do ! isppol
1608 : end if ! ll=<lmax
1609 : end do ! il
1610 : end do ! iat
1611 :
1612 1 : if(iall.ne.(natsph*inbfatbands)) then
1613 0 : ABI_ERROR("error1 ")
1614 : end if
1615 :
1616 : !-------------- WRITE FATBANDS IN FILES
1617 1 : if (pawfatbnd>0) then
1618 : ! Store eigenvalues with nkpt as first dimension for efficiency reasons
1619 5 : ABI_MALLOC(eigenvalues,(nkpt,mband,dtset%nsppol))
1620 2 : do isppol=1,dtset%nsppol
1621 18 : do ikpt=1,nkpt
1622 16 : nband_k=dtset%nband(ikpt+(isppol-1)*nkpt)
1623 497 : do iband=1,mband
1624 496 : eigenvalues(ikpt,iband,isppol)= ebands%eig(iband, ikpt, isppol) - ebands%fermie
1625 : end do
1626 : end do
1627 : end do
1628 1 : iall=0
1629 5 : do iat=1,natsph
1630 4 : lmax=(pawtab(dtset%typat(dtset%iatsph(iat)))%l_size-1)/2
1631 69 : do il=1,inbfatbands
1632 64 : iall=iall+1
1633 64 : ll=int(sqrt(float(il-1)))
1634 68 : if(ll.le.lmax) then
1635 72 : do isppol=1,dtset%nsppol
1636 1116 : do iband=1,mband
1637 1080 : write(unitfatbands_arr(iall,isppol),'(a,a,i8)') ch10,"# BAND number :",iband
1638 18396 : do ikpt=1,nkpt
1639 17280 : if(pawfatbnd==1) then
1640 0 : xfatband=0.d0
1641 0 : do ixfat=(il-1)**2+1,il**2
1642 0 : xfatband=xfatband+dos%fractions_m(ikpt,iband,isppol,(iat-1)*mbesslang**2+ixfat)
1643 : end do ! ixfat
1644 17280 : else if (pawfatbnd==2) then
1645 17280 : xfatband=dos%fractions_m(ikpt,iband,isppol,(iat-1)*mbesslang**2+il)
1646 : end if
1647 : write(unitfatbands_arr(iall,isppol),'(i5,e20.5,e20.5)')&
1648 18360 : ikpt-1,eigenvalues(ikpt,iband,isppol)*Ha_eV,xfatband
1649 : end do ! ikpt
1650 : end do !iband
1651 72 : write(unitfatbands_arr(iall,isppol),'(a)') '&'
1652 : !close(unitfatbands_arr(iall,isppol))
1653 : end do !isppol
1654 : end if
1655 : end do ! il
1656 : end do ! iat
1657 1 : ABI_FREE(eigenvalues)
1658 : end if
1659 :
1660 2 : do isppol=1,size(unitfatbands_arr, dim=2)
1661 66 : do iat=1,size(unitfatbands_arr, dim=1)
1662 65 : if (unitfatbands_arr(iat, isppol) /= -3) close (unitfatbands_arr(iat, isppol))
1663 : end do
1664 : end do
1665 :
1666 1 : ABI_FREE(unitfatbands_arr)
1667 :
1668 1 : call cwtime(cpu,wall,gflops,"stop")
1669 1 : write(msg,'(2(a,f8.2),a)')" epjdos_prtfatbands: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
1670 1 : call wrtout(std_out,msg,"PERS")
1671 :
1672 : DBG_EXIT("COLL")
1673 :
1674 2 : end subroutine epjdos_prtfatbands
1675 : !!***
1676 :
1677 : !----------------------------------------------------------------------
1678 :
1679 : !!****f* m_epjdos/epjdos_ncwrite
1680 : !! NAME
1681 : !! epjdos_ncwrite
1682 : !!
1683 : !! FUNCTION
1684 : !! Write PJDOS contributions to netcdf file.
1685 : !!
1686 : !! INPUTS
1687 : !! crystal<crystal_t>=Object defining the unit cell and its symmetries.
1688 : !! ebands<ebands_t>=Band structure data.
1689 : !! hdr<hdr_t>=Abinit header
1690 : !! dtset<dtset_type>=Dataset type
1691 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
1692 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
1693 : !! ncid=NC file handle.
1694 : !!
1695 : !! OUTPUT
1696 : !! Only writing
1697 : !!
1698 : !! SOURCE
1699 :
1700 28 : subroutine epjdos_ncwrite(dos, crystal, ebands, hdr, dtset, psps, pawtab, ncid)
1701 :
1702 : !Arguments ------------------------------------
1703 : !scalars
1704 : class(epjdos_t),intent(in) :: dos
1705 : integer,intent(in) :: ncid
1706 : type(crystal_t),intent(in) :: crystal
1707 : type(ebands_t),intent(in) :: ebands
1708 : type(hdr_type),intent(in) :: hdr
1709 : type(dataset_type),intent(in) :: dtset
1710 : type(pseudopotential_type),intent(in) :: psps
1711 : !arrays
1712 : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*dtset%usepaw)
1713 :
1714 : !Local variables-------------------------------
1715 : !scalars
1716 : integer :: itype,ncerr,fform
1717 : real(dp) :: cpu,wall,gflops
1718 : character(len=500) :: msg
1719 : !arrays
1720 56 : integer :: lmax_type(crystal%ntypat)
1721 : !*************************************************************************
1722 :
1723 28 : ABI_CHECK(dtset%natsph > 0, "natsph <=0")
1724 28 : call cwtime(cpu, wall, gflops, "start")
1725 :
1726 28 : fform = fform_from_ext("FATBANDS.nc")
1727 28 : ABI_CHECK(fform /= 0, "Cannot find fform associated to FATBANDS.nc")
1728 :
1729 : ! Write header, crystal structure and band energies.
1730 28 : NCF_CHECK(hdr%ncwrite(ncid, fform, nc_define=.True.))
1731 28 : NCF_CHECK(crystal%ncwrite(ncid))
1732 28 : NCF_CHECK(ebands%ncwrite(ncid))
1733 :
1734 : ! Add fatband-specific quantities
1735 : ncerr = nctk_def_dims(ncid, [ &
1736 : nctkdim_t("natsph", dtset%natsph), &
1737 84 : nctkdim_t("ndosfraction", dos%ndosfraction)], defmode=.True.)
1738 28 : NCF_CHECK(ncerr)
1739 :
1740 28 : if (dos%ndosfraction*dos%mbesslang > 0) then
1741 : ncerr = nctk_def_dims(ncid, [ &
1742 : nctkdim_t("mbesslang", dos%mbesslang), &
1743 72 : nctkdim_t("dos_fractions_m_lastsize", dos%ndosfraction*dos%mbesslang)])
1744 24 : NCF_CHECK(ncerr)
1745 : end if
1746 28 : if (dtset%natsph_extra /= 0) then
1747 2 : NCF_CHECK(nctk_def_dims(ncid, [nctkdim_t("natsph_extra", dtset%natsph_extra)]))
1748 : end if
1749 :
1750 112 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "prtdos", "pawprtdos", "prtdosm"])
1751 28 : NCF_CHECK(ncerr)
1752 56 : ncerr = nctk_def_dpscalars(ncid, [character(len=nctk_slen) :: "ratsph_extra"])
1753 28 : NCF_CHECK(ncerr)
1754 :
1755 : ncerr = nctk_def_arrays(ncid, [&
1756 : nctkarr_t("lmax_type", "int", "number_of_atom_species"), &
1757 : nctkarr_t("iatsph", "int", "natsph"), &
1758 : nctkarr_t("ratsph", "dp", "number_of_atom_species"), &
1759 : nctkarr_t("dos_fractions", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction") &
1760 140 : ])
1761 28 : NCF_CHECK(ncerr)
1762 :
1763 28 : if (allocated(dos%fractions_m)) then
1764 : ncerr = nctk_def_arrays(ncid, &
1765 : nctkarr_t("dos_fractions_m", "dp", &
1766 3 : "number_of_kpoints, max_number_of_states, number_of_spins, dos_fractions_m_lastsize"))
1767 3 : NCF_CHECK(ncerr)
1768 : end if
1769 :
1770 28 : if (allocated(dos%fractions_paw1)) then
1771 : ncerr = nctk_def_arrays(ncid, [&
1772 : nctkarr_t("dos_fractions_paw1", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction"), &
1773 : nctkarr_t("dos_fractions_pawt1", "dp", "number_of_kpoints, max_number_of_states, number_of_spins, ndosfraction") &
1774 24 : ])
1775 8 : NCF_CHECK(ncerr)
1776 : end if
1777 :
1778 28 : if (dtset%natsph_extra /= 0) then
1779 : ncerr = nctk_def_arrays(ncid, [&
1780 : nctkarr_t("xredsph_extra", "dp", "number_of_reduced_dimensions, natsph_extra") &
1781 2 : ])
1782 1 : NCF_CHECK(ncerr)
1783 : end if
1784 :
1785 : ! Write variables
1786 28 : NCF_CHECK(nctk_set_datamode(ncid))
1787 :
1788 : ! scalars
1789 28 : NCF_CHECK(nf90_put_var(ncid, vid("pawprtdos"), dtset%pawprtdos))
1790 28 : NCF_CHECK(nf90_put_var(ncid, vid("prtdos"), dos%prtdos))
1791 28 : NCF_CHECK(nf90_put_var(ncid, vid("prtdosm"), dos%prtdosm))
1792 :
1793 : ! arrays
1794 28 : if (dtset%usepaw == 1) then
1795 22 : lmax_type = (pawtab(:)%l_size - 1) / 2
1796 : else
1797 38 : do itype=1,crystal%ntypat
1798 121 : lmax_type(itype) = maxval(psps%indlmn(1, :, itype))
1799 : end do
1800 : end if
1801 28 : NCF_CHECK(nf90_put_var(ncid, vid("lmax_type"), lmax_type))
1802 28 : NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions"), dos%fractions))
1803 :
1804 28 : if (dos%prtdos == 3) then
1805 20 : NCF_CHECK(nf90_put_var(ncid, vid("iatsph"), dtset%iatsph(1:dtset%natsph)))
1806 20 : NCF_CHECK(nf90_put_var(ncid, vid("ratsph"), dtset%ratsph(1:dtset%ntypat)))
1807 20 : NCF_CHECK(nf90_put_var(ncid, vid("ratsph_extra"), dtset%ratsph_extra))
1808 20 : if (dtset%natsph_extra /= 0) then
1809 1 : NCF_CHECK(nf90_put_var(ncid, vid("xredsph_extra"), dtset%xredsph_extra(:, 1:dtset%natsph_extra)))
1810 : end if
1811 : end if
1812 :
1813 28 : if (allocated(dos%fractions_m)) then
1814 3 : NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_m"), dos%fractions_m))
1815 : end if
1816 28 : if (allocated(dos%fractions_paw1)) then
1817 8 : NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_paw1"), dos%fractions_paw1))
1818 8 : NCF_CHECK(nf90_put_var(ncid, vid("dos_fractions_pawt1"), dos%fractions_pawt1))
1819 : end if
1820 :
1821 28 : call cwtime(cpu,wall,gflops,"stop")
1822 28 : write(msg,'(2(a,f8.2),a)')" epjdos_ncwrite: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
1823 28 : call wrtout(std_out,msg,"PERS")
1824 :
1825 : contains
1826 220 : integer function vid(vname)
1827 : character(len=*),intent(in) :: vname
1828 220 : vid = nctk_idname(ncid, vname)
1829 : end function vid
1830 :
1831 : end subroutine epjdos_ncwrite
1832 : !!***
1833 :
1834 : !!****f* m_epjdos/partial_dos_fractions
1835 : !! NAME
1836 : !! partial_dos_fractions
1837 : !!
1838 : !! FUNCTION
1839 : !! calculate partial DOS fractions to feed to the tetrahedron method
1840 : !! 1: project states on angular momenta
1841 : !! 2: should be able to choose certain atoms or atom types, slabs of space...
1842 : !!
1843 : !! INPUTS
1844 : !! crystal<crystal_t>= data type gathering info on symmetries and unit cell
1845 : !! dtset<type(dataset_type)>=all input variables for this dataset
1846 : !! npwarr(nkpt)=number of planewaves in basis at this k point
1847 : !! kg(3,mpw*mkmem)=reduced planewave coordinates.
1848 : !! cg(2,mcg)=planewave coefficients of wavefunctions
1849 : !! mcg=size of wave-functions array (cg) =mpw*my_nspinor*mband*mkmem*nsppol
1850 : !! collect=1 if fractions should be MPI collected at the end, 0 otherwise.
1851 : !! mpi_enreg=information about MPI parallelization
1852 : !!
1853 : !! SIDE EFFECTS
1854 : !! dos%fractions(ikpt,iband,isppol,ndosfraction) = percentage of s, p, d..
1855 : !! character on each atom for the wavefunction # ikpt,iband, isppol
1856 : !! == if prtdosm /= 0
1857 : !! dos%fractions_m(ikpt,iband,isppol,ndosfraction*mbesslang) = percentage of s, p, d..
1858 : !! character on each atom for the wavefunction # ikpt,iband, isppol (m-resolved)
1859 : !!
1860 : !! NOTES
1861 : !!
1862 : !! psi(r) = (4pi/sqrt(ucvol)) \sum_{LMG} i**l u(G) e^{i(k+G).Ra} x Y_{LM}^*(k+G) Y_{LM}(r-Ra) j_L(|k+G||r-Ra|)
1863 : !!
1864 : !! int_(ratsph) |psi(r)|**2 = \sum_LM rho(LM)
1865 : !!
1866 : !! where
1867 : !!
1868 : !! rho_{LM} = 4pi \int_o^{rc} dr r**2 ||\sum_G u(G) Y_LM^*(k+G) e^{i(k+G).Ra} j_L(|k+G| r)||**2
1869 : !!
1870 : !! where S is a RSH. The final expression is:
1871 : !!
1872 : !! When k = G0/2, we have u_{G0/2}(G) = u_{G0/2}(-G-G0)^* and P can be rewritten as
1873 : !!
1874 : !! P = (4pi i^L}/sqrt(ucvol) \sum^'_G w(G) S_{LM}(k+G) \int_0^ratsph dr r^2 j_L(|k+G|r) x
1875 : !! 2 Re[u_k(G) e^{i(k+G).R_atom}] if L = 2n
1876 : !! 2 Im[u_k(G) e^{i(k+G).R_atom}] if L = 2n + 1
1877 : !!
1878 : !! where the sum over G is done on the reduced G-sphere and w(G) = 1/2 if G=G0 else 1.
1879 : !!
1880 : !! SOURCE
1881 :
1882 42 : subroutine partial_dos_fractions(dos,crystal,dtset,eigen,occ,npwarr,kg,cg,mcg,collect,mpi_enreg)
1883 :
1884 : !Arguments ------------------------------------
1885 : !scalars
1886 : class(epjdos_t),intent(inout) :: dos
1887 : integer,intent(in) :: mcg,collect
1888 : type(MPI_type),intent(inout) :: mpi_enreg
1889 : type(dataset_type),intent(in) :: dtset
1890 : type(crystal_t),intent(in) :: crystal
1891 : !arrays
1892 : integer,intent(in) :: kg(3,dtset%mpw*dtset%mkmem),npwarr(dtset%nkpt)
1893 : real(dp),intent(in) :: cg(2,mcg)
1894 : real(dp),intent(in) :: eigen(dtset%mband*dtset%nkpt*dtset%nsppol)
1895 : real(dp),intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol)
1896 :
1897 : !Local variables-------------------------------
1898 : !scalars
1899 : integer,parameter :: prtsphere0=0 ! do not output all the band by band details for projections.
1900 : integer :: shift_b,shift_sk,iat,iatom,iband,ierr,ikpt,ilang,ioffkg,is1, is2, isoff
1901 : integer :: ipw,isppol,ixint,mbess,mcg_disk,me_kpt,shift_cg
1902 : integer :: mgfft,my_nspinor,n1,n2,n3,natsph_tot,npw_k,nradintmax
1903 : integer :: rc_ylm,itypat,nband_k, abs_shift_b, unit_procar, ipauli
1904 : real(dp),parameter :: bessint_delta = 0.1_dp
1905 : real(dp) :: arg,bessarg,bessargmax,kpgmax,rmax, cpu,wall,gflops
1906 : character(len=500) :: msg
1907 : character(len=4) :: ikproc_str
1908 : character(len=fnlen) :: filename
1909 42 : type(jlspline_t) :: jlspl
1910 42 : type(MPI_type) :: mpi_enreg_seq
1911 : !arrays
1912 42 : integer :: iindex(dtset%mpw),nband_tmp(1),npwarr_tmp(1)
1913 42 : integer,allocatable :: iatsph(:),nradint(:),atindx(:),typat_extra(:),kg_k(:,:)
1914 : real(dp) :: kpoint(3),spin(3),ylmgr_dum(1)
1915 84 : real(dp) :: xfit(dtset%mpw),yfit(dtset%mpw)
1916 42 : real(dp),allocatable :: ylm_k(:,:), bess_fit(:,:,:)
1917 42 : real(dp),allocatable :: cg_1band(:,:),cg_1kpt(:,:),kpgnorm(:),ph1d(:,:)
1918 42 : real(dp),allocatable :: ph3d(:,:,:),ratsph(:),rint(:),sum_1ll_1atom(:,:,:)
1919 42 : real(dp),allocatable :: sum_1lm_1atom(:,:,:), cplx_1lm_1atom(:,:,:,:)
1920 42 : real(dp),allocatable :: xred_sph(:,:),znucl_sph(:),phkxred(:,:)
1921 : complex(dp) :: cgcmat(2,2)
1922 : !*************************************************************************
1923 :
1924 : ! write(std_out, '(a)') ' m_epjdos%partial_dos_fractions : enter '
1925 :
1926 42 : if(dtset%natsph==dtset%natom)then
1927 26 : write(msg, '(a)') ' Compute the partial DOS fractions. This can be time-consuming. Think using natsph and iatsph.'
1928 : else
1929 16 : write(msg, '(a)') ' Compute the partial DOS fractions.'
1930 : endif
1931 42 : call wrtout(std_out,msg)
1932 :
1933 : ! for the moment, only support projection on angular momenta
1934 42 : if (dos%partial_dos_flag /= 1 .and. dos%partial_dos_flag /= 2) then
1935 0 : write(std_out,*) 'Error: partial_dos_fractions only supports angular '
1936 0 : write(std_out,*) ' momentum projection and spinor components for the moment. return to outscfcv'
1937 0 : write(std_out,*) ' partial_dos = ', dos%partial_dos_flag
1938 0 : return
1939 : end if
1940 :
1941 : ! impose all kpoints have same number of bands
1942 89 : do isppol=1,dtset%nsppol
1943 443 : do ikpt=1,dtset%nkpt
1944 401 : if (dtset%nband((isppol-1)*dtset%nkpt + ikpt) /= dtset%mband) then
1945 0 : write(std_out,*) 'Error: partial_dos_fractions wants same number of bands at each kpoint'
1946 0 : write(std_out,*) ' isppol, ikpt = ', isppol,ikpt, dtset%nband((isppol-1)*dtset%nkpt + ikpt), dtset%mband
1947 0 : write(std_out,*) ' all nband = ', dtset%nband
1948 0 : return
1949 : end if
1950 : end do
1951 : end do
1952 :
1953 : ! Real or complex spherical harmonics?
1954 42 : rc_ylm = 2; if (dos%prtdosm == 2) rc_ylm = 1
1955 :
1956 42 : me_kpt = mpi_enreg%me_kpt
1957 42 : my_nspinor = max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
1958 :
1959 42 : n1 = dtset%ngfft(1); n2 = dtset%ngfft(2); n3 = dtset%ngfft(3)
1960 : mgfft = maxval(dtset%ngfft(1:3))
1961 :
1962 42 : call cwtime(cpu, wall, gflops, "start")
1963 :
1964 42 : if (dtset%prtprocar /= 0) then
1965 : ! open file for each proc, and print header for master node
1966 1 : call int2char4(me_kpt, ikproc_str)
1967 1 : filename = 'PROCAR_'//ikproc_str
1968 1 : if (open_file(filename, msg, newunit=unit_procar, form="formatted", action="write") /= 0) then
1969 0 : ABI_ERROR(msg)
1970 : end if
1971 1 : if(mpi_enreg%me==0) then
1972 1 : write (unit_procar,'(a)') 'PROCAR lm decomposed - need to merge files yourself in parallel case!!! Or use pyprocar package'
1973 1 : if (dtset%prtprocar == 2) then
1974 1 : write (unit_procar,'(a)') ' Requested complex output of PROCAR file (prtprocar 2)'
1975 : end if
1976 1 : write (unit_procar,'(a,I10,a,I10,a,I10,a)') '# of k-points: ', dtset%nkpt, &
1977 2 : ' # of bands:', dtset%mband, ' # of ions:', dtset%natom, ch10
1978 : end if
1979 : end if
1980 :
1981 : !##############################################################
1982 : !FIRST CASE: project on angular momenta to get dos parts
1983 : !##############################################################
1984 :
1985 42 : if (dos%partial_dos_flag == 1) then
1986 41 : natsph_tot = dtset%natsph + dtset%natsph_extra
1987 :
1988 123 : ABI_MALLOC(iatsph, (natsph_tot))
1989 82 : ABI_MALLOC(typat_extra, (natsph_tot))
1990 123 : ABI_MALLOC(ratsph, (natsph_tot))
1991 82 : ABI_MALLOC(znucl_sph, (natsph_tot))
1992 82 : ABI_MALLOC(nradint, (natsph_tot))
1993 82 : ABI_MALLOC(atindx, (natsph_tot))
1994 123 : ABI_MALLOC(phkxred, (2,natsph_tot))
1995 :
1996 : ! initialize atindx
1997 114 : do iatom=1,natsph_tot
1998 114 : atindx(iatom) = iatom
1999 : end do
2000 :
2001 102 : iatsph(1:dtset%natsph) = dtset%iatsph(1:dtset%natsph)
2002 102 : do iatom=1,dtset%natsph
2003 61 : itypat = dtset%typat(iatsph(iatom))
2004 61 : typat_extra(iatom) = itypat
2005 61 : ratsph(iatom) = dtset%ratsph(itypat)
2006 102 : znucl_sph(iatom) = dtset%znucl(itypat)
2007 : end do
2008 :
2009 : ! fictitious atoms are declared with
2010 : ! %natsph_extra, %ratsph_extra and %xredsph_extra(3, dtset%natsph_extra)
2011 : ! they have atom index (natom + ii) and itype = ntypat + 1
2012 53 : do iatom=1,dtset%natsph_extra
2013 12 : typat_extra(iatom+dtset%natsph) = dtset%ntypat + 1
2014 12 : ratsph(iatom+dtset%natsph) = dtset%ratsph_extra
2015 12 : znucl_sph(iatom+dtset%natsph) = zero
2016 53 : iatsph(iatom+dtset%natsph) = dtset%natom + iatom
2017 : end do
2018 :
2019 : ! init bessel function integral for recip_ylm max ang mom + 1
2020 205 : ABI_MALLOC(sum_1ll_1atom,(dtset%nspinor**2,dos%mbesslang,natsph_tot))
2021 205 : ABI_MALLOC(sum_1lm_1atom,(dtset%nspinor**2,dos%mbesslang**2,natsph_tot))
2022 205 : ABI_MALLOC(cplx_1lm_1atom,(2,dtset%nspinor,dos%mbesslang**2,natsph_tot))
2023 :
2024 : ! Note ecuteff instead of ecut.
2025 41 : kpgmax = sqrt(dtset%ecut * dtset%dilatmx**2)
2026 41 : rmax = zero; bessargmax = zero; nradintmax = 0
2027 114 : do iatom=1,natsph_tot
2028 73 : rmax = max(rmax, ratsph(iatom))
2029 73 : bessarg = ratsph(iatom)*two_pi*kpgmax
2030 73 : bessargmax = max(bessargmax, bessarg)
2031 73 : nradint(iatom) = int (bessarg / bessint_delta) + 1
2032 114 : nradintmax = max(nradintmax,nradint(iatom))
2033 : end do
2034 : !write(std_out,*)' partial_dos_fractions: rmax=', rmax,' nradintmax: ", nradintmax
2035 : ! use same number of grid points to calculate Bessel function and to do the integration later on r
2036 : ! and make sure bessargmax is a multiple of bessint_delta
2037 41 : mbess = nradintmax
2038 41 : bessargmax = bessint_delta*mbess
2039 :
2040 123 : ABI_MALLOC(rint,(nradintmax))
2041 205 : ABI_MALLOC(bess_fit,(dtset%mpw,nradintmax,dos%mbesslang))
2042 :
2043 : ! initialize general Bessel function array on uniform grid xx, from 0 to (2 \pi |k+G|_{max} |r_{max}|)
2044 41 : call jlspl%init(mbess, bessint_delta, dos%mbesslang)
2045 :
2046 123 : ABI_MALLOC(xred_sph, (3, natsph_tot))
2047 102 : do iatom=1,dtset%natsph
2048 285 : xred_sph(:,iatom) = crystal%xred(:,iatsph(iatom))
2049 : end do
2050 53 : do iatom=1,dtset%natsph_extra
2051 89 : xred_sph(:,dtset%natsph+iatom) = dtset%xredsph_extra(:, iatom)
2052 : end do
2053 :
2054 123 : ABI_MALLOC(ph1d,(2,(2*n1+1 + 2*n2+1 + 2*n3+1)*natsph_tot))
2055 41 : call getph(atindx,natsph_tot,n1,n2,n3,ph1d,xred_sph)
2056 :
2057 : ! Fake MPI data to be used for sequential call to initylmg.
2058 41 : call initmpi_seq(mpi_enreg_seq)
2059 41 : mpi_enreg_seq%my_natom = dtset%natom
2060 :
2061 41 : shift_sk = 0
2062 41 : abs_shift_b = 0 ! offset to allow for automatic update with +1 below
2063 87 : do isppol=1,dtset%nsppol
2064 46 : ioffkg = 0
2065 :
2066 431 : do ikpt=1,dtset%nkpt
2067 344 : nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt)
2068 344 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) then
2069 78 : abs_shift_b = abs_shift_b + nband_k ! jump the whole kpt in the eig and occ arrays
2070 78 : cycle
2071 : end if
2072 266 : npw_k = npwarr(ikpt)
2073 1064 : kpoint(:) = dtset%kpt(:,ikpt)
2074 :
2075 266 : if (dtset%prtprocar /= 0) then
2076 : write (unit_procar,'(a,I7,a,3F12.6,a,F12.6,a)') &
2077 8 : ' k-point ', ikpt, ' : ', kpoint(:), ' weight = ', dtset%wtk(ikpt), ch10
2078 : end if
2079 :
2080 : ! make phkred for all atoms
2081 677 : do iat=1,natsph_tot
2082 411 : arg=two_pi*( kpoint(1)*xred_sph(1,iat) + kpoint(2)*xred_sph(2,iat) + kpoint(3)*xred_sph(3,iat) )
2083 411 : phkxred(1,iat)=cos(arg)
2084 677 : phkxred(2,iat)=sin(arg)
2085 : end do
2086 :
2087 798 : ABI_MALLOC(kg_k, (3, npw_k))
2088 210416 : kg_k = kg(:,ioffkg+1:ioffkg+npw_k)
2089 :
2090 : ! kpgnorm contains norms only for kpoints used by this processor
2091 798 : ABI_MALLOC(kpgnorm, (npw_k))
2092 266 : call getkpgnorm(crystal%gprimd, kpoint, kg_k, kpgnorm, npw_k)
2093 :
2094 : ! Now get Ylm(k, G) factors: returns "real Ylms", which are real (+m) and
2095 : ! imaginary (-m) parts of actual complex Ylm. Yl-m = Ylm*
2096 : ! Call initylmg for a single k-point (mind mpi_enreg_seq).
2097 1064 : ABI_MALLOC(ylm_k, (npw_k, dos%mbesslang**2))
2098 266 : npwarr_tmp(1) = npw_k; nband_tmp(1) = nband_k
2099 : call initylmg(crystal%gprimd,kg_k,kpoint,1,mpi_enreg_seq,dos%mbesslang,&
2100 266 : npw_k,nband_tmp,1,npwarr_tmp,1,0,crystal%rprimd,ylm_k,ylmgr_dum)
2101 :
2102 : ! get phases exp (2 pi i (k+G).x_tau) in ph3d
2103 1064 : ABI_MALLOC(ph3d,(2,npw_k,natsph_tot))
2104 266 : call ph1d3d(1,natsph_tot,kg_k,natsph_tot,natsph_tot,npw_k,n1,n2,n3,phkxred,ph1d,ph3d)
2105 :
2106 : ! get Bessel function factors on array of |k+G|*r distances
2107 : ! since we need many r distances and have a large number of different
2108 : ! |k+G|, get j_l on uniform grid (above, in array gen_besj),
2109 : ! and spline it for each kpt Gvector set.
2110 : ! Note that we use the same step based on rmax, this can lead to (hopefully small)
2111 : ! inaccuracies when we integrate from 0 up to rmax(iatom)
2112 113667 : do ixint=1,nradintmax
2113 113401 : rint(ixint) = (ixint-1)*rmax / (nradintmax-1)
2114 28972190 : do ipw=1,npw_k
2115 28858789 : xfit(ipw) = two_pi * kpgnorm(ipw) * rint(ixint)
2116 28972190 : iindex(ipw) = ipw
2117 : end do
2118 :
2119 113401 : call sort_dp(npw_k,xfit,iindex,tol14)
2120 680672 : do ilang=1,dos%mbesslang
2121 567005 : call splint(mbess, jlspl%xx, jlspl%bess_spl(:,ilang), jlspl%bess_spl_der(:,ilang), npw_k, xfit, yfit)
2122 : ! re-order results for different G vectors
2123 144974351 : do ipw=1,npw_k
2124 144860950 : bess_fit(iindex(ipw),ixint,ilang) = yfit(ipw)
2125 : end do
2126 : end do
2127 : end do ! ixint
2128 :
2129 266 : shift_b = 0
2130 4510 : do iband=1,nband_k
2131 4244 : abs_shift_b = abs_shift_b + 1
2132 4244 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle
2133 : !write(std_out,*)"in band:",iband
2134 : ! TODO: eventually import eig and occ down to here - a pain, but printing outside would imply saving a huge array in memory
2135 4244 : if (dtset%prtprocar /= 0) then
2136 200 : write (unit_procar,'(a,I7,a,F12.6,a,F12.6,a)') 'band ', iband, ' # energy ', &
2137 400 : eigen(abs_shift_b), ' # occ. ', occ(abs_shift_b), ch10
2138 : end if
2139 :
2140 : ! Select wavefunction in cg array
2141 4244 : shift_cg = shift_sk + shift_b
2142 :
2143 : call recip_ylm(bess_fit, cg(:,shift_cg+1:shift_cg+my_nspinor*npw_k), dtset%istwfk(ikpt),&
2144 : & mpi_enreg, nradint, nradintmax, dos%mbesslang , dtset%mpw, natsph_tot, typat_extra, dos%mlang_type,&
2145 : & npw_k, dtset%nspinor, ph3d, prtsphere0, rint, ratsph, rc_ylm, sum_1ll_1atom, sum_1lm_1atom, cplx_1lm_1atom,&
2146 4244 : & crystal%ucvol, ylm_k, znucl_sph)
2147 : ! on exit the sum_1atom_* have both spinors counted
2148 :
2149 : ! Accumulate
2150 10625 : do iatom=1,natsph_tot
2151 42530 : do ilang=1,dos%mbesslang
2152 : dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) &
2153 : & = dos%fractions(ikpt,iband,isppol,dos%mbesslang*(iatom-1) + ilang) &
2154 38286 : & + sum_1ll_1atom(1,ilang,iatom)
2155 : end do
2156 : end do
2157 :
2158 4244 : if (dos%prtdosm /= 0) then
2159 409 : do iatom=1,natsph_tot
2160 5559 : do ilang=1,dos%mbesslang**2
2161 : dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) &
2162 : & = dos%fractions_m(ikpt,iband,isppol,dos%mbesslang**2*(iatom-1) + ilang) &
2163 5356 : & + sum_1lm_1atom(1,ilang,iatom)
2164 : end do
2165 : end do
2166 : end if
2167 :
2168 : ! Increment band, spinor shift
2169 : !shift_b = shift_b + npw_k
2170 4244 : shift_b = shift_b + my_nspinor*npw_k
2171 :
2172 : ! now we have both spinor components. The first option is for real values projections, eventually decomposed by Pauli spinor components
2173 4510 : if (dtset%prtprocar /= 0) then
2174 200 : write (unit_procar,'(a)') 'ion s py pz px dxy dyz dz2 dxz dx2 tot'
2175 400 : do ipauli= 1,dtset%nspinor**2
2176 : ! Contract with Pauli matrices to get projections for this k and band, all atoms and ilang
2177 400 : do iatom = 1, natsph_tot
2178 200 : write (unit_procar, '(1x,I5)', advance='no') iatom
2179 2000 : do ilang=1,min(dos%mbesslang**2,9)
2180 2000 : write (unit_procar, '(F7.3)',advance='no') sum_1lm_1atom(ipauli,ilang,iatom)
2181 : end do
2182 5400 : write (unit_procar, '(F7.3)',advance='yes') sum(sum_1lm_1atom(ipauli,:,iatom))
2183 : end do
2184 : ! final line with sum over atoms
2185 200 : write (unit_procar, '(a)', advance='no') 'tot '
2186 2000 : do ilang=1,min(dos%mbesslang**2,9)
2187 3800 : write (unit_procar, '(F7.3)',advance='no') sum(sum_1lm_1atom(ipauli,ilang,:))
2188 : end do
2189 5600 : write (unit_procar, '(F7.3)',advance='yes') sum(sum_1lm_1atom(ipauli,:,:))
2190 : end do
2191 :
2192 : ! second option is to also print the complex projection on the atomic like orbital: <psi_nk | Y_lm> in a sphere
2193 : ! Two blocks are printed, first real then imaginary part
2194 200 : if (dtset%prtprocar == 2) then
2195 200 : write (unit_procar,'(2a)') 'ion s py pz px',&
2196 400 : & ' dxy dyz dz2 dxz dx2 tot'
2197 400 : do is1= 1,dtset%nspinor
2198 : ! Contracted with Pauli matrices to get projections for this k and band, all atoms and ilang
2199 400 : do iatom = 1, natsph_tot
2200 200 : write (unit_procar, '(1x,I5)', advance='no') iatom
2201 2000 : do ilang=1,min(dos%mbesslang**2,9)
2202 2000 : write (unit_procar, '(2(F7.3,1x))',advance='no') cplx_1lm_1atom(:,is1,ilang,iatom)
2203 : end do
2204 5200 : write (unit_procar, '(2(F7.3,1x))',advance='yes') sum(cplx_1lm_1atom(1,is1,:,iatom)), &
2205 5600 : & sum(cplx_1lm_1atom(2,is1,:,iatom))
2206 : end do
2207 : ! final line with sum over atoms
2208 200 : write (unit_procar, '(a)', advance='no') 'charge'
2209 2000 : do ilang=1,min(dos%mbesslang**2,9)
2210 7400 : write (unit_procar, '(F7.3,9x)',advance='no') sum(cplx_1lm_1atom(:,is1,ilang,:)**2)
2211 : end do
2212 15600 : write (unit_procar, '(F7.3,9x)',advance='yes') sum(cplx_1lm_1atom(:,is1,:,:)**2)
2213 : end do
2214 200 : write (unit_procar,*)
2215 : end if
2216 : end if
2217 :
2218 : end do ! band
2219 :
2220 : ! Increment kpt and (spin, kpt) shifts
2221 266 : ioffkg = ioffkg + npw_k
2222 266 : shift_sk = shift_sk + nband_k*my_nspinor*npw_k
2223 :
2224 266 : ABI_FREE(kg_k)
2225 266 : ABI_FREE(kpgnorm)
2226 266 : ABI_FREE(ylm_k)
2227 312 : ABI_FREE(ph3d)
2228 : end do ! ikpt
2229 : end do ! isppol
2230 :
2231 : ! collect = 1 ==> gather all contributions from different processors
2232 41 : if (collect == 1) then
2233 41 : call xmpi_sum(dos%fractions,mpi_enreg%comm_kpt,ierr)
2234 41 : if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_kpt,ierr)
2235 :
2236 : ! this is now done inside recip_ylm
2237 : ! if (mpi_enreg%paral_spinor == 1)then
2238 : ! call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr)
2239 : ! if (dos%prtdosm /= 0) call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr)
2240 : ! end if
2241 : end if
2242 :
2243 41 : ABI_FREE(atindx)
2244 41 : ABI_FREE(bess_fit)
2245 41 : ABI_FREE(iatsph)
2246 41 : ABI_FREE(typat_extra)
2247 41 : ABI_FREE(nradint)
2248 41 : ABI_FREE(ph1d)
2249 41 : ABI_FREE(phkxred)
2250 41 : ABI_FREE(ratsph)
2251 41 : ABI_FREE(rint)
2252 41 : ABI_FREE(sum_1ll_1atom)
2253 41 : ABI_FREE(sum_1lm_1atom)
2254 41 : ABI_FREE(cplx_1lm_1atom)
2255 41 : ABI_FREE(xred_sph)
2256 41 : ABI_FREE(znucl_sph)
2257 :
2258 41 : call jlspl%free()
2259 41 : call destroy_mpi_enreg(mpi_enreg_seq)
2260 :
2261 : !##############################################################
2262 : !2ND CASE: project on spinors
2263 : !##############################################################
2264 :
2265 1 : else if (dos%partial_dos_flag == 2) then
2266 :
2267 1 : if (dtset%nsppol /= 1 .or. dtset%nspinor /= 2) then
2268 0 : ABI_WARNING("spinor projection is meaningless if nsppol==2 or nspinor/=2. Not calculating projections.")
2269 0 : return
2270 : end if
2271 1 : if (my_nspinor /= 2) then
2272 0 : ABI_WARNING("spinor projection with spinor parallelization is not coded. Not calculating projections.")
2273 0 : return
2274 : end if
2275 1 : ABI_CHECK(mpi_enreg%paral_spinor == 0, "prtdos 5 does not support spinor parallelism")
2276 :
2277 : ! FIXME: We should not allocate such a large chunk of memory!
2278 1 : mcg_disk = dtset%mpw*my_nspinor*dtset%mband
2279 3 : ABI_MALLOC(cg_1kpt,(2,mcg_disk))
2280 1 : shift_sk = 0
2281 1 : isppol = 1
2282 :
2283 11 : do ikpt=1,dtset%nkpt
2284 10 : nband_k = dtset%nband((isppol-1)*dtset%nkpt + ikpt)
2285 10 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle
2286 10 : npw_k = npwarr(ikpt)
2287 :
2288 793810 : cg_1kpt(:,:) = cg(:,shift_sk+1:shift_sk+mcg_disk)
2289 30 : ABI_MALLOC(cg_1band,(2,2*npw_k))
2290 10 : shift_b=0
2291 310 : do iband=1,nband_k
2292 300 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) cycle
2293 :
2294 : ! Select wavefunction in cg array
2295 : !shift_cg = shift_sk + shift_b
2296 780960 : cg_1band(:,:) = cg_1kpt(:,shift_b+1:shift_b+2*npw_k)
2297 300 : call cg_getspin(cg_1band, npw_k, spin, cgcmat=cgcmat)
2298 :
2299 : ! MG: TODO: imag part of off-diagonal terms is missing.
2300 : ! I will add them later on.
2301 900 : do is1=1,2
2302 2100 : do is2=1,2
2303 1200 : isoff = is2 + (is1-1)*2
2304 : dos%fractions(ikpt,iband,isppol,isoff) = dos%fractions(ikpt,iband,isppol,isoff) &
2305 1800 : & + real(cgcmat(is1,is2))
2306 : end do
2307 : end do
2308 :
2309 300 : dos%fractions(ikpt,iband,isppol,5) = dos%fractions(ikpt,iband,isppol,5) + spin(1)
2310 300 : dos%fractions(ikpt,iband,isppol,6) = dos%fractions(ikpt,iband,isppol,6) + spin(2)
2311 300 : dos%fractions(ikpt,iband,isppol,7) = dos%fractions(ikpt,iband,isppol,7) + spin(3)
2312 :
2313 310 : shift_b = shift_b + 2*npw_k
2314 : end do
2315 10 : ABI_FREE(cg_1band)
2316 11 : shift_sk = shift_sk + nband_k*2*npw_k
2317 : end do
2318 1 : ABI_FREE(cg_1kpt)
2319 :
2320 : ! Gather all contributions from different processors
2321 1 : if (collect == 1) then
2322 1 : call xmpi_sum(dos%fractions,mpi_enreg%comm_kpt,ierr)
2323 1 : call xmpi_sum(dos%fractions,mpi_enreg%comm_bandfft,ierr)
2324 : !below for future use - spinors should not be parallelized for the moment
2325 : !if (mpi_enreg%paral_spinor == 1)then
2326 : ! call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr)
2327 : !end if
2328 : end if
2329 :
2330 : else
2331 0 : ABI_WARNING('only partial_dos==1 or 2 is coded')
2332 : end if
2333 :
2334 42 : if (dtset%prtprocar /= 0) close(unit_procar)
2335 :
2336 42 : call cwtime(cpu,wall,gflops,"stop")
2337 42 : write(msg,'(2(a,f8.2),a)')" partial_dos_fractions: cpu_time: ",cpu,"[s], walltime: ",wall," [s]"
2338 42 : call wrtout(std_out,msg,"PERS")
2339 : !write(std_out,*) ' m_epjdos%partial_dos_fractions : exit '
2340 :
2341 42 : end subroutine partial_dos_fractions
2342 : !!***
2343 :
2344 : !!****f* m_epjdos/partial_dos_fractions_paw
2345 : !! NAME
2346 : !! partial_dos_fractions_paw
2347 : !!
2348 : !! FUNCTION
2349 : !! Calculate PAW contributions to the partial DOS fractions (tetrahedron method)
2350 : !!
2351 : !! INPUTS
2352 : !! cprj(natom,mcprj)= <p_lmn|Cnk> coefficients for each WF |Cnk> and each |p_lmn> non-local projector
2353 : !! dimcprj(natom)=array of dimensions of array cprj (not ordered)
2354 : !! dtset structured datatype, from which one uses :
2355 : !! iatsph(nasph)=number of atoms used to project dos
2356 : !! kpt(3,nkpt) =irreducible kpoints
2357 : !! mband =max number of bands per k-point
2358 : !! mkmem =number of kpoints in memory
2359 : !! natom =number of atoms in total
2360 : !! natsph =number of atoms ofor which the spherical decomposition must be done
2361 : !! nband =number of electronic bands for each kpoint
2362 : !! nkpt =number of irreducible kpoints
2363 : !! nspinor =number of spinor components
2364 : !! nsppol =1 or 2 spin polarization channels
2365 : !! fatbands_flag =1 if pawfatbnd=1 or 2
2366 : !! mbesslang=maximum angular momentum for Bessel function expansion
2367 : !! mpi_enreg=information about MPI parallelization
2368 : !! prtdosm=option for the m-contributions to the partial DOS
2369 : !! ndosfraction=natsph*mbesslang
2370 : !! paw_dos_flag=option for the PAW contributions to the partial DOS
2371 : !! pawrad(ntypat) <type(pawrad_type)>=paw radial mesh and related data:
2372 : !! pawtab(ntypat) <type(pawtab_type)>=paw tabulated starting data
2373 : !!
2374 : !! OUTPUT
2375 : !! === If paw_dos_flag==1:
2376 : !! dos%fractions_paw1(ikpt,iband,isppol,natom*mbesslang) = contribution to
2377 : !! dos fractions from the PAW partial waves (phi)
2378 : !! dos%fractions_pawt1(ikpt,iband,isppol,natom*mbesslang) = contribution to
2379 : !! dos fractions from the PAW pseudo partial waves (phi_tild)
2380 : !!
2381 : !! SIDE EFFECTS
2382 : !! dos%fractions(ikpt,iband,isppol,ndosfraction) = percentage of s, p, d..
2383 : !! character on each atom for the wavefunction # ikpt,iband, isppol
2384 : !! As input: contains only the pseudo contribution
2385 : !! As output: contains pseudo contribution + PAW corrections
2386 : !! == if prtdosm==1
2387 : !! dos%fractions_m(ikpt,iband,isppol,ndosfraction*mbesslang*prtdosm) =
2388 : !! m discretization of partial DOS fractions
2389 : !!
2390 : !! SOURCE
2391 :
2392 24 : subroutine partial_dos_fractions_paw(dos,cprj,dimcprj,dtset,mcprj,mkmem,mpi_enreg,pawrad,pawtab)
2393 :
2394 : !Arguments ------------------------------------
2395 : !scalars
2396 : class(epjdos_t),intent(inout) :: dos
2397 : integer,intent(in) :: mcprj,mkmem
2398 : type(MPI_type),intent(in) :: mpi_enreg
2399 : type(dataset_type),intent(in) :: dtset
2400 : !arrays
2401 : integer,intent(in) :: dimcprj(dtset%natom)
2402 : type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj)
2403 : type(pawrad_type),intent(in) :: pawrad(dtset%ntypat)
2404 : type(pawtab_type),target,intent(in) :: pawtab(dtset%ntypat)
2405 :
2406 : !Local variables-------------------------------
2407 : !scalars
2408 : integer :: bandpp,basis_size,comm_kptband,cplex,fatbands_flag,iat,iatom,iband,ibg,ibsp
2409 : integer :: ierr,ikpt,il,ilang,ilmn,iln,im,iorder_cprj,ispinor,isppol,itypat,j0lmn,j0ln
2410 : integer :: jl,jlmn,jln,jm,klmn,kln,lmn_size,mbesslang,me_band,me_kpt,my_nspinor
2411 : integer :: nband_cprj_k,nband_k,ndosfraction,nprocband,nproc_spkptband,paw_dos_flag,prtdosm
2412 : real(dp) :: cpij,one_over_nproc
2413 : !character(len=500) :: msg
2414 : !arrays
2415 24 : integer ,allocatable :: dimcprj_atsph(:)
2416 24 : integer,contiguous, pointer :: indlmn(:,:)
2417 : real(dp) :: tsec(2)
2418 24 : real(dp),allocatable :: int1(:,:),int2(:,:),int1m2(:,:)
2419 24 : type(pawcprj_type),allocatable :: cprj_k(:,:)
2420 : !******************************************************************************************
2421 :
2422 : DBG_ENTER("COLL")
2423 :
2424 0 : ABI_CHECK(mkmem/=0,"mkmem==0 not supported anymore!")
2425 :
2426 24 : fatbands_flag = dos%fatbands_flag
2427 24 : mbesslang = dos%mbesslang
2428 24 : prtdosm = dos%prtdosm
2429 24 : ndosfraction = dos%ndosfraction
2430 24 : paw_dos_flag = dos%paw_dos_flag
2431 :
2432 : !m-decomposed DOS not compatible with PAW-decomposed DOS
2433 24 : if(prtdosm>=1.and.paw_dos_flag==1) then
2434 0 : ABI_ERROR('m-decomposed DOS not compatible with PAW-decomposed DOS!')
2435 : end if
2436 :
2437 : !Prepare some useful integrals
2438 24 : basis_size=pawtab(1)%basis_size
2439 24 : if (dtset%ntypat>1) then
2440 30 : do itypat=1,dtset%ntypat
2441 30 : basis_size=max(basis_size,pawtab(itypat)%basis_size)
2442 : end do
2443 : end if
2444 96 : ABI_MALLOC(int1 ,(basis_size*(basis_size+1)/2,dtset%natsph))
2445 72 : ABI_MALLOC(int2,(basis_size*(basis_size+1)/2,dtset%natsph))
2446 72 : ABI_MALLOC(int1m2,(basis_size*(basis_size+1)/2,dtset%natsph))
2447 2181 : int1=zero;int2=zero;int1m2=zero
2448 62 : do iat=1,dtset%natsph
2449 38 : iatom=dtset%iatsph(iat)
2450 38 : itypat= dtset%typat(iatom)
2451 264 : do jln=1,pawtab(itypat)%basis_size
2452 202 : j0ln=jln*(jln-1)/2
2453 897 : do iln=1,jln
2454 657 : kln=j0ln+iln
2455 657 : call simp_gen(int1(kln,iat),pawtab(itypat)%phiphj(:,kln),pawrad(itypat))
2456 859 : if (dtset%pawprtdos<2) then
2457 549 : call simp_gen(int2(kln,iat),pawtab(itypat)%tphitphj(:,kln),pawrad(itypat))
2458 549 : int1m2(kln,iat)=int1(kln,iat)-int2(kln,iat)
2459 : else
2460 108 : int2(kln,iat)=zero;int1m2(kln,iat)=int1(kln,iat)
2461 : end if
2462 : end do !iln
2463 : end do !jln
2464 : end do
2465 :
2466 : !Antiferro case
2467 24 : if (dtset%nspden==2.and.dtset%nsppol==1.and.dtset%nspinor==1) then
2468 178 : int1m2(:,:)=half*int1m2(:,:)
2469 2 : if (paw_dos_flag==1.or.fatbands_flag==1.or.prtdosm==2) then
2470 178 : int1(:,:)=half*int1(:,:);int2(:,:)=half*int2(:,:)
2471 : end if
2472 : end if
2473 :
2474 : !Init parallelism
2475 24 : comm_kptband=mpi_enreg%comm_kptband
2476 24 : nproc_spkptband=xmpi_comm_size(comm_kptband)*mpi_enreg%nproc_spinor
2477 24 : me_kpt=mpi_enreg%me_kpt ; me_band=mpi_enreg%me_band
2478 24 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
2479 24 : bandpp=1;if (mpi_enreg%paral_kgb==1) bandpp=mpi_enreg%bandpp
2480 : !Check if cprj is distributed over bands
2481 24 : nprocband=my_nspinor*dtset%mband*dtset%mkmem*dtset%nsppol/mcprj
2482 24 : if (nprocband/=mpi_enreg%nproc_band) then
2483 0 : ABI_BUG('wrong mcprj/nproc_band!')
2484 : end if
2485 :
2486 : !Quick hack: in case of parallelism, dos_fractions have already
2487 : ! been reduced over MPI processes; they have to be prepared before
2488 : ! the next reduction (at the end of the following loop).
2489 24 : if (nproc_spkptband>1) then
2490 16 : one_over_nproc=one/real(nproc_spkptband,kind=dp)
2491 : !$OMP PARALLEL DO COLLAPSE(4) DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt)
2492 136 : do ilang=1,ndosfraction
2493 256 : do isppol=1,dtset%nsppol
2494 4720 : do iband=1,dtset%mband
2495 25400 : do ikpt=1,dtset%nkpt
2496 : dos%fractions(ikpt,iband,isppol,ilang)= &
2497 25280 : & one_over_nproc*dos%fractions(ikpt,iband,isppol,ilang)
2498 : end do
2499 : end do
2500 : end do
2501 : end do
2502 : !$OMP END PARALLEL DO
2503 :
2504 16 : if (fatbands_flag==1.or.prtdosm==1.or.prtdosm==2) then
2505 : !$OMP PARALLEL DO COLLAPSE(4) DEFAULT(SHARED) PRIVATE(ilang,isppol,iband,ikpt)
2506 0 : do ilang=1,ndosfraction*mbesslang
2507 16 : do isppol=1,dtset%nsppol
2508 0 : do iband=1,dtset%mband
2509 0 : do ikpt=1,dtset%nkpt
2510 : dos%fractions_m(ikpt,iband,isppol,ilang)= &
2511 0 : & one_over_nproc*dos%fractions_m(ikpt,iband,isppol,ilang)
2512 : end do
2513 : end do
2514 : end do
2515 : end do
2516 : !$OMP END PARALLEL DO
2517 : end if
2518 : end if
2519 :
2520 24 : iorder_cprj=0
2521 :
2522 : !LOOPS OVER SPINS,KPTS
2523 24 : ibg=0
2524 50 : do isppol =1,dtset%nsppol
2525 241 : do ikpt=1,dtset%nkpt
2526 :
2527 191 : nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt)
2528 191 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me_kpt)) cycle
2529 :
2530 133 : cplex=2;if (dtset%istwfk(ikpt)>1) cplex=1
2531 133 : nband_cprj_k=nband_k/nprocband
2532 10742 : ABI_MALLOC(cprj_k,(dtset%natsph,my_nspinor*nband_cprj_k))
2533 399 : ABI_MALLOC(dimcprj_atsph,(dtset%natsph))
2534 347 : do iat=1,dtset%natsph
2535 347 : dimcprj_atsph(iat)=dimcprj(dtset%iatsph(iat))
2536 : end do
2537 133 : call pawcprj_alloc(cprj_k,0,dimcprj_atsph)
2538 133 : ABI_FREE(dimcprj_atsph)
2539 :
2540 : ! Extract cprj for this k-point.
2541 133 : ibsp=0
2542 3183 : do iband=1,nband_cprj_k
2543 6633 : do ispinor=1,my_nspinor
2544 3450 : ibsp=ibsp+1
2545 13260 : do iat=1,dtset%natsph
2546 6760 : iatom=dtset%iatsph(iat)
2547 343030 : cprj_k(iat,ibsp)%cp(:,:)=cprj(iatom,ibsp+ibg)%cp(:,:)
2548 : end do
2549 : end do
2550 : end do
2551 :
2552 : ! LOOP OVER ATOMS (natsph_extra is not included on purpose)
2553 347 : do iat=1,dtset%natsph
2554 214 : iatom=dtset%iatsph(iat)
2555 214 : itypat= dtset%typat(iatom)
2556 214 : lmn_size=pawtab(itypat)%lmn_size
2557 214 : indlmn => pawtab(itypat)%indlmn
2558 :
2559 : ! LOOP OVER BANDS
2560 214 : ibsp=0
2561 6307 : do iband=1,nband_k
2562 :
2563 5960 : if (mod((iband-1)/bandpp,nprocband)/=me_band) cycle
2564 5960 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,me_kpt)) then
2565 0 : ibsp=ibsp+my_nspinor;cycle
2566 : end if
2567 :
2568 12934 : do ispinor=1,my_nspinor
2569 6760 : ibsp=ibsp+1
2570 :
2571 46520 : do ilang=1,mbesslang
2572 :
2573 595260 : do jlmn=1,lmn_size
2574 554700 : jl=indlmn(1,jlmn)
2575 554700 : jm=indlmn(2,jlmn)
2576 554700 : j0lmn=jlmn*(jlmn-1)/2
2577 5773700 : do ilmn=1,jlmn
2578 5185200 : il=indlmn(1,ilmn)
2579 5185200 : im=indlmn(2,ilmn)
2580 5185200 : klmn=j0lmn+ilmn
2581 5185200 : kln=pawtab(itypat)%indklmn(2,klmn)
2582 :
2583 5739900 : if (il==ilang-1.and.jl==ilang-1.and.im==jm) then
2584 :
2585 163560 : cpij=cprj_k(iat,ibsp)%cp(1,ilmn)*cprj_k(iat,ibsp)%cp(1,jlmn)
2586 163560 : if (cplex==2) cpij=cpij+cprj_k(iat,ibsp)%cp(2,ilmn)*cprj_k(iat,ibsp)%cp(2,jlmn)
2587 163560 : cpij=pawtab(itypat)%dltij(klmn)*cpij
2588 :
2589 : dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= &
2590 : & dos%fractions(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + &
2591 163560 : & cpij*int1m2(kln,iat)
2592 163560 : if (prtdosm==1) then
2593 : dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= &
2594 : & dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + &
2595 5400 : & cpij*int1m2(kln,iat)
2596 : end if
2597 163560 : if (fatbands_flag==1.or.prtdosm==2) then
2598 : dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im)= &
2599 : & dos%fractions_m(ikpt,iband,isppol,mbesslang**2*(iat-1)+il**2+il+1+im) + &
2600 42240 : & cpij*int1(kln,iat)
2601 : end if
2602 163560 : if (paw_dos_flag==1) then
2603 : dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= &
2604 : & dos%fractions_paw1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + &
2605 0 : & cpij*int1(kln,iat)
2606 : dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang)= &
2607 : & dos%fractions_pawt1(ikpt,iband,isppol,mbesslang*(iat-1)+ilang) + &
2608 0 : & cpij*int2(kln,iat)
2609 : end if
2610 :
2611 : end if
2612 :
2613 : end do !ilmn
2614 : end do !jlmn
2615 :
2616 : end do ! ilang
2617 : end do ! ispinor
2618 : end do ! iband
2619 :
2620 : end do !iatom
2621 :
2622 133 : if (mkmem/=0) ibg = ibg + my_nspinor*nband_cprj_k
2623 133 : call pawcprj_free(cprj_k)
2624 6977 : ABI_FREE(cprj_k)
2625 : end do ! ikpt
2626 : end do ! isppol
2627 :
2628 24 : ABI_FREE(int1)
2629 24 : ABI_FREE(int2)
2630 24 : ABI_FREE(int1m2)
2631 :
2632 : !Reduce data in case of parallelism
2633 24 : call timab(48,1,tsec)
2634 24 : call xmpi_sum(dos%fractions,comm_kptband,ierr)
2635 24 : if (prtdosm>=1.or.fatbands_flag==1) then
2636 2 : call xmpi_sum(dos%fractions_m,comm_kptband,ierr)
2637 : end if
2638 24 : if (paw_dos_flag==1) then
2639 0 : call xmpi_sum(dos%fractions_paw1,comm_kptband,ierr)
2640 0 : call xmpi_sum(dos%fractions_pawt1,comm_kptband,ierr)
2641 : end if
2642 24 : call timab(48,2,tsec)
2643 24 : if (mpi_enreg%paral_spinor==1) then
2644 4 : call xmpi_sum(dos%fractions,mpi_enreg%comm_spinor,ierr)
2645 4 : if (prtdosm>=1.or.fatbands_flag==1) then
2646 0 : call xmpi_sum(dos%fractions_m,mpi_enreg%comm_spinor,ierr)
2647 : end if
2648 4 : if (paw_dos_flag==1) then
2649 0 : call xmpi_sum(dos%fractions_paw1, mpi_enreg%comm_spinor,ierr)
2650 0 : call xmpi_sum(dos%fractions_pawt1, mpi_enreg%comm_spinor,ierr)
2651 : end if
2652 : end if
2653 :
2654 : !Averaging: A quick hack for m-decomposed LDOS:
2655 : !BA: not valid in presence of spin-orbit coupling !
2656 24 : if (prtdosm==1.and.fatbands_flag==0) then
2657 : ! if pawfatbnd is activated, one think in the cubic harmonics basis
2658 : ! whereas prtdosm=1 is in the spherical harmonics basis.
2659 : ! the following trick is done in order to have everything
2660 : ! in the complex spherical basis (not useful for pawfatbnd if we want to
2661 : ! have e.g t2g and eg d-orbitals).
2662 2 : do iat=1,dtset%natsph
2663 7 : do il = 0, mbesslang-1
2664 16 : do im = 1, il
2665 : dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) = &
2666 : (dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im) + &
2667 2530 : dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im))/2
2668 : dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1-im) = &
2669 2535 : dos%fractions_m(:,:,:,mbesslang**2*(iat-1)+il**2+il+1+im)
2670 : end do
2671 : end do
2672 : end do !iatom
2673 : end if
2674 :
2675 : DBG_EXIT("COLL")
2676 :
2677 72 : end subroutine partial_dos_fractions_paw
2678 : !!***
2679 :
2680 0 : end module m_epjdos
2681 : !!***
|
