Line data Source code
1 : !!****m* ABINIT/m_hamiltonian
2 : !! NAME
3 : !! m_hamiltonian
4 : !!
5 : !! FUNCTION
6 : !! This module provides the definition of the gs_hamiltonian_type and of the rf_hamiltonian_type
7 : !! datastructures used in the "getghc" and "getgh1c" routines to apply the Hamiltonian (or
8 : !! its derivative) on a wavefunction. Methods to initialize or destroy the objects are defined here.
9 : !!
10 : !! TODO
11 : !! All array pointers in H datatypes should be declared as contiguous for efficiency reasons.
12 : !! (well, here performance is critical). Client code should make sure they always point contiguous targets.
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2009-2026 ABINIT group (MG, MT)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !!
20 : !! SOURCE
21 :
22 : #if defined HAVE_CONFIG_H
23 : #include "config.h"
24 : #endif
25 :
26 : #include "abi_common.h"
27 :
28 : module m_hamiltonian
29 :
30 : use iso_fortran_env, only : int32, int64, real32, real64
31 :
32 : use defs_basis
33 : use m_abicore
34 : use m_errors
35 : use m_xmpi
36 :
37 : use m_fstrings, only : sjoin, itoa, ktoa, yesno
38 : use defs_datatypes, only : pseudopotential_type
39 : use defs_abitypes, only : MPI_type
40 : use m_dtset, only : dataset_type
41 : use m_copy, only : addr_copy
42 : use m_crystal, only : crystal_t
43 : use m_geometry, only : metric
44 : use m_pawtab, only : pawtab_type
45 : use m_pawfgr, only : pawfgr_type
46 : use m_fftcore, only : sphereboundary
47 : use m_fft, only : fftpac
48 : use m_fourier_interpol, only : transgrid
49 : use m_pawcprj, only : pawcprj_getdim
50 : use m_paw_ij, only : paw_ij_type
51 : use m_paral_atom, only : get_my_atmtab, free_my_atmtab
52 : use m_electronpositron, only : electronpositron_type, electronpositron_calctype
53 : use m_kg, only : ph1d3d, getph, mkkin, mkkpg
54 : use m_fock, only : fock_common_type, fock_BZ_type, fock_ACE_type, fock_type
55 : use m_mkffnl, only : mkffnl_objs
56 : use m_initylmg, only : initylmg_k
57 : use m_ompgpu_fourwf, only : free_ompgpu_fourwf
58 :
59 : #if defined HAVE_GPU_CUDA
60 : use m_manage_cuda
61 : #endif
62 :
63 : #if defined HAVE_FC_ISO_C_BINDING
64 : use, intrinsic :: iso_c_binding, only : c_ptr,c_loc,c_f_pointer,c_int32_t,c_int64_t,c_size_t
65 : #endif
66 :
67 : #if defined HAVE_GPU && defined HAVE_YAKL
68 : use gator_mod
69 : #endif
70 :
71 : implicit none
72 :
73 : private
74 :
75 : public :: pawdij2ekb
76 : public :: pawdij2e1kb
77 : public :: gspot_transgrid_and_pack ! Set up local potential vlocal on the coarse FFT mesh from vtrial on the fine mesh.
78 :
79 : ! These constants select how H_{k',k} is applied in reciprocal space
80 : integer,parameter,public :: KPRIME_H_K=1, K_H_KPRIME=2, K_H_K=3, KPRIME_H_KPRIME=4
81 : !!***
82 :
83 : !----------------------------------------------------------------------
84 :
85 : !!****t* m_hamiltonian/gs_hamiltonian_type
86 : !! NAME
87 : !! gs_hamiltonian_type
88 : !!
89 : !! FUNCTION
90 : !! This datastructure contains the information about one Hamiltonian,
91 : !! needed in the "getghc" routine, that applies the Hamiltonian on a wavefunction.
92 : !! The Hamiltonian is expressed in reciprocal space:
93 : !!
94 : !! H_k^prime,k = exp(-i.k^prime.r^prime) H exp(i.k.r)
95 : !!
96 : !! In most cases k = k^prime and the k^prime objects are simply pointers to k objects.
97 : !!
98 : !! SOURCE
99 :
100 : type,public :: gs_hamiltonian_type
101 :
102 : ! ===== Integer scalars
103 :
104 : integer :: dimekb1 = -1
105 : ! First dimension of Ekb
106 : ! Same as psps%dimekb
107 : ! ->Norm conserving : Max. number of Kleinman-Bylander energies
108 : ! for each atom type
109 : ! dimekb1=lnmax
110 : ! ->PAW : Max. number of Dij coefficients connecting projectors
111 : ! for each atom
112 : ! dimekb1=cplex_dij*lmnmax*(lmnmax+1)/2
113 :
114 : integer :: dimekb2 = -1
115 : ! Second dimension of Ekb
116 : ! ->Norm conserving psps: dimekb2=ntypat
117 : ! ->PAW : dimekb2=natom
118 :
119 : integer :: dimekbq = -1
120 : ! Fourth dimension of Ekb
121 : ! 2 if Ekb factors contain a exp(-iqR) phase, 1 otherwise
122 :
123 : integer :: ispin_gbt = -1
124 : ! 1, GBT spin up channel (uses k - q/2)
125 : ! 2, GBT spin down channel (uses k + q/2)
126 :
127 : integer :: istwf_k = -1
128 : ! option parameter that describes the storage of wfs at k
129 :
130 : integer :: istwf_kp = -1
131 : ! option parameter that describes the storage of wfs at k^prime
132 :
133 : integer :: lmnmax = -1
134 : ! Maximum number of different l,m,n components over all types of psps.
135 : ! same as dtset%lmnmax
136 :
137 : integer :: matblk = -1
138 : ! dimension of the array ph3d
139 :
140 : integer :: mgfft = -1
141 : ! maximum size for 1D FFTs (same as dtset%mgfft)
142 :
143 : integer :: mpsang = -1
144 : ! Highest angular momentum of non-local projectors over all type of psps.
145 : ! shifted by 1 : for all local psps, mpsang=0; for largest s, mpsang=1,
146 : ! for largest p, mpsang=2; for largest d, mpsang=3; for largest f, mpsang=4
147 : ! This gives also the number of non-local "channels"
148 : ! same as psps%mpsang
149 :
150 : integer :: mpssoang = -1
151 : ! Maximum number of channels, including those for treating the spin-orbit coupling
152 : ! For NC pseudopotentials only:
153 : ! when mpspso=1, mpssoang=mpsang
154 : ! when mpspso=2, mpssoang=2*mpsang-1
155 : ! For PAW: same as mpsang
156 : ! same as psps%mpssoang
157 :
158 : integer :: natom = -1
159 : ! The number of atoms for this dataset; same as dtset%natom
160 :
161 : integer :: nfft = -1
162 : ! number of FFT grid points same as dtset%nfft
163 :
164 : integer :: nfft_blocks = 1
165 : ! number of blocks fourwf computation is divided into
166 :
167 : integer :: npw_k = -1
168 : ! number of plane waves at k
169 : ! In case of band-FFT parallelism, npw_k is the number of plane waves
170 : ! processed by current proc
171 :
172 : integer :: npw_fft_k = -1
173 : ! number of plane waves at k used to apply Hamiltonian when band-FFT
174 : ! parallelism is activated (i.e. data are distributed in the "FFT" configuration)
175 :
176 : integer :: npw_kp = -1
177 : ! number of plane waves at k^prime
178 : ! In case of band-FFT parallelism, npw_kp is the number of plane waves
179 : ! processed by current proc
180 :
181 : integer :: npw_fft_kp = -1
182 : ! number of plane waves at k^prime used to apply Hamiltonian when band-FFT
183 : ! parallelism is activated (i.e. data are distributed in the "FFT" configuration)
184 :
185 : integer :: nspinor = -1
186 : ! Number of spinorial components
187 :
188 : integer :: nsppol = -1
189 : ! Total number of spin components (1=non-polarized, 2=polarized)
190 :
191 : integer :: ntypat = -1
192 : ! Number of types of pseudopotentials same as dtset%ntypat
193 :
194 : integer :: nvloc = -1
195 : ! Number of components of vloc
196 : ! usually, nvloc=1, except in the non-collinear magnetism case, where nvloc=4
197 :
198 : integer :: n4 = -1, n5 = -1, n6 = -1
199 : ! same as ngfft(4:6)
200 :
201 : integer :: gpu_option = -1
202 : ! Governs the choice of the GPU implementation:
203 : ! = 0 ==> do not use GPU
204 : ! > 0 ==> see defs_basis.F90 to have the list of possible GPU implementations
205 :
206 : integer :: usecprj = -1
207 : ! usecprj= 1 if cprj projected WF are stored in memory
208 : ! = 0 if they are to be computed on the fly
209 :
210 : integer :: usepaw = -1
211 : ! if usepaw=0 , use norm-conserving psps part of the code
212 : ! is usepaw=1 , use paw part of the code
213 :
214 : integer :: useylm = -1
215 : ! governs the way the nonlocal operator is to be applied:
216 : ! 1=using Ylm, 0=using Legendre polynomials
217 :
218 : integer :: use_gbt = 0
219 : ! 0, use normal non-collinear calculation
220 : ! 1, use spin spiral calculation
221 : ! 2, use spin spiral with z-component SOC
222 :
223 : integer :: zora = 0
224 : ! zora=0: no zora terms. zora=1: use available zora terms
225 : ! currently this is limited to nuclear dipole moment terms,
226 : ! although pawspnorb should eventually be included here as well
227 :
228 : ! ===== Integer arrays
229 :
230 : #if defined HAVE_GPU && defined HAVE_YAKL
231 : integer(c_int32_t), contiguous, pointer :: atindx(:) => null()
232 : #else
233 : integer, allocatable :: atindx(:)
234 : #endif
235 : ! atindx(natom)
236 : ! index table for atoms (see gstate.f)
237 :
238 : #if defined HAVE_GPU && defined HAVE_YAKL
239 : integer(c_int32_t), contiguous, pointer :: atindx1(:) => null()
240 : #else
241 : integer, allocatable :: atindx1(:)
242 : #endif
243 : ! atindx1(natom)
244 : ! index table for atoms, inverse of atindx (see gstate.f)
245 :
246 : integer, allocatable :: dimcprj(:)
247 : ! dimcprj(natom*usepaw)=dimensions of array cprj
248 : ! dimcprj(ia)=cprj(ia,:)%nlmn
249 : ! atoms are ordered by atom-type
250 :
251 : integer, allocatable :: gbound_k(:,:)
252 : ! gbound_k(2*mgfft+8,2)
253 : ! G sphere boundary, for each plane wave at k
254 :
255 : #if defined HAVE_GPU && defined HAVE_YAKL
256 : integer(c_int32_t), contiguous, pointer :: indlmn(:,:,:) => null()
257 : #else
258 : integer(c_int32_t), allocatable :: indlmn(:,:,:)
259 : #endif
260 : ! indlmn(6,lmnmax,ntypat)
261 : ! For each type of psp,
262 : ! array giving l,m,n,lm,ln,spin for i=ln (if useylm=0)
263 : ! or i=lmn (if useylm=1)
264 :
265 : #if defined HAVE_GPU && defined HAVE_YAKL
266 : integer(c_int32_t), contiguous, pointer :: nattyp(:) => null()
267 : #else
268 : integer, allocatable :: nattyp(:)
269 : #endif
270 : ! nattyp(ntypat)
271 : ! # of atoms of each type
272 :
273 : integer :: ngfft(18)
274 : ! ngfft(1:3)=integer fft box dimensions
275 : ! ngfft(4:6)=integer fft box dimensions, might be augmented for CPU speed
276 : ! ngfft(7)=fftalg
277 : ! ngfft(8)=fftalg
278 :
279 : integer :: nloalg(3)
280 : ! governs the choice of the algorithm for non-local operator same as dtset%nloalg
281 :
282 : integer, allocatable :: pspso(:)
283 : ! pspso(ntypat)
284 : ! For each type of psp, 1 if no spin-orbit component is taken
285 : ! into account, 2 if a spin-orbit component is used
286 : ! Relevant for NC-psps and PAW.
287 :
288 : #if defined HAVE_GPU && defined HAVE_YAKL
289 : integer(c_int32_t), contiguous, pointer :: typat(:) => null()
290 : #else
291 : integer, allocatable :: typat(:)
292 : #endif
293 : ! typat(natom)
294 : ! type of each atom
295 :
296 : ! integer, allocatable :: indpw_k(:,:)
297 : ! indpw_k(4,npw_fft_k)
298 : ! array which gives fft box index for given basis sphere
299 :
300 : ! Integer pointers
301 :
302 : integer, contiguous, pointer :: gbound_kp(:,:) => null()
303 : ! gbound_kp(2*mgfft+8,2)
304 : ! G sphere boundary, for each plane wave at k^prime
305 :
306 : #if defined HAVE_GPU && defined HAVE_YAKL
307 : integer(int32), contiguous, pointer :: kg_k(:,:) => null()
308 : #else
309 : integer, contiguous, pointer :: kg_k(:,:) => null()
310 : #endif
311 : ! kg_k(3,npw_fft_k)
312 : ! G vector coordinates with respect to reciprocal lattice translations
313 : ! at k
314 :
315 : integer, contiguous, pointer :: kg_kp(:,:) => null()
316 : ! kg_kp(3,npw_fft_kp)
317 : ! G vector coordinates with respect to reciprocal lattice translations
318 : ! at k^prime
319 :
320 : ! ===== Real scalars
321 :
322 : real(dp) :: ucvol = -one
323 : ! unit cell volume (Bohr**3)
324 :
325 : ! ===== Real arrays
326 :
327 : real(dp), allocatable :: ekb_spin(:,:,:,:,:)
328 : ! ekb_spin(dimekb1,dimekb2,nspinor**2,dimekbq,my_nsppol)
329 : ! Contains the values of ekb array for all spins treated by current process
330 : ! See ekb description ; ekb is pointer to ekb_spin(:,:,:,:,my_isppol)
331 :
332 : real(dp), allocatable :: sij(:,:)
333 : ! sij(dimekb1,ntypat*usepaw) = overlap matrix for paw calculation
334 :
335 : real(dp) :: gmet(3,3)
336 : ! reciprocal space metric tensor in Bohr**-2
337 :
338 : real(dp) :: gprimd(3,3)
339 : ! dimensional reciprocal space primitive translations (Bohr^-1)
340 :
341 : real(dp) :: kpt_k(3)
342 : ! dimensionless k point coordinates wrt reciprocal lattice vectors
343 :
344 : real(dp) :: kpt_kp(3)
345 : ! dimensionless k^prime point coordinates wrt reciprocal lattice vectors
346 :
347 : real(dp) :: spinaxis(3)
348 : ! spin quantization axis
349 :
350 : real(dp), allocatable :: nucdipmom(:,:)
351 : ! nucdipmom(3,natom)
352 : ! nuclear dipole moments at each atomic position
353 :
354 : #if defined HAVE_GPU && defined HAVE_YAKL
355 : real(c_double), contiguous, pointer :: ph1d(:,:) => null()
356 : #else
357 : real(dp), allocatable :: ph1d(:,:)
358 : #endif
359 : ! ph1d(2,3*(2*mgfft+1)*natom)
360 : ! 1-dim phase arrays for structure factor (see getph.f).
361 :
362 : real(dp), allocatable :: phkxred(:,:)
363 : ! phkxred(2,natom)
364 : ! phase factors exp(2 pi k.xred) at k
365 :
366 : ! ===== Real pointers
367 :
368 : real(dp), contiguous, pointer :: ekb(:,:,:,:) => null()
369 : ! ekb(dimekb1,dimekb2,nspinor**2,dimekbq)
370 : ! ->Norm conserving : (Real) Kleinman-Bylander energies (hartree)
371 : ! for number of basis functions (l,n) (lnmax)
372 : ! and number of atom types (ntypat)
373 : ! dimekb1=lnmax ; dimekb2=ntypat ; dimekbq=1
374 : ! ->PAW : (Real, symmetric) Frozen part of Dij coefficients
375 : ! to connect projectors
376 : ! for number of basis functions (l,m,n) (lmnmax)
377 : ! and number of atom (natom)
378 : ! dimekb1=lmnmax*(lmnmax+1)/2 ; dimekb2=natom ; dimekbq=1
379 : ! ekb is spin dependent in the case of PAW calculations.
380 : ! For each spin component, ekb points to ekb_spin(:,:,:,:,my_isppol)
381 : ! dimekbq=2 if Ekb factors contain a exp(-iqR) phase, dimekbq=1 otherwise
382 : ! About the non-local factors symmetry:
383 : ! - The lower triangular part of the Dij matrix can be deduced from the upper one
384 : ! with the following relation: D^s2s1_ji = (D^s1s2_ij)^*
385 : ! where s1,s2 are spinor components
386 :
387 : real(dp), contiguous, pointer :: ffnl_k(:,:,:,:) => null()
388 : ! ffnl_k(npw_fft_k,2,dimffnl_k,ntypat)
389 : ! nonlocal form factors at k
390 :
391 : real(dp), contiguous, pointer :: ffnl_kp(:,:,:,:) => null()
392 : ! ffnl_kp(npw_fft_kp,2,dimffnl_kp,ntypat)
393 : ! nonlocal form factors at k_prime
394 :
395 : real(dp), contiguous, pointer :: kinpw_k(:) => null()
396 : ! kinpw_k(npw_fft_k)
397 : ! (modified) kinetic energy for each plane wave at k
398 : ! CAVEAT: In band mode, this array is NOT EQUIVALENT to kinpw(npw_k)
399 :
400 : real(dp), contiguous, pointer :: kinpw_kp(:) => null()
401 : ! kinpw_kp(npw_fft_kp)
402 : ! (modified) kinetic energy for each plane wave at k^prime
403 :
404 : real(dp), contiguous, pointer :: kpg_k(:,:) => null()
405 : ! kpg_k(3,npw_fft_k)
406 : ! k+G vector coordinates at k
407 :
408 : real(dp), contiguous, pointer :: kpg_kp(:,:) => null()
409 : ! kpg_kp(3,npw_fft_kp)
410 : ! k^prime+G vector coordinates at k^prime
411 :
412 : real(dp), contiguous, pointer :: phkpxred(:,:) => null()
413 : ! phkpxred(2,natom)
414 : ! phase factors exp(2 pi k^prime.xred) at k^prime
415 :
416 : real(dp), contiguous, pointer :: ph3d_k(:,:,:) => null()
417 : ! ph3d_k(2,npw_fft_k,matblk)
418 : ! 3-dim structure factors, for each atom and plane wave at k
419 :
420 : real(dp), contiguous, pointer :: ph3d_kp(:,:,:) => null()
421 : ! ph3d_kp(2,npw_fft_kp,matblk)
422 : ! 3-dim structure factors, for each atom and plane wave at k^prime
423 :
424 : real(dp), contiguous, pointer :: vectornd(:,:,:,:,:) => null()
425 : ! vectornd(n4,n5,n6,nvloc,3)
426 : ! vector potential of nuclear magnetic dipoles
427 : ! in real space, on the augmented fft grid
428 :
429 : real(dp), contiguous, pointer :: vlocal(:,:,:,:) => null()
430 : ! vlocal(n4,n5,n6,nvloc)
431 : ! local potential in real space, on the augmented fft grid
432 :
433 : real(dp), contiguous, pointer :: vxctaulocal(:,:,:,:,:) => null()
434 : ! vxctaulocal(n4,n5,n6,nvloc,4)
435 : ! derivative of XC energy density with respect to kinetic energy density,
436 : ! in real space, on the augmented fft grid
437 :
438 : real(dp), contiguous, pointer :: xred(:,:) => null()
439 : ! xred(3,natom)
440 : ! reduced coordinates of atoms (dimensionless)
441 :
442 : real(dp), allocatable :: fofr_work(:,:,:,:)
443 : ! (2,n4,n5,n6,ndat)
444 : ! Buffer used in getgh1c when calling fourwf to compute <r|vlocal1|u_nk> for n =1, ndat.
445 : ! It is automatically allocated/reallocated by alloc_fofr according on the input ndat.
446 : ! In the case of gpu_option = 2, this is the buffer that is mapped to the GPU but only when ndat changes.
447 :
448 : ! ===== Structured datatype pointers
449 :
450 : type(fock_common_type), pointer :: fockcommon => null()
451 : ! common quantities needed to calculate Fock exact exchange
452 :
453 : type(fock_BZ_type), pointer :: fockbz => null()
454 : ! total brillouin zone quantities needed to calculate Fock exact exchange
455 :
456 : type(fock_ACE_type), pointer :: fockACE_k => null()
457 : ! ACE quantities needed to calculate Fock exact exchange in the ACE context
458 :
459 : contains
460 :
461 : procedure :: init => gsham_init
462 : ! Initialize the GS Hamiltonian
463 :
464 : procedure :: free => gsham_free
465 : ! Free the memory in the GS Hamiltonian
466 :
467 : procedure :: load_spin => gsham_load_spin
468 : ! Setup of the spin-dependent part of the GS Hamiltonian
469 :
470 : procedure :: load_k => gsham_load_k
471 : ! Setup of the k-dependent part of the GS Hamiltonian
472 :
473 : procedure :: load_kprime => gsham_load_kprime
474 : ! Setup of the k^prime-dependent part of the GS Hamiltonian
475 :
476 : procedure :: eph_setup_k => gsham_eph_setup_k
477 : ! Simplified interface to load either k or kprime in the case of e-ph calculations.
478 :
479 : procedure :: copy => gsham_copy
480 : ! Copy the object
481 :
482 : procedure :: print => gsham_print
483 : ! Print the object
484 :
485 : procedure :: alloc_fofr_work => gsham_alloc_fofr_work
486 : ! Allocate work space array before calling fourwf for ndat bands and map it to GPU
487 :
488 : end type gs_hamiltonian_type
489 : !!***
490 :
491 : !----------------------------------------------------------------------
492 :
493 : !!****t* m_hamiltonian/rf_hamiltonian_type
494 : !! NAME
495 : !! rf_hamiltonian_type
496 : !!
497 : !! FUNCTION
498 : !! This datastructure contains few data about one 1st-order Hamiltonian,
499 : !! needed in the "getgh1c" routine, that apply the 1st-order Hamiltonian
500 : !! on a wavefunction.
501 : !!
502 : !! SOURCE
503 :
504 : type,public :: rf_hamiltonian_type
505 :
506 : ! ===== Integer scalars
507 :
508 : integer :: cplex = -1
509 : ! if 1, real space 1-order functions on FFT grid are REAL; if 2, COMPLEX
510 :
511 : integer :: dime1kb1 = -1
512 : ! First dimension of E1kb, derivative of Ekb with respect to a perturbation
513 :
514 : integer :: dime1kb2 = -1
515 : ! Second dimension of E1kb, derivative of Ekb with respect to a perturbation
516 : ! NCPP: dime1kb2=ntypat, PAW: dime1kb2=natom
517 :
518 : integer :: npw_k = -1
519 : ! number of plane waves at k
520 :
521 : integer :: npw_kp = -1
522 : ! number of plane waves at k^prime
523 :
524 : integer:: nspinor = -1
525 : ! Number of spinorial components
526 :
527 : integer :: nsppol = -1
528 : ! Total number of spin components (1=non-polarized, 2=polarized)
529 :
530 : integer :: nvloc = -1
531 : ! Number of components of vloc
532 : ! usually, nvloc=1, except in the non-collinear magnetism case, where nvloc=4
533 :
534 : integer :: n4 = -1, n5 = -1, n6 = -1
535 : ! same as ngfft(4:6)
536 :
537 : integer :: zora = 0
538 : ! determines zora term use. 0 means nonrelativistic, 1 means use zora terms
539 : ! where available (currently only nuclear dipole moment terms)
540 :
541 : ! ===== Real arrays
542 :
543 : real(dp), allocatable :: e1kbfr_spin(:,:,:,:,:)
544 : ! e1kbfr_spin(dimekb1,dimekb2,nspinor**2,cplex,my_nsppol)
545 : ! Contains the values of e1kbfr array for all spins treated by current process
546 : ! See e1kbfr description; e1kbfr is pointer to e1kbfr_spin(:,:,:,:,isppol)
547 :
548 : real(dp), allocatable :: e1kbsc_spin(:,:,:,:,:)
549 : ! e1kbsc_spin(dimekb1,dimekb2,nspinor**2,cplex,my_nsppol)
550 : ! Contains the values of e1kbsc array for all spins treated by current process
551 : ! See e1kbsc description; e1kbsc is pointer to e1kbsc_spin(:,:,:,:,isppol)
552 :
553 : ! ===== Real pointers
554 :
555 : real(dp), contiguous, pointer :: dkinpw_k(:) => null()
556 : ! dkinpw_k(npw_k)
557 : ! 1st derivative of the (modified) kinetic energy for each plane wave at k
558 :
559 : real(dp), contiguous, pointer :: dkinpw_kp(:) => null()
560 : ! dkinpw_kp(npw_kp)
561 : ! 1st derivative of the (modified) kinetic energy for each plane wave at k^prime
562 :
563 : real(dp), contiguous, pointer :: ddkinpw_k(:) => null()
564 : ! ddkinpw_k(npw_k)
565 : ! 2nd derivative of the (modified) kinetic energy for each plane wave at k
566 :
567 : real(dp), contiguous, pointer :: ddkinpw_kp(:) => null()
568 : ! ddkinpw_kp(npw_kp)
569 : ! 2nd derivative of the (modified) kinetic energy for each plane wave at k^prime
570 :
571 : real(dp), contiguous, pointer :: e1kbfr(:,:,:,:) => null()
572 : ! Frozen part of 1st derivative of ekb for the considered perturbation
573 : ! (part not depending on VHxc^(1))
574 : ! e1kbfr(dime1kb1,dime1kb2,nspinor**2,cplex)
575 : ! For each spin component, e1kbfr points to e1kbfr_spin(:,:,:,:,my_isppol)
576 :
577 : real(dp), contiguous, pointer :: e1kbsc(:,:,:,:) => null()
578 : ! Self-consistent 1st derivative of ekb for the considered perturbation
579 : ! (part depending only on self-consistent VHxc^(1))
580 : ! e1kbsc(dime1kb1,dime1kb2,nspinor**2,cplex)
581 : ! For each spin component, e1kbfr points to e1kbfr_spin(:,:,:,:,my_isppol)
582 :
583 : real(dp), contiguous, pointer :: vectornd(:,:,:,:) => null()
584 : ! vectornd(n4,n5,n6,nvloc)
585 : ! vector potential of nuclear magnetic dipoles
586 : ! in real space, on the augmented fft grid, in direction idir
587 : ! (the ddk pert direction)
588 :
589 : real(dp), contiguous, pointer :: vlocal1(:,:,:,:) => null()
590 : ! vlocal1(cplex*n4,n5,n6,nvloc)
591 : ! 1st-order local potential in real space, on the augmented fft grid
592 :
593 : real(dp), contiguous, pointer :: vxctaulocal(:,:,:,:,:) => null()
594 : ! vxctaulocal(n4,n5,n6,nvloc,4)
595 : ! derivative of XC energy density with respect to kinetic energy density,
596 : ! in real space, on the augmented fft grid
597 :
598 : contains
599 : procedure :: init => rfham_init
600 : ! Initialize the RF Hamiltonian
601 :
602 : procedure :: free => rfham_free
603 : ! Free the memory in the RF Hamiltonian
604 :
605 : procedure :: load_spin => rfham_load_spin
606 : ! Setup of the spin-dependent part of the RF Hamiltonian.
607 :
608 : procedure :: load_k => rfham_load_k
609 : ! Setup of the k-dependent part of the RF Hamiltonian
610 :
611 : end type rf_hamiltonian_type
612 : !!***
613 :
614 : contains !===========================================================
615 :
616 : !----------------------------------------------------------------------
617 :
618 : !!****f* m_hamiltonian/gsham_free
619 : !! NAME
620 : !! gsham_free
621 : !!
622 : !! FUNCTION
623 : !! Clean and destroy gs_hamiltonian_type datastructure
624 : !!
625 : !! SOURCE
626 :
627 79904 : subroutine gsham_free(Ham)
628 :
629 : !Arguments ------------------------------------
630 : class(gs_hamiltonian_type),target,intent(inout) :: Ham
631 :
632 : !Local variables-------------------------------
633 79904 : real(dp), contiguous, pointer :: fofr_work_ptr(:,:,:,:)
634 : ! *************************************************************************
635 :
636 : DBG_ENTER("COLL")
637 :
638 : !@gs_hamiltonian_type
639 :
640 : ! Integer Pointers
641 79904 : if (associated(Ham%gbound_kp,Ham%gbound_k)) then
642 59736 : nullify(Ham%gbound_kp)
643 20168 : else if (associated(Ham%gbound_kp)) then
644 19771 : ABI_FREE(Ham%gbound_kp)
645 : end if
646 :
647 : ! Integer arrays
648 79904 : if(Ham%gpu_option == ABI_GPU_KOKKOS) then
649 : #if defined HAVE_GPU && defined HAVE_YAKL
650 : ABI_SFREE_MANAGED(Ham%atindx)
651 : ABI_SFREE_MANAGED(Ham%atindx1)
652 : ABI_SFREE_MANAGED(Ham%typat)
653 : ABI_SFREE_MANAGED(Ham%indlmn)
654 : ABI_SFREE_MANAGED(Ham%nattyp)
655 : #endif
656 : else
657 79904 : ABI_FREE(Ham%atindx)
658 79904 : ABI_FREE(Ham%atindx1)
659 79904 : ABI_FREE(Ham%typat)
660 79904 : ABI_FREE(Ham%indlmn)
661 79904 : ABI_FREE(Ham%nattyp)
662 : end if
663 79904 : ABI_SFREE(Ham%gbound_k)
664 79904 : ABI_SFREE(Ham%pspso)
665 79904 : ABI_SFREE(Ham%dimcprj)
666 :
667 : ! Real Pointers
668 79904 : if (associated(Ham%phkpxred,Ham%phkxred)) then
669 56541 : nullify(Ham%phkpxred)
670 23363 : else if (associated(Ham%phkpxred)) then
671 22966 : ABI_FREE(Ham%phkpxred)
672 : end if
673 79904 : ABI_SFREE(Ham%phkxred)
674 79904 : if (associated(Ham%ekb)) nullify(Ham%ekb)
675 79904 : if (associated(Ham%vectornd)) nullify(Ham%vectornd)
676 79904 : if (associated(Ham%vlocal)) nullify(Ham%vlocal)
677 79904 : if (associated(Ham%vxctaulocal)) nullify(Ham%vxctaulocal)
678 79904 : if (associated(Ham%xred)) nullify(Ham%xred)
679 79904 : if (associated(Ham%kinpw_k)) nullify(Ham%kinpw_k)
680 79904 : if (associated(Ham%kinpw_kp)) nullify(Ham%kinpw_kp)
681 79904 : if (associated(Ham%kg_k)) nullify(Ham%kg_k)
682 79904 : if (associated(Ham%kg_kp)) nullify(Ham%kg_kp)
683 79904 : if (associated(Ham%kpg_k)) nullify(Ham%kpg_k)
684 79904 : if (associated(Ham%kpg_kp)) nullify(Ham%kpg_kp)
685 79904 : if (associated(Ham%ffnl_k)) nullify(Ham%ffnl_k)
686 79904 : if (associated(Ham%ffnl_kp)) nullify(Ham%ffnl_kp)
687 79904 : if (associated(Ham%ph3d_k)) nullify(Ham%ph3d_k)
688 79904 : if (associated(Ham%ph3d_kp)) nullify(Ham%ph3d_kp)
689 :
690 :
691 : ! Real arrays
692 79904 : ABI_SFREE(Ham%ekb_spin)
693 79904 : ABI_SFREE(Ham%sij)
694 79904 : ABI_SFREE(Ham%nucdipmom)
695 79904 : if(Ham%gpu_option == ABI_GPU_KOKKOS) then
696 : #if defined HAVE_GPU && defined HAVE_YAKL
697 : ABI_SFREE_MANAGED(Ham%ph1d)
698 : #endif
699 : else
700 79904 : ABI_FREE(Ham%ph1d)
701 : end if
702 :
703 79904 : if(Ham%gpu_option==ABI_GPU_OPENMP) then
704 0 : call free_ompgpu_fourwf()
705 : end if
706 :
707 : ! Structured datatype pointers
708 79904 : if (associated(Ham%fockcommon)) nullify(Ham%fockcommon)
709 79904 : if (associated(Ham%fockACE_k)) nullify(Ham%fockACE_k)
710 79904 : if (associated(Ham%fockbz)) nullify(Ham%fockbz)
711 : #if defined HAVE_GPU_CUDA
712 : if(Ham%gpu_option==ABI_GPU_LEGACY .or. Ham%gpu_option==ABI_GPU_KOKKOS) then
713 : call gpu_finalize_ham_data()
714 : end if
715 : #endif
716 :
717 : if (Ham%gpu_option==ABI_GPU_OPENMP) then
718 : fofr_work_ptr => Ham%fofr_work
719 : #ifdef HAVE_OPENMP_OFFLOAD
720 : !$OMP TARGET EXIT DATA MAP(delete:fofr_work_ptr) IF (Ham%gpu_option==ABI_GPU_OPENMP)
721 : #endif
722 : end if
723 79904 : ABI_SFREE(Ham%fofr_work)
724 :
725 : DBG_EXIT("COLL")
726 :
727 79904 : end subroutine gsham_free
728 : !!***
729 :
730 : !----------------------------------------------------------------------
731 :
732 : !!****f* m_hamiltonian/gsham_init
733 : !! NAME
734 : !! gsham_init
735 : !!
736 : !! FUNCTION
737 : !! Creation method for the gs_hamiltonian_type structure.
738 : !! It allocates memory and initializes all quantities that do not depend on the k-point or spin.
739 : !!
740 : !! INPUTS
741 : !! [comm_atom]=optional, MPI communicator over atoms
742 : !! [fock <type(fock_type)>]= common quantities to calculate Fock exact exchange
743 : !! natom=Number of atoms in the unit cell.
744 : !! nfft=Number of FFT grid points (for this processors).
745 : !! nspinor=Number of spinorial components
746 : !! nsppol=1 for unpolarized, 2 for spin-polarized
747 : !! nspden=Number of spin density components.
748 : !! mgfft=Maximum size for 1D FFTs i.e., MAXVAL(ngfft(1:3))
749 : !! [mpi_atmtab(:)]=optional, indexes of the atoms treated by current proc
750 : !! [mpi_spintab(2)]=optional, flags defining the spin(s) treated be current process:
751 : !! mpi_spintab(1)=1 if non-polarized or spin-up treated
752 : !! mpi_spintab(2)=1 if polarized and spin-dn treated
753 : !! psps<pseudopotential_type>=structure datatype gathering data on the pseudopotentials.
754 : !! [electronpositron<electronpositron_type>]=Structured datatype storing data for the
755 : !! electron-positron two-component DFT (optional).
756 : !! ngfft(18)=integer array with FFT box dimensions and other information on FFTs, for the FINE rectangular grid.
757 : !! nloalg(3)=governs the choice of the algorithm for non-local operator
758 : !! [nucdipmom(3,natom)]= (optional) array of nuclear dipole moments at atomic sites
759 : !! [ph1d(2,3*(2*mgfft+1)*natom)]=1-dimensions phase arrays for structure factor (see getph.f).
760 : !! Optional, recalculated inside the routine if not present in input.
761 : !! rprimd(3,3)=Direct lattice vectors in Bohr.
762 : !! typat(natom)=Type of each atom.
763 : !! [usecprj]=flag use only for PAW; 1 if cprj datastructure is allocated
764 : !! [gpu_option] = GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
765 : !! xred(3,natom)=Reduced coordinates of the atoms.
766 : !! pawtab(ntypat*psps%usepaw)<pawtab_type>=PAW TABulated data initialized at start.
767 : !! [paw_ij(:) <type(paw_ij_type)>]=optional, paw arrays given on (i,j) channels
768 : !!
769 : !! SIDE EFFECTS
770 : !! Ham<gs_hamiltonian_type>=Structured datatype almost completely initialized:
771 : !! * Basic variables and dimensions are transferred to the structure.
772 : !! * All pointers are allocated with correct dimensions.
773 : !! * Quantities that do not depend on the k-point or spin are initialized.
774 : !!
775 : !! SOURCE
776 :
777 78546 : subroutine gsham_init(ham,Psps,pawtab,nspinor,nsppol,nspden,natom,typat,&
778 78546 : xred,nfft,mgfft,ngfft,rprimd,nloalg,&
779 78546 : ph1d,usecprj,comm_atom,mpi_atmtab,mpi_spintab,paw_ij,& ! optional
780 74533 : electronpositron,fock,nucdipmom,gpu_option,spinaxis,& ! optional
781 : use_gbt,zora,nfft_blocks) ! optional
782 :
783 : !Arguments ------------------------------------
784 : !scalars
785 : class(gs_hamiltonian_type),intent(inout),target :: ham
786 : integer,intent(in) :: nfft,natom,nspinor,nsppol,nspden,mgfft
787 : integer,optional,intent(in) :: comm_atom,usecprj,gpu_option,use_gbt,zora,nfft_blocks
788 : type(electronpositron_type),optional,pointer :: electronpositron
789 : type(fock_type),optional,pointer :: fock
790 : type(pseudopotential_type),intent(in) :: psps
791 : !arrays
792 : integer,intent(in) :: ngfft(18),nloalg(3),typat(natom)
793 : integer,optional,intent(in) :: mpi_atmtab(:),mpi_spintab(2)
794 : real(dp),intent(in) :: rprimd(3,3)
795 : real(dp),intent(in),target :: xred(3,natom)
796 : real(dp),optional,intent(in) :: nucdipmom(3,natom),ph1d(2,3*(2*mgfft+1)*natom),spinaxis(3)
797 : type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw)
798 : type(paw_ij_type),optional,intent(in) :: paw_ij(:)
799 :
800 : !Local variables-------------------------------
801 : !scalars
802 : integer :: my_comm_atom,my_nsppol,my_zora,itypat,iat,ilmn,indx,isp
803 : integer :: cplex_dij,jsp,l_gpu_option
804 : real(dp) :: ucvol
805 : !arrays
806 : integer :: my_spintab(2),ispin_gbt
807 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
808 78546 : real(dp),allocatable,target :: ekb_tmp(:,:,:,:)
809 : ! *************************************************************************
810 :
811 : DBG_ENTER("COLL")
812 :
813 : !@gs_hamiltonian_type
814 :
815 : !Manage optional parameters
816 78546 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
817 164394 : my_spintab=0;my_spintab(1:nsppol)=1;if (present(mpi_spintab)) my_spintab(1:2)=mpi_spintab(1:2)
818 235638 : my_nsppol=count(my_spintab==1)
819 78546 : l_gpu_option=ABI_GPU_DISABLED; if(present(gpu_option)) l_gpu_option=gpu_option
820 78546 : my_zora=0; if (present(zora)) my_zora=zora
821 :
822 78546 : ham%nfft_blocks=1;
823 78546 : if (present(nfft_blocks) .and. l_gpu_option==ABI_GPU_OPENMP) then
824 0 : ham%nfft_blocks=nfft_blocks
825 : end if
826 78546 : ham%use_gbt = 0; if (present(use_gbt)) ham%use_gbt = use_gbt
827 455163 : ham%spinaxis = zero; if (present(spinaxis)) ham%spinaxis = spinaxis
828 :
829 78546 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
830 :
831 314184 : ABI_CHECK_IEQ(mgfft, MAXVAL(ngfft(1:3)), "Wrong mgfft")
832 :
833 : !Allocate the arrays of the Hamiltonian whose dimensions do not depend on k
834 78546 : if(l_gpu_option == ABI_GPU_KOKKOS) then
835 : #if defined HAVE_GPU && defined HAVE_YAKL
836 : ABI_MALLOC_MANAGED(ham%atindx,(/natom/))
837 : ABI_MALLOC_MANAGED(ham%atindx1,(/natom/))
838 : ABI_MALLOC_MANAGED(ham%typat,(/natom/))
839 : ABI_MALLOC_MANAGED(ham%indlmn,(/6,psps%lmnmax,psps%ntypat/))
840 : ABI_MALLOC_MANAGED(ham%nattyp,(/psps%ntypat/))
841 : ABI_MALLOC_MANAGED(ham%ph1d,(/2,3*(2*mgfft+1)*natom/))
842 : #endif
843 : else
844 235638 : ABI_MALLOC(ham%atindx,(natom))
845 157092 : ABI_MALLOC(ham%atindx1,(natom))
846 157092 : ABI_MALLOC(ham%typat,(natom))
847 314184 : ABI_MALLOC(ham%indlmn,(6,psps%lmnmax,psps%ntypat))
848 235638 : ABI_MALLOC(ham%nattyp,(psps%ntypat))
849 235638 : ABI_MALLOC(ham%ph1d,(2,3*(2*mgfft+1)*natom))
850 : end if
851 :
852 235638 : ABI_MALLOC(ham%pspso,(psps%ntypat))
853 235638 : ABI_MALLOC(ham%nucdipmom,(3,natom))
854 :
855 : !Initialize most of the Hamiltonian
856 78546 : indx=1
857 192385 : do itypat=1,psps%ntypat
858 113839 : ham%nattyp(itypat)=0
859 495752 : do iat=1,natom
860 417206 : if (typat(iat)==itypat) then
861 202828 : ham%atindx (iat )=indx
862 202828 : ham%atindx1(indx)=iat
863 202828 : indx=indx+1
864 202828 : ham%nattyp(itypat)=ham%nattyp(itypat)+1
865 : end if
866 : end do
867 : end do
868 :
869 1021098 : ham%gmet(:,:) =gmet(:,:)
870 1021098 : ham%gprimd(:,:)=gprimd(:,:)
871 4606487 : ham%indlmn(:,:,:)=psps%indlmn(:,:,:)
872 78546 : ham%ispin_gbt =ispin_gbt
873 78546 : ham%lmnmax =psps%lmnmax
874 78546 : ham%mgfft =mgfft
875 78546 : ham%mpsang =psps%mpsang
876 78546 : ham%mpssoang =psps%mpssoang
877 78546 : ham%natom =natom
878 78546 : ham%nfft =nfft
879 1492374 : ham%ngfft(:) =ngfft(:)
880 314184 : ham%nloalg(:) =nloalg(:)
881 192385 : ham%matblk=min(NLO_MINCAT,maxval(ham%nattyp)); if (nloalg(2)>0) ham%matblk=natom
882 78546 : ham%nsppol =nsppol
883 78546 : ham%nspinor =nspinor
884 78546 : ham%ntypat =psps%ntypat
885 359920 : ham%typat =typat(1:natom)
886 78546 : ham%nvloc=1; if(nspden==4)ham%nvloc=4
887 78546 : ham%n4 =ngfft(4)
888 78546 : ham%n5 =ngfft(5)
889 78546 : ham%n6 =ngfft(6)
890 78546 : ham%usepaw =psps%usepaw
891 78546 : ham%ucvol =ucvol
892 78546 : ham%useylm =psps%useylm
893 78546 : ham%zora =my_zora
894 78546 : ham%gpu_option=ABI_GPU_DISABLED ; if(PRESENT(gpu_option)) ham%gpu_option=gpu_option
895 :
896 192385 : ham%pspso(:) =psps%pspso(1:psps%ntypat)
897 78546 : if (psps%usepaw==1) then
898 55581 : do itypat=1,psps%ntypat
899 55581 : ham%pspso(itypat)=1+pawtab(itypat)%usespnorb
900 : end do
901 : end if
902 :
903 78546 : if (present(nucdipmom)) then
904 851409 : ham%nucdipmom(:,:) = nucdipmom(:,:)
905 : else
906 38449 : ham%nucdipmom(:,:) = zero
907 : end if
908 :
909 78546 : ham%xred => xred
910 :
911 78546 : if (present(fock)) then
912 54521 : if (associated(fock)) then
913 2735 : ham%fockcommon => fock%fock_common
914 2735 : if (fock%fock_common%use_ACE==0) ham%fockbz=>fock%fock_BZ
915 : end if
916 : end if
917 :
918 78546 : if (present(ph1d)) then
919 90273459 : ham%ph1d(:,:) = ph1d(:,:)
920 : else ! Recalculate structure factor phases
921 252 : call getph(ham%atindx,natom,ngfft(1),ngfft(2),ngfft(3),ham%ph1d,xred)
922 : end if
923 :
924 78546 : if (ham%usepaw==1) then
925 21946 : ham%usecprj=0;if (present(usecprj)) ham%usecprj=usecprj
926 65838 : ABI_MALLOC(ham%dimcprj,(natom))
927 : !Be careful cprj are ordered by atom type (used in non-local operator)
928 21946 : call pawcprj_getdim(ham%dimcprj,natom,ham%nattyp,ham%ntypat,ham%typat,pawtab,'O')
929 : else
930 56600 : ham%usecprj=0
931 56600 : ABI_MALLOC(ham%dimcprj,(0))
932 : end if
933 :
934 : ! ===========================
935 : ! ==== Non-local factors ====
936 : ! ===========================
937 :
938 78546 : if (ham%usepaw==0) then
939 : ! Norm-conserving: use constant Kleimann-Bylander energies.
940 : ! nspinor ** 2 is a fake dimension here in the sense that
941 : ! the KB energies for the scalar part and the SOC part are packed in the firs dimension (dimekb).
942 : ! In nonlop_pl, ekb are accessed using ekb(iln,itypat,ispinor) where iln runs over all projects (scalar + SOC)
943 : ! The ispinor index is irrelevant as ekb(:,:,1) = ekb(:,:,2). See nonlop_pl
944 56600 : ham%dimekb1=psps%dimekb
945 56600 : ham%dimekb2=psps%ntypat
946 56600 : ham%dimekbq=1
947 283000 : ABI_MALLOC(ham%ekb_spin,(psps%dimekb,psps%ntypat,nspinor**2,1,1))
948 56600 : ham%ekb => ham%ekb_spin(:,:,:,:,1)
949 56600 : ABI_MALLOC(ham%sij,(0,0))
950 381396 : ham%ekb(:,:,1,1)=psps%ekb(:,:)
951 56600 : if (nspinor==2) then
952 36208 : ham%ekb(:,:,2,1)=psps%ekb(:,:)
953 75085 : ham%ekb(:,:,3:4,1)=zero
954 : end if
955 56600 : if (PRESENT(electronpositron)) then
956 36796 : if (electronpositron_calctype(electronpositron)==1) ham%ekb(:,:,:,:)=-ham%ekb(:,:,:,:)
957 : end if
958 :
959 : ! Update enl on GPU (will do it later for PAW)
960 : #if defined HAVE_GPU_CUDA
961 : if (ham%gpu_option==ABI_GPU_LEGACY .or. ham%gpu_option==ABI_GPU_KOKKOS) then
962 : call gpu_update_ham_data(&
963 : & ham%ekb(:,:,:,1), INT(size(ham%ekb), c_int64_t), &
964 : & ham%sij, INT(size(ham%sij), c_int64_t), &
965 : & ham%gprimd, INT(size(ham%gprimd),c_int64_t))
966 : end if
967 : #endif
968 :
969 : else
970 : ! PAW: store overlap coefficients (spin non dependent) and Dij coefficients (spin dependent)
971 21946 : cplex_dij=1
972 21946 : if (present(paw_ij)) then
973 21946 : if (size(paw_ij)>0) cplex_dij=paw_ij(1)%cplex_dij
974 : end if
975 264436 : if ((nspinor==2).or.any(abs(ham%nucdipmom)>tol8)) cplex_dij=2
976 21946 : ham%dimekb1=psps%dimekb*cplex_dij
977 21946 : ham%dimekb2=natom
978 21946 : ham%dimekbq=1
979 21946 : if (present(paw_ij)) then
980 21946 : if (size(paw_ij)>0) ham%dimekbq=paw_ij(1)%qphase
981 : end if
982 87784 : ABI_MALLOC(ham%sij,(ham%dimekb1,psps%ntypat))
983 55581 : do itypat=1,psps%ntypat
984 33635 : if (cplex_dij==1) then
985 1653135 : ham%sij(1:pawtab(itypat)%lmn2_size,itypat)=pawtab(itypat)%sij(:)
986 : else
987 403365 : do ilmn=1,pawtab(itypat)%lmn2_size
988 398470 : ham%sij(2*ilmn-1,itypat)=pawtab(itypat)%sij(ilmn)
989 403365 : ham%sij(2*ilmn ,itypat)=zero
990 : end do
991 : end if
992 55581 : if (cplex_dij*pawtab(itypat)%lmn2_size<ham%dimekb1) then
993 221815 : ham%sij(cplex_dij*pawtab(itypat)%lmn2_size+1:ham%dimekb1,itypat)=zero
994 : end if
995 : end do
996 : ! We preload here PAW non-local factors in order to avoid a communication over atoms
997 : ! inside the loop over spins.
998 153622 : ABI_MALLOC(ham%ekb_spin,(ham%dimekb1,ham%dimekb2,nspinor**2,ham%dimekbq,my_nsppol))
999 8080145 : ham%ekb_spin=zero
1000 21946 : if (present(paw_ij)) then
1001 21946 : if (my_nsppol<ham%nsppol) then
1002 3168 : ABI_MALLOC(ekb_tmp,(ham%dimekb1,ham%dimekb2,nspinor**2,ham%dimekbq))
1003 : end if
1004 21946 : jsp=0
1005 46594 : do isp=1,ham%nsppol
1006 24648 : if (my_spintab(isp)==1) then
1007 24092 : jsp=jsp+1 ; ham%ekb => ham%ekb_spin(:,:,:,:,jsp)
1008 : else
1009 556 : ham%ekb => ekb_tmp
1010 : end if
1011 46594 : if (present(mpi_atmtab)) then
1012 7556 : call pawdij2ekb(ham%ekb,paw_ij,isp,my_comm_atom,mpi_atmtab=mpi_atmtab)
1013 : else
1014 17092 : call pawdij2ekb(ham%ekb,paw_ij,isp,my_comm_atom)
1015 : end if
1016 : end do
1017 21946 : if (my_nsppol<ham%nsppol) then
1018 528 : ABI_FREE(ekb_tmp)
1019 : end if
1020 : end if
1021 21946 : nullify(ham%ekb)
1022 : end if
1023 :
1024 : DBG_EXIT("COLL")
1025 :
1026 231625 : end subroutine gsham_init
1027 : !!***
1028 :
1029 : !!****f* m_hamiltonian/gsham_load_k
1030 : !! NAME
1031 : !! gsham_load_k
1032 : !!
1033 : !! FUNCTION
1034 : !! Setup of the k-dependent part of the Hamiltonian H_k_k^prime
1035 : !!
1036 : !! INPUTS
1037 : !! [compute_gbound]=flag. if true, G sphere boundary is computed here
1038 : !! [compute_ph3d]=flag. if true, 3D structure factors are computed here (only if nloalg(1)>0)
1039 : !! [gbound_k]=G sphere boundary (not compatible with compute_gbound=TRUE)
1040 : !! [ffnl_k]=nonlocal form factors on basis sphere
1041 : !! [istwf_k]=parameter that describes the storage of wfs
1042 : !! [kinpw_k]=(modified) kinetic energy for each plane wave
1043 : !! [kg_k]=planewave reduced coordinates in basis sphere (g vectors)
1044 : !! [kpg_k]=(k+g) vectors in reciprocal space
1045 : !! [kpt_k]=k point coordinates
1046 : !! [npw_k]=number of plane waves (processed by current proc when band-FFT parallelism is on)
1047 : !! [npw_fft_k]=number of plane waves used to apply Hamiltonian (in the "FFT" configuration)
1048 : !! [ph3d_k]=3-dim structure factors, for each atom and plane wave
1049 : !!
1050 : !! SIDE EFFECTS
1051 : !! ham<gs_hamiltonian_type>=structured datatype completed with k-dependent quantities.
1052 : !! Quantities at k^prime are set equal to quantities at k.
1053 : !! k-dependent scalars and pointers associated
1054 : !! phkxred=exp(.k.xred) for each atom
1055 : !! [ham%gbound_k]=G sphere boundary, for each plane wave
1056 : !! [ham%ph3d_k]=3-dim structure factors, for each atom and plane wave
1057 : !!
1058 : !! SOURCE
1059 :
1060 9502854 : subroutine gsham_load_k(ham,ffnl_k,fockACE_k,gbound_k,istwf_k,kinpw_k,&
1061 6335236 : kg_k,kpg_k,kpt_k,npw_k,npw_fft_k,ph3d_k,&
1062 : compute_gbound,compute_ph3d)
1063 :
1064 : !Arguments ------------------------------------
1065 : !scalars
1066 : class(gs_hamiltonian_type),intent(inout),target :: ham
1067 : integer,intent(in),optional :: npw_k,npw_fft_k,istwf_k
1068 : logical,intent(in),optional :: compute_gbound,compute_ph3d
1069 : !arrays
1070 : integer,intent(in),optional,target :: gbound_k(:,:),kg_k(:,:)
1071 : real(dp),intent(in),optional :: kpt_k(3)
1072 : real(dp),intent(in),optional,target :: ffnl_k(:,:,:,:),kinpw_k(:),kpg_k(:,:),ph3d_k(:,:,:)
1073 : type(fock_ACE_type),intent(in),optional,target :: fockACE_k
1074 :
1075 : !Local variables-------------------------------
1076 : !scalars
1077 : integer :: iat,iatom
1078 : logical :: compute_gbound_
1079 : real(dp) :: arg
1080 : !character(len=500) :: msg
1081 : ! *************************************************************************
1082 :
1083 : DBG_ENTER("COLL")
1084 :
1085 : !@gs_hamiltonian_type
1086 :
1087 : !k-dependent scalars
1088 3167618 : if (present(kpt_k)) then
1089 7426540 : ham%kpt_k(:) = kpt_k(:)
1090 7426540 : ham%kpt_kp(:) = kpt_k(:)
1091 : end if
1092 3167618 : if (present(istwf_k)) then
1093 1856635 : ham%istwf_k = istwf_k
1094 1856635 : ham%istwf_kp = istwf_k
1095 : end if
1096 3167618 : if (present(npw_k)) then
1097 1856635 : ham%npw_k = npw_k
1098 1856635 : ham%npw_kp = npw_k
1099 : end if
1100 3167618 : if (present(npw_fft_k)) then
1101 46150 : ham%npw_fft_k = npw_fft_k
1102 46150 : ham%npw_fft_kp = npw_fft_k
1103 3121468 : else if (present(npw_k)) then
1104 1856635 : ham%npw_fft_k = npw_k
1105 1856635 : ham%npw_fft_kp = npw_k
1106 : end if
1107 :
1108 : !Pointers to k-dependent quantities
1109 3167618 : if (present(kinpw_k)) then
1110 493598 : ham%kinpw_k => kinpw_k
1111 493598 : ham%kinpw_kp => kinpw_k
1112 : end if
1113 3167618 : if (present(kg_k)) then
1114 1902793 : ham%kg_k => kg_k
1115 1902793 : ham%kg_kp => kg_k
1116 : end if
1117 3167618 : if (present(kpg_k)) then
1118 1892463 : ham%kpg_k => kpg_k
1119 1892463 : ham%kpg_kp => kpg_k
1120 : end if
1121 3167618 : if (present(ffnl_k)) then
1122 1961661 : ham%ffnl_k => ffnl_k
1123 1961661 : ham%ffnl_kp => ffnl_k
1124 : end if
1125 3167618 : if (present(ph3d_k)) then
1126 1782852 : ham%ph3d_k => ph3d_k
1127 1782852 : ham%ph3d_kp => ph3d_k
1128 : end if
1129 3167618 : if (present(fockACE_k)) then
1130 6063 : ham%fockACE_k => fockACE_k
1131 : end if
1132 : !Compute exp(i.k.R) for each atom
1133 3167618 : if (present(kpt_k)) then
1134 1856635 : if (associated(Ham%phkpxred).and.(.not.associated(Ham%phkpxred,Ham%phkxred))) then
1135 1355976 : ABI_FREE(Ham%phkpxred)
1136 : end if
1137 1856635 : ABI_SFREE(ham%phkxred)
1138 5569905 : ABI_MALLOC(ham%phkxred,(2,ham%natom))
1139 5440282 : do iat=1,ham%natom
1140 3583647 : iatom=ham%atindx(iat)
1141 14334588 : arg=two_pi*DOT_PRODUCT(kpt_k,ham%xred(:,iat))
1142 3583647 : ham%phkxred(1,iatom)=DCOS(arg)
1143 5440282 : ham%phkxred(2,iatom)=DSIN(arg)
1144 : end do
1145 1856635 : ham%phkpxred => ham%phkxred
1146 : end if
1147 :
1148 : ! Compute or copy G sphere boundary at k+g
1149 3167618 : compute_gbound_=.false.;if (present(compute_gbound)) compute_gbound_=compute_gbound
1150 3167618 : if (present(gbound_k)) compute_gbound_=.true.
1151 3121468 : if (compute_gbound_) then
1152 1681283 : if (associated(Ham%gbound_kp,Ham%gbound_k)) then
1153 365649 : nullify(Ham%gbound_kp)
1154 1315634 : else if (associated(Ham%gbound_kp)) then
1155 1253670 : ABI_FREE(Ham%gbound_kp)
1156 : end if
1157 1681283 : ABI_SFREE(ham%gbound_k)
1158 : end if
1159 3167618 : if (.not.allocated(ham%gbound_k)) then
1160 6789872 : ABI_MALLOC(ham%gbound_k,(2*ham%mgfft+8,2))
1161 138033904 : ham%gbound_k(:,:)=0
1162 1697468 : ham%gbound_kp => ham%gbound_k
1163 : end if
1164 3167618 : if (compute_gbound_) then
1165 1681283 : if (present(gbound_k)) then
1166 4357606 : ham%gbound_k(:,:)=gbound_k(:,:)
1167 : else
1168 1635133 : if (.not.associated(ham%kg_k)) then
1169 0 : ABI_BUG('Something is missing for gbound_k computation!')
1170 : end if
1171 : !write(std_out,*)"About to call sphereboundary"
1172 : !write(std_out,*)"size(kg_k), npw_k, mgfft",size(ham%kg_k, dim=2), ham%npw_k, ham%mgfft
1173 1635133 : call sphereboundary(ham%gbound_k,ham%istwf_k,ham%kg_k,ham%mgfft,ham%npw_k)
1174 : end if
1175 1681283 : ham%gbound_kp => ham%gbound_k
1176 : end if
1177 :
1178 : ! Compute 3D structure factors for each atom at k+g
1179 3167618 : if (present(compute_ph3d).and.present(ph3d_k)) then
1180 1736019 : if (compute_ph3d.and.ham%nloalg(2)>0) then
1181 1703155 : if ((.not.allocated(ham%phkxred)).or.(.not.associated(ham%kg_k)).or.&
1182 : (.not.associated(ham%ph3d_k))) then
1183 0 : ABI_BUG('Something is missing for ph3d_k computation!')
1184 : end if
1185 : call ph1d3d(1,ham%natom,ham%kg_k,ham%matblk,ham%natom,ham%npw_k,ham%ngfft(1),&
1186 1703155 : ham%ngfft(2),ham%ngfft(3),ham%phkxred,ham%ph1d,ham%ph3d_k)
1187 : end if
1188 : end if
1189 :
1190 : DBG_EXIT("COLL")
1191 :
1192 3167618 : end subroutine gsham_load_k
1193 : !!***
1194 :
1195 : !----------------------------------------------------------------------
1196 :
1197 : !!****f* m_hamiltonian/gsham_load_kprime
1198 : !! NAME
1199 : !! gsham_load_kprime
1200 : !!
1201 : !! FUNCTION
1202 : !! Setup of the k^prime-dependent part of the Hamiltonian H_k_k^prime
1203 : !!
1204 : !! INPUTS
1205 : !! [compute_gbound]=flag. if true, G sphere boundary is computed here
1206 : !! [compute_ph3d]=flag. if true, 3D structure factors are computed here (only if nloalg(2)>0)
1207 : !! [gbound_kp]=G sphere boundary (not compatible with compute_gbound=TRUE)
1208 : !! [ffnl_kp]=nonlocal form factors on basis sphere
1209 : !! [istwf_kp]=parameter that describes the storage of wfs
1210 : !! [kinpw_kp]=(modified) kinetic energy for each plane wave
1211 : !! [kg_kp]=planewave reduced coordinates in basis sphere (g vectors)
1212 : !! [kpg_kp]=(k+g) vectors in reciprocal space
1213 : !! [kpt_kp]=k point coordinates
1214 : !! [npw_kp]=number of plane waves (processed by current proc when band-FFT parallelism is on)
1215 : !! [npw_fft_kp]=number of plane waves used to apply Hamiltonian (in the "FFT" configuration)
1216 : !! [ph3d_kp]=3-dim structure factors, for each atom and plane wave
1217 : !!
1218 : !! SIDE EFFECTS
1219 : !! ham<gs_hamiltonian_type>=structured datatype completed with k^prime-dependent quantities.
1220 : !! k^prime-dependent scalars and pointers associated
1221 : !! phkpxred=exp(.k^prime.xred) for each atom
1222 : !! [ham%gbound_kp]=G sphere boundary, for each plane wave
1223 : !! [ham%ph3d_kp]=3-dim structure factors at k^prime at k, for each atom and plane wave
1224 : !!
1225 : !! SOURCE
1226 :
1227 6117102 : subroutine gsham_load_kprime(ham,ffnl_kp,gbound_kp,istwf_kp,kinpw_kp,&
1228 2039034 : kg_kp,kpg_kp,kpt_kp,npw_kp,npw_fft_kp,&
1229 2039034 : ph3d_kp,compute_gbound,compute_ph3d)
1230 :
1231 : !Arguments ------------------------------------
1232 : !scalars
1233 : class(gs_hamiltonian_type),intent(inout),target :: ham
1234 : integer,intent(in),optional :: npw_kp,npw_fft_kp,istwf_kp
1235 : logical,intent(in),optional :: compute_gbound,compute_ph3d
1236 : !arrays
1237 : integer,intent(in),optional,target :: gbound_kp(:,:),kg_kp(:,:)
1238 : real(dp),intent(in),optional :: kpt_kp(3)
1239 : real(dp),intent(in),optional,target :: ffnl_kp(:,:,:,:),kinpw_kp(:),kpg_kp(:,:),ph3d_kp(:,:,:)
1240 :
1241 : !Local variables-------------------------------
1242 : !scalars
1243 : integer :: iat,iatom
1244 : logical :: compute_gbound_
1245 : real(dp) :: arg
1246 : !character(len=500) :: msg
1247 : ! *************************************************************************
1248 :
1249 : DBG_ENTER("COLL")
1250 :
1251 : !@gs_hamiltonian_type
1252 :
1253 : !k-dependent scalars
1254 6863409 : if (present(kpt_kp)) ham%kpt_kp(:)= kpt_kp(:)
1255 2039034 : if (present(istwf_kp)) ham%istwf_kp = istwf_kp
1256 2039034 : if (present(npw_kp)) ham%npw_kp = npw_kp
1257 2039034 : if (present(npw_fft_kp)) then
1258 0 : ham%npw_fft_kp = npw_fft_kp
1259 2039034 : else if (present(npw_kp)) then
1260 1607285 : ham%npw_fft_kp = npw_kp
1261 : end if
1262 :
1263 : !Pointers to k-dependent quantities
1264 2039034 : if (present(kinpw_kp)) ham%kinpw_kp => kinpw_kp
1265 2039034 : if (present(kg_kp)) ham%kg_kp => kg_kp
1266 2039034 : if (present(kpg_kp)) ham%kpg_kp => kpg_kp
1267 2039034 : if (present(ffnl_kp)) ham%ffnl_kp => ffnl_kp
1268 2039034 : if (present(ph3d_kp)) ham%ph3d_kp => ph3d_kp
1269 :
1270 : !Compute exp(i.k^prime.R) for each atom
1271 2039034 : if (present(kpt_kp)) then
1272 1608125 : if (associated(ham%phkpxred,ham%phkxred)) then
1273 1378942 : nullify(ham%phkpxred)
1274 229183 : else if (associated(ham%phkpxred)) then
1275 229183 : ABI_FREE(ham%phkpxred)
1276 : end if
1277 4824375 : ABI_MALLOC(ham%phkpxred,(2,ham%natom))
1278 4740757 : do iat=1,ham%natom
1279 3132632 : iatom=ham%atindx(iat)
1280 12530528 : arg=two_pi*DOT_PRODUCT(kpt_kp,ham%xred(:,iat))
1281 3132632 : ham%phkpxred(1,iatom)=DCOS(arg)
1282 4740757 : ham%phkpxred(2,iatom)=DSIN(arg)
1283 : end do
1284 : end if
1285 :
1286 : !Compute or copy G sphere boundary at k^prime+g
1287 2039034 : compute_gbound_=.false.
1288 2039034 : if (present(kpt_kp).and.present(compute_gbound)) compute_gbound_=compute_gbound
1289 2039034 : if (present(gbound_kp)) compute_gbound_=.true.
1290 2039034 : if (compute_gbound_) then
1291 1277209 : if (associated(ham%gbound_kp,ham%gbound_k)) then
1292 1273441 : nullify(ham%gbound_kp)
1293 3768 : else if (associated(ham%gbound_kp)) then
1294 3768 : ABI_FREE(ham%gbound_kp)
1295 : end if
1296 1277209 : if (present(gbound_kp)) then
1297 0 : ham%gbound_kp(:,:)=gbound_kp(:,:)
1298 : else
1299 1277209 : if (.not.associated(ham%kg_kp)) then
1300 0 : ABI_BUG('Something is missing for gbound_kp computation!')
1301 : end if
1302 5108836 : ABI_MALLOC(ham%gbound_kp,(2*ham%mgfft+8,2))
1303 1277209 : call sphereboundary(ham%gbound_kp,ham%istwf_kp,ham%kg_kp,ham%mgfft,ham%npw_kp)
1304 : end if
1305 : end if
1306 :
1307 : !Compute 3D structure factors for each atom at k^prime+g
1308 2039034 : if (present(compute_ph3d).and.present(ph3d_kp)) then
1309 603246 : if (compute_ph3d.and.ham%nloalg(2)>0) then
1310 585876 : if ((.not.associated(ham%phkpxred)).or.(.not.associated(ham%kg_kp)).or.&
1311 : & (.not.associated(ham%ph3d_kp))) then
1312 0 : ABI_BUG('Something is missing for ph3d_kp computation!')
1313 : end if
1314 : call ph1d3d(1,ham%natom,ham%kg_kp,ham%matblk,ham%natom,ham%npw_kp,ham%ngfft(1),&
1315 585876 : & ham%ngfft(2),ham%ngfft(3),ham%phkpxred,ham%ph1d,ham%ph3d_kp)
1316 : end if
1317 : end if
1318 :
1319 : DBG_EXIT("COLL")
1320 :
1321 2039034 : end subroutine gsham_load_kprime
1322 : !!***
1323 :
1324 : !----------------------------------------------------------------------
1325 :
1326 : !!****f* m_hamiltonian/gsham_eph_setup_k
1327 : !! NAME
1328 : !! gsham_eph_setup_k
1329 : !!
1330 : !! FUNCTION
1331 : !! Simplified interface to load either k or kprime in the case of e-ph calculations.
1332 : !!
1333 : !! INPUTS
1334 : !! which_k= "k" to load k, "kq" to load kprime
1335 : !! See load_k or load_kprime for the meaning of arguments.
1336 : !!
1337 : !! SOURCE
1338 :
1339 157025 : subroutine gsham_eph_setup_k(gs_ham, which_k, kk, istwf_k, npw_k, kg_k, dtset, cryst, psps, & ! in
1340 : nkpg_k, kpg_k, ffnl_k, kinpw_k, ph3d_k, comm) ! out
1341 :
1342 : !Arguments ------------------------------------
1343 : !scalars
1344 : class(gs_hamiltonian_type),intent(inout) :: gs_ham
1345 : character(len=*),intent(in) :: which_k
1346 : type(dataset_type),intent(in) :: dtset
1347 : type(crystal_t),intent(in) :: cryst
1348 : type(pseudopotential_type),intent(in) :: psps
1349 : integer,intent(in) :: istwf_k, npw_k, comm
1350 : integer,intent(out) :: nkpg_k
1351 : !arrays
1352 : real(dp),intent(in) :: kk(3)
1353 : integer,intent(in) :: kg_k(3,npw_k)
1354 : real(dp),allocatable,intent(out) :: kpg_k(:,:), ffnl_k(:,:,:,:), kinpw_k(:), ph3d_k(:,:,:)
1355 :
1356 : !Local variables-------------------------------
1357 : !scalars
1358 : integer,parameter :: ider0 = 0, idir0 = 0, optder0 = 0
1359 : !arrays
1360 : real(dp) :: ylmgr_k_dum(1,1,1)
1361 : real(dp),allocatable :: ylm_k(:,:)
1362 : ! *************************************************************************
1363 :
1364 : ! Compute k+G vectors
1365 157025 : nkpg_k = 3 * dtset%nloalg(3)
1366 628100 : ABI_MALLOC(kpg_k, (npw_k, nkpg_k))
1367 157025 : if (nkpg_k > 0) call mkkpg(kg_k, kpg_k, kk, nkpg_k, npw_k)
1368 :
1369 : ! Spherical Harmonics at k for useylm == 1.
1370 628100 : ABI_MALLOC(ylm_k, (npw_k, psps%mpsang**2 * psps%useylm))
1371 157025 : if (psps%useylm == 1) call initylmg_k(npw_k, psps%mpsang, optder0, cryst%rprimd, cryst%gprimd, kk, kg_k, ylm_k, ylmgr_k_dum)
1372 :
1373 : ! Compute nonlocal form factors ffnl_k at (k+G)
1374 785125 : ABI_MALLOC(ffnl_k, (npw_k, 1, psps%lmnmax, psps%ntypat))
1375 157025 : call mkffnl_objs(cryst, psps, 1, ffnl_k, ider0, idir0, kg_k, kpg_k, kk, nkpg_k, npw_k, ylm_k, ylmgr_k_dum, comm=comm)
1376 157025 : ABI_FREE(ylm_k)
1377 :
1378 : ! Compute (1/2) (2 Pi)**2 (kG)**2:
1379 134657935 : ABI_CALLOC(kinpw_k, (npw_k))
1380 157025 : call mkkin(dtset%ecut, dtset%ecutsm, dtset%effmass_free, cryst%gmet, kg_k, kinpw_k, kk, npw_k, 0, 0)
1381 :
1382 628100 : ABI_MALLOC(ph3d_k, (2, npw_k, gs_ham%matblk))
1383 :
1384 : ! Load the k dependent parts of the Hamiltonian
1385 76649 : select case (which_k)
1386 : case ("k")
1387 : call gs_ham%load_k(kpt_k=kk, npw_k=npw_k, istwf_k=istwf_k, kg_k=kg_k, kpg_k=kpg_k, kinpw_k=kinpw_k, &
1388 76649 : ph3d_k=ph3d_k, ffnl_k=ffnl_k, compute_ph3d=.true., compute_gbound=.true.)
1389 : case ("kq")
1390 : call gs_ham%load_kprime(kpt_kp=kk, npw_kp=npw_k, istwf_kp=istwf_k, kg_kp=kg_k, kpg_kp=kpg_k, kinpw_kp=kinpw_k, &
1391 80376 : ph3d_kp=ph3d_k, ffnl_kp=ffnl_k, compute_ph3d=.true., compute_gbound=.true.)
1392 : case default
1393 157025 : ABI_ERROR(sjoin("Invalid value for which_k:", which_k))
1394 : end select
1395 :
1396 157025 : end subroutine gsham_eph_setup_k
1397 : !!***
1398 :
1399 : !----------------------------------------------------------------------
1400 :
1401 : !!****f* m_hamiltonian/gsham_copy
1402 : !! NAME
1403 : !! gsham_copy
1404 : !!
1405 : !! INPUTS
1406 : !! gs_hamk_in<gs_hamiltonian_type>=Structured datatype completely initialized, to be copied.
1407 : !!
1408 : !! FUNCTION
1409 : !! Copy a gs_hamiltonian_type variable (gs_hamk_in) in another (gs_hamk_out).
1410 : !! In contrast to an assignment statement (gs_hamk_out=gs_hamk_in), this
1411 : !! subroutine allocate memory space for the pointers contained in the data
1412 : !! structure (gs_hamk_out) and copy the content of the corresponding memory
1413 : !! space of gs_hamk_in in it. In contrast, the assignment statement would
1414 : !! only associate the pointers of gs_hamk_out to the same memory space than
1415 : !! the corresponding ones in gs_hamk_in. This can cause trouble if one data
1416 : !! structure is destroyed before a reading/writing statement for the other
1417 : !! structure, causing access to unallocated memory space (silently, without
1418 : !! segmentation fault being generated).
1419 : !!
1420 : !! OUTPUT
1421 : !! gs_hamk_out<gs_hamiltonian_type>=Structured datatype containing separate
1422 : !! copies of all data of gs_hamk_in upon exit.
1423 : !!
1424 : !! SOURCE
1425 :
1426 1358 : subroutine gsham_copy(gs_hamk_in, gs_hamk_out)
1427 :
1428 : !Arguments ------------------------------------
1429 : class(gs_hamiltonian_type),intent(in),target :: gs_hamk_in
1430 : class(gs_hamiltonian_type),intent(out),target :: gs_hamk_out
1431 :
1432 : !Local variables-------------------------------
1433 : integer :: tmp2i(5)
1434 : type(C_PTR) :: ham_ptr
1435 : ! *************************************************************************
1436 :
1437 : DBG_ENTER("COLL")
1438 :
1439 : !@gs_hamiltonian_type
1440 :
1441 1358 : gs_hamk_out%dimekb1 = gs_hamk_in%dimekb1
1442 1358 : gs_hamk_out%dimekb2 = gs_hamk_in%dimekb2
1443 1358 : gs_hamk_out%dimekbq = gs_hamk_in%dimekbq
1444 1358 : gs_hamk_out%istwf_k = gs_hamk_in%istwf_k
1445 1358 : gs_hamk_out%istwf_kp = gs_hamk_in%istwf_kp
1446 1358 : gs_hamk_out%lmnmax = gs_hamk_in%lmnmax
1447 1358 : gs_hamk_out%matblk = gs_hamk_in%matblk
1448 1358 : gs_hamk_out%mgfft = gs_hamk_in%mgfft
1449 1358 : gs_hamk_out%mpsang = gs_hamk_in%mpsang
1450 1358 : gs_hamk_out%mpssoang = gs_hamk_in%mpssoang
1451 1358 : gs_hamk_out%natom = gs_hamk_in%natom
1452 1358 : gs_hamk_out%nfft = gs_hamk_in%nfft
1453 1358 : gs_hamk_out%nfft_blocks = gs_hamk_in%nfft_blocks
1454 1358 : gs_hamk_out%npw_k = gs_hamk_in%npw_k
1455 1358 : gs_hamk_out%npw_kp = gs_hamk_in%npw_kp
1456 1358 : gs_hamk_out%npw_fft_k = gs_hamk_in%npw_fft_k
1457 1358 : gs_hamk_out%npw_fft_kp = gs_hamk_in%npw_fft_kp
1458 1358 : gs_hamk_out%nspinor = gs_hamk_in%nspinor
1459 1358 : gs_hamk_out%nsppol = gs_hamk_in%nsppol
1460 1358 : gs_hamk_out%ntypat = gs_hamk_in%ntypat
1461 1358 : gs_hamk_out%nvloc = gs_hamk_in%nvloc
1462 1358 : gs_hamk_out%n4 = gs_hamk_in%n4
1463 1358 : gs_hamk_out%n5 = gs_hamk_in%n5
1464 1358 : gs_hamk_out%n6 = gs_hamk_in%n6
1465 1358 : gs_hamk_out%gpu_option = gs_hamk_in%gpu_option
1466 9506 : gs_hamk_out%spinaxis = gs_hamk_in%spinaxis
1467 1358 : gs_hamk_out%usecprj = gs_hamk_in%usecprj
1468 1358 : gs_hamk_out%usepaw = gs_hamk_in%usepaw
1469 1358 : gs_hamk_out%useylm = gs_hamk_in%useylm
1470 1358 : gs_hamk_out%use_gbt = gs_hamk_in%use_gbt
1471 1358 : gs_hamk_out%zora = gs_hamk_in%zora
1472 50246 : gs_hamk_out%ngfft = gs_hamk_in%ngfft
1473 9506 : gs_hamk_out%nloalg = gs_hamk_in%nloalg
1474 1358 : gs_hamk_out%ucvol = gs_hamk_in%ucvol
1475 33950 : gs_hamk_out%gmet = gs_hamk_in%gmet
1476 33950 : gs_hamk_out%gprimd = gs_hamk_in%gprimd
1477 1358 : gs_hamk_out%ispin_gbt = gs_hamk_in%ispin_gbt
1478 9506 : gs_hamk_out%kpt_k = gs_hamk_in%kpt_k
1479 9506 : gs_hamk_out%kpt_kp = gs_hamk_in%kpt_kp
1480 :
1481 4074 : ABI_MALLOC(gs_hamk_out%atindx,(gs_hamk_out%natom))
1482 7160 : gs_hamk_out%atindx = gs_hamk_in%atindx
1483 4074 : ABI_MALLOC(gs_hamk_out%atindx1,(gs_hamk_out%natom))
1484 7160 : gs_hamk_out%atindx1 = gs_hamk_in%atindx1
1485 4074 : ABI_MALLOC(gs_hamk_out%dimcprj,(gs_hamk_out%natom*gs_hamk_out%usepaw))
1486 7072 : if (gs_hamk_out%usepaw==1) gs_hamk_out%dimcprj = gs_hamk_in%dimcprj
1487 4074 : ABI_MALLOC(gs_hamk_out%typat,(gs_hamk_out%natom))
1488 7160 : gs_hamk_out%typat = gs_hamk_in%typat
1489 5432 : ABI_MALLOC(gs_hamk_out%gbound_k,(2*gs_hamk_out%mgfft+8,2))
1490 227380 : gs_hamk_out%gbound_k = gs_hamk_in%gbound_k
1491 5432 : ABI_MALLOC(gs_hamk_out%indlmn,(6,gs_hamk_out%lmnmax,gs_hamk_out%ntypat))
1492 252168 : gs_hamk_out%indlmn = gs_hamk_in%indlmn
1493 4074 : ABI_MALLOC(gs_hamk_out%nattyp,(gs_hamk_out%ntypat))
1494 7112 : gs_hamk_out%nattyp = gs_hamk_in%nattyp
1495 4074 : ABI_MALLOC(gs_hamk_out%nucdipmom,(3,gs_hamk_out%natom))
1496 20492 : gs_hamk_out%nucdipmom = gs_hamk_in%nucdipmom
1497 4074 : ABI_MALLOC(gs_hamk_out%phkxred,(2,gs_hamk_out%natom))
1498 16048 : gs_hamk_out%phkxred = gs_hamk_in%phkxred
1499 4074 : ABI_MALLOC(gs_hamk_out%ph1d,(2,3*(2*gs_hamk_out%mgfft+1)*gs_hamk_out%natom))
1500 1333672 : gs_hamk_out%ph1d = gs_hamk_in%ph1d
1501 4074 : ABI_MALLOC(gs_hamk_out%pspso,(gs_hamk_out%ntypat))
1502 7112 : gs_hamk_out%pspso = gs_hamk_in%pspso
1503 8148 : tmp2i(1:5)=shape(gs_hamk_in%ekb_spin)
1504 9506 : ABI_MALLOC(gs_hamk_out%ekb_spin,(tmp2i(1),tmp2i(2),tmp2i(3),tmp2i(4),tmp2i(5)))
1505 593500 : gs_hamk_out%ekb_spin = gs_hamk_in%ekb_spin
1506 1358 : gs_hamk_out%ekb => gs_hamk_out%ekb_spin(:,:,:,:,1)
1507 4074 : tmp2i(1:2)=shape(gs_hamk_in%sij)
1508 5432 : ABI_MALLOC(gs_hamk_out%sij,(tmp2i(1),tmp2i(2)))
1509 321288 : gs_hamk_out%sij = gs_hamk_in%sij
1510 :
1511 1358 : if (associated(gs_hamk_in%gbound_kp,gs_hamk_in%gbound_k)) then
1512 1358 : gs_hamk_out%gbound_kp => gs_hamk_out%gbound_k
1513 : else
1514 0 : ABI_MALLOC(gs_hamk_out%gbound_kp,(2,gs_hamk_out%natom))
1515 0 : gs_hamk_out%gbound_kp = gs_hamk_in%gbound_kp
1516 : end if
1517 1358 : if (associated(gs_hamk_in%phkpxred,gs_hamk_in%phkxred)) then
1518 1358 : gs_hamk_out%phkpxred => gs_hamk_out%phkxred
1519 : else
1520 0 : ABI_MALLOC(gs_hamk_out%phkpxred,(2,gs_hamk_out%natom))
1521 0 : gs_hamk_out%phkpxred = gs_hamk_in%phkpxred
1522 : end if
1523 :
1524 1358 : call addr_copy(gs_hamk_in%xred,gs_hamk_out%xred)
1525 1358 : call addr_copy(gs_hamk_in%vectornd,gs_hamk_out%vectornd)
1526 1358 : call addr_copy(gs_hamk_in%vlocal,gs_hamk_out%vlocal)
1527 1358 : call addr_copy(gs_hamk_in%vxctaulocal,gs_hamk_out%vxctaulocal)
1528 1358 : call addr_copy(gs_hamk_in%kinpw_k,gs_hamk_out%kinpw_k)
1529 1358 : call addr_copy(gs_hamk_in%kinpw_kp,gs_hamk_out%kinpw_kp)
1530 1358 : call addr_copy(gs_hamk_in%kg_k,gs_hamk_out%kg_k)
1531 1358 : call addr_copy(gs_hamk_in%kg_kp,gs_hamk_out%kg_kp)
1532 1358 : call addr_copy(gs_hamk_in%kpg_k,gs_hamk_out%kpg_k)
1533 1358 : call addr_copy(gs_hamk_in%kpg_kp,gs_hamk_out%kpg_kp)
1534 1358 : call addr_copy(gs_hamk_in%ffnl_k,gs_hamk_out%ffnl_k)
1535 1358 : call addr_copy(gs_hamk_in%ffnl_kp,gs_hamk_out%ffnl_kp)
1536 1358 : call addr_copy(gs_hamk_in%ph3d_k,gs_hamk_out%ph3d_k)
1537 1358 : call addr_copy(gs_hamk_in%ph3d_kp,gs_hamk_out%ph3d_kp)
1538 :
1539 : !For pointers to structured datatypes, have to copy the address
1540 : !manually because there is no generic addr_copy function for that
1541 1358 : if (associated(gs_hamk_in%fockcommon)) then
1542 : #if defined HAVE_FC_ISO_C_BINDING
1543 0 : ham_ptr=c_loc(gs_hamk_in%fockcommon)
1544 0 : call c_f_pointer(ham_ptr,gs_hamk_out%fockcommon)
1545 : #else
1546 : gs_hamk_out%fockcommon=transfer(gs_hamk_in%fockcommon,gs_hamk_out%fockcommon)
1547 : #endif
1548 : else
1549 1358 : nullify(gs_hamk_out%fockcommon)
1550 : end if
1551 1358 : if (associated(gs_hamk_in%fockbz)) then
1552 : #if defined HAVE_FC_ISO_C_BINDING
1553 0 : ham_ptr=c_loc(gs_hamk_in%fockbz)
1554 0 : call c_f_pointer(ham_ptr,gs_hamk_out%fockbz)
1555 : #else
1556 : gs_hamk_out%fockbz=transfer(gs_hamk_in%fockbz,gs_hamk_out%fockbz)
1557 : #endif
1558 : else
1559 1358 : nullify(gs_hamk_out%fockbz)
1560 : end if
1561 1358 : if (associated(gs_hamk_in%fockACE_k)) then
1562 : #if defined HAVE_FC_ISO_C_BINDING
1563 0 : ham_ptr=c_loc(gs_hamk_in%fockACE_k)
1564 0 : call c_f_pointer(ham_ptr,gs_hamk_out%fockACE_k)
1565 : #else
1566 : gs_hamk_out%fockACE_k=transfer(gs_hamk_in%fockACE_k,gs_hamk_out%fockACE_k)
1567 : #endif
1568 : else
1569 1358 : nullify(gs_hamk_out%fockACE_k)
1570 : end if
1571 :
1572 : DBG_EXIT("COLL")
1573 :
1574 1358 : end subroutine gsham_copy
1575 : !!***
1576 :
1577 : !----------------------------------------------------------------------
1578 :
1579 : !!****f* m_hamiltonian/gsham_load_spin
1580 : !! NAME
1581 : !! gsham_load_spin
1582 : !!
1583 : !! INPUTS
1584 : !! isppol=index of current spin
1585 : !! [vectornd(n4,n5,n6,nvloc,3)]=optional, vector potential of nuclear magnetic dipoles in real space
1586 : !! [vlocal(n4,n5,n6,nvloc)]=optional, local potential in real space
1587 : !! [vxctaulocal(n4,n5,n6,nvloc,4)]=optional, derivative of XC energy density with respect
1588 : !! to kinetic energy density in real space
1589 : !! [with_nonlocal]=optional, true if non-local factors have to be loaded
1590 : !!
1591 : !! FUNCTION
1592 : !! Setup of the spin-dependent part of the GS Hamiltonian.
1593 : !!
1594 : !! SIDE EFFECTS
1595 : !! Ham<gs_hamiltonian_type>=Structured datatype initialization phase:
1596 : !! * Quantities that depend spin are initialized.
1597 : !!
1598 : !! SOURCE
1599 :
1600 198902 : subroutine gsham_load_spin(Ham,isppol,vectornd,vlocal,vxctaulocal,with_nonlocal)
1601 :
1602 : !Arguments ------------------------------------
1603 : !scalars
1604 : class(gs_hamiltonian_type),intent(inout),target :: Ham
1605 : integer,intent(in) :: isppol
1606 : logical,optional,intent(in) :: with_nonlocal
1607 : !arrays
1608 : real(dp),optional,intent(in),target :: vectornd(:,:,:,:,:)
1609 : real(dp),optional,intent(in),target :: vlocal(:,:,:,:),vxctaulocal(:,:,:,:,:)
1610 :
1611 : !Local variables-------------------------------
1612 : integer :: jsppol
1613 : ! *************************************************************************
1614 :
1615 : DBG_ENTER("COLL")
1616 :
1617 : !@gs_hamiltonian_type
1618 198902 : if (present(vlocal)) then
1619 920175 : ABI_CHECK_IEQ(size(vlocal), Ham%n4*Ham%n5*Ham%n6*Ham%nvloc, "Wrong vlocal")
1620 184035 : Ham%vlocal => vlocal
1621 : end if
1622 198902 : if (present(vxctaulocal)) then
1623 4560 : ABI_CHECK_IEQ(size(vxctaulocal), Ham%n4*Ham%n5*Ham%n6*Ham%nvloc*4, "Wrong vxctaulocal")
1624 760 : Ham%vxctaulocal => vxctaulocal
1625 : end if
1626 198902 : if (present(vectornd)) then
1627 5274 : ABI_CHECK_IEQ(size(vectornd), Ham%n4*Ham%n5*Ham%n6*Ham%nvloc*3, "Wrong vectornd")
1628 879 : Ham%vectornd => vectornd
1629 : end if
1630 :
1631 : ! Retrieve non-local factors for this spin component
1632 198902 : if (present(with_nonlocal)) then
1633 197263 : if (with_nonlocal) then
1634 197263 : jsppol=min(isppol,size(Ham%ekb_spin,5))
1635 197263 : if (jsppol>0) Ham%ekb => Ham%ekb_spin(:,:,:,:,jsppol)
1636 : end if
1637 : end if
1638 :
1639 : ! Update enl and sij on GPU
1640 : #if defined HAVE_GPU_CUDA
1641 : if (Ham%gpu_option==ABI_GPU_LEGACY .or. Ham%gpu_option==ABI_GPU_KOKKOS) then
1642 : call gpu_update_ham_data(&
1643 : & Ham%ekb(:,:,:,1), INT(size(Ham%ekb), c_int64_t), &
1644 : & Ham%sij, INT(size(Ham%sij), c_int64_t), &
1645 : & Ham%gprimd, INT(size(Ham%gprimd),c_int64_t))
1646 : end if
1647 : #endif
1648 :
1649 : DBG_EXIT("COLL")
1650 :
1651 198902 : end subroutine gsham_load_spin
1652 : !!***
1653 :
1654 : !!****f* m_hamiltonian/gsham_alloc_fofr_work
1655 : !! NAME
1656 : !! gsham_alloc_fofr_work
1657 : !!
1658 : !! FUNCTION
1659 : !!
1660 : !! INPUTS
1661 : !!
1662 : !! SOURCE
1663 :
1664 0 : subroutine gsham_alloc_fofr_work(gs_ham, ndat)
1665 :
1666 : !Arguments ------------------------------------
1667 : class(gs_hamiltonian_type),target,intent(inout) :: gs_ham
1668 : integer,intent(in) :: ndat
1669 :
1670 : !Local variables-------------------------------
1671 0 : real(dp), contiguous, pointer :: fofr_work_ptr(:,:,:,:)
1672 : ! *************************************************************************
1673 :
1674 0 : if (.not. allocated(gs_ham%fofr_work)) then
1675 : !print *, "first allocation"
1676 : ! First allocation on CPU and GPU.
1677 0 : ABI_MALLOC(gs_ham%fofr_work, (2, gs_ham%n4, gs_ham%n5, gs_ham%n6*ndat))
1678 0 : fofr_work_ptr => gs_ham%fofr_work
1679 : #ifdef HAVE_OPENMP_OFFLOAD
1680 : !$OMP TARGET ENTER DATA MAP(alloc:fofr_work_ptr) IF (gs_ham%gpu_option==ABI_GPU_OPENMP)
1681 : #endif
1682 : end if
1683 :
1684 : ! Realloc and remap if the buffer is not large enough.
1685 : !if (gs_ham%n6*ndat > size(gs_ham%fofr_work, dim=4)) then
1686 : ! Realloc and remap if buffer size changed.
1687 0 : if (gs_ham%n6*ndat /= size(gs_ham%fofr_work, dim=4)) then
1688 : !print *, "reallocating:", gs_ham%n6*ndat, size(gs_ham%fofr_work, dim=4)
1689 0 : fofr_work_ptr => gs_ham%fofr_work
1690 : #ifdef HAVE_OPENMP_OFFLOAD
1691 : !$OMP TARGET EXIT DATA MAP(delete:fofr_work_ptr) IF (gs_ham%gpu_option==ABI_GPU_OPENMP)
1692 : #endif
1693 0 : ABI_REMALLOC(gs_ham%fofr_work, (2, gs_ham%n4, gs_ham%n5, gs_ham%n6*ndat))
1694 0 : fofr_work_ptr => gs_ham%fofr_work
1695 : #ifdef HAVE_OPENMP_OFFLOAD
1696 : !$OMP TARGET ENTER DATA MAP(alloc:fofr_work_ptr) IF (gs_ham%gpu_option==ABI_GPU_OPENMP)
1697 : #endif
1698 : end if ! realloc condition.
1699 :
1700 0 : end subroutine gsham_alloc_fofr_work
1701 : !!***
1702 :
1703 : !!****f* m_hamiltonian/gsham_print
1704 : !! NAME
1705 : !! gsham_print
1706 : !!
1707 : !! FUNCTION
1708 : !! Print info on the object.
1709 : !!
1710 : !! INPUTS
1711 : !! units=Unit numbers for output.
1712 : !! header=Header string
1713 : !! prtvol=verbosity level
1714 : !!
1715 : !! SOURCE
1716 :
1717 0 : subroutine gsham_print(gs_ham, units, header, prtvol)
1718 :
1719 : !Arguments ------------------------------------
1720 : class(gs_hamiltonian_type),intent(in) :: gs_ham
1721 : integer,intent(in) :: units(:), prtvol
1722 : character(len=*),intent(in) :: header
1723 :
1724 : !Local variables-------------------------------
1725 : character(len=5000) :: msg
1726 : ! *************************************************************************
1727 :
1728 0 : if (len_trim(header) == 0) then
1729 0 : msg = ' ==== Info on the gs_hamiltonian_type object ==== '
1730 : else
1731 0 : msg = ' ==== '//trim(adjustl(header))//' ==== '
1732 : end if
1733 0 : call wrtout(units, msg)
1734 :
1735 0 : call wrtout(units, sjoin(" kpt_k:", ktoa(gs_ham%kpt_k)))
1736 0 : call wrtout(units, sjoin(" kpt_kp:", ktoa(gs_ham%kpt_kp)))
1737 0 : call wrtout(units, sjoin(" npw_k:", itoa(gs_ham%npw_k)))
1738 0 : call wrtout(units, sjoin(" npw_kp:", itoa(gs_ham%npw_kp)))
1739 0 : call wrtout(units, sjoin(" associated(vlocal):", yesno(associated(gs_ham%vlocal))))
1740 :
1741 : if (prtvol > 0) then
1742 : end if
1743 :
1744 0 : end subroutine gsham_print
1745 : !!***
1746 :
1747 : !----------------------------------------------------------------------
1748 :
1749 : !!****f* m_hamiltonian/rfham_free
1750 : !! NAME
1751 : !! rfham_free
1752 : !!
1753 : !! FUNCTION
1754 : !! Clean and destroy rf_hamiltonian_type datastructure
1755 : !!
1756 : !! SOURCE
1757 :
1758 2155923 : subroutine rfham_free(rf_Ham)
1759 :
1760 : !Arguments ------------------------------------
1761 : class(rf_hamiltonian_type),intent(inout) :: rf_Ham
1762 : ! *************************************************************************
1763 :
1764 : DBG_ENTER("COLL")
1765 :
1766 : !@rf_hamiltonian_type
1767 :
1768 : ! Real arrays
1769 2155923 : ABI_SFREE(rf_Ham%e1kbfr_spin)
1770 2155923 : ABI_SFREE(rf_Ham%e1kbsc_spin)
1771 :
1772 : ! Real pointers
1773 2155923 : if (associated(rf_Ham%dkinpw_k)) nullify(rf_Ham%dkinpw_k)
1774 2155923 : if (associated(rf_Ham%dkinpw_kp)) nullify(rf_Ham%dkinpw_kp)
1775 2155923 : if (associated(rf_Ham%ddkinpw_k)) nullify(rf_Ham%ddkinpw_k)
1776 2155923 : if (associated(rf_Ham%ddkinpw_kp)) nullify(rf_Ham%ddkinpw_kp)
1777 2155923 : if (associated(rf_Ham%vectornd)) nullify(rf_Ham%vectornd)
1778 2155923 : if (associated(rf_Ham%vlocal1)) nullify(rf_Ham%vlocal1)
1779 2155923 : if (associated(rf_Ham%vxctaulocal)) nullify(rf_Ham%vxctaulocal)
1780 2155923 : if (associated(rf_Ham%e1kbfr)) nullify(rf_Ham%e1kbfr)
1781 2155923 : if (associated(rf_Ham%e1kbsc)) nullify(rf_Ham%e1kbsc)
1782 :
1783 : DBG_EXIT("COLL")
1784 :
1785 2155923 : end subroutine rfham_free
1786 : !!***
1787 :
1788 : !----------------------------------------------------------------------
1789 :
1790 : !!****f* m_hamiltonian/rfham_init
1791 : !! NAME
1792 : !! rfham_init
1793 : !!
1794 : !! FUNCTION
1795 : !! Creation method for the rf_hamiltonian_type structure.
1796 : !! It allocates memory and initializes all quantities that do not depend on the k-point or spin.
1797 : !!
1798 : !! INPUTS
1799 : !! [comm_atom]=optional, MPI communicator over atoms
1800 : !! cplex_paw=1 if all on-site PAW quantities are real (GS), 2 if they are complex (RF)
1801 : !! gs_Ham<gs_hamiltonian_type>=Structured datatype containing data for ground-state Hamiltonian at (k+q)
1802 : !! [has_e1kbsc]=optional, true if rf_Ham%e1kbsc has to be initialized.
1803 : !! e1kbsc contains the self-consistent 1st-order PAW Dij coefficients (depending on VHxc^(1))
1804 : !! ipert=index of perturbation
1805 : !! [mpi_atmtab(:)]=optional, indexes of the atoms treated by current proc
1806 : !! [mpi_spintab(2)]=optional, flags defining the spin(s) treated be current process:
1807 : !! mpi_spintab(1)=1 if non-polarized or spin-up treated
1808 : !! mpi_spintab(2)=1 if polarized and spin-dn treated
1809 : !! [paw_ij1(:)<paw_ij_type>]=Various 1st-order arrays given on (i,j) (partial waves)
1810 : !! channels (paw_ij1%dij and paw_ij1%difr only used here).
1811 : !!
1812 : !! SIDE EFFECTS
1813 : !! rf_Ham<rf_hamiltonian_type>=Structured datatype almost completely initialized:
1814 : !! * Basic variables and dimensions are transferred to the structure.
1815 : !! * All pointers are allocated with correct dimensions.
1816 : !! * Quantities that do not depend on the k-point or spin are initialized.
1817 : !!
1818 : !! SOURCE
1819 :
1820 2229623 : subroutine rfham_init(rf_ham, cplex, gs_Ham, ipert,&
1821 2229623 : comm_atom, mpi_atmtab, mpi_spintab, paw_ij1, has_e1kbsc) ! optional arguments
1822 :
1823 : !Arguments ------------------------------------
1824 : !scalars
1825 : class(rf_hamiltonian_type),intent(inout),target :: rf_Ham
1826 : integer,intent(in) :: cplex,ipert
1827 : integer,intent(in),optional :: comm_atom
1828 : logical,intent(in),optional :: has_e1kbsc
1829 : type(gs_hamiltonian_type),intent(in) :: gs_Ham
1830 : !arrays
1831 : integer,optional,intent(in) :: mpi_atmtab(:),mpi_spintab(2)
1832 : type(paw_ij_type),optional,intent(in) :: paw_ij1(:)
1833 :
1834 : !Local variables-------------------------------
1835 : !scalars
1836 : integer :: cplex_dij1,isp,jsp,my_comm_atom,my_nsppol
1837 : logical :: has_e1kbsc_
1838 : !arrays
1839 : integer :: my_spintab(2)
1840 2229623 : real(dp),allocatable,target :: e1kb_tmp(:,:,:,:)
1841 : ! *************************************************************************
1842 :
1843 : DBG_ENTER("COLL")
1844 :
1845 : !@rf_hamiltonian_type
1846 :
1847 : !Manage optional parameters
1848 2229623 : has_e1kbsc_=.false.;if (present(has_e1kbsc)) has_e1kbsc_=has_e1kbsc
1849 2229623 : my_comm_atom=xmpi_comm_self;if (present(comm_atom)) my_comm_atom=comm_atom
1850 4474770 : my_spintab=0;my_spintab(1:gs_Ham%nsppol)=1;if(present(mpi_spintab)) my_spintab=mpi_spintab
1851 6688869 : my_nsppol=count(my_spintab==1)
1852 :
1853 2229623 : rf_Ham%cplex =cplex
1854 2229623 : rf_Ham%n4 =gs_Ham%n4
1855 2229623 : rf_Ham%n5 =gs_Ham%n5
1856 2229623 : rf_Ham%n6 =gs_Ham%n6
1857 2229623 : rf_Ham%nvloc =gs_Ham%nvloc
1858 2229623 : rf_Ham%nsppol =gs_Ham%nsppol
1859 2229623 : rf_Ham%nspinor =gs_Ham%nspinor
1860 2229623 : rf_Ham%zora =gs_Ham%zora
1861 :
1862 2229623 : rf_Ham%dime1kb1=0
1863 2229623 : rf_Ham%dime1kb2=gs_Ham%dimekb2
1864 2229623 : if (gs_Ham%usepaw==1.and.ipert/=gs_Ham%natom+1.and.ipert/=gs_Ham%natom+10) then
1865 61684 : cplex_dij1=1;if ((gs_Ham%nspinor==2).or.any(abs(gs_Ham%nucdipmom)>tol8)) cplex_dij1=2
1866 4144 : rf_Ham%dime1kb1=cplex_dij1*(gs_Ham%lmnmax*(gs_Ham%lmnmax+1))/2
1867 : end if
1868 :
1869 : ! Allocate the arrays of the 1st-order Hamiltonian
1870 : ! We preload here 1st-order non-local factors in order to avoid
1871 : ! a communication over atoms inside the loop over spins.
1872 2229623 : if (gs_Ham%usepaw==1.and.rf_Ham%dime1kb1>0) then
1873 : if ((ipert>=1.and.ipert<=gs_Ham%natom).or.ipert==gs_Ham%natom+2.or.&
1874 4144 : ipert==gs_Ham%natom+3.or.ipert==gs_Ham%natom+4.or.ipert==gs_Ham%natom+11) then
1875 :
1876 29008 : ABI_MALLOC(rf_Ham%e1kbfr_spin,(rf_Ham%dime1kb1,rf_Ham%dime1kb2,rf_Ham%nspinor**2,cplex,my_nsppol))
1877 484083 : rf_Ham%e1kbfr_spin=zero
1878 4144 : if (has_e1kbsc_) then
1879 20671 : ABI_MALLOC(rf_Ham%e1kbsc_spin,(rf_Ham%dime1kb1,rf_Ham%dime1kb2,rf_Ham%nspinor**2,cplex,my_nsppol))
1880 467155 : rf_Ham%e1kbsc_spin=zero
1881 : end if
1882 :
1883 4144 : if (present(paw_ij1)) then
1884 :
1885 3068 : if (my_nsppol<rf_Ham%nsppol) then
1886 0 : ABI_MALLOC(e1kb_tmp,(rf_Ham%dime1kb1,rf_Ham%dime1kb2,rf_Ham%nspinor**2,cplex))
1887 : end if
1888 :
1889 : ! === Frozen term
1890 3068 : jsp=0
1891 6816 : do isp=1,rf_Ham%nsppol
1892 3748 : if (my_spintab(isp)==1) then
1893 3748 : jsp=jsp+1 ; rf_Ham%e1kbfr => rf_Ham%e1kbfr_spin(:,:,:,:,jsp)
1894 : else
1895 0 : rf_Ham%e1kbfr => e1kb_tmp
1896 : end if
1897 6816 : if (present(mpi_atmtab)) then
1898 1314 : call pawdij2e1kb(paw_ij1,isp,my_comm_atom,e1kbfr=rf_Ham%e1kbfr,mpi_atmtab=mpi_atmtab)
1899 : else
1900 2434 : call pawdij2e1kb(paw_ij1,isp,my_comm_atom,e1kbfr=rf_Ham%e1kbfr)
1901 : end if
1902 : end do
1903 :
1904 : ! === Self-consistent term
1905 3068 : if (has_e1kbsc_) then
1906 2953 : jsp=0
1907 6586 : do isp=1,rf_Ham%nsppol
1908 3633 : if (my_spintab(isp)==1) then
1909 3633 : jsp=jsp+1 ; rf_Ham%e1kbsc => rf_Ham%e1kbsc_spin(:,:,:,:,jsp)
1910 : else
1911 0 : rf_Ham%e1kbsc => e1kb_tmp
1912 : end if
1913 6586 : if (present(mpi_atmtab)) then
1914 1254 : call pawdij2e1kb(paw_ij1,isp,my_comm_atom,e1kbsc=rf_Ham%e1kbsc,mpi_atmtab=mpi_atmtab)
1915 : else
1916 2379 : call pawdij2e1kb(paw_ij1,isp,my_comm_atom,e1kbsc=rf_Ham%e1kbsc)
1917 : end if
1918 : end do
1919 : end if
1920 :
1921 3068 : if (my_nsppol<rf_Ham%nsppol) then
1922 0 : ABI_FREE(e1kb_tmp)
1923 : end if
1924 :
1925 : end if
1926 : end if
1927 : end if
1928 :
1929 2229623 : if (.not.allocated(rf_Ham%e1kbfr_spin)) then
1930 2151779 : ABI_MALLOC(rf_Ham%e1kbfr_spin,(0,0,0,0,0))
1931 : end if
1932 2229623 : if (.not.allocated(rf_Ham%e1kbsc_spin)) then
1933 2152970 : ABI_MALLOC(rf_Ham%e1kbsc_spin,(0,0,0,0,0))
1934 : end if
1935 2229623 : nullify(rf_Ham%e1kbfr)
1936 2229623 : nullify(rf_Ham%e1kbsc)
1937 :
1938 : DBG_EXIT("COLL")
1939 :
1940 2229623 : end subroutine rfham_init
1941 : !!***
1942 :
1943 : !----------------------------------------------------------------------
1944 :
1945 : !!****f* m_hamiltonian/rfham_load_spin
1946 : !! NAME
1947 : !! rfham_load_spin
1948 : !!
1949 : !! FUNCTION
1950 : !! Setup of the spin-dependent part of the 1st- and 2nd- order Hamiltonian.
1951 : !!
1952 : !! INPUTS
1953 : !! isppol=index of current spin
1954 : !! [vectornd(n4,n5,n6,nvloc)]=optional, vector potential of nuclear magnetic dipoles in real space in
1955 : !! ddk direction idir
1956 : !! [vlocal1(cplex*n4,n5,n6,nvloc)]=optional, 1st-order local potential in real space
1957 : !! [vxctaulocal(n4,n5,n6,nvloc,4)]=optional, deriv of e_XC wrt kin energy, for mGGA
1958 : !! [with_nonlocal]=optional, true if non-local factors have to be loaded
1959 : !!
1960 : !! SIDE EFFECTS
1961 : !! rf_Ham<rf_hamiltonian_type>=Structured datatype initialization phase:
1962 : !! * Quantities that depend on spin are initialized.
1963 : !!
1964 : !! SOURCE
1965 :
1966 903676 : subroutine rfham_load_spin(rf_Ham, isppol, &
1967 903676 : vectornd, vlocal1, vxctaulocal, with_nonlocal) ! optional
1968 :
1969 : !Arguments ------------------------------------
1970 : !scalars
1971 : class(rf_hamiltonian_type),intent(inout),target :: rf_Ham
1972 : integer,intent(in) :: isppol
1973 : logical,optional,intent(in) :: with_nonlocal
1974 : !arrays
1975 : real(dp),optional,target,intent(in) :: vlocal1(:,:,:,:)
1976 : real(dp),optional,target,intent(in) :: vectornd(:,:,:,:)
1977 : real(dp),optional,target,intent(in) :: vxctaulocal(:,:,:,:,:)
1978 :
1979 : !Local variables-------------------------------
1980 : integer :: jsppol
1981 : ! *************************************************************************
1982 :
1983 : DBG_ENTER("COLL")
1984 :
1985 : !@rf_hamiltonian_type
1986 :
1987 903676 : if (present(vlocal1)) then
1988 4508645 : ABI_CHECK_IEQ(size(vlocal1), rf_Ham%cplex*rf_Ham%n4*rf_Ham%n5*rf_Ham%n6*rf_Ham%nvloc, "Wrong vlocal1")
1989 901729 : rf_Ham%vlocal1 => vlocal1
1990 : end if
1991 :
1992 903676 : if (present(vectornd)) then
1993 2190 : ABI_CHECK_IEQ(size(vectornd), rf_Ham%cplex*rf_Ham%n4*rf_Ham%n5*rf_Ham%n6*rf_Ham%nvloc, "Wrong vectornd")
1994 438 : rf_Ham%vectornd => vectornd
1995 : end if
1996 :
1997 903676 : if (present(vxctaulocal)) then
1998 288 : ABI_CHECK_IEQ(size(vxctaulocal), rf_Ham%n4*rf_Ham%n5*rf_Ham%n6*rf_Ham%nvloc*4, "Wrong vxctaulocal")
1999 48 : rf_Ham%vxctaulocal => vxctaulocal
2000 : end if
2001 :
2002 : ! Retrieve non-local factors for this spin component
2003 903676 : if (present(with_nonlocal)) then
2004 902943 : if (with_nonlocal) then
2005 5417658 : if (size(rf_Ham%e1kbfr_spin)>0) then
2006 3748 : jsppol=min(isppol,size(rf_Ham%e1kbfr_spin,5))
2007 3748 : if (jsppol>0) rf_Ham%e1kbfr => rf_Ham%e1kbfr_spin(:,:,:,:,jsppol)
2008 : end if
2009 5417658 : if (size(rf_Ham%e1kbsc_spin)>0) then
2010 3633 : jsppol=min(isppol,size(rf_Ham%e1kbsc_spin,5))
2011 3633 : if (jsppol>0) rf_Ham%e1kbsc => rf_Ham%e1kbsc_spin(:,:,:,:,jsppol)
2012 : end if
2013 : end if
2014 : end if
2015 :
2016 : DBG_EXIT("COLL")
2017 :
2018 903676 : end subroutine rfham_load_spin
2019 : !!***
2020 :
2021 : !----------------------------------------------------------------------
2022 :
2023 : !!****f* m_hamiltonian/rfham_load_k
2024 : !! NAME
2025 : !! rfham_load_k
2026 : !!
2027 : !! FUNCTION
2028 : !! Setup of the k-dependent part of the 1st- and 2nd- order Hamiltonian
2029 : !!
2030 : !! INPUTS
2031 : !! [dkinpw_k]=1st derivative of the (modified) kinetic energy for each plane wave
2032 : !! [ddkinpw_k]=2nd derivative of the (modified) kinetic energy for each plane wave
2033 : !! [npw_k]=number of plane waves
2034 : !!
2035 : !! SIDE EFFECTS
2036 : !! rf_Ham<rf_hamiltonian_type>=structured datatype completed with k-dependent quantities.
2037 : !! Quantities at k^prime are set equal to quantities at k.
2038 : !!
2039 : !! SOURCE
2040 :
2041 1323496 : subroutine rfham_load_k(rf_Ham, &
2042 1323496 : dkinpw_k, ddkinpw_k, npw_k) ! optional.
2043 :
2044 : !Arguments ------------------------------------
2045 : !scalars
2046 : class(rf_hamiltonian_type),intent(inout),target :: rf_Ham
2047 : integer,intent(in),optional :: npw_k
2048 : !arrays
2049 : real(dp),intent(in),optional,target :: dkinpw_k(:),ddkinpw_k(:)
2050 : ! *************************************************************************
2051 :
2052 : DBG_ENTER("COLL")
2053 :
2054 : !@gs_hamiltonian_type
2055 :
2056 : !k-dependent scalars
2057 1323496 : if (present(npw_k)) then
2058 1305280 : rf_Ham%npw_k = npw_k
2059 1305280 : rf_Ham%npw_kp = npw_k
2060 : end if
2061 :
2062 : !Pointers to k-dependent quantities
2063 1323496 : if (present(dkinpw_k)) then
2064 1312336 : rf_Ham%dkinpw_k => dkinpw_k
2065 1312336 : rf_Ham%dkinpw_kp => dkinpw_k
2066 : end if
2067 1323496 : if (present(ddkinpw_k)) then
2068 18216 : rf_Ham%ddkinpw_k => ddkinpw_k
2069 18216 : rf_Ham%ddkinpw_kp => ddkinpw_k
2070 : end if
2071 :
2072 : DBG_EXIT("COLL")
2073 :
2074 1323496 : end subroutine rfham_load_k
2075 : !!***
2076 :
2077 : !----------------------------------------------------------------------
2078 :
2079 : !!****f* m_hamiltonian/pawdij2ekb
2080 : !! NAME
2081 : !! pawdij2ekb
2082 : !!
2083 : !! FUNCTION
2084 : !! Transfer PAW Dij (on-site GS Hamiltonian) values
2085 : !! from paw_ij datastructure to ekb array
2086 : !!
2087 : !! INPUTS
2088 : !!
2089 : !! OUTPUT
2090 : !!
2091 : !! SOURCE
2092 :
2093 24648 : subroutine pawdij2ekb(ekb,paw_ij,isppol,comm_atom,mpi_atmtab)
2094 :
2095 : !Arguments ------------------------------------
2096 : !scalars
2097 : integer,intent(in) :: isppol,comm_atom
2098 : !arrays
2099 : integer,intent(in),optional,target :: mpi_atmtab(:)
2100 : real(dp),intent(out) :: ekb(:,:,:,:)
2101 : type(paw_ij_type),intent(in) :: paw_ij(:)
2102 :
2103 : !Local variables-------------------------------
2104 : !scalars
2105 : integer :: dimdij,dimekb1,dimekb3,dimekb4,iatom,iatom_tot,ierr,ii,isp,ispden,my_natom,natom,qphase
2106 : logical :: my_atmtab_allocated,paral_atom
2107 : !arrays
2108 24648 : integer,pointer :: my_atmtab(:)
2109 : ! *************************************************************************
2110 :
2111 : DBG_ENTER("COLL")
2112 :
2113 8281691 : ekb=zero
2114 :
2115 : !Set up parallelism over atoms
2116 24648 : natom=size(ekb,2); my_natom=size(paw_ij)
2117 24648 : paral_atom=(xmpi_comm_size(comm_atom)>1)
2118 24648 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2119 24648 : call get_my_atmtab(comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2120 :
2121 : !Retrieve PAW Dij coefficients for this spin component
2122 24648 : if (my_natom>0) then
2123 22662 : if (allocated(paw_ij(1)%dij)) then
2124 22662 : dimekb1=size(ekb,1) ; dimekb3=size(ekb,3) ; dimekb4=size(ekb,4)
2125 22662 : qphase=paw_ij(1)%qphase
2126 22662 : ABI_CHECK(qphase <= dimekb4, 'paw_ij%qphase>dimekb4!')
2127 45324 : do ii=1,qphase
2128 74916 : do ispden=1,dimekb3
2129 29592 : isp=isppol; if (dimekb3==4) isp=ispden
2130 115946 : do iatom=1,my_natom
2131 63692 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
2132 63692 : dimdij=paw_ij(iatom)%cplex_dij*paw_ij(iatom)%lmn2_size
2133 63692 : ABI_CHECK(dimdij <= dimekb1, 'Size of paw_ij%dij>dimekb1!')
2134 6019649 : ekb(1:dimdij,iatom_tot,ispden,ii)=paw_ij(iatom)%dij(1+(ii-1)*dimdij:ii*dimdij,isp)
2135 : end do
2136 : end do
2137 : end do
2138 : end if
2139 : end if
2140 :
2141 : ! Communication in case of distribution over atomic sites
2142 24648 : if (paral_atom) call xmpi_sum(ekb,comm_atom,ierr)
2143 :
2144 : ! Destroy atom table used for parallelism
2145 24648 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2146 :
2147 : DBG_EXIT("COLL")
2148 :
2149 24648 : end subroutine pawdij2ekb
2150 : !!***
2151 :
2152 : !----------------------------------------------------------------------
2153 :
2154 : !!****f* m_hamiltonian/pawdij2e1kb
2155 : !! NAME
2156 : !! pawdij2e1kb
2157 : !!
2158 : !! FUNCTION
2159 : !! Transfer PAW Dij (on-site RF Hamiltonian) values
2160 : !! from paw_ij datastructure to e1kb array
2161 : !!
2162 : !! INPUTS
2163 : !!
2164 : !! OUTPUT
2165 : !!
2166 : !! SOURCE
2167 :
2168 10984 : subroutine pawdij2e1kb(paw_ij1,isppol,comm_atom,mpi_atmtab,e1kbfr,e1kbsc)
2169 :
2170 : !Arguments ------------------------------------
2171 : !scalars
2172 : integer,intent(in) :: isppol,comm_atom
2173 : !arrays
2174 : integer,intent(in),optional,target :: mpi_atmtab(:)
2175 : real(dp),optional,intent(out) :: e1kbfr(:,:,:,:),e1kbsc(:,:,:,:)
2176 : type(paw_ij_type),intent(in) :: paw_ij1(:)
2177 :
2178 : !Local variables-------------------------------
2179 : !scalars
2180 : integer :: dimdij1,dime1kb1,dime1kb3,dime1kb4,iatom,iatom_tot,ierr,isp,ispden
2181 : integer :: my_natom,natom,qphase
2182 : logical :: my_atmtab_allocated,paral_atom
2183 : !arrays
2184 10984 : integer,pointer :: my_atmtab(:)
2185 : ! *************************************************************************
2186 :
2187 : DBG_ENTER("COLL")
2188 :
2189 0 : if ((.not.present(e1kbfr)).and.(.not.present(e1kbsc))) return
2190 10984 : if (present(e1kbfr)) then
2191 1041260 : e1kbfr=zero ; natom=size(e1kbfr,2)
2192 : end if
2193 10984 : if (present(e1kbsc)) then
2194 464202 : e1kbsc=zero ; natom=size(e1kbsc,2)
2195 : end if
2196 :
2197 : ! Set up parallelism over atoms
2198 10984 : my_natom=size(paw_ij1) ; paral_atom=(xmpi_comm_size(comm_atom)>1)
2199 10984 : nullify(my_atmtab);if (present(mpi_atmtab)) my_atmtab => mpi_atmtab
2200 10984 : call get_my_atmtab(comm_atom,my_atmtab,my_atmtab_allocated,paral_atom,natom,my_natom_ref=my_natom)
2201 :
2202 : ! Retrieve 1st-order PAW Dij coefficients for this spin component (frozen)
2203 10984 : if (my_natom>0.and.present(e1kbfr)) then
2204 7289 : if (allocated(paw_ij1(1)%dijfr)) then
2205 7289 : dime1kb1=size(e1kbfr,1) ; dime1kb3=size(e1kbfr,3) ; dime1kb4=size(e1kbfr,4)
2206 7289 : ABI_CHECK_IEQ(paw_ij1(1)%qphase, dime1kb4,'BUG in pawdij2e1kb (1)!')
2207 14578 : do ispden=1,dime1kb3
2208 7289 : isp=isppol;if (dime1kb3==4) isp=ispden
2209 29700 : do iatom=1,my_natom
2210 15122 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
2211 15122 : qphase=paw_ij1(iatom)%qphase
2212 15122 : dimdij1=paw_ij1(iatom)%cplex_dij*paw_ij1(iatom)%lmn2_size
2213 15122 : ABI_CHECK(dimdij1<=dime1kb1,'BUG: size of paw_ij1%dij>dime1kb1!')
2214 587176 : e1kbfr(1:dimdij1,iatom_tot,ispden,1)=paw_ij1(iatom)%dijfr(1:dimdij1,isp)
2215 60859 : if (qphase==2) e1kbfr(1:dimdij1,iatom_tot,ispden,2)=paw_ij1(iatom)%dijfr(dimdij1+1:2*dimdij1,isp)
2216 : end do
2217 : end do
2218 : end if
2219 : end if
2220 :
2221 : ! Retrieve 1st-order PAW Dij coefficients for this spin component (self-consistent)
2222 10984 : if (my_natom>0.and.present(e1kbsc)) then
2223 3613 : if (allocated(paw_ij1(1)%dijfr).and.allocated(paw_ij1(1)%dij)) then
2224 3613 : dime1kb1=size(e1kbsc,1) ; dime1kb3=size(e1kbsc,3) ; dime1kb4=size(e1kbsc,4)
2225 3613 : ABI_CHECK_IEQ(paw_ij1(1)%qphase, dime1kb4, 'BUG in pawdij2e1kb (1)!')
2226 7226 : do ispden=1,dime1kb3
2227 3613 : isp=isppol;if (dime1kb3==4) isp=ispden
2228 14577 : do iatom=1,my_natom
2229 7351 : iatom_tot=iatom;if (paral_atom) iatom_tot=my_atmtab(iatom)
2230 7351 : qphase=paw_ij1(iatom)%qphase
2231 7351 : dimdij1=paw_ij1(iatom)%cplex_dij*paw_ij1(iatom)%lmn2_size
2232 7351 : ABI_CHECK(dimdij1<=dime1kb1, 'BUG: size of paw_ij1%dij>dime1kb1!')
2233 : e1kbsc(1:dimdij1,iatom_tot,ispden,1)=paw_ij1(iatom)%dij (1:dimdij1,isp) &
2234 281835 : & -paw_ij1(iatom)%dijfr(1:dimdij1,isp)
2235 7351 : if (qphase==2) e1kbsc(1:dimdij1,iatom_tot,ispden,2)=paw_ij1(iatom)%dij (dimdij1+1:2*dimdij1,isp) &
2236 28165 : & -paw_ij1(iatom)%dijfr(dimdij1+1:2*dimdij1,isp)
2237 : end do
2238 : end do
2239 : end if
2240 : end if
2241 :
2242 : ! Communication in case of distribution over atomic sites
2243 10984 : if (paral_atom) then
2244 4494 : if (present(e1kbfr)) call xmpi_sum(e1kbfr,comm_atom,ierr)
2245 4494 : if (present(e1kbsc)) call xmpi_sum(e1kbsc,comm_atom,ierr)
2246 : end if
2247 :
2248 : ! Destroy atom table used for parallelism
2249 10984 : call free_my_atmtab(my_atmtab,my_atmtab_allocated)
2250 :
2251 : DBG_EXIT("COLL")
2252 :
2253 10984 : end subroutine pawdij2e1kb
2254 : !!***
2255 :
2256 : !!****f* ABINIT/gspot_transgrid_and_pack
2257 : !! NAME
2258 : !! gspot_transgrid_and_pack
2259 : !!
2260 : !! FUNCTION
2261 : !! Set up local potential vlocal on the coarse FFT mesh with proper dimensioning from vtrial given on the fine mesh.
2262 : !! Also take into account the spin.
2263 : !!
2264 : !! INPUTS
2265 : !! isppol=Spin polarization.
2266 : !! usepaw=1 if PAW
2267 : !! paral_kgb: 1 if paral_kgb
2268 : !! nfft=(effective) number of FFT grid points on the coarse mesh (for this processor)
2269 : !! ngfft(18)contain all needed information about 3D FFT, for the coarse FFT mesh. see ~abinit/doc/variables/vargs.htm#ngfft
2270 : !! nfftf=(effective) number of FFT grid points on the fine mesh (for this processor)
2271 : !! nvloc==1 if nspden <=2, nvloc==4 for nspden==4,
2272 : !! nspden=Number of spin density components.
2273 : !! ncomp=Number of extra components in vtrial and vlocal (e.g. 1 if LDA/GGA pot, 4 for Meta-GGA, etc).
2274 : !! pawfgr <type(pawfgr_type)>=fine grid parameters and related data
2275 : !! vtrial(nfftf,nspden)=INPUT potential Vtrial(r).
2276 : !! mpi_enreg=information about MPI parallelization
2277 : !!
2278 : !! OUTPUT
2279 : !! vlocal(n4,n5,n6,nvloc,ncomp): Potential on the coarse grid.
2280 : !!
2281 : !! SOURCE
2282 :
2283 54924 : subroutine gspot_transgrid_and_pack(isppol, usepaw, paral_kgb, nfft, ngfft, nfftf, &
2284 54924 : nspden, nvloc, ncomp, pawfgr, mpi_enreg, vtrial, vlocal)
2285 :
2286 : !Arguments -------------------------------
2287 : integer,intent(in) :: isppol, nspden, ncomp, usepaw, paral_kgb, nfft, nfftf, nvloc
2288 : type(pawfgr_type), intent(in) :: pawfgr
2289 : type(MPI_type), intent(in) :: mpi_enreg
2290 : !arrays
2291 : integer,intent(in) :: ngfft(18)
2292 : real(dp),intent(inout) :: vtrial(nfftf, nspden, ncomp)
2293 : real(dp),intent(out) :: vlocal(ngfft(4), ngfft(5), ngfft(6), nvloc, ncomp)
2294 :
2295 : !Local variables-------------------------------
2296 : !scalars
2297 : integer :: n1,n2,n3,n4,n5,n6,ispden,ic
2298 : real(dp) :: rhodum(1)
2299 54924 : real(dp),allocatable :: cgrvtrial(:,:), vlocal_tmp(:,:,:)
2300 : ! *************************************************************************
2301 :
2302 : ! Coarse mesh.
2303 54924 : n1=ngfft(1); n2=ngfft(2); n3=ngfft(3)
2304 54924 : n4=ngfft(4); n5=ngfft(5); n6=ngfft(6)
2305 :
2306 : ! Set up local potential vlocal with proper dimensioning, from vtrial
2307 : ! Also take into account the spin.
2308 54924 : if (nspden /= 4) then
2309 52458 : if (usepaw == 0 .or. pawfgr%usefinegrid == 0) then
2310 : ! Fine mesh == Coarse mesh.
2311 70169 : do ic=1,ncomp
2312 70169 : call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial(:,:,ic),vlocal(:,:,:,:,ic),2)
2313 : end do
2314 : else
2315 : ! Transfer from fine mesh to coarse and then pack data
2316 70876 : ABI_MALLOC(cgrvtrial,(nfft,nspden))
2317 38737 : do ic=1,ncomp
2318 21018 : call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr, rhodum,rhodum,cgrvtrial,vtrial(:,:,ic))
2319 38737 : call fftpac(isppol,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft, cgrvtrial,vlocal(:,:,:,:,ic),2)
2320 : end do
2321 17719 : ABI_FREE(cgrvtrial)
2322 : end if
2323 : else
2324 : ! nspden == 4. replace isppol by loop over ispden.
2325 12330 : ABI_MALLOC(vlocal_tmp, (n4,n5,n6))
2326 2466 : if (usepaw == 0 .or. pawfgr%usefinegrid == 0) then
2327 : ! Fine mesh == Coarse mesh.
2328 2768 : do ic=1,ncomp
2329 8304 : do ispden=1,nspden
2330 5536 : call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,vtrial(:,:,ic),vlocal_tmp,2)
2331 71371728 : vlocal(:,:,:,ispden,ic) = vlocal_tmp(:,:,:)
2332 : end do
2333 : end do
2334 : else
2335 : ! Transfer from fine mesh to coarse and then pack data
2336 4328 : ABI_MALLOC(cgrvtrial,(nfft,nspden))
2337 2164 : do ic=1,ncomp
2338 1082 : call transgrid(1,mpi_enreg,nspden,-1,0,0,paral_kgb,pawfgr,rhodum,rhodum,cgrvtrial,vtrial(:,:,ic))
2339 6492 : do ispden=1,nspden
2340 4328 : call fftpac(ispden,mpi_enreg,nspden,n1,n2,n3,n4,n5,n6,ngfft,cgrvtrial,vlocal_tmp,2)
2341 21768666 : vlocal(:,:,:,ispden,ic) = vlocal_tmp(:,:,:)
2342 : end do
2343 : end do
2344 1082 : ABI_FREE(cgrvtrial)
2345 : end if
2346 2466 : ABI_FREE(vlocal_tmp)
2347 : end if ! nspden
2348 :
2349 54924 : end subroutine gspot_transgrid_and_pack
2350 : !!***
2351 :
2352 4074 : end module m_hamiltonian
2353 : !!***
|
