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