Line data Source code
1 : !!****m* ABINIT/m_paw_atomorb
2 : !! NAME
3 : !! m_paw_atomorb
4 : !!
5 : !! FUNCTION
6 : !! This module provides the definition of the atomorb_type used
7 : !! to store atomic orbitals on a radial mesh as well
8 : !! as methods to operate on it.
9 : !!
10 : !! Copyright (C) 2008-2026 ABINIT group (MG)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #include "libpaw.h"
18 :
19 : MODULE m_paw_atomorb
20 :
21 : USE_DEFS
22 : USE_MSG_HANDLING
23 : USE_MPI_WRAPPERS
24 : USE_MEMORY_PROFILING
25 :
26 : use m_paw_numeric
27 : use m_libpaw_tools, only : libpaw_basename, libpaw_get_free_unit
28 : use m_pawrad, only : pawrad_type, pawrad_init, bound_deriv,pawrad_copy, &
29 : & pawrad_free, pawrad_print, pawrad_isame, pawrad_ifromr, simp_gen,pawrad_deducer0
30 :
31 : implicit none
32 :
33 : private
34 : !!***
35 :
36 : !!****t* m_paw_atomorb/atomorb_type
37 : !! NAME
38 : !!
39 : !! FUNCTION
40 : !! Defines the atomorb_type datastructure type.
41 : !! It contains the atomic orbitals on a radial mesh for a given type of atom.
42 : !!
43 : !! NOTES
44 : !! * Should the radial mesh included in the object or not?
45 : !! * We might have different meshes, useful for deep states in heavy atoms!
46 : !! * Methods to be added: corekin
47 : !!
48 : !! SOURCE
49 :
50 : type, public :: atomorb_type
51 :
52 : !scalars
53 : integer :: ixc
54 : ! Exchange and correlation functional used to generate the orbitals
55 :
56 : integer :: method
57 : ! 1 for restricted, compatible only with nsppol=1.
58 : ! 2 for spin unrestricted, compatible only with nsppol=2.
59 :
60 : integer :: nspden
61 : ! Number of spin-density components.
62 :
63 : integer :: nsppol
64 : ! Number of independent spin-components.
65 : ! FIXME: here a lot of quantities might depend on nsppol in the
66 : ! case of magnetic atoms!
67 :
68 : integer :: nspinor
69 : ! Number of spinorial components
70 : ! TODO this is a quite delicate issue, for S.O. one should use J = L+S instead of L!
71 : ! If we use scalar relativistic then...
72 :
73 : integer :: l_max
74 : ! Maximum value of angular momentum l+1
75 :
76 : integer :: l_size
77 : ! Maximum value of l+1 leading to non zero Gaunt coeffs
78 : ! l_size=2*l_max-1
79 :
80 : integer :: ln_size
81 : ! Number of (l,n) components.
82 :
83 : integer :: ln2_size
84 : ! ln2_size=ln_size*(ln_size+1)/2
85 : ! where ln_size is the number of (l,n) elements for core orbitals.
86 :
87 : integer :: lmn_size
88 : ! Number of (l,m,n) elements.
89 :
90 : integer :: lmn2_size
91 : ! lmn2_size=lmn_size*(lmn_size+1)/2
92 : ! where lmn_size is the number of (l,m,n) elements for core orbitals.
93 :
94 : integer :: mesh_size
95 : ! Dimension of the radial mesh.
96 :
97 : integer :: mult
98 : ! Number of atoms of the same typat
99 :
100 : logical :: dirac
101 : ! Dirac relativism or not
102 :
103 : logical :: nc_conv
104 : ! nc has converged ?
105 :
106 : logical :: zcore_conv
107 : ! zcore has converged ?
108 :
109 : real(dp) :: edcc
110 : ! DC core energy
111 :
112 : real(dp) :: eeigc
113 : ! Core eigenvalue energy contribution
114 :
115 : real(dp) :: ehnzc
116 : ! Core Hartree nc+Z energy
117 :
118 : real(dp) :: ekinc
119 : ! Core kinetic energy
120 :
121 : real(dp) :: eigshift
122 : ! Eigenenergy shift due to truncation of local potentiel to the sphere in rcpaw
123 :
124 : real(dp) :: min_eigv
125 : ! Minimal eigenvalue of the valence orbitals
126 :
127 : real(dp) :: nresid_c
128 : ! Residual error on core density
129 :
130 : real(dp) :: rcore
131 : ! Radius of the sphere used to describe core electrons.
132 : ! It should be <= rpaw
133 :
134 : real(dp) :: zion
135 : ! zionpsp
136 : ! The ionic pseudo-charge, (giving raise to a long-range coulomb potential)
137 :
138 : real(dp) :: zcore
139 : ! Number of core electrons
140 :
141 : real(dp) :: zcore_orig
142 : ! original zcore
143 : ! This is used in RCPAW
144 :
145 : ! TODO alchemy?
146 : !real(dp) :: ziontypat
147 : ! ziontypat
148 : ! For each type of atom (might be alchemy wrt psps), the ionic pseudo-charge
149 : ! (giving raise to a long-range coulomb potential)
150 :
151 : real(dp) :: znucl
152 : ! The atomic number of the atom.
153 :
154 : ! TODO alchemy?
155 : !real(dp) :: znucltypat
156 : ! znucltypat
157 : ! The atomic number of each type of atom (might be alchemy wrt psps)
158 :
159 : character(len=fnlen) :: fname
160 : ! The filename for temporary storage.
161 :
162 : type(pawrad_type) :: radmesh
163 : ! Radial mesh
164 :
165 : !arrays
166 : integer, allocatable :: indlmn(:,:)
167 : ! indlmn(6,lmn_size)
168 : ! Array giving l,m,n,lm,ln,spin for i=lmn.
169 :
170 : integer, allocatable :: indln(:,:)
171 : ! indln(2,ln_size)
172 : ! Array giving l and n for i=ln
173 :
174 : integer, allocatable :: indklmn(:,:)
175 : ! indklmn(8,lmn2_size)
176 : ! Array giving klm, kln, abs(il-jl), (il+jl), ilm and jlm, ilmn and jlmn for each klmn=(ilmn,jlmn)
177 : ! Note: ilmn=(il,im,in) and ilmn<=jlmn
178 :
179 : !integer, allocatable :: klm2lm TODO add
180 : ! klm2lm(6,lm2_size)=Table giving il, jl ,im, jm, ilm, jlm for each klm=(ilm,jlm)
181 : ! where ilm=(il,im) and ilm<=jlm. NB: klm2lm is an application and not a bijection.
182 :
183 : integer, allocatable :: klm_diag(:)
184 : ! klm_diag(lmn2_size)
185 : ! 1 il==jl and im==jm, 0 otherwise.
186 :
187 : integer, allocatable :: klmntomn(:,:)
188 : ! klmntomn(4,lmn2_size)
189 : ! Array giving im, jm ,in, and jn for each klmn=(ilmn,jlmn)
190 : ! Note: ilmn=(il,im,in) and ilmn<=jlmn
191 : ! NB: klmntomn is an application and not a bijection
192 :
193 : integer, allocatable :: kln2ln(:,:)
194 : ! kln2ln(6,ln2_size)
195 : ! Table giving il, jl ,in, jn, iln, jln for each kln=(iln,jln)
196 : ! where iln=(il,in) and iln<=jln. NB: kln2ln is an application and not a bijection
197 :
198 : integer, allocatable :: kappa(:)
199 : ! Kappa for dirac relativism
200 :
201 : integer, allocatable :: mode(:,:,:)
202 : ! mode(ln_size,nsppol,2)
203 : ! Flag defining how the orbital is treated.
204 : ! During the pseudopotential generation we can have: ORB_FROZEN or ORB_VALENCE
205 : ! For calculations in extended systems we can have: ORB_FROZEN or ORB_RELAXED_CORE
206 : ! Namely different treatment depending of the degree of localization.
207 : ! For example the 1s in plutonium might be treated as ORB_FROZEN during
208 : ! a relaxed core calculation.
209 : ! Last index: 1 for current status, 2 for input status
210 : ! TODO define function to test the type, much safer!
211 :
212 : real(dp), allocatable :: eig(:,:)
213 : ! eig(ln_size,nsppol)
214 : ! Eigenvalues for each ln channel and spin.
215 :
216 : real(dp), allocatable :: max_occ(:,:)
217 : ! max_occ(ln_size,nsppol)
218 : ! Maximal occupancy for each, used in RCPAW
219 :
220 : real(dp), allocatable :: occ(:,:)
221 : ! occ(ln_size,nsppol)
222 : ! Occupation for each ln channel and spin.
223 :
224 : real(dp), allocatable :: occ_res(:,:)
225 : ! occ_res(ln_size,nsppol)
226 : ! Occupation residue for each ln channel and spin, used in RCPAW
227 :
228 : real(dp), allocatable :: occ_respc(:,:)
229 : ! occ_respc(ln_size,nsppol)
230 : ! Occupation preconditionned residue for each ln channel and spin, used in RCPAW.
231 :
232 : real(dp), allocatable :: phi(:,:,:)
233 : ! phi(mesh_size,ln_size,nsppol)
234 : ! Here we might have different meshes, useful for deep states in heavy atoms!
235 :
236 : ! this might be retrieved with a method get_atomden
237 : !real(dp), allocatable :: density(:)
238 : ! density(mesh_size,nspden)
239 : ! Gives the core density of the atom for each spin channel
240 : ! Total charge in first dimension,up component in second one (if present)
241 :
242 : real(dp), allocatable :: vhtnzc_orig(:)
243 : ! vhtnzc_orig(size(pawtab(itypat)%vhtnzc))
244 : ! Original vhtnzc, used in RCPAW
245 :
246 : end type atomorb_type
247 :
248 : ! public procedures.
249 : public :: destroy_atomorb
250 : public :: copy_atomorb
251 : public :: print_atomorb
252 : public :: get_overlap
253 : !!***
254 :
255 : !----------------------------------------------------------------------
256 :
257 : integer,public,parameter :: ORB_FROZEN =0
258 : integer,public,parameter :: ORB_RELAXED_CORE =1
259 : integer,public,parameter :: ORB_VALENCE =2
260 :
261 :
262 : CONTAINS !=========================================================================================================================
263 : !!***
264 :
265 : !----------------------------------------------------------------------
266 :
267 : !!****f* m_paw_atomorb/destroy_atomorb
268 : !! NAME
269 : !! destroy_atomorb
270 : !!
271 : !! FUNCTION
272 : !! Free the dynamic memory allocated in a structure of type atomorb_type.
273 : !!
274 : !! SIDE EFFECTS
275 : !! Atm <type(atomorb_type)>=datastructure containing atomic orbitals for a given type of atom.
276 : !!
277 : !! SOURCE
278 :
279 11 : subroutine destroy_atomorb(Atm)
280 :
281 : !Arguments ------------------------------------
282 : !scalars
283 : type(atomorb_type),intent(inout) :: Atm
284 :
285 : !************************************************************************
286 :
287 : !@atomorb_type
288 :
289 : ! integers
290 11 : if (allocated(Atm%indlmn)) then
291 11 : LIBPAW_DEALLOCATE(Atm%indlmn)
292 : end if
293 11 : if (allocated(Atm%indln)) then
294 11 : LIBPAW_DEALLOCATE(Atm%indln)
295 : end if
296 11 : if (allocated(Atm%indklmn)) then
297 11 : LIBPAW_DEALLOCATE(Atm%indklmn)
298 : end if
299 11 : if (allocated(Atm%klm_diag)) then
300 11 : LIBPAW_DEALLOCATE(Atm%klm_diag)
301 : end if
302 11 : if (allocated(Atm%klmntomn)) then
303 11 : LIBPAW_DEALLOCATE(Atm%klmntomn)
304 : end if
305 11 : if (allocated(Atm%kln2ln)) then
306 0 : LIBPAW_DEALLOCATE(Atm%kln2ln)
307 : end if
308 11 : if (allocated(Atm%kappa)) then
309 3 : LIBPAW_DEALLOCATE(Atm%kappa)
310 : end if
311 11 : if (allocated(Atm%mode)) then
312 11 : LIBPAW_DEALLOCATE(Atm%mode)
313 : end if
314 :
315 : !real
316 11 : if (allocated(Atm%eig)) then
317 11 : LIBPAW_DEALLOCATE(Atm%eig)
318 : end if
319 11 : if (allocated(Atm%max_occ)) then
320 11 : LIBPAW_DEALLOCATE(Atm%max_occ)
321 : end if
322 11 : if (allocated(Atm%occ)) then
323 11 : LIBPAW_DEALLOCATE(Atm%occ)
324 : end if
325 11 : if (allocated(Atm%occ_res)) then
326 11 : LIBPAW_DEALLOCATE(Atm%occ_res)
327 : end if
328 11 : if (allocated(Atm%occ_respc)) then
329 11 : LIBPAW_DEALLOCATE(Atm%occ_respc)
330 : end if
331 11 : if (allocated(Atm%phi)) then
332 11 : LIBPAW_DEALLOCATE(Atm%phi)
333 : end if
334 11 : if(allocated(atm%vhtnzc_orig)) then
335 5 : LIBPAW_DEALLOCATE(atm%vhtnzc_orig)
336 : endif
337 11 : call pawrad_free(atm%radmesh)
338 :
339 11 : end subroutine destroy_atomorb
340 : !!***
341 :
342 :
343 : !!****f* m_paw_atomorb/copy_atomorb
344 : !! NAME
345 : !! copy_atomorb
346 : !!
347 : !! FUNCTION
348 : !!
349 : !! SIDE EFFECTS
350 : !!
351 : !! SOURCE
352 :
353 1 : subroutine copy_atomorb(Atm_in,Atm_out)
354 :
355 : !Arguments ------------------------------------
356 : !scalars
357 : type(atomorb_type),intent(in) :: Atm_in
358 : type(atomorb_type),intent(inout) :: Atm_out
359 : !Local ------------------------------------
360 : !scalars
361 : integer :: s1,s2,s3
362 :
363 : !************************************************************************
364 1 : Atm_out%ixc=Atm_in%ixc
365 1 : Atm_out%method=Atm_in%method
366 1 : Atm_out%nspden=Atm_in%nspden
367 1 : Atm_out%nsppol=Atm_in%nsppol
368 1 : Atm_out%nspinor=Atm_in%nspinor
369 1 : Atm_out%l_max=Atm_in%l_max
370 1 : Atm_out%l_size=Atm_in%l_size
371 1 : Atm_out%ln_size=Atm_in%ln_size
372 1 : Atm_out%ln2_size=Atm_in%ln2_size
373 1 : Atm_out%lmn_size=Atm_in%lmn_size
374 1 : Atm_out%lmn2_size=Atm_in%lmn2_size
375 1 : Atm_out%mesh_size=Atm_in%mesh_size
376 1 : Atm_out%mult=Atm_in%mult
377 1 : Atm_out%dirac=Atm_in%dirac
378 1 : Atm_out%nc_conv=Atm_in%nc_conv
379 1 : Atm_out%zcore_conv=Atm_in%zcore_conv
380 1 : Atm_out%edcc=Atm_in%edcc
381 1 : Atm_out%eeigc=Atm_in%eeigc
382 1 : Atm_out%ehnzc=Atm_in%ehnzc
383 1 : Atm_out%ekinc=Atm_in%ekinc
384 1 : Atm_out%eigshift=Atm_in%eigshift
385 1 : Atm_out%min_eigv=Atm_in%min_eigv
386 1 : Atm_out%nresid_c=Atm_in%nresid_c
387 1 : Atm_out%rcore=Atm_in%rcore
388 1 : Atm_out%zion=Atm_in%zion
389 1 : Atm_out%zcore=Atm_in%zcore
390 1 : Atm_out%zcore_orig=Atm_in%zcore_orig
391 1 : Atm_out%znucl=Atm_in%znucl
392 1 : Atm_out%fname=Atm_in%fname
393 1 : call pawrad_copy(Atm_in%radmesh,Atm_out%radmesh)
394 :
395 :
396 1 : if(allocated(Atm_in%indlmn)) then
397 1 : s1=size(Atm_in%indlmn(1,:))
398 1 : s2=size(Atm_in%indlmn(:,1))
399 4 : LIBPAW_ALLOCATE(Atm_out%indlmn,(s2,s1))
400 202 : Atm_out%indlmn=Atm_in%indlmn
401 : endif
402 :
403 1 : if(allocated(Atm_in%indln)) then
404 1 : s1=size(Atm_in%indln(1,:))
405 1 : s2=size(Atm_in%indln(:,1))
406 4 : LIBPAW_ALLOCATE(Atm_out%indln,(s2,s1))
407 14 : Atm_out%indln=Atm_in%indln
408 : endif
409 :
410 1 : if(allocated(Atm_in%indklmn)) then
411 1 : s1=size(Atm_in%indklmn(1,:))
412 1 : s2=size(Atm_in%indklmn(:,1))
413 4 : LIBPAW_ALLOCATE(Atm_out%indklmn,(s2,s1))
414 1892 : Atm_out%indklmn=Atm_in%indklmn
415 : endif
416 :
417 1 : if(allocated(Atm_in%klm_diag)) then
418 1 : s1=size(Atm_in%klm_diag(:))
419 3 : LIBPAW_ALLOCATE(Atm_out%klm_diag,(s1))
420 212 : Atm_out%klm_diag=Atm_in%klm_diag
421 : endif
422 :
423 1 : if(allocated(Atm_in%klmntomn)) then
424 1 : s1=size(Atm_in%klmntomn(1,:))
425 1 : s2=size(Atm_in%klmntomn(:,1))
426 4 : LIBPAW_ALLOCATE(Atm_out%klmntomn,(s2,s1))
427 1052 : Atm_out%klmntomn=Atm_in%klmntomn
428 : endif
429 :
430 1 : if(allocated(Atm_in%kln2ln)) then
431 0 : s1=size(Atm_in%kln2ln(1,:))
432 0 : s2=size(Atm_in%kln2ln(:,1))
433 0 : LIBPAW_ALLOCATE(Atm_out%kln2ln,(s2,s1))
434 0 : Atm_out%kln2ln=Atm_in%kln2ln
435 : endif
436 :
437 1 : if(allocated(Atm_in%kappa)) then
438 1 : s1=size(Atm_in%kappa(:))
439 3 : LIBPAW_ALLOCATE(Atm_out%kappa,(s1))
440 6 : Atm_out%kappa=Atm_in%kappa
441 : endif
442 :
443 1 : if(allocated(Atm_in%mode)) then
444 1 : s1=size(Atm_in%mode(1,1,:))
445 1 : s2=size(Atm_in%mode(1,:,1))
446 1 : s3=size(Atm_in%mode(:,1,1))
447 5 : LIBPAW_ALLOCATE(Atm_out%mode,(s3,s2,s1))
448 14 : Atm_out%mode=Atm_in%mode
449 : endif
450 :
451 1 : if(allocated(Atm_in%eig)) then
452 1 : s1=size(Atm_in%eig(1,:))
453 1 : s2=size(Atm_in%eig(:,1))
454 4 : LIBPAW_ALLOCATE(Atm_out%eig,(s2,s1))
455 7 : Atm_out%eig=Atm_in%eig
456 : endif
457 :
458 1 : if(allocated(Atm_in%max_occ)) then
459 1 : s1=size(Atm_in%max_occ(1,:))
460 1 : s2=size(Atm_in%max_occ(:,1))
461 4 : LIBPAW_ALLOCATE(Atm_out%max_occ,(s2,s1))
462 7 : Atm_out%max_occ=Atm_in%max_occ
463 : endif
464 :
465 1 : if(allocated(Atm_in%occ)) then
466 1 : s1=size(Atm_in%occ(1,:))
467 1 : s2=size(Atm_in%occ(:,1))
468 4 : LIBPAW_ALLOCATE(Atm_out%occ,(s2,s1))
469 7 : Atm_out%occ=Atm_in%occ
470 : endif
471 :
472 1 : if(allocated(Atm_in%occ_res)) then
473 1 : s1=size(Atm_in%occ_res(1,:))
474 1 : s2=size(Atm_in%occ_res(:,1))
475 4 : LIBPAW_ALLOCATE(Atm_out%occ_res,(s2,s1))
476 7 : Atm_out%occ_res=Atm_in%occ_res
477 : endif
478 :
479 1 : if(allocated(Atm_in%occ_respc)) then
480 1 : s1=size(Atm_in%occ_respc(1,:))
481 1 : s2=size(Atm_in%occ_respc(:,1))
482 4 : LIBPAW_ALLOCATE(Atm_out%occ_respc,(s2,s1))
483 7 : Atm_out%occ_respc=Atm_in%occ_respc
484 : endif
485 :
486 1 : if(allocated(Atm_in%phi)) then
487 1 : s1=size(Atm_in%phi(1,1,:))
488 1 : s2=size(Atm_in%phi(1,:,1))
489 1 : s3=size(Atm_in%phi(:,1,1))
490 5 : LIBPAW_ALLOCATE(Atm_out%phi,(s3,s2,s1))
491 8011 : Atm_out%phi=Atm_in%phi
492 : endif
493 :
494 1 : if(allocated(atm_in%vhtnzc_orig)) then
495 1 : s1=size(Atm_in%vhtnzc_orig(:))
496 3 : LIBPAW_ALLOCATE(atm_out%vhtnzc_orig,(s1))
497 2003 : atm_out%vhtnzc_orig=atm_in%vhtnzc_orig
498 : endif
499 1 : end subroutine copy_atomorb
500 : !!***
501 :
502 :
503 : !----------------------------------------------------------------------
504 :
505 :
506 : !!****f* m_paw_atomorb/get_atomorb_charge
507 : !! NAME
508 : !! get_atomorb_charge
509 : !!
510 : !! FUNCTION
511 : !! Get core charge from a structure of type atomorb_type
512 : !! and optionally core density.
513 : !!
514 : !! INPUTS
515 : !! Atm<atomorb_type>=Structure defining the set of core orbitals.
516 : !! Radmesh<pawrad_type>=Info oh the Radial mesh used for core electrons.
517 : !!
518 : !! OUTPUT
519 : !! nele=core charge
520 : !! raddens(mesh_size)=core density (optional)
521 : !!
522 : !! SOURCE
523 :
524 : subroutine get_atomorb_charge(Atm,Radmesh,nele,radens)
525 :
526 : !Arguments ------------------------------------
527 : !scalars
528 : real(dp),intent(out) :: nele
529 : type(atomorb_type),intent(in) :: Atm
530 : type(pawrad_type),intent(in) :: Radmesh
531 : !arrays
532 : real(dp),optional,intent(out) :: radens(Atm%mesh_size,Atm%nspden)
533 :
534 : !Local variables-------------------------------
535 : !scalars
536 : integer :: iln,isppol
537 : real(dp) :: intg,focc
538 : real(dp),allocatable :: phi2nl(:)
539 :
540 : !************************************************************************
541 :
542 : if (Atm%nsppol==2) then
543 : LIBPAW_ERROR("nsppol==2 is Working in progress")
544 : end if
545 :
546 : LIBPAW_ALLOCATE(phi2nl,(Atm%mesh_size))
547 : if (PRESENT(radens)) radens = zero
548 :
549 : nele = zero
550 : do isppol=1,Atm%nsppol
551 : do iln=1,Atm%ln_size
552 : !Atm%mode(iln,isppol) TODO add option to select particular states
553 : focc = Atm%occ(iln,isppol)
554 : if (ABS(focc) > tol16) then
555 : phi2nl = Atm%phi(:,iln,isppol)**2
556 : call simp_gen(intg,phi2nl,Radmesh)
557 : nele = nele + focc*intg
558 : ! if (PRESENT(radens)) then !FIXME maybe it is better to rr**2 radens
559 : ! radens(2:Atm%mesh_size) = radens(2:Atm%mesh_size) &
560 : !& + focc * phi2nl(2:Atm%mesh_size)/(four_pi*Radmesh%rad(2:Atm%mesh_size)**2)
561 : ! end if
562 : end if
563 : end do
564 : end do
565 :
566 : LIBPAW_DEALLOCATE(phi2nl)
567 :
568 : end subroutine get_atomorb_charge
569 : !!***
570 :
571 : !----------------------------------------------------------------------
572 :
573 :
574 : !!****f* m_paw_atomorb/get_overlap
575 : !! NAME
576 : !! get_overlap
577 : !!
578 : !! FUNCTION
579 : !! Get overlap between core and valence states
580 : !!
581 : !! INPUTS
582 : !! Atm<atomorb_type>=Structure defining the set of core states
583 : !! Atmesh<pawrad_type>=Info oh the Radial mesh used for core states
584 : !! isppol=index for spin component
585 : !! nphi=number of core states
586 : !! phi(Radmesh2%mesh_size,nphi)=valence states
587 : !! phi_indln(nphi)=Array giving l and and n for i=1,nphi
588 : !! Radmesh2<pawrad_type>=Info oh the Radial mesh used for valence states
589 : !!
590 : !! OUTPUT
591 : !! overlap(ln_size,nphi)=core-valence overlap matrix
592 : !!
593 : !! SOURCE
594 :
595 1 : subroutine get_overlap(Atm,Atmesh,Radmesh2,isppol,nphi,phi,phi_indln,overlap)
596 :
597 : !Arguments ------------------------------------
598 : !scalars
599 : integer,intent(in) :: nphi,isppol
600 : type(atomorb_type),intent(in) :: Atm
601 : type(pawrad_type),target,intent(in) :: Atmesh,Radmesh2
602 : !arrays
603 : integer,intent(in) :: phi_indln(2,nphi)
604 : real(dp),target,intent(in) :: phi(Radmesh2%mesh_size,nphi)
605 : real(dp),intent(out) :: overlap(Atm%ln_size,nphi)
606 :
607 : !Local variables-------------------------------
608 : !scalars
609 : integer :: iln_atm,iphi,ll_phi,ll_atm,do_spline,iln
610 : integer :: whichdenser,size4spl,my_mesh_size
611 : real(dp) :: ybcbeg,ybcend,intg
612 : logical :: hasameq
613 : !arrays
614 1 : real(dp),pointer :: ff_spl(:,:)
615 1 : real(dp),allocatable :: der(:),ypp(:),func(:)
616 1 : real(dp),pointer :: rad4spl(:),my_pts(:)
617 :
618 : !************************************************************************
619 :
620 1 : if(isppol<=0.or.isppol>Atm%nsppol) LIBPAW_ERROR("Wrong isppol")
621 :
622 1 : call pawrad_isame(Atmesh,Radmesh2,hasameq,whichdenser)
623 :
624 1 : do_spline= 0; if (.not.hasameq) do_spline=1
625 :
626 1 : my_mesh_size = MIN(Atmesh%mesh_size,Radmesh2%mesh_size)
627 1 : ff_spl => phi
628 :
629 : ! === Spline valence onto Atom mesh (natural spline) ===
630 1 : if (do_spline==1) then
631 1 : LIBPAW_COMMENT("Splining in overlap")
632 :
633 1 : my_mesh_size = Atmesh%mesh_size
634 1 : my_pts => Atmesh%rad(1:my_mesh_size)
635 4 : LIBPAW_ALLOCATE(ff_spl,(my_mesh_size,nphi))
636 :
637 1 : size4spl = Radmesh2%mesh_size
638 1 : rad4spl => Radmesh2%rad
639 3 : LIBPAW_ALLOCATE(der,(size4spl))
640 2 : LIBPAW_ALLOCATE(ypp,(size4spl))
641 :
642 5 : do iln=1,nphi
643 2484 : ypp(:) = zero; ybcbeg = zero; ybcend = zero
644 4 : call paw_spline(rad4spl,phi(:,iln),size4spl,ybcbeg,ybcend,ypp)
645 :
646 5 : call paw_splint(size4spl,rad4spl,phi(:,iln),ypp,my_mesh_size,my_pts,ff_spl(:,iln))
647 : end do
648 :
649 1 : LIBPAW_DEALLOCATE(der)
650 1 : LIBPAW_DEALLOCATE(ypp)
651 : end if
652 :
653 3 : LIBPAW_ALLOCATE(func,(my_mesh_size))
654 17 : overlap = zero
655 :
656 5 : do iphi=1,nphi
657 4 : ll_phi = phi_indln(1,iphi)
658 17 : do iln_atm=1,Atm%ln_size
659 12 : ll_atm = Atm%indln(1,iln_atm)
660 :
661 16 : if (ll_atm == ll_phi) then ! selection rule on l
662 8298 : func(:) = Atm%phi(1:my_mesh_size,iln_atm,isppol) * ff_spl(1:my_mesh_size,iphi)
663 6 : call simp_gen(intg,func,Atmesh)
664 6 : overlap(iln_atm,iphi)=intg
665 6 : write(std_out,*)"overlap <phic_i|phi_j> for ll_phi",ll_phi,"ll_phic",ll_atm,"=",intg
666 : end if
667 :
668 : end do
669 : end do
670 1 : LIBPAW_DEALLOCATE(func)
671 :
672 1 : if (do_spline==1) then
673 1 : LIBPAW_DEALLOCATE(ff_spl)
674 : end if
675 :
676 2 : end subroutine get_overlap
677 : !!***
678 :
679 : !----------------------------------------------------------------------
680 :
681 : !!****f* m_paw_atomorb/print_atomorb
682 : !! NAME
683 : !! print_atomorb
684 : !!
685 : !! FUNCTION
686 : !! Reports info on a structure of type atomorb_type.
687 : !!
688 : !! INPUTS
689 : !! Atm <type(atomorb_type)>=datastructure containing atomic orbitals for a given type of atom.
690 : !!
691 : !! OUTPUT
692 : !!
693 : !! SOURCE
694 :
695 1 : subroutine print_atomorb(Atm,header,unit,prtvol,mode_paral)
696 :
697 : !Arguments ------------------------------------
698 : !scalars
699 : type(atomorb_type),intent(in) :: Atm
700 : integer,optional,intent(in) :: prtvol,unit
701 : character(len=*),optional,intent(in) :: header
702 : character(len=4),optional,intent(in) :: mode_paral
703 :
704 : !Local variables-------------------------------
705 : integer :: my_unt,my_prtvol,iln,ll,nn,isppol
706 : character(len=4) :: my_mode
707 : character(len=500) :: msg
708 : ! ************************************************************************
709 :
710 : !@atomorb_type
711 1 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
712 1 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
713 1 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
714 :
715 1 : msg=' ==== Info on the atomorb_type ==== '
716 1 : if (PRESENT(header)) msg=header
717 1 : call wrtout(my_unt,msg,my_mode)
718 :
719 2 : select case (Atm%method)
720 : case (1)
721 1 : msg = " Spin restricted"
722 : case(2)
723 0 : msg = " Spin unrestricted"
724 : case default
725 0 : write(msg,'(a,i3)')" Wrong method= ",Atm%method
726 1 : LIBPAW_BUG(msg)
727 : end select
728 1 : call wrtout(my_unt,msg,my_mode)
729 :
730 : write(msg,'(7(a,i5,a),(a,f8.5,a))')&
731 1 : & ' Number of spinorial components ...... ',Atm%nspinor,ch10,&
732 1 : & ' Number of ind. spin polarizations ... ',Atm%nsppol,ch10,&
733 1 : & ' Number of spin-density components ... ',Atm%nspden,ch10,&
734 1 : & ' Maximum angular momentum + 1 ........ ',Atm%l_max,ch10,&
735 1 : & ' Number of (l,n) orbitals ........... ',Atm%ln_size,ch10,&
736 1 : & ' Number of (l,m,n) orbitals ......... ',Atm%lmn_size,ch10,&
737 1 : & ' Dimensions of radial mesh ........... ',Atm%mesh_size,ch10,&
738 2 : & ' Core Radius ........................ ',Atm%rcore,ch10
739 1 : call wrtout(my_unt,msg,my_mode)
740 :
741 : write(msg,'(2(a,f8.5,a))')&
742 1 : & ' Ionic charge ........................ ',Atm%zion,ch10,&
743 2 : & ' Atomic number ....................... ',Atm%znucl,ch10
744 1 : call wrtout(my_unt,msg,my_mode)
745 :
746 2 : do isppol=1,Atm%nsppol
747 5 : do iln=1,Atm%ln_size
748 3 : ll = Atm%indln(1,iln)
749 3 : nn = Atm%indln(2,iln)
750 : write(msg,'(" n=",i2,", l=",i2,", spin=",i2,", nocc=",f15.7,", energy=",f15.7,2x,"(",a,")")')&
751 3 : & nn,ll,isppol,Atm%occ(iln,isppol),Atm%eig(iln,isppol),TRIM(my_mode2str(Atm%mode(iln,isppol,1)))
752 4 : call wrtout(my_unt,msg,my_mode)
753 : end do
754 : end do
755 :
756 1 : end subroutine print_atomorb
757 : !!***
758 :
759 : !----------------------------------------------------------------------
760 :
761 : !!****f* m_paw_atomorb/my_mode2str
762 : !! NAME
763 : !! my_mode2str
764 : !!
765 : !! FUNCTION
766 : !! Converts an integer flags defining the way an orbital is treated to a string.
767 : !!
768 : !! INPUTS
769 : !! mode=Integer
770 : !!
771 : !! OUTPUT
772 : !! str=mode. Either "Frozen", "Relazed Core", "Valence"
773 : !!
774 : !! SOURCE
775 :
776 3 : function my_mode2str(mode) result(str)
777 :
778 : !Arguments ------------------------------------
779 : !scalars
780 : integer,intent(in) :: mode
781 : character(len=50) :: str
782 :
783 : !Local variables
784 : character(len=500) :: msg
785 :
786 : !************************************************************************
787 :
788 6 : select case (mode)
789 : case (ORB_FROZEN)
790 3 : str="Frozen Orbital"
791 : case (ORB_RELAXED_CORE)
792 0 : str="Relaxed Core Orbital"
793 : case (ORB_VALENCE)
794 0 : str="Valence Orbital"
795 : case default
796 0 : write(msg,'(a,i3)')" Wrong mode= ",mode
797 3 : ABI_BUG(msg)
798 : end select
799 :
800 3 : end function my_mode2str
801 : !!***
802 :
803 0 : END MODULE m_paw_atomorb
804 : !!***
|
