Line data Source code
1 : !!****m* ABINIT/m_rcpaw
2 : !! NAME
3 : !! m_rcpaw
4 : !!
5 : !! FUNCTION
6 : !! This module contains types and subroutines linked to the PAW core relaxation
7 : !! approach
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2019-2026 ABINIT group (NBrouwer,MT, JBoust)
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 : !! PARENTS
16 : !!
17 : !! SOURCE
18 : !!
19 :
20 : #if defined HAVE_CONFIG_H
21 : #include "config.h"
22 : #endif
23 :
24 : #include "abi_common.h"
25 :
26 : module m_rcpaw
27 : use defs_basis
28 : use defs_abitypes
29 : use m_dtset
30 : use m_pawtab
31 : use m_pawrad
32 : use m_xmpi
33 : use m_abicore
34 : use m_errors
35 : use m_paw_atomorb
36 : use m_paw_atom
37 : use m_paral_atom
38 : use m_paw_atom_solve
39 : use m_pawpsp, only : pawpsp_init_core
40 : use m_extfpmd, only : extfpmd_type
41 : use defs_datatypes, only : pseudopotential_type
42 : use m_pawang, only : pawang_type
43 : use m_pawrhoij, only : pawrhoij_type
44 : use m_paw_an, only : paw_an_type
45 : use m_pawfgrtab, only : pawfgrtab_type
46 : use m_paw_finegrid, only : pawrfgd_fft
47 :
48 : #ifdef HAVE_MPI2
49 : use mpi
50 : #endif
51 :
52 :
53 : implicit none
54 :
55 : private
56 : !!***
57 :
58 : !----------------------------------------------------------------------
59 :
60 : !!****t* m_rcpaw/valdens_type
61 : !! NAME
62 : !! valdens_type
63 : !!
64 : !! FUNCTION
65 : !!
66 : !! SOURCE
67 : type,public :: valdens_type
68 : logical :: has_dens
69 : real(dp) :: compch_sph
70 : real(dp), allocatable :: rho1(:,:,:)
71 : real(dp), allocatable :: trho1(:,:,:)
72 : real(dp), allocatable :: nhat1(:,:,:)
73 : end type valdens_type
74 : !!***
75 :
76 : !----------------------------------------------------------------------
77 :
78 : !!****t* m_rcpaw/rcpaw_type
79 : !! NAME
80 : !! rcpaw_type
81 : !!
82 : !! FUNCTION
83 : !!
84 : !! SOURCE
85 : type,public :: rcpaw_type
86 : integer :: ntypat
87 : integer :: istep
88 : integer :: updatepaw(2)
89 : integer :: updateocc
90 : integer :: updatetnc
91 : logical :: frocc
92 : logical :: all_atoms_relaxed
93 : real(dp) :: nelect_core
94 : real(dp) :: nelect_core_orig
95 : real(dp) :: ehnzc
96 : real(dp) :: ekinc
97 : real(dp) :: edcc
98 : real(dp) :: eeigc
99 : real(dp) :: entropy
100 : real(dp) :: tolnc
101 : logical, allocatable :: eijkl_is_sym(:)
102 : type(atomorb_type),allocatable :: atm(:)
103 : type(atompaw_type),allocatable :: atp(:)
104 : type(valdens_type),allocatable :: val(:)
105 : end type rcpaw_type
106 : !!***
107 :
108 : !----------------------------------------------------------------------
109 :
110 : public :: rcpaw_destroy ! Destroy RCPAW
111 : public :: rcpaw_init ! Initialize RCPAW
112 : public :: rcpaw_reinit ! Re-Initialize RCPAW
113 : public :: rcpaw_core_eig ! Compute core eigenenergies
114 : public :: rcpaw_core_energies ! Compute total energy contributions from the core
115 : !!***
116 :
117 :
118 : CONTAINS !===========================================================
119 : !!***
120 :
121 : !----------------------------------------------------------------------
122 :
123 : !!****f* m_rcpaw/rcpaw_destroy
124 : !! NAME
125 : !! rcpaw_destroy
126 : !!
127 : !! FUNCTION
128 : !! Destroy RCPAW object
129 : !!
130 : !! INPUTS
131 : !!
132 : !!
133 : !! OUTPUT
134 : !!
135 : !!
136 : !! SOURCE
137 :
138 4 : subroutine rcpaw_destroy(rcpaw)
139 : !Arguments ------------------------------------
140 : !scalars
141 : integer :: ii
142 : type(rcpaw_type), pointer,intent(inout) :: rcpaw
143 :
144 : !******************************************************************************************
145 :
146 4 : if(allocated(rcpaw%atm)) then
147 8 : do ii=1,size(rcpaw%atm)
148 8 : call destroy_atomorb(rcpaw%atm(ii))
149 : enddo
150 8 : ABI_FREE(rcpaw%atm)
151 : endif
152 4 : if(allocated(rcpaw%atp)) then
153 8 : do ii=1,size(rcpaw%atp)
154 8 : call atompaw_destroy(rcpaw%atp(ii))
155 : enddo
156 8 : ABI_FREE(rcpaw%atp)
157 : endif
158 4 : if(allocated(rcpaw%val)) then
159 10 : do ii=1,size(rcpaw%val)
160 10 : call destroy_valdens(rcpaw%val(ii))
161 : enddo
162 10 : ABI_FREE(rcpaw%val)
163 : endif
164 4 : ABI_SFREE(rcpaw%eijkl_is_sym)
165 :
166 4 : end subroutine rcpaw_destroy
167 : !!***
168 :
169 :
170 : !----------------------------------------------------------------------
171 :
172 : !!****f* m_rcpaw/destroy_valdens
173 : !! NAME
174 : !! destroy_valdens
175 : !!
176 : !! FUNCTION
177 : !! Destroy valdens object
178 : !!
179 : !! INPUTS
180 : !!
181 : !!
182 : !! OUTPUT
183 : !!
184 : !!
185 : !! SOURCE
186 :
187 6 : subroutine destroy_valdens(val)
188 : !Arguments ------------------------------------
189 : !scalars
190 : type(valdens_type), intent(inout) :: val
191 :
192 : !******************************************************************************************
193 :
194 6 : val%compch_sph=zero
195 6 : val%has_dens=.false.
196 6 : ABI_SFREE(val%rho1)
197 6 : ABI_SFREE(val%trho1)
198 6 : ABI_SFREE(val%nhat1)
199 :
200 6 : end subroutine destroy_valdens
201 : !!***
202 :
203 :
204 : !----------------------------------------------------------------------
205 :
206 : !!****f* m_rcpaw/rcpaw_reinit
207 : !! NAME
208 : !! rcpaw_reinit
209 : !!
210 : !! FUNCTION
211 : !! Reinitialize rcpaw object
212 : !!
213 : !! INPUTS
214 : !!
215 : !!
216 : !! OUTPUT
217 : !!
218 : !!
219 : !! SOURCE
220 :
221 3 : subroutine rcpaw_reinit(rcpaw)
222 : !Arguments ------------------------------------
223 : !scalars
224 : integer :: itypat
225 : type(rcpaw_type), pointer,intent(inout) :: rcpaw
226 :
227 : !******************************************************************************************
228 :
229 3 : rcpaw%all_atoms_relaxed=.true.
230 6 : do itypat=1,size(rcpaw%atm)
231 3 : rcpaw%atm(itypat)%nresid_c=one
232 3 : rcpaw%atm(itypat)%nc_conv=.false.
233 3 : rcpaw%atm(itypat)%mode(1,1,1)=rcpaw%atm(itypat)%mode(1,1,2)
234 6 : if(rcpaw%atm(itypat)%mode(1,1,1)==orb_relaxed_core) then
235 3 : rcpaw%all_atoms_relaxed=.false.
236 : else
237 0 : rcpaw%atm(itypat)%nc_conv=.true.
238 : endif
239 : enddo
240 :
241 3 : end subroutine rcpaw_reinit
242 : !!***
243 :
244 :
245 :
246 : !----------------------------------------------------------------------
247 :
248 : !!****f* m_rcpaw/rcpaw_init
249 : !! NAME
250 : !! rcpaw_init
251 : !!
252 : !! FUNCTION
253 : !! Initialize the RCPAW functionality
254 : !!
255 : !! INPUTS
256 : !!
257 : !!
258 : !! OUTPUT
259 : !!
260 : !!
261 : !! SOURCE
262 :
263 4 : subroutine rcpaw_init(rcpaw,dtset,filpsp,pawrad,pawtab,ntypat,cplex,dirac,my_natom,comm_atom,mpi_atmtab)
264 : !Arguments ------------------------------------
265 : !scalars
266 : integer, intent(in) :: ntypat,my_natom,cplex
267 : integer,optional,intent(in) :: comm_atom
268 : logical, intent(in) :: dirac
269 : type(rcpaw_type), pointer, intent(inout) :: rcpaw
270 : type(dataset_type), intent(in) :: dtset
271 : !arrays
272 : integer,optional,target,intent(in) :: mpi_atmtab(:)
273 : character(len=fnlen), intent(in) :: filpsp(ntypat)
274 : type(pawrad_type), intent(in) :: pawrad(ntypat)
275 : type(pawtab_type), intent(inout) :: pawtab(ntypat)
276 :
277 : !Local variables-------------------------------
278 : !scalars
279 : integer :: itypat,iatom,lm_size,mesh_size,my_comm_atom,iat
280 : logical :: my_atmtab_allocated,paral_atom
281 : !arrays
282 8 : integer :: mult(ntypat)
283 4 : integer,pointer :: my_atmtab(:)
284 :
285 : !******************************************************************************************
286 :
287 4 : write(std_out, * ) 'RCPAW initialization'
288 :
289 : !Set up parallelism over atoms
290 4 : paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom))
291 4 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
292 4 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
293 4 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom)
294 :
295 : ! Set up eijkl_is_sym
296 12 : ABI_MALLOC(rcpaw%eijkl_is_sym,(dtset%ntypat))
297 8 : rcpaw%eijkl_is_sym(:)=.true.
298 :
299 : ! Set up multiplicity of atoms
300 4 : rcpaw%istep=0
301 8 : mult=0
302 10 : do iatom=1,dtset%natom
303 6 : itypat=dtset%typat(iatom)
304 10 : mult(itypat)=mult(itypat)+1
305 : enddo
306 :
307 : ! Allocate arrays
308 4 : if(.not.allocated(rcpaw%val)) then
309 18 : ABI_MALLOC(rcpaw%val,(my_natom))
310 : endif
311 4 : if(.not.allocated(rcpaw%atm)) then
312 16 : ABI_MALLOC(rcpaw%atm,(ntypat))
313 : endif
314 4 : if(.not.allocated(rcpaw%atp)) then
315 24 : ABI_MALLOC(rcpaw%atp,(ntypat))
316 : endif
317 :
318 : ! Init atm and atp
319 4 : rcpaw%all_atoms_relaxed=.true.
320 8 : do itypat=1,ntypat
321 4 : rcpaw%atp(itypat)%ixc=dtset%ixc
322 4 : rcpaw%atp(itypat)%xclevel=dtset%xclevel
323 4 : rcpaw%atp(itypat)%electrons=dtset%nelect
324 : call atompaw_init(pawtab(itypat),pawrad(itypat),rcpaw%atp(itypat),&
325 : & rcpaw%atm(itypat),dtset%rcpaw_sc(itypat),&
326 4 : & dtset%rcpaw_elin,dtset%rcpaw_vhtnzc,dtset%rcpaw_tpaw,dirac,filpsp(itypat),dtset%rcpaw_prtpaw)
327 12 : ABI_MALLOC(rcpaw%atm(itypat)%vhtnzc_orig,(size(pawtab(itypat)%vhtnzc)))
328 46 : rcpaw%atm(itypat)%mode=dtset%rcpaw_rctypat(itypat)
329 8012 : rcpaw%atm(itypat)%vhtnzc_orig=pawtab(itypat)%vhtnzc
330 4 : rcpaw%atm(itypat)%mult=mult(itypat)
331 4 : rcpaw%atm(itypat)%nspden=dtset%nspden
332 4 : rcpaw%atm(itypat)%eigshift=zero
333 8 : if(rcpaw%atm(itypat)%mode(1,1,1)==orb_relaxed_core) then
334 4 : rcpaw%all_atoms_relaxed=.false.
335 : else
336 0 : rcpaw%atm(itypat)%nc_conv=.true.
337 : endif
338 : enddo
339 :
340 : ! Init val
341 10 : do iat=1,my_natom
342 6 : iatom=iat;if (paral_atom) iatom=my_atmtab(iat)
343 6 : itypat=dtset%typat(iatom)
344 6 : lm_size=Pawtab(itypat)%lcut_size**2
345 6 : mesh_size=pawtab(itypat)%mesh_size
346 30 : ABI_MALLOC(rcpaw%val(iat)%nhat1,(mesh_size*cplex,lm_size,dtset%nspden))
347 24 : ABI_MALLOC(rcpaw%val(iat)%rho1,(mesh_size*cplex,lm_size,dtset%nspden))
348 24 : ABI_MALLOC(rcpaw%val(iat)%trho1,(mesh_size*cplex,lm_size,dtset%nspden))
349 347530 : rcpaw%val(iat)%nhat1=zero
350 347530 : rcpaw%val(iat)%rho1=zero
351 347530 : rcpaw%val(iat)%trho1=zero
352 6 : rcpaw%val(iat)%compch_sph=zero
353 10 : rcpaw%val(iat)%has_dens=.false.
354 : enddo
355 4 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
356 :
357 : ! Init non arrays
358 4 : rcpaw%edcc=zero
359 4 : rcpaw%eeigc=zero
360 4 : rcpaw%ehnzc=zero
361 4 : rcpaw%ekinc=zero
362 4 : rcpaw%entropy=zero
363 4 : rcpaw%tolnc=dtset%rcpaw_tolnc
364 4 : rcpaw%ntypat=ntypat
365 4 : rcpaw%nelect_core=zero
366 8 : do itypat=1,ntypat
367 8 : rcpaw%nelect_core=rcpaw%nelect_core+rcpaw%atm(itypat)%zcore*rcpaw%atm(itypat)%mult
368 : enddo
369 4 : rcpaw%nelect_core_orig=rcpaw%nelect_core
370 :
371 4 : if(dtset%rcpaw_frocc==1) then
372 0 : ABI_ERROR('rcpaw_frocc=1 is work in progress')
373 0 : rcpaw%frocc=.true.
374 : else
375 4 : rcpaw%frocc=.false.
376 : endif
377 12 : rcpaw%updatepaw=dtset%rcpaw_updatepaw
378 4 : if(rcpaw%frocc.and.rcpaw%updatepaw(2)>0) then
379 0 : rcpaw%updateocc=rcpaw%updatepaw(2)
380 : else
381 4 : rcpaw%updateocc=dtset%nstep
382 : endif
383 4 : rcpaw%updatetnc=dtset%rcpaw_updatetnc
384 :
385 : ! Init core energies
386 4 : call rcpaw_core_energies(rcpaw,ntypat)
387 4 : end subroutine rcpaw_init
388 : !!***
389 :
390 :
391 : !----------------------------------------------------------------------
392 :
393 : !!****f* m_rcpaw/rcpaw_core_eig
394 : !! NAME
395 : !! rcpaw_core_eig
396 : !!
397 : !! FUNCTION
398 : !! Computes the core eigenenergies
399 : !!
400 : !! INPUTS
401 : !!
402 : !!
403 : !! OUTPUT
404 : !!
405 : !!
406 : !! SOURCE
407 :
408 36 : subroutine rcpaw_core_eig(pawtab,pawrad,ntypat,rcpaw,dtset,&
409 18 : & nfft,vtrial,cplex,ucvol,&
410 18 : & gmet,rprimd,xred,ngfft,my_natom,&
411 18 : & distribfft,comm_fft,mpi_atmtab,comm_atom)
412 : !Arguments ------------------------------------
413 : !scalars
414 : integer, intent(in) :: ntypat,cplex
415 : integer,intent(in) :: nfft,my_natom
416 : integer,optional,intent(in) :: comm_atom
417 : integer,optional,intent(in) :: comm_fft
418 : real(dp), intent(in) :: ucvol
419 : type(distribfft_type),optional,target,intent(in) :: distribfft
420 : type(rcpaw_type), intent(inout) :: rcpaw
421 : type(dataset_type), intent(in) :: dtset
422 : !arrays
423 : integer,optional,target,intent(in) :: mpi_atmtab(:)
424 : integer,intent(in) :: ngfft(18)
425 : real(dp),intent(in) :: gmet(3,3)
426 : real(dp), intent(in) :: rprimd(3,3)
427 : real(dp), intent(in) :: xred(3,dtset%natom)
428 : real(dp),intent(in),target :: vtrial(cplex*nfft)
429 : type(pawtab_type), target,intent(inout) :: pawtab(ntypat)
430 : type(pawrad_type), intent(in) :: pawrad(ntypat)
431 :
432 : !Local variables-------------------------------
433 : !scalars
434 : integer :: me_fft,iatom
435 18 : integer, ABI_CONTIGUOUS pointer :: fftn3_distrib(:),ffti3_local(:)
436 : integer :: ii,itypat,ifft_old
437 : integer :: mesh_size,ind1,ind2
438 : integer :: il,nfgd,ifft,iln
439 : integer :: n1,n2,n3,i3,ispden
440 : integer :: my_comm_atom,iat,ierr
441 : character(len=500) :: message
442 : logical :: my_atmtab_allocated,paral_atom,grid_found
443 : real(dp) :: eigshift,r1,r2,est_err,vh1,vh2
444 : !arrays
445 18 : integer,pointer :: my_atmtab(:)
446 18 : integer,allocatable :: ifftsph(:)
447 18 : real(dp), allocatable :: nt1hat0(:)
448 18 : real(dp),allocatable :: rfgd(:,:)
449 18 : real(dp),allocatable :: vh_sph(:)
450 : integer :: local_arr(2),global_arr(2)
451 :
452 : !******************************************************************************************
453 :
454 : ! FFT grid
455 18 : if(.not.rcpaw%all_atoms_relaxed) then
456 18 : if(cplex.ne.1) then
457 0 : ABI_ERROR('cplex not 1')
458 : endif
459 18 : paral_atom=(present(comm_atom).and.(my_natom/=dtset%natom))
460 18 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
461 18 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
462 18 : call get_my_atmtab(my_comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,dtset%natom,my_natom_ref=my_natom)
463 18 : me_fft=0
464 18 : if (present(comm_fft)) then
465 18 : me_fft=xmpi_comm_rank(comm_fft)
466 : end if
467 18 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
468 18 : if (present(distribfft)) then
469 18 : grid_found=.false.
470 18 : if (n2 == distribfft%n2_coarse) then
471 0 : if (n3== size(distribfft%tab_fftdp3_distrib)) then
472 0 : fftn3_distrib => distribfft%tab_fftdp3_distrib
473 0 : ffti3_local => distribfft%tab_fftdp3_local
474 0 : grid_found=.true.
475 : end if
476 : end if
477 18 : if (n2 == distribfft%n2_fine) then
478 18 : if (n3 == size(distribfft%tab_fftdp3dg_distrib)) then
479 18 : fftn3_distrib => distribfft%tab_fftdp3dg_distrib
480 18 : ffti3_local => distribfft%tab_fftdp3dg_local
481 : grid_found = .true.
482 : end if
483 : end if
484 0 : if (.not.(grid_found)) then
485 0 : ABI_BUG('Unable to find an allocated distrib for this fft grid!')
486 : end if
487 : else
488 0 : ABI_MALLOC(fftn3_distrib,(n3))
489 0 : ABI_MALLOC(ffti3_local,(n3))
490 0 : fftn3_distrib=0;ffti3_local=(/(i3,i3=1,n3)/)
491 : end if
492 : endif
493 :
494 : ! Loop on typat
495 36 : do itypat=1,dtset%ntypat
496 36 : if(.not.rcpaw%atm(itypat)%nc_conv) then
497 18 : eigshift=zero
498 18 : est_err=zero
499 48 : do iat=1,my_natom
500 30 : iatom=iat;if (paral_atom) iatom=my_atmtab(iat)
501 48 : if(dtset%typat(iatom)==itypat) then ! Average on atoms of same type
502 30 : mesh_size=pawtab(itypat)%mesh_size
503 90 : ABI_MALLOC(nt1hat0,(pawtab(itypat)%mesh_size)) ! spherical part of nt1+nhat
504 42378 : nt1hat0=zero
505 60 : do ispden=1,dtset%nspden
506 : nt1hat0(1:pawtab(itypat)%mesh_size)=nt1hat0(1:pawtab(itypat)%mesh_size)+&
507 : & rcpaw%val(iat)%trho1(1:pawtab(itypat)%mesh_size,1,ispden)*sqrt(four*pi)*&
508 : & pawrad(itypat)%rad(1:pawtab(itypat)%mesh_size)**2+&
509 : & rcpaw%val(iat)%nhat1(1:pawtab(itypat)%mesh_size,1,ispden)*sqrt(four*pi)*&
510 42408 : & pawrad(itypat)%rad(1:pawtab(itypat)%mesh_size)**2
511 : end do
512 60 : ABI_MALLOC(vh_sph,(mesh_size))
513 30 : call poisson(nt1hat0,0,pawrad(itypat),vh_sph)
514 42348 : do il=2,mesh_size
515 42348 : vh_sph(il)=vh_sph(il)/pawrad(itypat)%rad(il)
516 : enddo
517 30 : call pawrad_deducer0(vh_sph,mesh_size,pawrad(itypat))
518 30 : ABI_FREE(nt1hat0)
519 : call pawrfgd_fft(ifftsph,gmet,n1,n2,n3,nfgd,0.5_dp,rfgd,rprimd,ucvol,xred(:,iatom),&
520 30 : & fftn3_distrib,ffti3_local,me_fft)
521 30 : r1=0.6_dp
522 30 : r2=r1
523 30 : ifft=1
524 30 : ifft_old=1
525 6996 : do ii=1,nfgd
526 27894 : if(norm2(rfgd(:,ii))<r1) then
527 414 : ifft_old=ifft
528 414 : ifft=ifftsph(ii)
529 414 : r2=r1
530 1656 : r1=norm2(rfgd(:,ii))
531 : endif
532 : enddo
533 30 : ABI_FREE(ifftsph)
534 30 : ABI_FREE(rfgd)
535 30 : global_arr(1)=r1
536 30 : global_arr(2)=0
537 30 : if(present(comm_fft)) then
538 30 : local_arr(1)=r1
539 30 : local_arr(2)=me_fft
540 : ! Check which processor is closest to minimum
541 : #if defined HAVE_MPI
542 30 : call MPI_ALLREDUCE(local_arr,global_arr,1,MPI_2INT,MPI_MINLOC,comm_fft,ierr)
543 : #endif
544 : endif
545 30 : if(me_fft==global_arr(2)) then
546 30 : ind1=pawrad_ifromr(pawrad(itypat),r1)
547 30 : ind2=pawrad_ifromr(pawrad(itypat),r2)
548 : vh1=vh_sph(ind1)+pawtab(itypat)%vhtnzc(ind1)+(r1-pawrad(itypat)%rad(ind1))*&
549 : & (pawtab(itypat)%vhtnzc(ind1+1)+vh_sph(ind1+1)-vh_sph(ind1)-pawtab(itypat)%vhtnzc(ind1))/&
550 30 : & (pawrad(itypat)%rad(ind1+1)-pawrad(itypat)%rad(ind1))
551 : vh2=vh_sph(ind2)+pawtab(itypat)%vhtnzc(ind2)+(r2-pawrad(itypat)%rad(ind2))*&
552 : & (pawtab(itypat)%vhtnzc(ind2+1)+vh_sph(ind2+1)-vh_sph(ind2)-pawtab(itypat)%vhtnzc(ind2))/&
553 30 : & (pawrad(itypat)%rad(ind2+1)-pawrad(itypat)%rad(ind2))
554 30 : eigshift=eigshift+vtrial(ifft)-vh1
555 30 : if(r1>zero) then
556 : est_err=est_err+abs((vtrial(ifft)-vh1)-&
557 12 : & (vtrial(ifft_old)-vh2))
558 : endif
559 : endif
560 30 : if(present(comm_fft)) then
561 30 : call xmpi_bcast(eigshift,global_arr(2),comm_fft,ierr)
562 30 : call xmpi_bcast(est_err,global_arr(2),comm_fft,ierr)
563 : endif
564 30 : ABI_FREE(vh_sph)
565 : endif
566 : enddo
567 : ! mpi reduction
568 18 : if(paral_atom) then
569 0 : call xmpi_sum(eigshift,my_comm_atom,ierr)
570 0 : call xmpi_bcast(eigshift,0,my_comm_atom,ierr)
571 0 : call xmpi_sum(est_err,my_comm_atom,ierr)
572 0 : call xmpi_bcast(est_err,0,my_comm_atom,ierr)
573 : endif
574 18 : rcpaw%atm(itypat)%eigshift=eigshift/rcpaw%atm(itypat)%mult
575 18 : write(std_out,*) 'ESTIMATED ERROR ON CORE EIGS OF TYPAT',itypat,' = ',est_err*Ha_eV, ' eV'
576 90 : if(allocated(rcpaw%atm(itypat)%eig)) rcpaw%atm(itypat)%eig=rcpaw%atm(itypat)%eig+eigshift/rcpaw%atm(itypat)%mult ! Average on atoms of same type
577 18 : if(rcpaw%atm(itypat)%nresid_c<rcpaw%tolnc.and.rcpaw%istep>rcpaw%updatepaw(2).and.rcpaw%updatepaw(2)/=0.and.&
578 : & (dtset%rcpaw_vhtnzc/=2.or.dtset%rcpaw_frocc==1))then
579 0 : rcpaw%atm(itypat)%nc_conv=.true.
580 0 : write(message,'(a,i5,a)') 'RCPAW: core for typat ',itypat, ' converged'
581 0 : call wrtout(ab_out,message)
582 : endif
583 : endif
584 : enddo
585 :
586 : !Destroy atom table used for parallelism
587 18 : if(.not.rcpaw%all_atoms_relaxed) then
588 18 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
589 18 : if (.not.present(distribfft)) then
590 0 : ABI_FREE(fftn3_distrib)
591 0 : ABI_FREE(ffti3_local)
592 : end if
593 : endif
594 :
595 : ! Update convergence status of cores
596 18 : rcpaw%all_atoms_relaxed=.true.
597 36 : do itypat=1,dtset%ntypat
598 36 : if(rcpaw%atm(itypat)%zcore_conv.and.rcpaw%atm(itypat)%nc_conv) then
599 0 : rcpaw%atm(itypat)%mode(:,:,1)=ORB_FROZEN
600 : else
601 18 : rcpaw%all_atoms_relaxed=.false.
602 : endif
603 : enddo
604 :
605 : ! Print core eigenenergies and occupations
606 36 : do itypat=1,dtset%ntypat
607 18 : write(std_out,*) 'RCPAW core eigenergies (Ha) and occupations for typat ',itypat
608 90 : do iln=1,rcpaw%atm(itypat)%ln_size
609 72 : write(std_out,*) rcpaw%atm(itypat)%eig(iln,1),rcpaw%atm(itypat)%occ(iln,1)
610 : enddo
611 : enddo
612 :
613 : ! Print in abo
614 18 : if(rcpaw%istep>=rcpaw%updatepaw(1).and.rcpaw%updatepaw(2)/=0.and.rcpaw%istep<=rcpaw%updatepaw(2))then
615 2 : write(message,'(a)') 'RCPAW: updated PAW transform'
616 2 : call wrtout(ab_out,message)
617 2 : write(std_out,*) 'RCPAW: updated PAW transform'
618 : endif
619 18 : if(rcpaw%istep==rcpaw%updateocc.and.rcpaw%frocc) then
620 0 : write(message,'(a)') 'RCPAW: freezing core occupations'
621 0 : call wrtout(ab_out,message)
622 0 : write(std_out,*) 'RCPAW: freezing core occupations'
623 : endif
624 :
625 36 : end subroutine rcpaw_core_eig
626 : !!***
627 :
628 :
629 :
630 : !----------------------------------------------------------------------
631 :
632 : !!****f* m_rcpaw/rcpaw_core_energies
633 : !! NAME
634 : !! rcpaw_core_energies
635 : !!
636 : !! FUNCTION
637 : !!
638 : !! INPUTS
639 : !!
640 : !!
641 : !! OUTPUT
642 : !!
643 : !!
644 : !! SOURCE
645 :
646 23 : subroutine rcpaw_core_energies(rcpaw,ntypat)
647 : !Arguments ------------------------------------
648 : !scalars
649 : integer, intent(in) :: ntypat
650 : type(rcpaw_type), pointer, intent(inout) :: rcpaw
651 : !arrays
652 :
653 : !Local variables-------------------------------
654 : !scalars
655 : integer :: itypat
656 :
657 : !******************************************************************************************
658 :
659 23 : rcpaw%ehnzc=zero
660 23 : rcpaw%edcc=zero
661 23 : rcpaw%eeigc=zero
662 23 : rcpaw%ekinc=zero
663 46 : do itypat=1,ntypat
664 46 : if(rcpaw%atm(itypat)%zcore_orig>zero) then
665 23 : rcpaw%edcc=rcpaw%edcc+rcpaw%atm(itypat)%edcc*rcpaw%atm(itypat)%mult
666 23 : rcpaw%ekinc=rcpaw%ekinc+rcpaw%atm(itypat)%ekinc*rcpaw%atm(itypat)%mult
667 23 : rcpaw%eeigc=rcpaw%eeigc+rcpaw%atm(itypat)%eeigc*rcpaw%atm(itypat)%mult
668 23 : rcpaw%ehnzc=rcpaw%ehnzc+rcpaw%atm(itypat)%ehnzc*rcpaw%atm(itypat)%mult
669 : endif
670 : enddo
671 23 : write(std_out,*)'RCPAW energies at step ',rcpaw%istep,':', ' ekinc = ',rcpaw%ekinc,' eeig = ',&
672 46 : & rcpaw%eeigc,' edcc = ',rcpaw%edcc,' ehnzc = ',rcpaw%ehnzc
673 :
674 23 : end subroutine rcpaw_core_energies
675 : !!***
676 :
677 :
678 0 : end module m_rcpaw
679 : !!***
680 :
|
