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