Branch data 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 : : !!***
|
