Line data Source code
1 : !!****m* ABINIT/m_xg_nonlop
2 : !! NAME
3 : !! m_xg_nonlop
4 : !!
5 : !! FUNCTION
6 : !! This module provides functions to compute the nonlocal operator by means of the BLAS GEMM
7 : !! routine. By treating ndat simultaneous wavefunctions, it is able to exploit BLAS3 routines,
8 : !! which leads to excellent CPU efficiency and OpenMP scalability.
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2022-2026 ABINIT group (LB)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_xg_nonlop
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_xmpi
30 : ! use m_abi_linalg
31 : use m_xg
32 : use m_xomp
33 : #ifdef HAVE_OPENMP
34 : use omp_lib
35 : #endif
36 :
37 : use m_time, only : timab
38 : use defs_abitypes, only : MPI_type
39 :
40 : use m_paw_ij, only : paw_ij_type
41 : use m_pawtab, only : pawtab_type
42 :
43 : implicit none
44 :
45 : private
46 :
47 : double precision, parameter :: inv_sqrt2 = 1/sqrt2
48 :
49 : ! Independent timers of xg_nonlop :
50 : integer, parameter :: tim_getcprj = 2101
51 : integer, parameter :: tim_apply_prj = 2102
52 : integer, parameter :: tim_apply_Aij = 2103
53 : integer, parameter :: tim_mult_cprj = 2104
54 : integer, parameter :: tim_make_k = 2105
55 : integer, parameter :: tim_make_Dij = 2106
56 : integer, parameter :: tim_make_Sij = 2107
57 : integer, parameter :: tim_make_ekb = 2108
58 : integer, parameter :: tim_apply_diag = 2109
59 : integer, parameter :: tim_init = 2110
60 :
61 : ! Timers that depend on other xg_nonlop timers :
62 : integer, parameter :: tim_getXSY = 2120
63 : integer, parameter :: tim_getXHY = 2121
64 : integer, parameter :: tim_getHmeSX = 2122
65 : integer, parameter :: tim_iter_refinement = 2123
66 :
67 : integer, parameter :: tim_getcprj_gemm = 2130
68 : integer, parameter :: tim_getcprj_copy = 2131
69 : integer, parameter :: tim_getcprj_mpi = 2132
70 : integer, parameter :: tim_getcprj_otf = 2133
71 :
72 : integer, parameter :: tim_apply_prj_gemm = 2135
73 : integer, parameter :: tim_apply_prj_copy = 2136
74 : integer, parameter :: tim_apply_prj_mpi = 2137
75 : integer, parameter :: tim_apply_prj_otf = 2138
76 :
77 : integer, parameter :: tim_mult_cprj_gemm = 2140
78 : integer, parameter :: tim_mult_cprj_copy = 2141
79 : integer, parameter :: tim_mult_cprj_mpi = 2142
80 :
81 : integer, parameter :: tim_forces_stress = 2150
82 : integer, parameter :: tim_fst_start = 2151
83 : integer, parameter :: tim_fst_cprj_deriv_f = 2152
84 : integer, parameter :: tim_fst_cprj_deriv_str = 2153
85 : integer, parameter :: tim_fst_mult_cprj_f = 2154
86 : integer, parameter :: tim_fst_mult_cprj_str = 2155
87 : integer, parameter :: tim_fst_work_str = 2156
88 :
89 : integer, parameter, public :: DERIV_ATOM = 1
90 : integer, parameter, public :: DERIV_STRESS = 2
91 :
92 : type,public :: xg_nonlop_t
93 :
94 : integer :: cplex
95 : integer :: cplex_alldij
96 : integer :: cprjdim
97 : integer :: comm_atom
98 : integer :: comm_band
99 : integer :: npw_k
100 : integer :: total_npw_k
101 : integer :: max_npw_k
102 : integer :: me_band
103 : integer :: my_natom
104 : integer :: natom
105 : integer :: mkmem
106 : integer :: nlmn_max
107 : integer :: ntypat
108 : integer :: nspinor
109 : integer :: space_pw
110 : integer :: space_cprj
111 : integer :: space_Dij
112 : logical :: paw
113 : integer :: option
114 : real(dp) :: weight
115 :
116 : integer, pointer :: mpi_atmtab(:)
117 : integer, pointer :: indlmn(:,:,:)
118 : integer, pointer :: nattyp(:)
119 :
120 : integer,allocatable :: l_npw_k(:)
121 : integer,allocatable :: l_shift_npw_k(:)
122 :
123 : real(dp), pointer :: sij_triangular_mat(:,:)
124 :
125 : integer, allocatable :: nlmn_natom(:)
126 : integer, allocatable :: nlmn_ntypat(:)
127 :
128 : real(dp), pointer :: ffnl_k(:,:,:,:)
129 : real(dp), pointer :: ph3d_k(:,:,:)
130 :
131 : real(dp), pointer :: kpg_k(:,:)
132 :
133 : type(xg_t),pointer :: projectors(:)
134 : type(xg_t),pointer :: projectors_k
135 :
136 : type(xg_t),pointer :: ffnl_gather(:)
137 : type(xg_t),pointer :: ffnl_gather_k
138 :
139 : type(xg_t),pointer :: ph3d_gather(:)
140 : type(xg_t),pointer :: ph3d_gather_k
141 :
142 : ! non paw only:
143 : type(xg_t) :: ekb
144 : ! end non paw only
145 :
146 : ! paw only:
147 : type(xg_t),pointer :: gram_proj(:)
148 : type(xg_t),pointer :: gram_proj_k
149 :
150 : type(xg_t) :: Dij
151 : type(xgBlock_t) :: Dij_spin
152 : type(xg_t) :: Sij
153 :
154 : type(xg_t) :: Sijm1
155 : type(xg_t), pointer :: invSij_approx(:)
156 : type(xg_t), pointer :: invSij_approx_k
157 : ! end paw only
158 :
159 : end type xg_nonlop_t
160 : !!***
161 :
162 : ! To initialize/make/destroy xg_nonlop object
163 : public :: xg_nonlop_init
164 : public :: xg_nonlop_make_k
165 : public :: xg_nonlop_destroy
166 : public :: xg_nonlop_make_ekb ! non paw only
167 : public :: xg_nonlop_destroy_ekb ! non paw only
168 : public :: xg_nonlop_update_weight
169 : public :: xg_nonlop_init_cplex_alldij ! paw only
170 : public :: xg_nonlop_make_Dij ! paw only
171 : public :: xg_nonlop_set_Dij_spin! paw only
172 : public :: xg_nonlop_make_Sij ! paw only
173 : public :: xg_nonlop_destroy_Dij ! paw only
174 : public :: xg_nonlop_destroy_Sij ! paw only
175 : ! Generic operations
176 : public :: xg_nonlop_getcprj
177 : public :: xg_nonlop_apply_Aij
178 : public :: xg_nonlop_precond_iterative_refinement
179 : public :: xg_nonlop_mult_cprj
180 : public :: xg_nonlop_apply_prj
181 : public :: xg_nonlop_colwiseXAX
182 : public :: xg_nonlop_colwiseXDX
183 : public :: xg_nonlop_getXAY
184 : public :: xg_nonlop_getXDY
185 : public :: xg_nonlop_getAX
186 : public :: xg_nonlop_getDX
187 : ! Specific operations (using Sij/Dij or ekb arrays)
188 : public :: xg_nonlop_getXHY
189 : public :: xg_nonlop_getHX
190 : public :: xg_nonlop_colwiseXHX
191 : public :: xg_nonlop_getXSY ! paw only
192 : public :: xg_nonlop_getSX ! paw only
193 : public :: xg_nonlop_getSm1X ! paw only
194 : public :: xg_nonlop_getHmeSX ! paw only
195 : public :: xg_nonlop_forces_stress
196 :
197 : contains
198 : !!***
199 :
200 : !!****f* m_xg_nonlop/xg_nonlop_init
201 : !! NAME
202 : !! xg_nonlop_init
203 : !!
204 : !! FUNCTION
205 : !! Initalization of the xg_nonlop_kpt array
206 : !!
207 : !! INPUTS
208 : !!
209 : !! SOURCE
210 :
211 162 : subroutine xg_nonlop_init(xg_nonlop,indlmn,my_natom,nattyp,mkmem,ntypat,nspinor,ucvol,usepaw,&
212 162 : xg_nonlop_option,me_band,comm_band,comm_atom,mpi_atmtab)
213 :
214 : integer ,intent(in) :: me_band,comm_band,comm_atom
215 : integer ,intent(in) :: my_natom
216 : integer ,intent(in) :: mkmem
217 : integer ,intent(in) :: ntypat
218 : integer ,intent(in) :: nspinor
219 : integer ,intent(in) :: usepaw
220 : integer ,intent(in) :: xg_nonlop_option
221 : real(dp),intent(in) :: ucvol
222 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
223 :
224 : integer,optional,intent(in),target :: mpi_atmtab(:)
225 : integer,intent(in),target :: indlmn(:,:,:)
226 : integer,intent(in),target :: nattyp(:)
227 : real(dp) :: tsec(2)
228 :
229 : integer :: itypat,cprjdim,nlmn,nlmn_max,natom,shift,nmpi,nattyp_i
230 :
231 162 : call timab(tim_init,1,tsec)
232 :
233 162 : xg_nonlop%mkmem=mkmem
234 162 : xg_nonlop%nspinor=nspinor
235 162 : xg_nonlop%comm_atom=comm_atom
236 162 : if (xmpi_comm_size(comm_atom)>1) then
237 108 : if (present(mpi_atmtab)) then
238 108 : xg_nonlop%mpi_atmtab=>mpi_atmtab
239 : else
240 0 : ABI_ERROR("mpi_atmtab must be present")
241 : end if
242 : else
243 54 : xg_nonlop%mpi_atmtab=>null()
244 : end if
245 162 : xg_nonlop%my_natom=my_natom
246 162 : xg_nonlop%me_band=me_band
247 162 : xg_nonlop%comm_band=comm_band
248 :
249 162 : if (xg_nonlop_option==0.or.xg_nonlop_option==1) then
250 162 : xg_nonlop%option = xg_nonlop_option
251 : else
252 0 : ABI_ERROR('Wrong value of xg_nonlop_option')
253 : end if
254 :
255 162 : xg_nonlop%paw=usepaw==1
256 :
257 162 : xg_nonlop%space_pw=0
258 162 : xg_nonlop%space_cprj=0
259 162 : xg_nonlop%space_Dij=0
260 :
261 162 : xg_nonlop%weight=four_pi/sqrt(ucvol)
262 162 : xg_nonlop%ntypat=ntypat
263 :
264 162 : xg_nonlop%nattyp=>nattyp
265 162 : xg_nonlop%indlmn=>indlmn
266 :
267 162 : natom=0
268 476 : do itypat=1,ntypat
269 476 : natom = natom + nattyp(itypat)
270 : end do
271 :
272 486 : ABI_MALLOC(xg_nonlop%nlmn_ntypat,(ntypat))
273 486 : ABI_MALLOC(xg_nonlop%nlmn_natom,(natom))
274 :
275 162 : cprjdim=0
276 162 : shift=0
277 162 : nlmn_max=0
278 476 : do itypat=1,ntypat
279 4842 : nlmn = count(indlmn(3,:,itypat)>0)
280 314 : if (nlmn>nlmn_max) nlmn_max=nlmn
281 314 : xg_nonlop%nlmn_ntypat(itypat) = nlmn
282 314 : nattyp_i = nattyp(itypat)
283 476 : if (nattyp_i>0) then
284 310 : cprjdim = cprjdim + nattyp_i*nlmn
285 674 : xg_nonlop%nlmn_natom(1+shift:nattyp_i+shift) = nlmn
286 : shift = shift + nattyp_i
287 : end if
288 : end do
289 162 : xg_nonlop%nlmn_max = nlmn_max
290 162 : xg_nonlop%natom = natom
291 162 : xg_nonlop%cprjdim = cprjdim
292 :
293 1585 : ABI_MALLOC(xg_nonlop%projectors,(mkmem))
294 1423 : ABI_MALLOC(xg_nonlop%ffnl_gather,(mkmem))
295 1423 : ABI_MALLOC(xg_nonlop%ph3d_gather,(mkmem))
296 162 : if (xg_nonlop%paw) then
297 1183 : ABI_MALLOC(xg_nonlop%gram_proj,(mkmem))
298 1183 : ABI_MALLOC(xg_nonlop%invSij_approx,(mkmem))
299 : end if
300 :
301 162 : nmpi = xmpi_comm_size(xg_nonlop%comm_band)
302 486 : ABI_MALLOC(xg_nonlop%l_npw_k,(nmpi))
303 324 : ABI_MALLOC(xg_nonlop%l_shift_npw_k,(nmpi))
304 :
305 162 : call timab(tim_init,2,tsec)
306 :
307 162 : end subroutine xg_nonlop_init
308 : !!***
309 :
310 : !!****f* m_xg_nonlop/xg_nonlop_update_weight
311 : !! NAME
312 : !! xg_nonlop_init
313 : !!
314 : !! FUNCTION
315 : !! Compute new weights with respect to ucvol
316 : !!
317 : !! INPUTS
318 : !!
319 : !! SOURCE
320 2572 : subroutine xg_nonlop_update_weight(xg_nonlop,ucvol)
321 :
322 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
323 : real(dp),intent(in) :: ucvol
324 :
325 2572 : xg_nonlop%weight=four_pi/sqrt(ucvol)
326 :
327 2572 : end subroutine xg_nonlop_update_weight
328 : !!***
329 :
330 : !!****f* m_xg_nonlop/xg_nonlop_init_cplex_alldij
331 : !! NAME
332 : !! xg_nonlop_init_cplex_alldij
333 : !!
334 : !! FUNCTION
335 : !!
336 : !! INPUTS
337 : !!
338 : !! SOURCE
339 142 : subroutine xg_nonlop_init_cplex_alldij(xg_nonlop,paw_ij)
340 :
341 : type(paw_ij_type),intent(in) :: paw_ij(:)
342 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
343 :
344 : logical :: paral_atom
345 : integer :: iatom,iatom_tot,ierr
346 : integer,allocatable :: l_cplex(:)
347 : real(dp) :: tsec(2)
348 :
349 142 : call timab(tim_init,1,tsec)
350 :
351 142 : paral_atom=(xmpi_comm_size(xg_nonlop%comm_atom)>1)
352 :
353 426 : ABI_MALLOC(l_cplex,(xg_nonlop%natom))
354 466 : l_cplex=0
355 280 : do iatom=1,xg_nonlop%my_natom ! loop over atoms treated by this proc
356 138 : iatom_tot=iatom;if (paral_atom) iatom_tot=xg_nonlop%mpi_atmtab(iatom)
357 280 : l_cplex(iatom_tot)=paw_ij(iatom)%cplex_dij
358 : end do
359 142 : call xmpi_sum(l_cplex,xg_nonlop%comm_atom,ierr)
360 142 : xg_nonlop%cplex_alldij = 0
361 466 : do iatom=1,xg_nonlop%natom ! loop over all atoms
362 466 : if (xg_nonlop%cplex_alldij<l_cplex(iatom)) xg_nonlop%cplex_alldij=l_cplex(iatom)
363 : end do
364 142 : ABI_FREE(l_cplex)
365 :
366 142 : call timab(tim_init,2,tsec)
367 :
368 142 : end subroutine xg_nonlop_init_cplex_alldij
369 : !!***
370 :
371 162 : subroutine xg_nonlop_destroy(xg_nonlop)
372 :
373 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
374 :
375 : integer :: ikpt
376 :
377 162 : if (xg_nonlop%paw) then
378 138 : call xg_nonlop_destroy_Sij(xg_nonlop) ! Can be destroyed before
379 138 : call xg_nonlop_destroy_Dij(xg_nonlop) ! Can be destroyed before
380 : else
381 24 : call xg_nonlop_destroy_ekb(xg_nonlop) ! Can be destroyed before
382 : end if
383 :
384 162 : ABI_FREE(xg_nonlop%nlmn_ntypat)
385 162 : ABI_FREE(xg_nonlop%nlmn_natom)
386 :
387 162 : ABI_FREE(xg_nonlop%l_npw_k)
388 162 : ABI_FREE(xg_nonlop%l_shift_npw_k)
389 :
390 1261 : do ikpt=1,xg_nonlop%mkmem
391 1099 : call xg_free(xg_nonlop%projectors(ikpt))
392 1099 : call xg_free(xg_nonlop%ffnl_gather(ikpt))
393 1099 : call xg_free(xg_nonlop%ph3d_gather(ikpt))
394 1261 : if (xg_nonlop%paw) then
395 907 : call xg_free(xg_nonlop%gram_proj(ikpt))
396 907 : call xg_free(xg_nonlop%invSij_approx(ikpt))
397 : end if
398 : end do
399 162 : ABI_FREE(xg_nonlop%projectors)
400 162 : ABI_FREE(xg_nonlop%ffnl_gather)
401 162 : ABI_FREE(xg_nonlop%ph3d_gather)
402 162 : if (xg_nonlop%paw) then
403 138 : ABI_FREE(xg_nonlop%gram_proj)
404 138 : ABI_FREE(xg_nonlop%invSij_approx)
405 : end if
406 :
407 162 : end subroutine xg_nonlop_destroy
408 : !!***
409 :
410 24 : subroutine xg_nonlop_make_ekb(xg_nonlop,ekb)
411 :
412 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
413 : real(dp), intent(in) :: ekb(:,:)
414 :
415 : integer :: itypat, ilmn, iln, nlmn, nlmn_max, ntypat
416 24 : real(dp),pointer :: ekb_(:)
417 : real(dp) :: tsec(2)
418 : type(xgBlock_t) :: ekb_itypat
419 :
420 : ! *************************************************************************
421 :
422 24 : call timab(tim_make_ekb,1,tsec)
423 :
424 24 : if (xg_nonlop%paw) then
425 0 : ABI_ERROR('Not implemented with paw=True.')
426 : end if
427 :
428 24 : ntypat = xg_nonlop%ntypat
429 24 : nlmn_max = xg_nonlop%nlmn_max
430 :
431 24 : call xg_init(xg_nonlop%ekb,SPACE_R,nlmn_max,ntypat,xmpi_comm_self)
432 :
433 72 : do itypat=1,ntypat
434 :
435 48 : nlmn=xg_nonlop%nlmn_ntypat(itypat)
436 48 : call xg_setBlock(xg_nonlop%ekb,ekb_itypat,nlmn,1,fcol=itypat)
437 48 : call xgBlock_reverseMap_1D(ekb_itypat,ekb_)
438 1584 : do ilmn=1,nlmn
439 1512 : iln=xg_nonlop%indlmn(5,ilmn,itypat)
440 1560 : ekb_(ilmn) = ekb(iln,itypat)
441 : end do
442 :
443 : end do
444 :
445 24 : call timab(tim_make_ekb,2,tsec)
446 :
447 24 : end subroutine xg_nonlop_make_ekb
448 : !!***
449 :
450 2188 : subroutine xg_nonlop_make_Dij(xg_nonlop,paw_ij,nsppol,atindx)
451 :
452 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
453 : type(paw_ij_type),intent(in) :: paw_ij(:)
454 : integer,intent(in) :: nsppol
455 : integer,intent(in) :: atindx(:)
456 :
457 : logical :: paral_atom
458 : integer :: isppol, iatom, iatom_input, iatom_type, nlmn, nlmn_max, natom, nspinor, shift
459 : integer :: ilmn, jlmn, j0lmn, jjlmn, ijlmn
460 : integer :: cplex_alldij,cplex_dij,isp,isps,jsp,jsps,ijsp
461 2188 : real(dp),pointer :: Dij_iatom_(:,:)
462 : real(dp) :: tsec(2)
463 : type(xgBlock_t) :: Dij_iatom
464 :
465 : ! *************************************************************************
466 :
467 2188 : call timab(tim_make_Dij,1,tsec)
468 :
469 2188 : if (.not.xg_nonlop%paw) then
470 0 : ABI_ERROR('Not implemented with paw=False.')
471 : end if
472 :
473 2188 : nspinor = xg_nonlop%nspinor
474 2188 : natom = xg_nonlop%natom
475 2188 : nlmn_max = xg_nonlop%nlmn_max
476 :
477 2188 : paral_atom=(xmpi_comm_size(xg_nonlop%comm_atom)>1)
478 :
479 2188 : cplex_alldij = xg_nonlop%cplex_alldij
480 2188 : if (cplex_alldij==1) then
481 1645 : xg_nonlop%space_Dij=SPACE_R
482 543 : else if (cplex_alldij==2) then
483 543 : xg_nonlop%space_Dij=SPACE_C
484 : else
485 0 : ABI_ERROR('Bad cplex_alldij')
486 : end if
487 :
488 2188 : call xg_init(xg_nonlop%Dij,xg_nonlop%space_Dij,nspinor*nlmn_max,nspinor*nlmn_max*natom*nsppol,xmpi_comm_null)
489 :
490 4440 : do isppol=1,nsppol
491 :
492 6518 : do iatom=1,xg_nonlop%my_natom ! loop over atoms treated by this proc
493 2078 : iatom_input = iatom
494 2078 : if (paral_atom) iatom_input=xg_nonlop%mpi_atmtab(iatom) ! mpi_atmtab(iatom) has the ordering of the input file
495 2078 : iatom_type = atindx(iatom_input) ! convert iatom from input file to iatom ordered by type
496 2078 : nlmn=xg_nonlop%nlmn_natom(iatom_type)
497 2078 : shift=1+(iatom_type-1)*nspinor*nlmn_max+(isppol-1)*natom*nspinor*nlmn_max
498 2078 : call xg_setBlock(xg_nonlop%Dij,Dij_iatom,nspinor*nlmn,nspinor*nlmn,fcol=shift)
499 2078 : cplex_dij=paw_ij(iatom)%cplex_dij
500 2078 : call xgBlock_reverseMap(Dij_iatom,Dij_iatom_)
501 6822 : do jsp=1,nspinor
502 2492 : jsps = (jsp-1)*nlmn
503 27246 : do jlmn=1,nlmn
504 22676 : j0lmn=jlmn*(jlmn-1)/2
505 22676 : jjlmn=j0lmn+jlmn
506 55708 : do isp=1,nspinor
507 30540 : isps = (isp-1)*nlmn
508 30540 : if (nspinor==1) then
509 : ijsp = isppol
510 : else
511 15728 : if (isp==jsp) then
512 : ijsp = isp
513 7864 : else if (isp==1) then
514 : ijsp = 3 ! up/down
515 : else
516 3932 : ijsp = 4 ! down/up
517 : end if
518 : end if
519 : ! see m_hamiltonian:pawdij2ekb
520 30540 : if (cplex_dij==1) then
521 14812 : Dij_iatom_(cplex_alldij*(jlmn+isps-1)+1,jlmn+jsps) = paw_ij(iatom)%dij(jjlmn,ijsp)
522 14812 : if (cplex_alldij==2) Dij_iatom_(2*(jlmn+isps),jlmn+jsps) = zero
523 : else
524 15728 : Dij_iatom_(2*(jlmn+isps-1)+1,jlmn+jsps) = paw_ij(iatom)%dij(2*jjlmn-1,ijsp)
525 15728 : Dij_iatom_(2*(jlmn+isps) ,jlmn+jsps) = paw_ij(iatom)%dij(2*jjlmn ,ijsp)
526 : end if
527 195926 : do ilmn=1,jlmn-1
528 : ! see m_hamiltonian:pawdij2ekb and opernlc_ylm
529 142710 : ijlmn=j0lmn+ilmn
530 173250 : if (cplex_dij==1) then
531 65342 : Dij_iatom_(cplex_alldij*(ilmn+isps-1)+1,jlmn+jsps) = paw_ij(iatom)%dij(ijlmn,ijsp)
532 65342 : if (cplex_alldij==2) Dij_iatom_(2*(ilmn+isps),jlmn+jsps) = zero
533 65342 : Dij_iatom_(cplex_alldij*(jlmn+jsps-1)+1,ilmn+isps) = paw_ij(iatom)%dij(ijlmn,ijsp)
534 65342 : if (cplex_alldij==2) Dij_iatom_(2*(jlmn+jsps),ilmn+isps) = zero
535 : else
536 77368 : Dij_iatom_(2*(ilmn+isps-1)+1,jlmn+jsps) = paw_ij(iatom)%dij(2*ijlmn-1,ijsp)
537 77368 : Dij_iatom_(2*(ilmn+isps) ,jlmn+jsps) = paw_ij(iatom)%dij(2*ijlmn ,ijsp)
538 77368 : Dij_iatom_(2*(jlmn+jsps-1)+1,ilmn+isps) = paw_ij(iatom)%dij(2*ijlmn-1,ijsp)
539 77368 : Dij_iatom_(2*(jlmn+jsps) ,ilmn+isps) =-paw_ij(iatom)%dij(2*ijlmn ,ijsp)
540 : end if
541 : end do
542 : end do
543 : end do
544 : end do
545 : end do
546 :
547 : end do
548 :
549 : ! Communication in case of distribution over atomic sites
550 2188 : if (paral_atom) then
551 1692 : call xgBlock_mpi_sum(xg_nonlop%Dij%self,comm=xg_nonlop%comm_atom)
552 : end if
553 :
554 2188 : call timab(tim_make_Dij,2,tsec)
555 :
556 2188 : end subroutine xg_nonlop_make_Dij
557 : !!***
558 :
559 2252 : subroutine xg_nonlop_set_Dij_spin(xg_nonlop,isppol)
560 :
561 : integer,intent(in) :: isppol
562 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
563 :
564 : integer :: shift,nlmn_max,nspinor,natom
565 :
566 2252 : nspinor = xg_nonlop%nspinor
567 2252 : natom = xg_nonlop%natom
568 2252 : nlmn_max = xg_nonlop%nlmn_max
569 :
570 2252 : if (isppol/=1) then ! isppol must be 1 or 2 if nspinor==1, and must be 1 of nspinor==2
571 64 : if (isppol/=2.or.nspinor/=1) then
572 0 : ABI_ERROR('wrong isppol')
573 : end if
574 : end if
575 :
576 2252 : shift=1+(isppol-1)*natom*nspinor*nlmn_max
577 2252 : call xg_setBlock(xg_nonlop%Dij,xg_nonlop%Dij_spin,nspinor*nlmn_max,nspinor*nlmn_max*natom,fcol=shift)
578 :
579 2252 : end subroutine xg_nonlop_set_Dij_spin
580 : !!***
581 :
582 138 : subroutine xg_nonlop_make_Sij(xg_nonlop,pawtab,inv_sij)
583 :
584 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
585 : type(pawtab_type),intent(in) :: pawtab(:)
586 : logical,optional,intent(in) :: inv_sij
587 :
588 : logical :: inv_sij_
589 : integer :: itypat, nlmn, nlmn_max, ntypat, shift
590 : integer :: ilmn, jlmn, j0lmn, jjlmn, ijlmn
591 138 : real(dp),pointer :: Sij_itypat_(:,:)
592 : real(dp) :: tsec(2)
593 : type(xg_t) :: work
594 : type(xgBlock_t) :: Sij_itypat,Sijm1_itypat
595 :
596 : ! *************************************************************************
597 :
598 138 : call timab(tim_make_Sij,1,tsec)
599 :
600 138 : if (.not.xg_nonlop%paw) then
601 0 : ABI_ERROR('Not implemented with paw=False.')
602 : end if
603 :
604 138 : ntypat = xg_nonlop%ntypat
605 138 : nlmn_max = xg_nonlop%nlmn_max
606 :
607 138 : call xg_init(xg_nonlop%Sij,SPACE_R,nlmn_max,nlmn_max*ntypat,xmpi_comm_self)
608 :
609 138 : inv_sij_ = .false.
610 138 : if (present(inv_sij)) inv_sij_ = inv_sij
611 :
612 138 : if (inv_sij_) then
613 102 : call xg_init(xg_nonlop%Sijm1,SPACE_R,nlmn_max,nlmn_max*ntypat,xmpi_comm_self)
614 : end if
615 :
616 404 : do itypat=1,ntypat
617 :
618 266 : nlmn=xg_nonlop%nlmn_ntypat(itypat)
619 :
620 266 : shift=1+(itypat-1)*nlmn_max
621 266 : call xg_setBlock(xg_nonlop%Sij,Sij_itypat,nlmn,nlmn,fcol=shift)
622 266 : call xgBlock_reverseMap(Sij_itypat,Sij_itypat_)
623 2584 : do jlmn=1,nlmn
624 2318 : j0lmn=jlmn*(jlmn-1)/2
625 2318 : jjlmn=j0lmn+jlmn
626 2318 : Sij_itypat_(jlmn,jlmn) = pawtab(itypat)%sij(jjlmn)
627 12407 : do ilmn=1,jlmn-1
628 9823 : ijlmn=j0lmn+ilmn
629 9823 : Sij_itypat_(ilmn,jlmn) = pawtab(itypat)%sij(ijlmn)
630 12141 : Sij_itypat_(jlmn,ilmn) = pawtab(itypat)%sij(ijlmn)
631 : end do
632 : end do
633 404 : if (inv_sij_) then
634 207 : call xg_init(work,SPACE_R,nlmn,nlmn,xmpi_comm_self)
635 207 : call xg_setBlock(xg_nonlop%Sijm1,Sijm1_itypat,nlmn,nlmn,fcol=shift)
636 207 : call xgBlock_invert_sy(Sijm1_itypat,work%self,xg_input=Sij_itypat)
637 207 : call xg_free(work)
638 : end if
639 :
640 : end do
641 :
642 138 : call timab(tim_make_Sij,2,tsec)
643 :
644 138 : end subroutine xg_nonlop_make_Sij
645 : !!***
646 :
647 24 : subroutine xg_nonlop_destroy_ekb(xg_nonlop)
648 :
649 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
650 :
651 : ! *************************************************************************
652 :
653 24 : if (xg_nonlop%paw) then
654 0 : ABI_ERROR('Not implemented with paw=True.')
655 : end if
656 :
657 24 : call xg_free(xg_nonlop%ekb)
658 :
659 24 : end subroutine xg_nonlop_destroy_ekb
660 : !!***
661 :
662 2326 : subroutine xg_nonlop_destroy_Dij(xg_nonlop)
663 :
664 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
665 :
666 : ! *************************************************************************
667 :
668 2326 : if (.not.xg_nonlop%paw) then
669 0 : ABI_ERROR('Not implemented with paw=False.')
670 : end if
671 :
672 2326 : call xg_free(xg_nonlop%Dij)
673 :
674 2326 : end subroutine xg_nonlop_destroy_Dij
675 : !!***
676 :
677 138 : subroutine xg_nonlop_destroy_Sij(xg_nonlop)
678 :
679 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
680 :
681 : ! *************************************************************************
682 :
683 138 : if (.not.xg_nonlop%paw) then
684 0 : ABI_ERROR('Not implemented with paw=False.')
685 : end if
686 :
687 138 : call xg_free(xg_nonlop%Sij)
688 138 : call xg_free(xg_nonlop%Sijm1)
689 :
690 138 : end subroutine xg_nonlop_destroy_Sij
691 : !!***
692 :
693 2320 : subroutine xg_nonlop_compute_projs(xg_nonlop)
694 :
695 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
696 :
697 2320 : complex(dp),pointer :: projectors_k_(:,:)
698 2320 : real(dp),pointer :: projectors_k_real(:,:)
699 2320 : real(dp),pointer :: ffnl_gather_k_(:,:)
700 2320 : complex(dp),pointer :: ph3d_gather_k_(:,:)
701 2320 : real(dp),pointer :: ph3d_gather_k_real(:,:)
702 :
703 : integer :: shift_itypat,shift_itypat_nlmn,ntypat,nattyp_i,shift_ipw
704 : integer :: icol,ilmn,nlmn,il,ipw,ia,iatom,itypat
705 : real(dp) :: ffnl_ipw
706 : complex(dp) :: cil(4),ph3d_ipw,ctmp
707 : logical :: compute_gather
708 :
709 2320 : ntypat = xg_nonlop%ntypat
710 :
711 : ! 4pi/sqrt(ucvol) * (-i)^l
712 2320 : cil(1) = ( 1.0_DP, 0.0_DP) * xg_nonlop%weight
713 2320 : cil(2) = ( 0.0_DP,-1.0_DP) * xg_nonlop%weight
714 2320 : cil(3) = (-1.0_DP, 0.0_DP) * xg_nonlop%weight
715 2320 : cil(4) = ( 0.0_DP, 1.0_DP) * xg_nonlop%weight
716 :
717 2320 : compute_gather = xg_nonlop%option==0
718 :
719 2320 : if (compute_gather) then
720 2192 : call xgBlock_zero(xg_nonlop%ffnl_gather_k%self)
721 2192 : call xgBlock_zero(xg_nonlop%ph3d_gather_k%self)
722 2192 : call xgBlock_reverseMap(xg_nonlop%ffnl_gather_k%self,ffnl_gather_k_)
723 : end if
724 :
725 2320 : shift_ipw = xg_nonlop%l_shift_npw_k(xg_nonlop%me_band+1)
726 :
727 3654 : select case(xg_nonlop%space_pw)
728 :
729 : case (SPACE_C)
730 :
731 1334 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_)
732 1334 : if (compute_gather) then
733 1206 : call xgBlock_reverseMap(xg_nonlop%ph3d_gather_k%self,ph3d_gather_k_)
734 : end if
735 : shift_itypat=0
736 : shift_itypat_nlmn=0
737 : !$omp parallel default (none) &
738 : !$omp& shared(xg_nonlop,ph3d_gather_k_,ffnl_gather_k_,projectors_k_), &
739 : !$omp& firstprivate(compute_gather,ntypat,shift_ipw,shift_itypat,shift_itypat_nlmn,cil), &
740 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,ipw,iatom,il,ph3d_ipw,ffnl_ipw,icol)
741 4048 : do itypat = 1, ntypat
742 2714 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
743 2714 : nattyp_i = xg_nonlop%nattyp(itypat)
744 4048 : if (nattyp_i>0) then
745 : !! projectors = 4pi/sqrt(ucvol) * (-i)^l * conj(ph3d) * ffnl
746 : !$omp do collapse(3)
747 5350 : do ia = 1, nattyp_i
748 30478 : do ilmn=1,nlmn
749 2192660 : do ipw=1,xg_nonlop%npw_k
750 2164832 : iatom = ia + shift_itypat
751 2164832 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
752 2164832 : ph3d_ipw = cmplx( xg_nonlop%ph3d_k(1,ipw,iatom), xg_nonlop%ph3d_k(2,ipw,iatom), kind=DP)
753 2164832 : ffnl_ipw = xg_nonlop%ffnl_k(ipw, 1, ilmn, itypat)
754 : !
755 2164832 : if (compute_gather) then
756 2096544 : ph3d_gather_k_(ipw+shift_ipw,iatom) = ph3d_ipw
757 2096544 : ffnl_gather_k_(ipw+shift_ipw,ilmn+(itypat-1)*xg_nonlop%nlmn_max) = ffnl_ipw
758 : end if
759 : !
760 2164832 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
761 2189960 : projectors_k_(ipw,icol) = cil(il) * conjg(ph3d_ipw) * ffnl_ipw
762 : end do
763 : end do
764 : end do
765 : !$omp end do
766 2650 : shift_itypat = shift_itypat + nattyp_i
767 2650 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
768 : end if
769 : end do
770 : !$omp end parallel
771 :
772 : case (SPACE_CR)
773 :
774 986 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_real)
775 986 : if (compute_gather) then
776 986 : call xgBlock_reverseMap(xg_nonlop%ph3d_gather_k%self,ph3d_gather_k_real)
777 : end if
778 : shift_itypat=0
779 : shift_itypat_nlmn=0
780 : !$omp parallel default (none) &
781 : !$omp& shared(xg_nonlop,ph3d_gather_k_real,ffnl_gather_k_,projectors_k_real), &
782 : !$omp& firstprivate(compute_gather,ntypat,shift_ipw,shift_itypat,shift_itypat_nlmn,cil), &
783 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,ipw,iatom,il,ph3d_ipw,ffnl_ipw,icol,ctmp)
784 2948 : do itypat = 1, ntypat
785 1962 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
786 1962 : nattyp_i = xg_nonlop%nattyp(itypat)
787 2948 : if (nattyp_i>0) then
788 : !! projectors = 4pi/sqrt(ucvol) * (-i)^l * conj(ph3d) * ffnl
789 : !$omp do collapse(3)
790 3934 : do ia = 1, nattyp_i
791 34950 : do ilmn=1,nlmn
792 952300 : do ipw=1,xg_nonlop%npw_k
793 919312 : iatom = ia + shift_itypat
794 919312 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
795 919312 : ph3d_ipw = cmplx( xg_nonlop%ph3d_k(1,ipw,iatom), xg_nonlop%ph3d_k(2,ipw,iatom), kind=DP)
796 919312 : ffnl_ipw = xg_nonlop%ffnl_k(ipw, 1, ilmn, itypat)
797 : !
798 919312 : if (compute_gather) then
799 919312 : ph3d_gather_k_real(2*(ipw+shift_ipw)-1,iatom) = dble(ph3d_ipw)
800 919312 : ph3d_gather_k_real(2*(ipw+shift_ipw) ,iatom) = dimag(ph3d_ipw)
801 919312 : ffnl_gather_k_(ipw+shift_ipw,ilmn+(itypat-1)*xg_nonlop%nlmn_max) = ffnl_ipw
802 : end if
803 : !
804 919312 : ctmp = cil(il) * conjg(ph3d_ipw) * ffnl_ipw
805 919312 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
806 919312 : projectors_k_real(2*ipw-1,icol) = dble(ctmp)
807 950328 : projectors_k_real(2*ipw ,icol) = dimag(ctmp)
808 : end do
809 : end do
810 : end do
811 : !$omp end do
812 1962 : shift_itypat = shift_itypat + nattyp_i
813 1962 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
814 : end if
815 : end do
816 : !$omp end parallel
817 :
818 : case default
819 2320 : ABI_ERROR("Wrong space")
820 :
821 : end select
822 :
823 2320 : if (compute_gather) then
824 2192 : call xgBlock_mpi_sum(xg_nonlop%ffnl_gather_k%self,comm=xg_nonlop%comm_band)
825 2192 : call xgBlock_mpi_sum(xg_nonlop%ph3d_gather_k%self,comm=xg_nonlop%comm_band)
826 : end if
827 :
828 2320 : end subroutine xg_nonlop_compute_projs
829 : !!***
830 :
831 18300 : subroutine xg_nonlop_compute_projs_deriv_atom(xg_nonlop,projs_deriv_atom)
832 :
833 : type(xgBlock_t),intent(inout) :: projs_deriv_atom
834 : type(xg_nonlop_t),intent(in) :: xg_nonlop
835 :
836 18300 : complex(dp),pointer :: projectors_k_(:,:)
837 18300 : complex(dp),pointer :: projectors_deriv_atom_k_(:,:)
838 18300 : real(dp),pointer :: projectors_k_real(:,:)
839 18300 : real(dp),pointer :: projectors_deriv_atom_k_real(:,:)
840 :
841 : integer :: shift_itypat_nlmn,shift_itypat_3nlmn,ntypat,nattyp_i
842 : integer :: icol,icol_deriv,ilmn,nlmn,ipw,ia,itypat,idir
843 : complex(dp) :: ctmp
844 : real(dp) :: tmp,proj_deriv_ipw_re,proj_deriv_ipw_im
845 :
846 18300 : if (.not.associated(xg_nonlop%projectors_k)) then
847 0 : ABI_ERROR('projectors_k should be associated')
848 : end if
849 :
850 18300 : if (rows(projs_deriv_atom)/=xg_nonlop%npw_k) then
851 0 : ABI_ERROR('rows(projs_deriv_atom)/=npw_k')
852 : end if
853 18300 : if (cols(projs_deriv_atom)/=3*xg_nonlop%cprjdim) then
854 0 : ABI_ERROR('cols(projs_deriv_atom)/=3*cprjdim')
855 : end if
856 :
857 18300 : ntypat = xg_nonlop%ntypat
858 :
859 28903 : select case(xg_nonlop%space_pw)
860 :
861 : case (SPACE_C)
862 :
863 10603 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_)
864 10603 : call xgBlock_reverseMap(projs_deriv_atom,projectors_deriv_atom_k_)
865 10603 : shift_itypat_nlmn=0
866 10603 : shift_itypat_3nlmn=0
867 : !$omp parallel default (none) &
868 : !$omp& shared(xg_nonlop,projectors_k_,projectors_deriv_atom_k_), &
869 : !$omp& firstprivate(ntypat,shift_itypat_nlmn,shift_itypat_3nlmn), &
870 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,ipw,idir,icol,icol_deriv,ctmp)
871 32376 : do itypat = 1, ntypat
872 21773 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
873 21773 : nattyp_i = xg_nonlop%nattyp(itypat)
874 32376 : if (nattyp_i>0) then
875 : !! projectors_deriv_atom(k+G) = -i * 2pi * (k+G)_idir * projectors(k+G)
876 : !$omp do collapse(4)
877 42379 : do ia = 1, nattyp_i
878 230887 : do ilmn=1,nlmn
879 16228673 : do ipw=1,xg_nonlop%npw_k
880 64264280 : do idir=1,3
881 48056829 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
882 48056829 : ctmp = ( 0.0_DP, -1.0_DP) * two_pi * xg_nonlop%kpg_k(ipw,idir)
883 48056829 : icol_deriv = ilmn + (idir-1)*nlmn + (ia-1)*3*nlmn + shift_itypat_3nlmn
884 64075772 : projectors_deriv_atom_k_(ipw,icol_deriv) = ctmp * projectors_k_(ipw,icol)
885 : end do
886 : end do
887 : end do
888 : end do
889 : !$omp end do
890 21157 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
891 21157 : shift_itypat_3nlmn = shift_itypat_3nlmn + nattyp_i*3*nlmn
892 : end if
893 : end do
894 : !$omp end parallel
895 :
896 : case (SPACE_CR)
897 :
898 7697 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_real)
899 7697 : call xgBlock_reverseMap(projs_deriv_atom,projectors_deriv_atom_k_real)
900 7697 : shift_itypat_nlmn=0
901 7697 : shift_itypat_3nlmn=0
902 : !$omp parallel default (none) &
903 : !$omp& shared(xg_nonlop,projectors_k_real,projectors_deriv_atom_k_real), &
904 : !$omp& firstprivate(ntypat,shift_itypat_nlmn,shift_itypat_3nlmn), &
905 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,ipw,idir,icol,icol_deriv,tmp), &
906 : !$omp& private(proj_deriv_ipw_re,proj_deriv_ipw_im)
907 23074 : do itypat = 1, ntypat
908 15377 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
909 15377 : nattyp_i = xg_nonlop%nattyp(itypat)
910 23074 : if (nattyp_i>0) then
911 : !! projectors_deriv_atom(k+G) = -i * 2pi * (k+G)_idir * projectors(k+G)
912 : !$omp do collapse(4)
913 30771 : do ia = 1, nattyp_i
914 220063 : do ilmn=1,nlmn
915 6104064 : do ipw=1,xg_nonlop%npw_k
916 23786804 : do idir=1,3
917 17698134 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
918 17698134 : tmp = - two_pi * xg_nonlop%kpg_k(ipw,idir)
919 17698134 : icol_deriv = ilmn + (idir-1)*nlmn + (ia-1)*3*nlmn + shift_itypat_3nlmn
920 : !! Re(projectors_deriv_atom) = 2pi * (k+G)_idir * Im(projectors)
921 : !! Im(projectors_deriv_atom) = -2pi * (k+G)_idir * Re(projectors)
922 17698134 : proj_deriv_ipw_re = - tmp * projectors_k_real(2*ipw ,icol)
923 17698134 : proj_deriv_ipw_im = tmp * projectors_k_real(2*ipw-1,icol)
924 17698134 : projectors_deriv_atom_k_real(2*ipw-1,icol_deriv) = proj_deriv_ipw_re
925 23597512 : projectors_deriv_atom_k_real(2*ipw ,icol_deriv) = proj_deriv_ipw_im
926 : end do
927 : end do
928 : end do
929 : end do
930 : !$omp end do
931 15377 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
932 15377 : shift_itypat_3nlmn = shift_itypat_3nlmn + nattyp_i*3*nlmn
933 : end if
934 : end do
935 : !$omp end parallel
936 :
937 : case default
938 18300 : ABI_ERROR("Wrong space")
939 :
940 : end select
941 :
942 18300 : end subroutine xg_nonlop_compute_projs_deriv_atom
943 : !!***
944 :
945 1140 : subroutine xg_nonlop_compute_projs_deriv_stress(xg_nonlop,projs_deriv_stress)
946 :
947 : type(xgBlock_t),intent(inout) :: projs_deriv_stress
948 : type(xg_nonlop_t),intent(in) :: xg_nonlop
949 :
950 1140 : complex(dp),pointer :: projectors_k_(:,:)
951 1140 : complex(dp),pointer :: projectors_deriv_stress_k_(:,:)
952 1140 : real(dp),pointer :: projectors_k_real(:,:)
953 1140 : real(dp),pointer :: projectors_deriv_stress_k_real(:,:)
954 :
955 : integer :: shift_itypat,shift_itypat_nlmn,shift_itypat_6nlmn,ntypat,nattyp_i
956 : integer :: iatom,icol_shift,icol_deriv,ilmn,nlmn,ipw,ia,itypat,idir,il
957 : complex(dp) :: ctmp(3),cil(4),ph3d_ipw,cipw
958 : real(dp) :: ffnl_ipw(3)
959 :
960 1140 : if (rows(projs_deriv_stress)/=xg_nonlop%npw_k) then
961 0 : ABI_ERROR('rows(projs_deriv_atom)/=npw_k')
962 : end if
963 1140 : if (cols(projs_deriv_stress)/=6*xg_nonlop%cprjdim) then
964 0 : ABI_ERROR('cols(projs_deriv_stress)/=6*cprjdim')
965 : end if
966 :
967 1140 : ntypat = xg_nonlop%ntypat
968 :
969 : ! 4pi/sqrt(ucvol) * (-i)^l
970 1140 : cil(1) = ( 1.0_DP, 0.0_DP) * xg_nonlop%weight
971 1140 : cil(2) = ( 0.0_DP,-1.0_DP) * xg_nonlop%weight
972 1140 : cil(3) = (-1.0_DP, 0.0_DP) * xg_nonlop%weight
973 1140 : cil(4) = ( 0.0_DP, 1.0_DP) * xg_nonlop%weight
974 :
975 1871 : select case(xg_nonlop%space_pw)
976 :
977 : case (SPACE_C)
978 :
979 731 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_)
980 731 : call xgBlock_reverseMap(projs_deriv_stress,projectors_deriv_stress_k_)
981 731 : shift_itypat=0
982 731 : shift_itypat_nlmn=0
983 731 : shift_itypat_6nlmn=0
984 : !$omp parallel default (none) &
985 : !$omp& shared(xg_nonlop,projectors_k_,projectors_deriv_stress_k_), &
986 : !$omp& firstprivate(cil,ntypat,shift_itypat,shift_itypat_nlmn,shift_itypat_6nlmn), &
987 : !$omp& private(il,iatom,itypat,nattyp_i,nlmn,ia,ilmn,ipw,idir,icol_shift,icol_deriv), &
988 : !$omp& private(ctmp,ph3d_ipw,ffnl_ipw)
989 2176 : do itypat = 1, ntypat
990 1445 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
991 1445 : nattyp_i = xg_nonlop%nattyp(itypat)
992 2176 : if (nattyp_i>0) then
993 : !! projectors_deriv_stress(k+G)_ab = 4pi/sqrt(ucvol) * (-i)^l * conj(ph3d) * (-(k+G)_b) * d/d(K_a)[ffnl_deriv(k+G)]
994 : !$omp do collapse(3)
995 2891 : do ia = 1, nattyp_i
996 17279 : do ilmn=1,nlmn
997 1381698 : do ipw=1,xg_nonlop%npw_k
998 1365832 : iatom = ia + shift_itypat
999 1365832 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
1000 1365832 : ph3d_ipw = cmplx( xg_nonlop%ph3d_k(1,ipw,iatom), xg_nonlop%ph3d_k(2,ipw,iatom), kind=DP)
1001 5463328 : ffnl_ipw(:) = xg_nonlop%ffnl_k(ipw, 2:4, ilmn, itypat)
1002 5463328 : ctmp(:) = - cil(il) * conjg(ph3d_ipw) * xg_nonlop%kpg_k(ipw,:)
1003 1365832 : icol_shift = ilmn + (ia-1)*6*nlmn + shift_itypat_6nlmn
1004 : ! diagonal part
1005 5463328 : do idir=1,3
1006 4097496 : icol_deriv = icol_shift + (idir-1)*nlmn
1007 5463328 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(idir) * ffnl_ipw(idir)
1008 : end do
1009 : ! off-diagonal part (which is symmetric)
1010 5463328 : ctmp(:) = half*ctmp(:)
1011 1365832 : icol_deriv = icol_shift + (4-1)*nlmn
1012 1365832 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(2) * ffnl_ipw(3) + ctmp(3) * ffnl_ipw(2)
1013 1365832 : icol_deriv = icol_shift + (5-1)*nlmn
1014 1365832 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(1) * ffnl_ipw(3) + ctmp(3) * ffnl_ipw(1)
1015 1365832 : icol_deriv = icol_shift + (6-1)*nlmn
1016 1380220 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(1) * ffnl_ipw(2) + ctmp(2) * ffnl_ipw(1)
1017 : end do
1018 : end do
1019 : end do
1020 : !$omp end do
1021 1413 : shift_itypat = shift_itypat + nattyp_i
1022 1413 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
1023 1413 : shift_itypat_6nlmn = shift_itypat_6nlmn + nattyp_i*6*nlmn
1024 : end if
1025 : end do
1026 : !$omp end parallel
1027 :
1028 : case (SPACE_CR)
1029 :
1030 409 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_real)
1031 409 : call xgBlock_reverseMap(projs_deriv_stress,projectors_deriv_stress_k_real)
1032 409 : shift_itypat=0
1033 409 : shift_itypat_nlmn=0
1034 409 : shift_itypat_6nlmn=0
1035 : !$omp parallel default (none) &
1036 : !$omp& shared(xg_nonlop,projectors_k_real,projectors_deriv_stress_k_real), &
1037 : !$omp& firstprivate(cil,ntypat,shift_itypat,shift_itypat_nlmn,shift_itypat_6nlmn), &
1038 : !$omp& private(il,iatom,itypat,nattyp_i,nlmn,ia,ilmn,ipw,idir,icol_shift,icol_deriv), &
1039 : !$omp& private(ctmp,cipw,ph3d_ipw,ffnl_ipw)
1040 1210 : do itypat = 1, ntypat
1041 801 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
1042 801 : nattyp_i = xg_nonlop%nattyp(itypat)
1043 1210 : if (nattyp_i>0) then
1044 : !! projectors_deriv_stress(k+G)_ab = 4pi/sqrt(ucvol) * (-i)^l * conj(ph3d) * (-(k+G)_b) * d/d(K_a)[ffnl_deriv(k+G)]
1045 : !$omp do collapse(3)
1046 1619 : do ia = 1, nattyp_i
1047 14519 : do ilmn=1,nlmn
1048 454776 : do ipw=1,xg_nonlop%npw_k
1049 441058 : iatom = ia + shift_itypat
1050 441058 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
1051 441058 : ph3d_ipw = cmplx( xg_nonlop%ph3d_k(1,ipw,iatom), xg_nonlop%ph3d_k(2,ipw,iatom), kind=DP)
1052 1764232 : ffnl_ipw(:) = xg_nonlop%ffnl_k(ipw, 2:4, ilmn, itypat)
1053 1764232 : ctmp(:) = - cil(il) * conjg(ph3d_ipw) * xg_nonlop%kpg_k(ipw,:)
1054 441058 : icol_shift = ilmn + (ia-1)*6*nlmn + shift_itypat_6nlmn
1055 : ! diagonal part
1056 1764232 : do idir=1,3
1057 1323174 : icol_deriv = icol_shift + (idir-1)*nlmn
1058 1323174 : projectors_deriv_stress_k_real(2*ipw-1,icol_deriv) = dble(ctmp(idir)) * ffnl_ipw(idir)
1059 1764232 : projectors_deriv_stress_k_real(2*ipw ,icol_deriv) = dimag(ctmp(idir)) * ffnl_ipw(idir)
1060 : end do
1061 : ! off-diagonal part (which is symmetric)
1062 1764232 : ctmp(:) = half*ctmp(:)
1063 :
1064 441058 : icol_deriv = icol_shift + (4-1)*nlmn
1065 441058 : cipw = ctmp(2) * ffnl_ipw(3) + ctmp(3) * ffnl_ipw(2)
1066 441058 : projectors_deriv_stress_k_real(2*ipw-1,icol_deriv) = dble(cipw)
1067 441058 : projectors_deriv_stress_k_real(2*ipw ,icol_deriv) = dimag(cipw)
1068 :
1069 441058 : icol_deriv = icol_shift + (5-1)*nlmn
1070 441058 : cipw = ctmp(1) * ffnl_ipw(3) + ctmp(3) * ffnl_ipw(1)
1071 441058 : projectors_deriv_stress_k_real(2*ipw-1,icol_deriv) = dble(cipw)
1072 441058 : projectors_deriv_stress_k_real(2*ipw ,icol_deriv) = dimag(cipw)
1073 :
1074 441058 : icol_deriv = icol_shift + (6-1)*nlmn
1075 441058 : cipw = ctmp(1) * ffnl_ipw(2) + ctmp(2) * ffnl_ipw(1)
1076 441058 : projectors_deriv_stress_k_real(2*ipw-1,icol_deriv) = dble(cipw)
1077 453958 : projectors_deriv_stress_k_real(2*ipw ,icol_deriv) = dimag(cipw)
1078 :
1079 : end do
1080 : end do
1081 : end do
1082 : !$omp end do
1083 801 : shift_itypat = shift_itypat + nattyp_i
1084 801 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
1085 801 : shift_itypat_6nlmn = shift_itypat_6nlmn + nattyp_i*6*nlmn
1086 : end if
1087 : end do
1088 : !$omp end parallel
1089 :
1090 : case default
1091 1140 : ABI_ERROR("Wrong space")
1092 :
1093 : end select
1094 :
1095 1140 : end subroutine xg_nonlop_compute_projs_deriv_stress
1096 : !!***
1097 :
1098 865236 : subroutine xg_nonlop_compute_projs_otf(xg_nonlop,projs_otf,index_mpi)
1099 :
1100 : integer,intent(in) :: index_mpi
1101 : type(xg_nonlop_t),intent(in) :: xg_nonlop
1102 : type(xgBlock_t),intent(inout) :: projs_otf
1103 :
1104 865236 : complex(dp),pointer :: projectors_k_(:,:)
1105 865236 : real(dp),pointer :: projectors_k_real(:,:)
1106 865236 : real(dp),pointer :: ffnl_gather_k_(:,:)
1107 865236 : complex(dp),pointer :: ph3d_gather_k_(:,:)
1108 865236 : real(dp),pointer :: ph3d_gather_k_real(:,:)
1109 :
1110 : integer :: nmpi,npw_k
1111 : integer :: shift_itypat,shift_itypat_nlmn,ntypat,nattyp_i,shift_ipw
1112 : integer :: icol,ilmn,nlmn,il,ipw,ia,iatom,itypat
1113 : real(dp) :: ffnl_ipw,ph3d_ipw_r(2)
1114 : complex(dp) :: cil(4),ph3d_ipw,ctmp
1115 :
1116 865236 : if (xg_nonlop%option/=0) then
1117 0 : ABI_ERROR('xg_nonlop%option/=0')
1118 : end if
1119 :
1120 865236 : nmpi = xmpi_comm_size(xg_nonlop%comm_band)
1121 865236 : if (index_mpi<0.or.index_mpi>nmpi-1) then
1122 0 : ABI_ERROR('index_mpi should be between 0 and size(comm_band)-1')
1123 : end if
1124 :
1125 865236 : npw_k = xg_nonlop%l_npw_k(index_mpi+1)
1126 865236 : shift_ipw = xg_nonlop%l_shift_npw_k(index_mpi+1)
1127 :
1128 865236 : if (rows(projs_otf)/=npw_k) then
1129 0 : ABI_ERROR('rows(projs_otf)/=npw_k!')
1130 : end if
1131 865236 : if (cols(projs_otf)/=xg_nonlop%cprjdim) then
1132 0 : ABI_ERROR('cols(projs_otf)/=cprjdim!')
1133 : end if
1134 :
1135 865236 : ntypat = xg_nonlop%ntypat
1136 :
1137 : ! 4pi/sqrt(ucvol) * (-i)^l
1138 865236 : cil(1) = ( 1.0_DP, 0.0_DP) * xg_nonlop%weight
1139 865236 : cil(2) = ( 0.0_DP,-1.0_DP) * xg_nonlop%weight
1140 865236 : cil(3) = (-1.0_DP, 0.0_DP) * xg_nonlop%weight
1141 865236 : cil(4) = ( 0.0_DP, 1.0_DP) * xg_nonlop%weight
1142 :
1143 865236 : call xgBlock_reverseMap(xg_nonlop%ffnl_gather_k%self,ffnl_gather_k_)
1144 :
1145 1273380 : select case(xg_nonlop%space_pw)
1146 :
1147 : case (SPACE_C)
1148 :
1149 408144 : call xgBlock_reverseMap(projs_otf,projectors_k_)
1150 408144 : call xgBlock_reverseMap(xg_nonlop%ph3d_gather_k%self,ph3d_gather_k_)
1151 408144 : shift_itypat=0
1152 408144 : shift_itypat_nlmn=0
1153 : !$omp parallel default (none) &
1154 : !$omp& shared(xg_nonlop,ph3d_gather_k_,ffnl_gather_k_,projectors_k_), &
1155 : !$omp& firstprivate(ntypat,npw_k,shift_ipw,shift_itypat,shift_itypat_nlmn,cil), &
1156 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,ipw,iatom,il,ph3d_ipw,ffnl_ipw,icol)
1157 1224432 : do itypat = 1, ntypat
1158 816288 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
1159 816288 : nattyp_i = xg_nonlop%nattyp(itypat)
1160 1224432 : if (nattyp_i>0) then
1161 : !! projectors = 4pi/sqrt(ucvol)* conj(ph3d) * ffnl * (-i)^l
1162 : !$omp do collapse(3)
1163 1632576 : do ia = 1, nattyp_i
1164 8162880 : do ilmn=1,nlmn
1165 239950512 : do ipw=1,npw_k
1166 232603920 : iatom = ia + shift_itypat
1167 232603920 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
1168 232603920 : ph3d_ipw = ph3d_gather_k_(ipw+shift_ipw,iatom)
1169 232603920 : ffnl_ipw = ffnl_gather_k_(ipw+shift_ipw,ilmn+(itypat-1)*xg_nonlop%nlmn_max)
1170 : !
1171 232603920 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
1172 239134224 : projectors_k_(ipw,icol) = cil(il) * conjg(ph3d_ipw) * ffnl_ipw
1173 : end do
1174 : end do
1175 : end do
1176 : !$omp end do
1177 816288 : shift_itypat = shift_itypat + nattyp_i
1178 816288 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
1179 : end if
1180 : end do
1181 : !$omp end parallel
1182 :
1183 : case(SPACE_CR)
1184 :
1185 457092 : call xgBlock_reverseMap(projs_otf,projectors_k_real)
1186 457092 : call xgBlock_reverseMap(xg_nonlop%ph3d_gather_k%self,ph3d_gather_k_real)
1187 457092 : shift_itypat=0
1188 457092 : shift_itypat_nlmn=0
1189 : !$omp parallel default (none) &
1190 : !$omp& shared(xg_nonlop,ph3d_gather_k_real,ffnl_gather_k_,projectors_k_real), &
1191 : !$omp& firstprivate(ntypat,npw_k,shift_ipw,shift_itypat,shift_itypat_nlmn,cil), &
1192 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,ipw,iatom,il,ph3d_ipw_r,ffnl_ipw,icol,ctmp)
1193 1371276 : do itypat = 1, ntypat
1194 914184 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
1195 914184 : nattyp_i = xg_nonlop%nattyp(itypat)
1196 1371276 : if (nattyp_i>0) then
1197 : !! projectors = 4pi/sqrt(ucvol)* conj(ph3d) * ffnl * (-i)^l
1198 : !$omp do collapse(3)
1199 1828368 : do ia = 1, nattyp_i
1200 14977548 : do ilmn=1,nlmn
1201 247695219 : do ipw=1,npw_k
1202 233631855 : iatom = ia + shift_itypat
1203 233631855 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
1204 233631855 : ph3d_ipw_r(1) = ph3d_gather_k_real(2*(ipw+shift_ipw)-1,iatom)
1205 233631855 : ph3d_ipw_r(2) = ph3d_gather_k_real(2*(ipw+shift_ipw) ,iatom)
1206 233631855 : ffnl_ipw = ffnl_gather_k_(ipw+shift_ipw,ilmn+(itypat-1)*xg_nonlop%nlmn_max)
1207 : !
1208 233631855 : ctmp = cmplx( ph3d_ipw_r(1), ph3d_ipw_r(2), kind=DP)
1209 233631855 : ctmp = cil(il) * conjg(ctmp) * ffnl_ipw
1210 233631855 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
1211 233631855 : projectors_k_real(2*ipw-1,icol) = dble(ctmp)
1212 246781035 : projectors_k_real(2*ipw ,icol) = dimag(ctmp)
1213 : end do
1214 : end do
1215 : end do
1216 : !$omp end do
1217 914184 : shift_itypat = shift_itypat + nattyp_i
1218 914184 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
1219 : end if
1220 : end do
1221 : !$omp end parallel
1222 :
1223 : case default
1224 865236 : ABI_ERROR("Wrong space")
1225 :
1226 : end select
1227 :
1228 865236 : end subroutine xg_nonlop_compute_projs_otf
1229 : !!***
1230 :
1231 18244 : subroutine xg_nonlop_make_k(xg_nonlop,ikpt,istwf_k,me_g0,me_g0_fft,npw_k,ffnl_k,ph3d_k,kpg_k,compute_proj,&
1232 : compute_invS_approx,compute_gram)
1233 :
1234 : type(xg_nonlop_t),intent(inout) :: xg_nonlop
1235 :
1236 : logical ,intent(in) :: compute_proj
1237 : integer ,intent(in) :: ikpt
1238 : integer ,intent(in) :: istwf_k
1239 : integer, intent(in) :: me_g0
1240 : integer, intent(in) :: me_g0_fft
1241 : integer, intent(in) :: npw_k
1242 : real(dp), intent(in), target :: ffnl_k(:,:,:,:)
1243 : real(dp), intent(in), target :: ph3d_k(:,:,:)
1244 : real(dp), intent(in), target :: kpg_k(:,:)
1245 : logical ,optional,intent(in) :: compute_invS_approx
1246 : logical ,optional,intent(in) :: compute_gram
1247 :
1248 : logical :: compute_gram_,compute_invS_approx_
1249 18244 : real(dp),pointer :: gram_proj_k_(:,:),Sijm1_(:,:)
1250 : integer :: ierr, iblock, shift, shiftc, shift_sij, shift_itypat, itypat, ilmn, jlmn, nlmn, nlmn_max, ia
1251 : integer :: cplex,nattyp_i,ntypat,ncols,nmpi,me_g0_loc,me_g0_fft_loc,space_cprj
1252 : !integer :: cols,cond
1253 : real(dp) :: tsec(2)
1254 : type(xg_t) :: work
1255 : type(xgBlock_t) :: projs,invSij_approx_k_itypat,Sijm1_itypat
1256 :
1257 : ! *************************************************************************
1258 :
1259 18244 : call timab(tim_make_k,1,tsec)
1260 :
1261 18244 : me_g0_loc = -1
1262 18244 : me_g0_fft_loc = -1
1263 18244 : if (istwf_k==1) then
1264 10357 : xg_nonlop%cplex=2
1265 10357 : xg_nonlop%space_pw=SPACE_C
1266 10357 : xg_nonlop%space_cprj=SPACE_C
1267 : else ! istwf_k>1
1268 7887 : xg_nonlop%cplex=1
1269 7887 : xg_nonlop%space_pw=SPACE_CR
1270 7887 : xg_nonlop%space_cprj=SPACE_R
1271 7887 : me_g0_loc = 0
1272 7887 : me_g0_fft_loc = 0
1273 7887 : if (istwf_k==2.and.me_g0==1) then
1274 431 : me_g0_loc = me_g0
1275 431 : me_g0_fft_loc = me_g0_fft
1276 : end if
1277 : end if
1278 :
1279 18244 : cplex = xg_nonlop%cplex
1280 18244 : space_cprj = xg_nonlop%space_cprj
1281 :
1282 18244 : xg_nonlop%npw_k = npw_k
1283 :
1284 73416 : xg_nonlop%l_npw_k(:) = 0
1285 18244 : xg_nonlop%l_npw_k(xg_nonlop%me_band+1) = npw_k
1286 18244 : call xmpi_sum(xg_nonlop%l_npw_k,xg_nonlop%comm_band,ierr)
1287 :
1288 73416 : xg_nonlop%total_npw_k = sum(xg_nonlop%l_npw_k)
1289 73416 : xg_nonlop%max_npw_k = maxval(xg_nonlop%l_npw_k)
1290 :
1291 18244 : xg_nonlop%l_shift_npw_k(1) = 0
1292 18244 : nmpi = xmpi_comm_size(xg_nonlop%comm_band)
1293 55172 : do iblock=2,nmpi
1294 55172 : xg_nonlop%l_shift_npw_k(iblock) = xg_nonlop%l_shift_npw_k(iblock-1) + xg_nonlop%l_npw_k(iblock-1)
1295 : end do
1296 :
1297 18244 : ntypat = xg_nonlop%ntypat
1298 18244 : nlmn_max = xg_nonlop%nlmn_max
1299 :
1300 18244 : xg_nonlop%ph3d_k => ph3d_k
1301 18244 : xg_nonlop%ffnl_k => ffnl_k
1302 :
1303 18244 : xg_nonlop%kpg_k => kpg_k
1304 :
1305 18244 : xg_nonlop%projectors_k => xg_nonlop%projectors(ikpt)
1306 18244 : if (xg_nonlop%option==0) then
1307 17220 : xg_nonlop%ffnl_gather_k => xg_nonlop%ffnl_gather(ikpt)
1308 17220 : xg_nonlop%ph3d_gather_k => xg_nonlop%ph3d_gather(ikpt)
1309 : end if
1310 :
1311 18244 : if (xg_nonlop%paw) then
1312 15172 : xg_nonlop%gram_proj_k => xg_nonlop%gram_proj(ikpt)
1313 15172 : xg_nonlop%invSij_approx_k => xg_nonlop%invSij_approx(ikpt)
1314 : end if
1315 :
1316 18244 : if (compute_proj) then
1317 :
1318 : call xg_init(xg_nonlop%projectors_k,xg_nonlop%space_pw,npw_k,xg_nonlop%cprjdim,&
1319 2320 : xg_nonlop%comm_band,me_g0=me_g0_loc)
1320 :
1321 2320 : if (xg_nonlop%option==0) then
1322 2192 : call xg_init(xg_nonlop%ffnl_gather_k,SPACE_R,xg_nonlop%total_npw_k,nlmn_max*ntypat,xmpi_comm_null)
1323 : call xg_init(xg_nonlop%ph3d_gather_k,xg_nonlop%space_pw,xg_nonlop%total_npw_k,xg_nonlop%natom,&
1324 2192 : xmpi_comm_null,me_g0=me_g0_fft_loc)
1325 : end if
1326 :
1327 2320 : call xg_nonlop_compute_projs(xg_nonlop)
1328 :
1329 2320 : compute_invS_approx_=.false.
1330 2320 : if (present(compute_invS_approx)) compute_invS_approx_ = compute_invS_approx
1331 1160 : if (compute_invS_approx_) then
1332 :
1333 770 : if (.not.xg_nonlop%paw) then
1334 0 : ABI_ERROR('Not implemented with paw=False.')
1335 : end if
1336 :
1337 : !invSij_approx_k is allocated here as space_pw depends on k-point
1338 770 : call xg_init(xg_nonlop%invSij_approx_k,space_cprj,nlmn_max,nlmn_max*ntypat,xmpi_comm_self)
1339 :
1340 770 : shift_itypat=1
1341 2334 : do itypat = 1, ntypat
1342 1564 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
1343 1564 : nattyp_i = xg_nonlop%nattyp(itypat)
1344 2334 : if (nattyp_i>0) then
1345 1540 : shift_sij = 1+(itypat-1)*nlmn_max
1346 1540 : call xgBlock_setBlock(xg_nonlop%projectors_k%self,projs,npw_k,nlmn,fcol=shift_itypat)
1347 1540 : call xg_setBlock(xg_nonlop%invSij_approx_k,invSij_approx_k_itypat,nlmn,nlmn,fcol=shift_sij)
1348 1540 : call xg_setBlock(xg_nonlop%Sijm1,Sijm1_itypat,nlmn,nlmn,fcol=shift_sij)
1349 1540 : if (space_cprj==SPACE_R) then
1350 496 : call xgBlock_copy(Sijm1_itypat,invSij_approx_k_itypat)
1351 : else
1352 1044 : call xgBlock_r2c(Sijm1_itypat,invSij_approx_k_itypat,1)
1353 : end if
1354 1540 : call xg_init(work,space_cprj,nlmn,nlmn,xmpi_comm_self)
1355 1540 : call xgBlock_gemm('t','n',1.0d0,projs,projs,0.0d0,work%self,comm=xg_nonlop%comm_band)
1356 1540 : call xgBlock_add(invSij_approx_k_itypat,work%self)
1357 1540 : call xgBlock_invert_sy(invSij_approx_k_itypat,work%self)
1358 1540 : call xg_free(work)
1359 1540 : shift_itypat = shift_itypat + nlmn*nattyp_i
1360 : end if
1361 : end do
1362 :
1363 : end if
1364 :
1365 2320 : compute_gram_=.false.
1366 2320 : if (present(compute_gram)) compute_gram_ = compute_gram
1367 1160 : if (compute_gram_) then
1368 :
1369 770 : if (.not.xg_nonlop%paw) then
1370 0 : ABI_ERROR('Not implemented with paw=False.')
1371 : end if
1372 770 : cplex=xg_nonlop%cplex
1373 770 : ncols = xg_nonlop%cprjdim
1374 770 : call xg_init(xg_nonlop%gram_proj_k,space_cprj,ncols,ncols,xmpi_comm_self)
1375 770 : projs = xg_nonlop%projectors_k%self
1376 770 : call xgBlock_gemm('t','n',1.0d0,projs,projs,0.0d0,xg_nonlop%gram_proj_k%self,comm=xg_nonlop%comm_band)
1377 770 : call xgBlock_reverseMap(xg_nonlop%gram_proj_k%self,gram_proj_k_)
1378 770 : shift=0
1379 770 : shiftc=0
1380 2334 : do itypat = 1, ntypat
1381 1564 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
1382 1564 : nattyp_i = xg_nonlop%nattyp(itypat)
1383 2334 : if (nattyp_i>0) then
1384 1540 : shift_sij = 1+(itypat-1)*nlmn_max
1385 1540 : call xg_setBlock(xg_nonlop%Sijm1,Sijm1_itypat,nlmn,nlmn,fcol=shift_sij)
1386 1540 : call xgBlock_reverseMap(Sijm1_itypat,Sijm1_)
1387 3096 : do ia = 1, nattyp_i
1388 14084 : do jlmn=1,nlmn
1389 115748 : do ilmn=1,nlmn
1390 : gram_proj_k_(shiftc+cplex*(ilmn-1)+1,shift+jlmn) = gram_proj_k_(shiftc+cplex*(ilmn-1)+1,shift+jlmn) &
1391 114192 : & + Sijm1_(ilmn,jlmn)
1392 : end do
1393 : end do
1394 1556 : shift = shift + nlmn
1395 3096 : shiftc = shiftc + cplex*nlmn
1396 : end do
1397 : end if
1398 : end do
1399 :
1400 : !ITEST
1401 : !write(903,*) 'make_cprj, compute gram', cols, space_cprj==SPACE_C
1402 : !call xgBlock_hermitian_pd_cond(xg_nonlop%gram_proj_k%self, cols, cond)
1403 : !write(903,*) 'cond(B)=', cond
1404 : !flush(903)
1405 : !ITEST
1406 :
1407 : end if
1408 :
1409 : end if ! compute_proj
1410 :
1411 18244 : call timab(tim_make_k,2,tsec)
1412 :
1413 18244 : end subroutine xg_nonlop_make_k
1414 : !!***
1415 :
1416 679140 : subroutine xg_nonlop_set_nmpi(xg_nonlop,X,cprjX,nmpi,blocksize,fft_representation)
1417 :
1418 : type(xg_nonlop_t), intent(in) :: xg_nonlop
1419 : type(xgBlock_t), intent(in) :: X
1420 : type(xgBlock_t), intent(in) :: cprjX
1421 : integer, intent(out) :: blocksize,nmpi
1422 : logical, intent(out) :: fft_representation
1423 :
1424 : integer :: nspinor,npw
1425 : type(xgBlock_t) :: projs
1426 :
1427 679140 : if (.not.associated(xg_nonlop%projectors_k)) then
1428 0 : ABI_ERROR('projectors_k should be associated')
1429 : end if
1430 :
1431 679140 : projs = xg_nonlop%projectors_k%self
1432 :
1433 679140 : npw = xg_nonlop%npw_k
1434 679140 : nspinor = xg_nonlop%nspinor
1435 679140 : blocksize = cols(cprjX)
1436 :
1437 : ! Check projector sizes
1438 679140 : if (npw/=rows(projs)) then
1439 0 : ABI_ERROR('npw/=rows(projs)')
1440 : end if
1441 679140 : if (xg_nonlop%cprjdim/=cols(projs)) then
1442 0 : ABI_ERROR('cols(projs)/=cprjdim')
1443 : end if
1444 : ! Check cprj sizes
1445 679140 : if (rows(cprjX)/=xg_nonlop%cprjdim) then
1446 0 : ABI_ERROR('rows(cprjX)/=cprjdim')
1447 : end if
1448 : ! rows(projs), cols(projs) and rows(cprj) are checked.
1449 : ! Now we check : rows(X),cols(X) and cols(cprj) depending on mpi
1450 679140 : nmpi = xmpi_comm_size(comm(cprjX))
1451 679140 : if (nmpi==1) then ! sequential
1452 196864 : if (xmpi_comm_size(comm(X))/=1) then
1453 0 : ABI_ERROR('size(comm(X))/=1')
1454 : end if
1455 196864 : if (rows(X)/=npw*nspinor) then
1456 0 : ABI_ERROR('rows(X)/=npw*nspinor')
1457 : end if
1458 196864 : if (cols(X)*nspinor/=blocksize) then
1459 0 : ABI_ERROR('cols(cprjX)/=cols(X)*nspinor')
1460 : end if
1461 : else ! MPI
1462 : ! FFT representation (X have all rows, cols are distributed)
1463 482276 : if (rows(X)==xg_nonlop%total_npw_k*nspinor) then
1464 343900 : fft_representation = .True.
1465 343900 : if (xmpi_comm_size(comm(X))/=1) then
1466 0 : ABI_ERROR('size(comm(X))/=1')
1467 : end if
1468 343900 : if (cols(X)*nspinor/=blocksize) then
1469 0 : ABI_ERROR('cols(cprjX)/=cols(X)*nspinor')
1470 : end if
1471 : ! Linalg representation (X have all cols, rows are distributed)
1472 : else
1473 138376 : fft_representation = .False.
1474 138376 : if (comm(X)/=comm(cprjX)) then
1475 0 : ABI_ERROR('comm(X)/=comm(cprjX)')
1476 : end if
1477 138376 : if (rows(X)/=npw*nspinor) then
1478 0 : ABI_ERROR('rows(X)/=npw*nspinor')
1479 : end if
1480 138376 : if (cols(X)*nspinor/=blocksize*nmpi) then
1481 0 : ABI_ERROR('cols(cprjX)*nmpi/=cols(X)*nspinor')
1482 : end if
1483 : end if
1484 : end if
1485 :
1486 679140 : end subroutine xg_nonlop_set_nmpi
1487 :
1488 : !!****f* m_xg_nonlop/xg_nonlop_getcprj
1489 : !! NAME
1490 : !!
1491 : !! FUNCTION
1492 : !!
1493 : !! INPUTS
1494 : !!
1495 728616 : subroutine xg_nonlop_getcprj(xg_nonlop,X,cprjX,work_mpi)
1496 :
1497 : type(xg_nonlop_t), intent(in) :: xg_nonlop
1498 : type(xgBlock_t) , intent(in ) :: X
1499 : type(xgBlock_t) , intent(inout) :: cprjX,work_mpi
1500 :
1501 : real(dp) :: tsec(2)
1502 : integer :: ierr,iblock,nmpi,npw,npw_max,blocksize,nspinor,shift,shift_row
1503 : integer :: source,dest,tag,request,me_band,me_g0_loc
1504 : logical :: fft_representation
1505 : type(xgBlock_t) :: X_block,X_spinor,projs,work_mpi_npw
1506 : type(xg_t) :: work_npw
1507 :
1508 364308 : call timab(tim_getcprj,1,tsec)
1509 :
1510 364308 : call xg_nonlop_set_nmpi(xg_nonlop,X,cprjX,nmpi,blocksize,fft_representation)
1511 :
1512 364308 : nspinor = xg_nonlop%nspinor
1513 364308 : projs = xg_nonlop%projectors_k%self
1514 :
1515 364308 : call xgBlock_reshape_spinor(X,X_spinor,nspinor,ROWS2COLS)
1516 :
1517 364308 : if (nmpi==1) then
1518 :
1519 106394 : call timab(tim_getcprj_gemm,1,tsec)
1520 106394 : call xgBlock_gemm('t','n',1.0d0,projs,X_spinor,0.d0,cprjX)
1521 106394 : call timab(tim_getcprj_gemm,2,tsec)
1522 :
1523 : else
1524 :
1525 257914 : me_band=xg_nonlop%me_band
1526 :
1527 257914 : if (fft_representation) then ! FFT representation (X have all rows, cols are distributed)
1528 :
1529 177488 : npw_max = xg_nonlop%max_npw_k
1530 :
1531 177488 : if (rows(work_mpi)/=npw_max) then
1532 0 : ABI_ERROR('rows(work)/=npw_max')
1533 : end if
1534 177488 : if (cols(work_mpi)/=xg_nonlop%cprjdim) then
1535 0 : ABI_ERROR('rows(work)/=cprjdim')
1536 : end if
1537 :
1538 : ! npw_max=npw_tot/Nmpi and blocksize=nband/Nmpi so size(work_npw) ~ 1/Nmpi^2
1539 177488 : call xg_init(work_npw,xg_nonlop%space_pw,npw_max,blocksize,xmpi_comm_null,me_g0=me_g0(X))
1540 :
1541 859696 : do iblock=1,nmpi
1542 :
1543 859696 : if (iblock==1) then
1544 177488 : npw = xg_nonlop%l_npw_k(me_band+1)
1545 177488 : shift_row = xg_nonlop%l_shift_npw_k(me_band+1)
1546 177488 : call xg_setBlock(work_npw,X_block,npw,blocksize)
1547 :
1548 177488 : call timab(tim_getcprj_copy,1,tsec)
1549 177488 : call xgBlock_partialcopy(X_spinor,X_block,shift_row,0,BIG2SMALL)
1550 177488 : call timab(tim_getcprj_copy,2,tsec)
1551 :
1552 177488 : call timab(tim_getcprj_gemm,1,tsec)
1553 177488 : call xgBlock_gemm('t','n',1.0d0,projs,X_block,0.d0,cprjX)
1554 177488 : call timab(tim_getcprj_gemm,2,tsec)
1555 :
1556 : else
1557 504720 : tag = iblock
1558 504720 : dest = mod(me_band-(iblock-1),nmpi)
1559 504720 : if (dest<0) dest=dest+nmpi
1560 :
1561 504720 : if (xg_nonlop%option==1) then
1562 55488 : call timab(tim_getcprj_mpi,1,tsec)
1563 55488 : call xgBlock_mpi_isend(projs,dest,tag,request,comm=xg_nonlop%comm_band)
1564 55488 : call timab(tim_getcprj_mpi,2,tsec)
1565 : end if
1566 :
1567 504720 : source = mod(me_band+(iblock-1),nmpi)
1568 504720 : npw = xg_nonlop%l_npw_k(source+1)
1569 504720 : shift_row = xg_nonlop%l_shift_npw_k(source+1)
1570 504720 : me_g0_loc = -1
1571 504720 : if (xg_nonlop%space_pw==SPACE_CR) then
1572 247608 : me_g0_loc = me_g0(X)
1573 247608 : if (shift_row>0) me_g0_loc = 0
1574 : end if
1575 504720 : call xg_setBlock(work_npw,X_block,npw,blocksize)
1576 :
1577 504720 : call timab(tim_getcprj_copy,1,tsec)
1578 504720 : call xgBlock_partialcopy(X_spinor,X_block,shift_row,0,BIG2SMALL)
1579 504720 : call timab(tim_getcprj_copy,2,tsec)
1580 :
1581 504720 : call xgBlock_setBlock(work_mpi,work_mpi_npw,npw_max,xg_nonlop%cprjdim)
1582 504720 : call xgBlock_free_reshape(work_mpi_npw,npw,xg_nonlop%cprjdim,new_me_g0=me_g0_loc)
1583 :
1584 504720 : if (xg_nonlop%option==1) then
1585 55488 : call timab(tim_getcprj_mpi,1,tsec)
1586 55488 : call xgBlock_mpi_recv(work_mpi_npw,source,tag,comm=xg_nonlop%comm_band)
1587 55488 : call timab(tim_getcprj_mpi,2,tsec)
1588 449232 : else if (xg_nonlop%option==0) then
1589 449232 : call timab(tim_getcprj_otf,1,tsec)
1590 449232 : call xg_nonlop_compute_projs_otf(xg_nonlop,work_mpi_npw,source)
1591 449232 : call timab(tim_getcprj_otf,2,tsec)
1592 : else
1593 0 : ABI_ERROR("Wrong xg_nonlop%option")
1594 : end if
1595 :
1596 504720 : call timab(tim_getcprj_gemm,1,tsec)
1597 504720 : call xgBlock_gemm('t','n',1.0d0,work_mpi_npw,X_block,1.d0,cprjX)
1598 504720 : call timab(tim_getcprj_gemm,2,tsec)
1599 :
1600 504720 : if (xg_nonlop%option==1) then
1601 55488 : call timab(tim_getcprj_mpi,1,tsec)
1602 55488 : call xmpi_wait(request,ierr)
1603 55488 : call timab(tim_getcprj_mpi,2,tsec)
1604 : end if
1605 :
1606 : end if
1607 :
1608 : end do
1609 :
1610 177488 : call xg_free(work_npw)
1611 :
1612 : else ! Linalg representation (X have all cols, rows are distributed)
1613 :
1614 80426 : call xgBlock_check(cprjX,work_mpi)
1615 :
1616 387078 : do iblock=1,nmpi
1617 306652 : shift=1+(iblock-1)*blocksize
1618 306652 : call xgBlock_setBlock(X_spinor,X_block,rows(X_spinor),blocksize,fcol=shift)
1619 :
1620 306652 : call timab(tim_getcprj_gemm,1,tsec)
1621 306652 : call xgBlock_gemm('t','n',1.0d0,projs,X_block,0.d0,work_mpi)
1622 306652 : call timab(tim_getcprj_gemm,2,tsec)
1623 : ! We do the mpi sum outside xgBlock_gemm just to include the timing in tim_getcprj_mpi,
1624 : ! (instead of tim_gemm_mpi).
1625 306652 : call timab(tim_getcprj_mpi,1,tsec)
1626 306652 : call xgBlock_mpi_sum(work_mpi,comm=xg_nonlop%comm_band)
1627 306652 : call timab(tim_getcprj_mpi,2,tsec)
1628 :
1629 387078 : if (me_band==iblock-1) then
1630 80426 : call timab(tim_getcprj_copy,1,tsec)
1631 80426 : call xgBlock_copy(work_mpi,cprjX)
1632 80426 : call timab(tim_getcprj_copy,2,tsec)
1633 : end if
1634 : end do
1635 :
1636 : end if
1637 :
1638 : end if
1639 :
1640 364308 : call timab(tim_getcprj,2,tsec)
1641 :
1642 364308 : end subroutine xg_nonlop_getcprj
1643 : !!***
1644 :
1645 : !!****f* m_xg_nonlop/xg_nonlop_getcprj_deriv
1646 : !! NAME
1647 : !!
1648 : !! FUNCTION
1649 : !!
1650 : !! INPUTS
1651 : !!
1652 19440 : subroutine xg_nonlop_getcprj_deriv(xg_nonlop,X,cprjX,work_mpi,option)
1653 :
1654 : type(xg_nonlop_t), intent(in) :: xg_nonlop
1655 : type(xgBlock_t) , intent(in) :: X
1656 : integer , intent(in) :: option
1657 : type(xgBlock_t) , intent(inout) :: cprjX,work_mpi
1658 :
1659 : ! real(dp) :: tsec(2)
1660 : integer :: iblock,nmpi,npw,blocksize,nspinor,shift
1661 : integer :: me_band,proj_size
1662 : type(xgBlock_t) :: X_block,X_spinor
1663 : type(xg_t) :: projs_deriv
1664 :
1665 : ! call timab(tim_getcprj,1,tsec)
1666 :
1667 19440 : npw = xg_nonlop%npw_k
1668 19440 : nspinor = xg_nonlop%nspinor
1669 19440 : blocksize = cols(cprjX)
1670 :
1671 19440 : if (option==DERIV_ATOM) then
1672 18300 : proj_size = 3*xg_nonlop%cprjdim
1673 1140 : else if (option==DERIV_STRESS) then
1674 1140 : proj_size = 6*xg_nonlop%cprjdim
1675 : else
1676 0 : ABI_ERROR('Bad option')
1677 : end if
1678 : ! Check cprj sizes
1679 19440 : if (rows(cprjX)/=proj_size) then
1680 0 : ABI_ERROR('rows(cprjX)/=proj_size')
1681 : end if
1682 : ! Now we check : rows(X),cols(X) and cols(cprj) depending on mpi
1683 19440 : nmpi = xmpi_comm_size(comm(cprjX))
1684 19440 : if (nmpi==1) then ! sequential
1685 6496 : if (xmpi_comm_size(comm(X))/=1) then
1686 0 : ABI_ERROR('size(comm(X))/=1')
1687 : end if
1688 6496 : if (rows(X)/=npw*nspinor) then
1689 0 : ABI_ERROR('rows(X)/=npw*nspinor')
1690 : end if
1691 6496 : if (cols(X)*nspinor/=blocksize) then
1692 0 : ABI_ERROR('cols(cprjX)/=cols(X)*nspinor')
1693 : end if
1694 : else ! MPI
1695 : ! Linalg representation (X have all cols, rows are distributed)
1696 12944 : if (comm(X)/=comm(cprjX)) then
1697 0 : ABI_ERROR('comm(X)/=comm(cprjX)')
1698 : end if
1699 12944 : if (rows(X)/=npw*nspinor) then
1700 0 : ABI_ERROR('rows(X)/=npw*nspinor')
1701 : end if
1702 12944 : if (cols(X)*nspinor/=blocksize*nmpi) then
1703 0 : ABI_ERROR('cols(cprjX)*nmpi/=cols(X)*nspinor')
1704 : end if
1705 : end if
1706 :
1707 19440 : call xg_init(projs_deriv,xg_nonlop%space_pw,npw,proj_size,comm=xg_nonlop%comm_band,me_g0=me_g0(X))
1708 :
1709 18300 : select case (option)
1710 : case (DERIV_ATOM)
1711 18300 : call xg_nonlop_compute_projs_deriv_atom(xg_nonlop,projs_deriv%self)
1712 : case (DERIV_STRESS)
1713 19440 : call xg_nonlop_compute_projs_deriv_stress(xg_nonlop,projs_deriv%self)
1714 : end select
1715 :
1716 19440 : call xgBlock_reshape_spinor(X,X_spinor,nspinor,ROWS2COLS)
1717 :
1718 19440 : if (nmpi==1) then
1719 :
1720 : ! call timab(tim_getcprj_gemm,1,tsec)
1721 6496 : call xgBlock_gemm('t','n',1.0d0,projs_deriv%self,X_spinor,0.d0,cprjX)
1722 : ! call timab(tim_getcprj_gemm,2,tsec)
1723 :
1724 : else
1725 :
1726 12944 : me_band=xg_nonlop%me_band
1727 :
1728 12944 : call xgBlock_check(cprjX,work_mpi)
1729 :
1730 62128 : do iblock=1,nmpi
1731 49184 : shift=1+(iblock-1)*blocksize
1732 49184 : call xgBlock_setBlock(X_spinor,X_block,rows(X_spinor),blocksize,fcol=shift)
1733 :
1734 : ! call timab(tim_getcprj_gemm,1,tsec)
1735 49184 : call xgBlock_gemm('t','n',1.0d0,projs_deriv%self,X_block,0.d0,work_mpi)
1736 : ! call timab(tim_getcprj_gemm,2,tsec)
1737 : ! We do the mpi sum outside xgBlock_gemm just to include the timing in tim_getcprj_mpi,
1738 : ! (instead of tim_gemm_mpi).
1739 : ! call timab(tim_getcprj_mpi,1,tsec)
1740 49184 : call xgBlock_mpi_sum(work_mpi,comm=xg_nonlop%comm_band)
1741 : ! call timab(tim_getcprj_mpi,2,tsec)
1742 :
1743 62128 : if (me_band==iblock-1) then
1744 : ! call timab(tim_getcprj_copy,1,tsec)
1745 12944 : call xgBlock_copy(work_mpi,cprjX)
1746 : ! call timab(tim_getcprj_copy,2,tsec)
1747 : end if
1748 : end do
1749 :
1750 : end if
1751 :
1752 19440 : call xg_free(projs_deriv)
1753 : ! call timab(tim_getcprj,2,tsec)
1754 :
1755 19440 : end subroutine xg_nonlop_getcprj_deriv
1756 : !!***
1757 :
1758 629664 : subroutine xg_nonlop_apply_prj(xg_nonlop,cprjX,X,work_mpi)
1759 :
1760 : type(xg_nonlop_t), intent(in) :: xg_nonlop
1761 : type(xgBlock_t), intent(in) :: cprjX
1762 : type(xgBlock_t), intent(inout) :: X,work_mpi
1763 :
1764 : logical :: fft_representation
1765 : integer :: nmpi,npw,blocksize,iblock,shift_col,nspinor!,nband,npw
1766 : integer :: me_band,shift_row,npw_max
1767 : integer :: ierr,source,dest,tag,request
1768 : real(dp) :: tsec(2)
1769 : type(xgBlock_t) :: projs
1770 : type(xgBlock_t) :: X_spinor,X_block,work_mpi_npw
1771 : type(xg_t) :: work_npw
1772 :
1773 314832 : call timab(tim_apply_prj,1,tsec)
1774 :
1775 314832 : call xg_nonlop_set_nmpi(xg_nonlop,X,cprjX,nmpi,blocksize,fft_representation)
1776 :
1777 314832 : projs = xg_nonlop%projectors_k%self
1778 :
1779 314832 : nspinor = xg_nonlop%nspinor
1780 314832 : call xgBlock_reshape_spinor(X,X_spinor,nspinor,ROWS2COLS)
1781 :
1782 314832 : if (nmpi==1) then
1783 :
1784 90470 : call timab(tim_apply_prj_gemm,1,tsec)
1785 90470 : call xgBlock_gemm('n','n',1.0d0,projs,cprjX,1.d0,X_spinor)
1786 90470 : call timab(tim_apply_prj_gemm,2,tsec)
1787 :
1788 : else
1789 :
1790 224362 : me_band = xg_nonlop%me_band
1791 :
1792 224362 : if (fft_representation) then ! FFT representation (X have all rows, cols are distributed)
1793 :
1794 166412 : npw_max = xg_nonlop%max_npw_k
1795 :
1796 166412 : if (rows(work_mpi)/=npw_max) then
1797 0 : ABI_ERROR('rows(work)/=npw_max')
1798 : end if
1799 166412 : if (cols(work_mpi)/=xg_nonlop%cprjdim) then
1800 0 : ABI_ERROR('rows(work)/=cprjdim')
1801 : end if
1802 :
1803 : ! npw_max=npw_tot/Nmpi and blocksize=nband/Nmpi so size(work_npw) ~ 1/Nmpi^2
1804 166412 : call xg_init(work_npw,xg_nonlop%space_pw,npw_max,blocksize,xmpi_comm_null)
1805 :
1806 804316 : do iblock=1,nmpi
1807 :
1808 804316 : if (iblock==1) then
1809 166412 : npw = xg_nonlop%l_npw_k(me_band+1)
1810 166412 : call xg_setBlock(work_npw,X_block,npw,blocksize)
1811 166412 : shift_row = xg_nonlop%l_shift_npw_k(me_band+1)
1812 :
1813 166412 : call timab(tim_apply_prj_copy,1,tsec)
1814 166412 : call xgBlock_partialcopy(X_spinor,X_block,shift_row,0,BIG2SMALL)
1815 166412 : call timab(tim_apply_prj_copy,2,tsec)
1816 :
1817 166412 : call timab(tim_apply_prj_gemm,1,tsec)
1818 166412 : call xgBlock_gemm('n','n',1.0d0,projs,cprjX,1.d0,X_block)
1819 166412 : call timab(tim_apply_prj_gemm,2,tsec)
1820 :
1821 166412 : call timab(tim_apply_prj_copy,1,tsec)
1822 166412 : call xgBlock_partialcopy(X_block,X_spinor,shift_row,0,SMALL2BIG)
1823 166412 : call timab(tim_apply_prj_copy,2,tsec)
1824 : else
1825 471492 : tag = iblock
1826 471492 : dest = mod(me_band-(iblock-1),nmpi)
1827 471492 : if (dest<0) dest=dest+nmpi
1828 :
1829 471492 : if (xg_nonlop%option==1) then
1830 55488 : call timab(tim_apply_prj_mpi,1,tsec)
1831 55488 : call xgBlock_mpi_isend(projs,dest,tag,request,comm=xg_nonlop%comm_band)
1832 55488 : call timab(tim_apply_prj_mpi,2,tsec)
1833 : end if
1834 :
1835 471492 : source = mod(me_band+(iblock-1),nmpi)
1836 471492 : npw = xg_nonlop%l_npw_k(source+1)
1837 471492 : call xg_setBlock(work_npw,X_block,npw,blocksize)
1838 471492 : shift_row = xg_nonlop%l_shift_npw_k(source+1)
1839 :
1840 471492 : call timab(tim_apply_prj_copy,1,tsec)
1841 471492 : call xgBlock_partialcopy(X_spinor,X_block,shift_row,0,BIG2SMALL)
1842 471492 : call timab(tim_apply_prj_copy,2,tsec)
1843 :
1844 471492 : call xgBlock_setBlock(work_mpi,work_mpi_npw,npw_max,xg_nonlop%cprjdim)
1845 471492 : call xgBlock_free_reshape(work_mpi_npw,npw,xg_nonlop%cprjdim)
1846 :
1847 471492 : if (xg_nonlop%option==1) then
1848 55488 : call timab(tim_apply_prj_mpi,1,tsec)
1849 55488 : call xgBlock_mpi_recv(work_mpi_npw,source,tag,comm=xg_nonlop%comm_band)
1850 55488 : call timab(tim_apply_prj_mpi,2,tsec)
1851 416004 : else if (xg_nonlop%option==0) then
1852 416004 : call timab(tim_apply_prj_otf,1,tsec)
1853 416004 : call xg_nonlop_compute_projs_otf(xg_nonlop,work_mpi_npw,source)
1854 416004 : call timab(tim_apply_prj_otf,2,tsec)
1855 : else
1856 0 : ABI_ERROR("Wrong xg_nonlop%option")
1857 : end if
1858 :
1859 471492 : call timab(tim_apply_prj_gemm,1,tsec)
1860 471492 : call xgBlock_gemm('n','n',1.0d0,work_mpi_npw,cprjX,1.d0,X_block)
1861 471492 : call timab(tim_apply_prj_gemm,2,tsec)
1862 :
1863 471492 : call timab(tim_apply_prj_copy,1,tsec)
1864 471492 : call xgBlock_partialcopy(X_block,X_spinor,shift_row,0,SMALL2BIG)
1865 471492 : call timab(tim_apply_prj_copy,2,tsec)
1866 :
1867 471492 : if (xg_nonlop%option==1) then
1868 55488 : call timab(tim_apply_prj_mpi,1,tsec)
1869 55488 : call xmpi_wait(request,ierr)
1870 55488 : call timab(tim_apply_prj_mpi,2,tsec)
1871 : end if
1872 : end if
1873 :
1874 : end do
1875 :
1876 166412 : call xg_free(work_npw)
1877 :
1878 : else ! Linalg representation (X have all cols, rows are distributed)
1879 :
1880 57950 : call xgBlock_check(cprjX,work_mpi)
1881 :
1882 278762 : do iblock=1,nmpi
1883 :
1884 220812 : shift_col = 1 + mod(me_band+iblock-1,nmpi) * cols(cprjX)
1885 :
1886 220812 : call xgBlock_setBlock(X_spinor,X_block,rows(X_spinor),cols(cprjX),fcol=shift_col)
1887 :
1888 278762 : if (iblock==1) then
1889 :
1890 57950 : call timab(tim_apply_prj_gemm,1,tsec)
1891 57950 : call xgBlock_gemm('n','n',1.0d0,projs,cprjX,1.d0,X_block)
1892 57950 : call timab(tim_apply_prj_gemm,2,tsec)
1893 :
1894 : else
1895 :
1896 162862 : call timab(tim_apply_prj_mpi,1,tsec)
1897 162862 : tag = iblock
1898 162862 : dest = mod(me_band-(iblock-1),nmpi)
1899 162862 : if (dest<0) dest=dest+nmpi
1900 162862 : call xgBlock_mpi_isend(cprjX,dest,tag,request)
1901 162862 : source = mod(me_band+(iblock-1),nmpi)
1902 162862 : call xgBlock_mpi_recv(work_mpi,source,tag)
1903 162862 : call timab(tim_apply_prj_mpi,2,tsec)
1904 :
1905 162862 : call timab(tim_apply_prj_gemm,1,tsec)
1906 162862 : call xgBlock_gemm('n','n',1.0d0,projs,work_mpi,1.d0,X_block)
1907 162862 : call timab(tim_apply_prj_gemm,2,tsec)
1908 :
1909 162862 : call timab(tim_apply_prj_mpi,1,tsec)
1910 162862 : call xmpi_wait(request,ierr)
1911 162862 : call timab(tim_apply_prj_mpi,2,tsec)
1912 :
1913 : end if
1914 :
1915 : end do
1916 :
1917 : end if
1918 :
1919 : end if
1920 :
1921 314832 : call timab(tim_apply_prj,2,tsec)
1922 :
1923 314832 : end subroutine xg_nonlop_apply_prj
1924 : !!***
1925 :
1926 142520 : subroutine xg_nonlop_apply_diag(xg_nonlop,diag_op,cprjin,cprjout)
1927 :
1928 : type(xg_nonlop_t), intent(in) :: xg_nonlop
1929 : type(xgBlock_t), intent(in) :: diag_op
1930 : type(xgBlock_t), intent(in) :: cprjin
1931 : type(xgBlock_t), intent(inout) :: cprjout
1932 :
1933 : logical :: loop_over_atoms
1934 : integer :: ia, iband, cprjdim, shift_itypat, iatom, itypat, nattyp_i, nlmn, shift
1935 : integer :: space_cprj, cplex, nlmn_max
1936 : integer :: nspinor, nrows, ncols
1937 :
1938 : type(xg_t) :: cprjin_nlmn_max,cprjout_nlmn_max
1939 : type(xgBlock_t) :: cprjin_nlmn,cprjout_nlmn,diag_op_iatom
1940 71260 : real(dp),pointer :: cprjin_nlmn_(:,:),cprjout_nlmn_(:,:)
1941 71260 : real(dp),pointer :: cprjin_(:,:),cprjout_(:,:)
1942 : real(dp) :: tsec(2)
1943 :
1944 71260 : call timab(tim_apply_diag,1,tsec)
1945 :
1946 71260 : call xgBlock_getsize(cprjin,nrows,ncols)
1947 :
1948 71260 : cprjdim = xg_nonlop%cprjdim
1949 71260 : if (nrows/=cprjdim) then
1950 0 : ABI_ERROR('nrows/=cprjdim')
1951 : end if
1952 71260 : call xgBlock_check(cprjin,cprjout)
1953 :
1954 71260 : space_cprj = space(cprjin)
1955 71260 : if (space_cprj==SPACE_C) then
1956 : cplex=2
1957 : else
1958 59308 : cplex=1
1959 : end if
1960 :
1961 71260 : nspinor=xg_nonlop%nspinor
1962 71260 : nlmn_max=xg_nonlop%nlmn_max
1963 :
1964 : ! Create work spaces for cprj of ONE atom for ALL bands
1965 71260 : call xg_init(cprjin_nlmn_max ,space_cprj,nlmn_max,ncols)
1966 71260 : call xg_init(cprjout_nlmn_max,space_cprj,nlmn_max,ncols)
1967 :
1968 71260 : call xgBlock_reverseMap(cprjin ,cprjin_ )
1969 71260 : call xgBlock_reverseMap(cprjout,cprjout_)
1970 :
1971 71260 : if (cols(diag_op) == xg_nonlop%ntypat) then
1972 : loop_over_atoms = .false.
1973 0 : else if (cols(diag_op) == xg_nonlop%natom) then
1974 : loop_over_atoms = .true.
1975 : else
1976 0 : ABI_ERROR('wrong cols for diag_op!')
1977 : end if
1978 71260 : shift = 0
1979 71260 : shift_itypat = 0
1980 213780 : do itypat=1,xg_nonlop%ntypat
1981 142520 : nlmn=xg_nonlop%nlmn_ntypat(itypat)
1982 142520 : nattyp_i=xg_nonlop%nattyp(itypat)
1983 213780 : if (nattyp_i>0) then
1984 142520 : call xg_setBlock(cprjin_nlmn_max ,cprjin_nlmn ,nlmn,ncols)
1985 142520 : call xg_setBlock(cprjout_nlmn_max,cprjout_nlmn,nlmn,ncols)
1986 142520 : call xgBlock_reverseMap(cprjin_nlmn ,cprjin_nlmn_ )
1987 142520 : call xgBlock_reverseMap(cprjout_nlmn,cprjout_nlmn_)
1988 142520 : if (.not.loop_over_atoms) then
1989 142520 : call xgBlock_setBlock(diag_op,diag_op_iatom,nlmn,1,fcol=itypat)
1990 : end if
1991 285040 : do ia=1,nattyp_i
1992 142520 : if (loop_over_atoms) then
1993 0 : iatom = ia + shift_itypat
1994 0 : call xgBlock_setBlock(diag_op,diag_op_iatom,nlmn,1,fcol=iatom)
1995 : end if
1996 : ! Copy cprj of ONE atom for ALL bands from cprjin to cprin_nlmn
1997 546560 : do iband=1,ncols
1998 32411960 : cprjin_nlmn_(1:cplex*nlmn,iband) = cprjin_(1+shift:cplex*nlmn+shift,iband)
1999 : end do
2000 142520 : call xgBlock_apply_diag(cprjin_nlmn,diag_op_iatom,1,Y=cprjout_nlmn)
2001 546560 : do iband=1,ncols
2002 : cprjout_(1+shift:cplex*nlmn+shift,iband) = cprjout_(1+shift:cplex*nlmn+shift,iband) &
2003 32411960 : & + cprjout_nlmn_(1:cplex*nlmn,iband)
2004 : end do
2005 285040 : shift=shift+cplex*nlmn
2006 : end do
2007 142520 : shift_itypat = shift_itypat + nattyp_i
2008 : end if
2009 :
2010 : end do
2011 :
2012 71260 : call xg_free(cprjin_nlmn_max)
2013 71260 : call xg_free(cprjout_nlmn_max)
2014 :
2015 71260 : call timab(tim_apply_diag,2,tsec)
2016 :
2017 71260 : end subroutine xg_nonlop_apply_diag
2018 : !!***
2019 :
2020 2827252 : subroutine xg_nonlop_apply_Aij(xg_nonlop,Aij,cprjin,cprjout,A_with_spin)
2021 :
2022 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2023 : type(xgBlock_t), intent(in) :: Aij
2024 : type(xgBlock_t), intent(in) :: cprjin
2025 : type(xgBlock_t), intent(inout) :: cprjout
2026 : logical,optional,intent(in) :: A_with_spin
2027 :
2028 : logical :: loop_over_atoms
2029 : integer :: ia, iband, cprjdim, shift_itypat, iatom, itypat, nattyp_i, nlmn, shift
2030 : integer :: space_aij, space_cprj, cplex, nlmn_max
2031 : integer :: nspinor, nrows, ncols, nrows_A, ncols_A
2032 : integer :: nlmn_1atom,nlmn_max_1atom,ncols_1atom
2033 :
2034 : type(xg_t) :: cprjin_nlmn_max,cprjout_nlmn_max
2035 : type(xg_t),target :: Aij_complex
2036 : type(xgBlock_t) :: cprjin_nlmn,cprjout_nlmn
2037 : type(xgBlock_t) :: Aij_iatom,Aij_iatom_
2038 1413626 : real(dp),pointer :: cprjin_nlmn_(:,:),cprjout_nlmn_(:,:)
2039 1413626 : real(dp),pointer :: cprjin_(:,:),cprjout_(:,:)
2040 : real(dp) :: tsec(2)
2041 : logical :: aij_r2c,A_with_spin_
2042 :
2043 1413626 : call timab(tim_apply_Aij,1,tsec)
2044 :
2045 1413626 : call xgBlock_getsize(cprjin,nrows,ncols)
2046 :
2047 1413626 : cprjdim = xg_nonlop%cprjdim
2048 1413626 : if (nrows/=cprjdim) then
2049 0 : ABI_ERROR('nrows/=cprjdim')
2050 : end if
2051 1413626 : call xgBlock_check(cprjin,cprjout)
2052 :
2053 1413626 : space_cprj = space(cprjin)
2054 1413626 : if (space_cprj==SPACE_C) then
2055 : cplex=2
2056 : else
2057 522161 : cplex=1
2058 : end if
2059 1413626 : space_aij = space(Aij)
2060 1413626 : aij_r2c = .false.
2061 1413626 : if (space_cprj/=space_aij) then
2062 345425 : if (space_aij==SPACE_R.and.space_cprj==SPACE_C) then
2063 : aij_r2c = .true.
2064 : else
2065 0 : ABI_ERROR('space_aij and space_cprj are not compatible')
2066 : end if
2067 : end if
2068 :
2069 1413626 : nspinor=xg_nonlop%nspinor
2070 1413626 : nlmn_max=xg_nonlop%nlmn_max
2071 :
2072 1413626 : A_with_spin_=.true.
2073 1413626 : if (present(A_with_spin)) then
2074 708025 : A_with_spin_=A_with_spin
2075 : end if
2076 708025 : if (A_with_spin_) then
2077 705601 : nlmn_max_1atom = nlmn_max * nspinor
2078 705601 : ncols_1atom = ncols / nspinor
2079 : else
2080 708025 : nlmn_max_1atom = nlmn_max
2081 708025 : ncols_1atom = ncols
2082 : end if
2083 :
2084 1413626 : if (aij_r2c) then
2085 345425 : if (rows(Aij) /= nlmn_max) then
2086 0 : ABI_ERROR('wrong rows for Aij')
2087 : end if
2088 345425 : if (cols(Aij) == xg_nonlop%ntypat*nlmn_max) then
2089 : loop_over_atoms = .false.
2090 18421 : else if (cols(Aij) == xg_nonlop%natom*nlmn_max) then
2091 : loop_over_atoms = .true.
2092 : else
2093 0 : ABI_ERROR('wrong cols for Aij (aij_r2c)')
2094 : end if
2095 : else
2096 1068201 : if (rows(Aij) /= nlmn_max_1atom) then
2097 0 : ABI_ERROR('wrong rows for Aij')
2098 : end if
2099 1068201 : if (cols(Aij) == xg_nonlop%ntypat*nlmn_max_1atom) then
2100 : loop_over_atoms = .false.
2101 18070 : else if (cols(Aij) == xg_nonlop%natom*nlmn_max_1atom) then
2102 : loop_over_atoms = .true.
2103 : else
2104 0 : ABI_ERROR('wrong cols for Aij')
2105 : end if
2106 : end if
2107 :
2108 : ! Create work spaces for cprj of ONE atom for ALL bands
2109 1413626 : call xg_init(cprjin_nlmn_max ,space_cprj,nlmn_max_1atom,ncols_1atom)
2110 1413626 : call xg_init(cprjout_nlmn_max,space_cprj,nlmn_max_1atom,ncols_1atom)
2111 :
2112 1413626 : call xgBlock_reverseMap(cprjin ,cprjin_ )
2113 1413626 : call xgBlock_reverseMap(cprjout,cprjout_)
2114 :
2115 1413626 : shift = 0
2116 1413626 : shift_itypat = 0
2117 4231761 : do itypat=1,xg_nonlop%ntypat
2118 2818135 : nlmn=xg_nonlop%nlmn_ntypat(itypat)
2119 2818135 : if (A_with_spin_) then
2120 1377521 : nlmn_1atom = nlmn*nspinor
2121 : else
2122 1440614 : nlmn_1atom = nlmn
2123 : end if
2124 2818135 : nattyp_i=xg_nonlop%nattyp(itypat)
2125 4231761 : if (nattyp_i>0) then
2126 2768036 : call xg_setBlock(cprjin_nlmn_max ,cprjin_nlmn ,nlmn_1atom,ncols_1atom)
2127 2768036 : call xg_setBlock(cprjout_nlmn_max,cprjout_nlmn,nlmn_1atom,ncols_1atom)
2128 2768036 : call xgBlock_reverseMap(cprjin_nlmn ,cprjin_nlmn_ )
2129 2768036 : call xgBlock_reverseMap(cprjout_nlmn,cprjout_nlmn_)
2130 2768036 : if (aij_r2c) call xg_init(Aij_complex,SPACE_C,nlmn_1atom,nlmn_1atom)
2131 2768036 : if (.not.loop_over_atoms) then
2132 2718988 : if (aij_r2c) then
2133 627647 : call xgBlock_setBlock(Aij,Aij_iatom,nlmn,nlmn,fcol=1+(itypat-1)*nlmn_max)
2134 : else
2135 2091341 : call xgBlock_setBlock(Aij,Aij_iatom,nlmn_1atom,nlmn_1atom,fcol=1+(itypat-1)*nlmn_max_1atom)
2136 : end if
2137 : end if
2138 5609800 : do ia=1,nattyp_i
2139 2841764 : if (loop_over_atoms) then
2140 73894 : iatom = ia + shift_itypat
2141 73894 : if (aij_r2c) then
2142 37754 : call xgBlock_setBlock(Aij,Aij_iatom,nlmn,nlmn,fcol=1+(iatom-1)*nlmn_max)
2143 : else
2144 36140 : call xgBlock_setBlock(Aij,Aij_iatom,nlmn_1atom,nlmn_1atom,fcol=1+(iatom-1)*nlmn_max_1atom)
2145 : end if
2146 : end if
2147 2841764 : call xgBlock_getsize(Aij_iatom,nrows_A,ncols_A)
2148 2841764 : if (.not.aij_r2c) then
2149 2149602 : if (nrows_A/=nlmn_1atom) then
2150 0 : ABI_ERROR('nrows_A/=nlmn_1atom')
2151 : end if
2152 2149602 : if (ncols_A/=nlmn_1atom) then
2153 0 : ABI_ERROR('ncols_A/=nlmn_1atom')
2154 : end if
2155 : else
2156 692162 : if (nrows_A/=nlmn) then
2157 0 : ABI_ERROR('nrows_A/=nlmn')
2158 : end if
2159 692162 : if (ncols_A/=nlmn) then
2160 0 : ABI_ERROR('ncols_A/=nlmn')
2161 : end if
2162 : end if
2163 : ! if needed, transfer real matrix to a complex one
2164 2841764 : if (aij_r2c) then
2165 692162 : call xgBlock_r2c(Aij_iatom,Aij_complex%self,nspinor)
2166 692162 : Aij_iatom_ = Aij_complex%self
2167 : else
2168 2149602 : Aij_iatom_ = Aij_iatom
2169 : end if
2170 : ! Copy cprj of ONE atom for ALL bands from cprjin to cprin_nlmn
2171 2841764 : if (A_with_spin_) then
2172 5332134 : do iband=1,ncols_1atom
2173 121269702 : cprjin_nlmn_(1:cplex*nlmn,iband) = cprjin_(1+shift:cplex*nlmn+shift,1+nspinor*(iband-1))
2174 : end do
2175 1412514 : if (nspinor==2) then
2176 1367358 : do iband=1,ncols_1atom
2177 39094014 : cprjin_nlmn_(1+cplex*nlmn:2*cplex*nlmn,iband) = cprjin_(1+shift:cplex*nlmn+shift,nspinor*iband)
2178 : end do
2179 : end if
2180 : else
2181 12099738 : do iband=1,ncols_1atom
2182 330920090 : cprjin_nlmn_(1:cplex*nlmn,iband) = cprjin_(1+shift:cplex*nlmn+shift,iband)
2183 : end do
2184 : end if
2185 :
2186 2841764 : call xgBlock_gemm('n','n',1.0d0,Aij_iatom_,cprjin_nlmn,0.d0,cprjout_nlmn,timing=.false.)
2187 :
2188 2841764 : if (A_with_spin_) then
2189 5332134 : do iband=1,ncols_1atom
2190 : cprjout_(1+shift:cplex*nlmn+shift,1+nspinor*(iband-1)) = cprjout_(1+shift:cplex*nlmn+shift,1+nspinor*(iband-1)) &
2191 121269702 : & + cprjout_nlmn_(1:cplex*nlmn,iband)
2192 : end do
2193 1412514 : if (nspinor==2) then
2194 1367358 : do iband=1,ncols_1atom
2195 : cprjout_(1+shift:cplex*nlmn+shift,nspinor*iband) = cprjout_(1+shift:cplex*nlmn+shift,nspinor*iband) &
2196 39094014 : & + cprjout_nlmn_(1+cplex*nlmn:2*cplex*nlmn,iband)
2197 : end do
2198 : end if
2199 : else
2200 12099738 : do iband=1,ncols_1atom
2201 : cprjout_(1+shift:cplex*nlmn+shift,iband) = cprjout_(1+shift:cplex*nlmn+shift,iband) &
2202 330920090 : & + cprjout_nlmn_(1:cplex*nlmn,iband)
2203 : end do
2204 : end if
2205 8451564 : shift=shift+cplex*nlmn
2206 : end do
2207 :
2208 2768036 : if (aij_r2c) call xg_free(Aij_complex)
2209 2768036 : shift_itypat = shift_itypat + nattyp_i
2210 : end if
2211 :
2212 : end do
2213 :
2214 1413626 : call xg_free(cprjin_nlmn_max)
2215 1413626 : call xg_free(cprjout_nlmn_max)
2216 :
2217 1413626 : call timab(tim_apply_Aij,2,tsec)
2218 :
2219 1413626 : end subroutine xg_nonlop_apply_Aij
2220 : !!***
2221 :
2222 189900 : subroutine xg_nonlop_precond_iterative_refinement(xg_nonlop,A,precond,cprj_in,cprj_out,cprj_work)
2223 :
2224 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2225 : type(xgBlock_t), intent(in ) :: A,precond
2226 : type(xgBlock_t), intent(in ) :: cprj_in
2227 : type(xgBlock_t), intent(inout) :: cprj_out,cprj_work
2228 :
2229 : integer :: iter,cprjdim,ncols,additional_steps_to_take
2230 : !integer :: rows_A
2231 : real(dp), parameter :: tolerance = 1e-14 ! maximum relative error. TODO: use tolwfr ?
2232 : type(xg_t) :: err
2233 : !real(dp) :: cond
2234 : real(dp) :: norm,max_err,previous_max_err,convergence_rate,tsec(2)
2235 :
2236 94950 : call timab(tim_iter_refinement,1,tsec)
2237 :
2238 : ! Note that precond is block-diagonal whereas A is not
2239 :
2240 94950 : cprjdim = xg_nonlop%cprjdim
2241 94950 : ncols = cols(cprj_in)
2242 94950 : if (cprjdim/=rows(cprj_in)) then
2243 0 : ABI_ERROR('Wrong size for cprj_in')
2244 : end if
2245 94950 : call xgBlock_check(cprj_in,cprj_out)
2246 94950 : call xgBlock_check(cprj_in,cprj_work)
2247 94950 : if (rows(A)/=cprjdim.or.cols(A)/=cprjdim) then
2248 0 : ABI_ERROR('Wrong size for A')
2249 : end if
2250 :
2251 94950 : call xg_init(err,SPACE_R,ncols,1,xmpi_comm_self)
2252 :
2253 94950 : call xgBlock_colwiseNorm2(cprj_in,err%self,max_val=norm)
2254 :
2255 : ! Y_0 = PX (with P block diagonal and "close" to A^-1)
2256 94950 : call xgBlock_zero(cprj_out)
2257 94950 : call xg_nonlop_apply_Aij(xg_nonlop,precond,cprj_in,cprj_out,A_with_spin=.false.)
2258 :
2259 94950 : additional_steps_to_take = -1
2260 708025 : do iter=1,30
2261 : ! compute AY_i
2262 708025 : call xgBlock_gemm('n','n',1.0d0,A,cprj_out,0.0d0,cprj_work)
2263 : ! ITEST
2264 : !write(903,*) 'Apply getBm1X to spd matrix of size', rows(A), cols(A)
2265 : !rows_A = rows(A)
2266 : !call xgBlock_hermitian_pd_cond(A, rows_A, cond)
2267 : !write(903,*) 'cond(B)=', cond
2268 : !flush(903)
2269 : ! ITEST
2270 : ! RES = AY_i - X
2271 708025 : call xgBlock_saxpy(cprj_work,-1.0d0,cprj_in)
2272 708025 : call xgBlock_colwiseNorm2(cprj_work,err%self,max_val=max_err)
2273 708025 : max_err = sqrt(max_err / norm)
2274 708025 : if(max_err < tolerance .or. additional_steps_to_take == 1) then
2275 : exit
2276 : ! We might stall and never get to the specified precision because of machine errors.
2277 : ! If we got to 1e-10, extrapolate convergence rate and determine the number of additional
2278 : ! steps to take to reach precision
2279 613085 : else if(max_err < 1e-10 .and. additional_steps_to_take == -1) then
2280 94950 : convergence_rate = -LOG(1e-10) / iter
2281 94950 : additional_steps_to_take = CEILING(-LOG(tolerance/1e-10)/convergence_rate) + 1
2282 518135 : else if(additional_steps_to_take > 0) then
2283 94555 : if(previous_max_err<max_err)exit
2284 94545 : additional_steps_to_take = additional_steps_to_take - 1
2285 : end if
2286 613075 : previous_max_err=max_err
2287 : ! RES = X - AY_i
2288 613075 : call xgBlock_scale(cprj_work,-1.0d0,1)
2289 : ! Y_(i+1) = Y_i + P RES
2290 1416050 : call xg_nonlop_apply_Aij(xg_nonlop,precond,cprj_work,cprj_out,A_with_spin=.false.)
2291 : end do
2292 :
2293 94950 : call xg_free(err)
2294 :
2295 94950 : call timab(tim_iter_refinement,2,tsec)
2296 :
2297 94950 : end subroutine xg_nonlop_precond_iterative_refinement
2298 : !!***
2299 :
2300 442258 : subroutine xg_nonlop_mult_cprj(xg_nonlop,cprj_left,cprj_right,res,blocksize)
2301 :
2302 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2303 : type(xgBlock_t), intent(in) :: cprj_left,cprj_right
2304 : type(xgBlock_t), intent(inout) :: res
2305 : integer, intent(in),optional :: blocksize
2306 :
2307 : integer :: space_res
2308 : integer :: blocksize_,blocksize_spinor,iblock_mpi,nblocks_mpi,shift_row,shift_col,shift_col_mpi
2309 : integer :: iblock_left,iblock_right,nblocks_left,nblocks_right
2310 : integer :: res_nrows,res_ncols,cprjdim
2311 : integer :: nrows_r,nrows_l,ncols_r,ncols_l,nspinor
2312 : integer :: comm_cprj,tag,request,ierr,me_band,source,dest
2313 : logical :: multiblock
2314 : real(dp) :: tsec(2)
2315 : type(xg_t) :: res_work_mpi,res_block,res_block_mpi,cprj_work_mpi
2316 : type(xgBlock_t) :: cprj_left_spinor,cprj_right_spinor
2317 :
2318 442258 : call timab(tim_mult_cprj,1,tsec)
2319 :
2320 442258 : comm_cprj = comm(cprj_right)
2321 442258 : nblocks_mpi = xmpi_comm_size(comm_cprj)
2322 :
2323 442258 : cprjdim = xg_nonlop%cprjdim
2324 442258 : nspinor = xg_nonlop%nspinor
2325 :
2326 442258 : call xgBlock_reshape_spinor(cprj_right,cprj_right_spinor,nspinor,COLS2ROWS)
2327 442258 : call xgBlock_reshape_spinor(cprj_left ,cprj_left_spinor ,nspinor,COLS2ROWS)
2328 :
2329 442258 : nrows_r = rows(cprj_right_spinor)
2330 442258 : nrows_l = rows(cprj_left_spinor)
2331 :
2332 442258 : ncols_r = cols(cprj_right_spinor)
2333 442258 : ncols_l = cols(cprj_left_spinor)
2334 :
2335 442258 : me_band = xg_nonlop%me_band
2336 :
2337 442258 : if (nrows_r/=nspinor*cprjdim) then
2338 0 : ABI_ERROR("rows(cprj_right)/=nspinor*cprjdim")
2339 : end if
2340 442258 : if (nrows_l/=nspinor*cprjdim) then
2341 0 : ABI_ERROR("rows(cprj_left)/=nspinor*cprjdim")
2342 : end if
2343 442258 : if (rows(res)/=nblocks_mpi*ncols_l) then
2344 0 : ABI_ERROR("rows(res)/=nblocks_mpi*cols(cprj_left)")
2345 : end if
2346 442258 : if (cols(res)/=nblocks_mpi*ncols_r) then
2347 0 : ABI_ERROR("cols(res)/=nblocks_mpi*cols(cprj_right)")
2348 : end if
2349 :
2350 442258 : blocksize_ = ncols_r
2351 442258 : if (present(blocksize)) then
2352 442258 : if (mod(blocksize,nspinor)/=0) then
2353 0 : ABI_ERROR("wrong blocksize (nspinor)")
2354 : end if
2355 442258 : if (mod(ncols_l,blocksize / nspinor)/=0) then
2356 0 : ABI_ERROR("wrong blocksize")
2357 : end if
2358 : blocksize_ = blocksize
2359 : end if
2360 :
2361 442258 : if (nblocks_mpi==1) then
2362 :
2363 159474 : call timab(tim_mult_cprj_gemm,1,tsec)
2364 159474 : call xgBlock_gemm('t','n',1.0d0,cprj_left_spinor,cprj_right_spinor,1.d0,res)
2365 159474 : call timab(tim_mult_cprj_gemm,2,tsec)
2366 :
2367 : else
2368 :
2369 282784 : blocksize_spinor = blocksize_ / nspinor
2370 282784 : nblocks_right = ncols_r / blocksize_spinor
2371 282784 : nblocks_left = ncols_l / blocksize_spinor
2372 :
2373 282784 : multiblock = .false.
2374 282784 : if (nblocks_right>1.or.nblocks_left>1) then
2375 164454 : multiblock = .true.
2376 : end if
2377 :
2378 282784 : res_nrows = rows(res)
2379 282784 : res_ncols = cols(res)
2380 :
2381 282784 : space_res = space(cprj_right)
2382 :
2383 282784 : call xg_init(res_block_mpi,space_res,ncols_l,ncols_r,xmpi_comm_null)
2384 282784 : if (multiblock) then
2385 164454 : call xg_init(res_block,space_res,blocksize_spinor,blocksize_spinor,xmpi_comm_null)
2386 : end if
2387 282784 : call xg_init(res_work_mpi,space_res,res_nrows,res_ncols,xmpi_comm_null)
2388 282784 : call xgBlock_zero(res_work_mpi%self)
2389 :
2390 282784 : call xg_init(cprj_work_mpi,space_res,nspinor*cprjdim,ncols_r,comm_cprj)
2391 :
2392 1353528 : do iblock_mpi=1,nblocks_mpi
2393 :
2394 1070744 : if (iblock_mpi==1) then
2395 :
2396 282784 : call timab(tim_mult_cprj_gemm,1,tsec)
2397 282784 : call xgBlock_gemm('t','n',1.0d0,cprj_left_spinor,cprj_right_spinor,0.d0,res_block_mpi%self)
2398 282784 : call timab(tim_mult_cprj_gemm,2,tsec)
2399 :
2400 : else
2401 :
2402 787960 : call timab(tim_mult_cprj_mpi,1,tsec)
2403 787960 : tag = iblock_mpi
2404 787960 : dest = mod(me_band-(iblock_mpi-1),nblocks_mpi)
2405 787960 : if (dest<0) dest=dest+nblocks_mpi
2406 787960 : call xgBlock_mpi_isend(cprj_right_spinor,dest,tag,request)
2407 787960 : source = mod(me_band+(iblock_mpi-1),nblocks_mpi)
2408 787960 : call xgBlock_mpi_recv(cprj_work_mpi%self,source,tag)
2409 787960 : call timab(tim_mult_cprj_mpi,2,tsec)
2410 :
2411 787960 : call timab(tim_mult_cprj_gemm,1,tsec)
2412 787960 : call xgBlock_gemm('t','n',1.0d0,cprj_left_spinor,cprj_work_mpi%self,0.d0,res_block_mpi%self)
2413 787960 : call timab(tim_mult_cprj_gemm,2,tsec)
2414 :
2415 : end if
2416 :
2417 1070744 : if (.not.multiblock) then
2418 447764 : shift_row = me_band*ncols_l
2419 447764 : shift_col = mod(me_band+iblock_mpi-1,res_ncols)*ncols_r
2420 447764 : if (shift_col>=res_ncols) shift_col = shift_col - res_ncols
2421 :
2422 447764 : call timab(tim_mult_cprj_copy,1,tsec)
2423 447764 : call xgBlock_partialcopy(res_block_mpi%self,res_work_mpi%self,shift_row,shift_col,SMALL2BIG)
2424 447764 : call timab(tim_mult_cprj_copy,2,tsec)
2425 :
2426 : else
2427 1458812 : do iblock_right=1,nblocks_right
2428 3689348 : do iblock_left=1,nblocks_left
2429 2230536 : shift_row = (iblock_left-1)*blocksize_spinor
2430 2230536 : shift_col = (iblock_right-1)*blocksize_spinor
2431 :
2432 2230536 : call timab(tim_mult_cprj_copy,1,tsec)
2433 2230536 : call xgBlock_partialcopy(res_block_mpi%self,res_block%self,shift_row,shift_col,BIG2SMALL)
2434 2230536 : call timab(tim_mult_cprj_copy,2,tsec)
2435 :
2436 2230536 : shift_row = (me_band*blocksize_spinor) + (iblock_left-1)*res_nrows/nblocks_left
2437 2230536 : shift_col_mpi = mod(me_band+iblock_mpi-1,res_ncols/nblocks_right)*blocksize_spinor
2438 2230536 : if (shift_col_mpi>=res_ncols/nblocks_right) shift_col_mpi = shift_col_mpi - res_ncols/nblocks_right
2439 2230536 : shift_col = shift_col_mpi + (iblock_right-1)*res_ncols/nblocks_right
2440 :
2441 2230536 : call timab(tim_mult_cprj_copy,1,tsec)
2442 2230536 : call xgBlock_partialcopy(res_block%self,res_work_mpi%self,shift_row,shift_col,SMALL2BIG)
2443 3066368 : call timab(tim_mult_cprj_copy,2,tsec)
2444 :
2445 : end do
2446 : end do
2447 : end if
2448 :
2449 1353528 : if (iblock_mpi>1) then
2450 787960 : call timab(tim_mult_cprj_mpi,1,tsec)
2451 787960 : call xmpi_wait(request,ierr)
2452 787960 : call timab(tim_mult_cprj_mpi,2,tsec)
2453 : end if
2454 :
2455 : end do
2456 :
2457 282784 : call timab(tim_mult_cprj_mpi,1,tsec)
2458 282784 : call xgBlock_mpi_sum(res_work_mpi%self,comm=comm_cprj)
2459 282784 : call timab(tim_mult_cprj_mpi,2,tsec)
2460 :
2461 282784 : call xgBlock_add(res,res_work_mpi%self)
2462 :
2463 282784 : if (multiblock) then
2464 164454 : call xg_free(res_block)
2465 : end if
2466 282784 : call xg_free(res_block_mpi)
2467 282784 : call xg_free(res_work_mpi)
2468 282784 : call xg_free(cprj_work_mpi)
2469 :
2470 : end if
2471 :
2472 442258 : call timab(tim_mult_cprj,2,tsec)
2473 :
2474 442258 : end subroutine xg_nonlop_mult_cprj
2475 : !!***
2476 :
2477 12920 : subroutine xg_nonlop_colwiseXAX(xg_nonlop,Aij,cprj,cprj_work,res)
2478 :
2479 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2480 : type(xgBlock_t), intent(in) :: cprj,Aij
2481 : type(xgBlock_t), intent(inout) :: cprj_work,res
2482 :
2483 : integer :: ncols
2484 : type(xgBlock_t) :: cprj_spinor,cprj_work_spinor
2485 :
2486 12920 : call xgBlock_check(cprj,cprj_work)
2487 12920 : ncols = cols(cprj)
2488 12920 : if (ncols/=xg_nonlop%nspinor*rows(res)) then
2489 0 : ABI_ERROR('Wrong cols for cprj or res.')
2490 : end if
2491 :
2492 12920 : call xgBlock_zero(cprj_work)
2493 12920 : call xg_nonlop_apply_Aij(xg_nonlop,Aij,cprj,cprj_work)
2494 :
2495 12920 : call xgBlock_reshape_spinor(cprj ,cprj_spinor ,xg_nonlop%nspinor,COLS2ROWS)
2496 12920 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,xg_nonlop%nspinor,COLS2ROWS)
2497 12920 : call xgBlock_colwiseDotProduct(cprj_spinor,cprj_work_spinor,res,comm_loc=xmpi_comm_null)
2498 :
2499 12920 : end subroutine xg_nonlop_colwiseXAX
2500 : !!***
2501 :
2502 3264 : subroutine xg_nonlop_colwiseXDX(xg_nonlop,diag,cprj,cprj_work,res)
2503 :
2504 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2505 : type(xgBlock_t), intent(in) :: cprj,diag
2506 : type(xgBlock_t), intent(inout) :: cprj_work,res
2507 :
2508 : integer :: ncols,space_diag,space_res
2509 : type(xgBlock_t) :: cprj_spinor,cprj_work_spinor
2510 : type(xg_t) :: res_complex
2511 :
2512 3264 : call xgBlock_check(cprj,cprj_work)
2513 3264 : ncols = cols(cprj)
2514 3264 : if (ncols/=xg_nonlop%nspinor*rows(res)) then
2515 0 : ABI_ERROR('Wrong cols for cprj or res.')
2516 : end if
2517 3264 : space_diag = space(diag)
2518 3264 : space_res = space(res)
2519 3264 : if (space_diag==SPACE_C) then
2520 0 : if (space_res/=SPACE_C) then
2521 0 : ABI_ERROR('space(res) should be SPACE_C.')
2522 : end if
2523 3264 : else if (space_diag/=SPACE_R) then
2524 0 : ABI_ERROR('space(diag) should be SPACE_C or SPACE_R.')
2525 : end if
2526 :
2527 3264 : call xgBlock_zero(cprj_work)
2528 3264 : call xg_nonlop_apply_diag(xg_nonlop,diag,cprj,cprj_work)
2529 :
2530 3264 : call xgBlock_reshape_spinor(cprj ,cprj_spinor ,xg_nonlop%nspinor,COLS2ROWS)
2531 3264 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,xg_nonlop%nspinor,COLS2ROWS)
2532 :
2533 : ! If space_diag==SPACE_R, the result is actually real and can be stored in a xgBlock with SPACE_R
2534 3264 : if ( space_diag==SPACE_R .and. space_res==SPACE_R .and. space(cprj)==SPACE_C) then
2535 544 : call xg_init(res_complex,SPACE_C,rows(res),cols(res))
2536 544 : call xgBlock_colwiseDotProduct(cprj_spinor,cprj_work_spinor,res_complex%self,comm_loc=xmpi_comm_null)
2537 544 : call xgBlock_c2r(res_complex%self,res)
2538 544 : call xg_free(res_complex)
2539 : else
2540 2720 : call xgBlock_colwiseDotProduct(cprj_spinor,cprj_work_spinor,res,comm_loc=xmpi_comm_null)
2541 : end if
2542 :
2543 3264 : end subroutine xg_nonlop_colwiseXDX
2544 : !!***
2545 :
2546 3264 : subroutine xg_nonlop_colwiseXHX(xg_nonlop,cprj,cprj_work,res)
2547 :
2548 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2549 : type(xgBlock_t), intent(in) :: cprj
2550 : type(xgBlock_t), intent(inout) :: cprj_work,res
2551 :
2552 : integer :: nmpi,ncols,nres,shift
2553 : ! type(xg_t) :: res_mpi
2554 : type(xgBlock_t) :: res_mpi,res_tmp
2555 :
2556 3264 : if (cols(res)/=1) then
2557 0 : ABI_ERROR('cols(res)/=1')
2558 : end if
2559 3264 : nmpi = xmpi_comm_size(comm(cprj))
2560 3264 : ncols = cols(cprj)
2561 3264 : nres = rows(res)
2562 3264 : if (nres*xg_nonlop%nspinor/=nmpi*ncols) then
2563 0 : ABI_ERROR('rows(res)*nspinor/=nmpi*cols(cprj))')
2564 : end if
2565 :
2566 3264 : call xgBlock_zero(res)
2567 :
2568 3264 : if (nmpi==1) then
2569 1088 : res_mpi = res
2570 : else
2571 2176 : call xgBlock_setBlock(res,res_tmp,nres,1)
2572 2176 : call xgBlock_reshape(res_tmp,1,nres)
2573 2176 : shift = xg_nonlop%me_band*ncols
2574 2176 : call xgBlock_setBlock(res_tmp,res_mpi,1,ncols,fcol=1+shift)
2575 2176 : call xgBlock_reshape(res_mpi,ncols,1)
2576 : end if
2577 :
2578 3264 : if (xg_nonlop%paw) then
2579 0 : call xg_nonlop_colwiseXAX(xg_nonlop,xg_nonlop%Dij_spin,cprj,cprj_work,res_mpi)
2580 : else
2581 3264 : call xg_nonlop_colwiseXDX(xg_nonlop,xg_nonlop%ekb%self,cprj,cprj_work,res_mpi)
2582 : end if
2583 :
2584 3264 : call xgBlock_mpi_sum(res,comm=xg_nonlop%comm_band)
2585 :
2586 3264 : end subroutine xg_nonlop_colwiseXHX
2587 : !!***
2588 :
2589 410705 : subroutine xg_nonlop_getXAY(xg_nonlop,Aij,cprj_left,cprj_right,cprj_work,res,blocksize)
2590 :
2591 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2592 : type(xgBlock_t), intent(in) :: cprj_left,cprj_right,Aij
2593 : type(xgBlock_t), intent(inout) :: cprj_work,res
2594 : integer,intent(in),optional :: blocksize
2595 :
2596 : integer :: blocksize_
2597 :
2598 410705 : call xgBlock_zero(cprj_work)
2599 410705 : call xg_nonlop_apply_Aij(xg_nonlop,Aij,cprj_right,cprj_work)
2600 :
2601 410705 : blocksize_ = cols(cprj_right)
2602 410705 : if (present(blocksize)) then
2603 410705 : blocksize_ = blocksize
2604 : end if
2605 410705 : call xg_nonlop_mult_cprj(xg_nonlop,cprj_left,cprj_work,res,blocksize=blocksize_)
2606 :
2607 410705 : end subroutine xg_nonlop_getXAY
2608 : !!***
2609 :
2610 31553 : subroutine xg_nonlop_getXDY(xg_nonlop,diag,cprj_left,cprj_right,cprj_work,res,blocksize)
2611 :
2612 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2613 : type(xgBlock_t), intent(in) :: cprj_left,cprj_right,diag
2614 : type(xgBlock_t), intent(inout) :: cprj_work,res
2615 : integer,intent(in),optional :: blocksize
2616 :
2617 : integer :: blocksize_
2618 :
2619 31553 : call xgBlock_zero(cprj_work)
2620 31553 : call xg_nonlop_apply_diag(xg_nonlop,diag,cprj_right,cprj_work)
2621 :
2622 31553 : blocksize_ = cols(cprj_right)
2623 31553 : if (present(blocksize)) then
2624 31553 : blocksize_ = blocksize
2625 : end if
2626 31553 : call xg_nonlop_mult_cprj(xg_nonlop,cprj_left,cprj_work,res,blocksize=blocksize_)
2627 :
2628 31553 : end subroutine xg_nonlop_getXDY
2629 : !!***
2630 :
2631 255088 : subroutine xg_nonlop_getXSY(xg_nonlop,cprj_left,cprj_right,cprj_work,res,blocksize)
2632 :
2633 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2634 : type(xgBlock_t), intent(inout) :: cprj_left,cprj_right,cprj_work,res
2635 : integer,intent(in),optional :: blocksize
2636 :
2637 : integer :: blocksize_
2638 :
2639 : real(dp) :: tsec(2)
2640 :
2641 255088 : call timab(tim_getXSY,1,tsec)
2642 :
2643 255088 : if (.not.xg_nonlop%paw) then
2644 0 : ABI_ERROR('Not implemented with paw=False.')
2645 : end if
2646 :
2647 255088 : blocksize_ = cols(cprj_right)
2648 255088 : if (present(blocksize)) then
2649 217434 : blocksize_ = blocksize
2650 : end if
2651 255088 : call xg_nonlop_getXAY(xg_nonlop,xg_nonlop%Sij%self,cprj_left,cprj_right,cprj_work,res,blocksize=blocksize_)
2652 :
2653 255088 : call timab(tim_getXSY,2,tsec)
2654 :
2655 255088 : end subroutine xg_nonlop_getXSY
2656 : !!***
2657 :
2658 187170 : subroutine xg_nonlop_getXHY(xg_nonlop,cprj_left,cprj_right,cprj_work,res,blocksize)
2659 :
2660 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2661 : type(xgBlock_t), intent(inout) :: cprj_left,cprj_right,cprj_work,res
2662 : integer,intent(in),optional :: blocksize
2663 :
2664 : integer :: blocksize_
2665 :
2666 : real(dp) :: tsec(2)
2667 :
2668 187170 : call timab(tim_getXHY,1,tsec)
2669 :
2670 187170 : blocksize_ = cols(cprj_right)
2671 187170 : if (present(blocksize)) then
2672 187170 : blocksize_ = blocksize
2673 : end if
2674 187170 : if (xg_nonlop%paw) then
2675 155617 : call xg_nonlop_getXAY(xg_nonlop,xg_nonlop%Dij_spin,cprj_left,cprj_right,cprj_work,res,blocksize=blocksize_)
2676 : else
2677 31553 : call xg_nonlop_getXDY(xg_nonlop,xg_nonlop%ekb%self,cprj_left,cprj_right,cprj_work,res,blocksize=blocksize_)
2678 : end if
2679 :
2680 187170 : call timab(tim_getXHY,2,tsec)
2681 :
2682 187170 : end subroutine xg_nonlop_getXHY
2683 : !!***
2684 :
2685 330546 : subroutine xg_nonlop_getAX(xg_nonlop,Aij,Xin,cprjin,cprj_work,work_mpi,Xout)
2686 :
2687 : use iso_c_binding
2688 :
2689 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2690 : type(xgBlock_t), intent(in) :: cprjin,Aij
2691 : type(xgBlock_t), intent(inout) :: Xin,cprj_work,work_mpi
2692 : type(xgBlock_t), optional, intent(inout) :: Xout
2693 :
2694 : integer :: nblocks,nspinor
2695 : integer :: nrows,nrows_cprj
2696 : integer :: ncols,ncols_cprj
2697 :
2698 110182 : nblocks = xmpi_comm_size(comm(cprjin))
2699 :
2700 110182 : nspinor = xg_nonlop%nspinor
2701 :
2702 : ! check sizes
2703 110182 : if (present(Xout)) then
2704 0 : call xgBlock_check(Xin,Xout)
2705 : end if
2706 110182 : call xgBlock_check(cprjin,cprj_work)
2707 110182 : call xgBlock_getsize(Xin,nrows,ncols)
2708 110182 : call xgBlock_getsize(cprjin,nrows_cprj,ncols_cprj)
2709 :
2710 : ! cprj_work = sum_j Saij cprjin
2711 110182 : call xgBlock_zero(cprj_work)
2712 110182 : call xg_nonlop_apply_Aij(xg_nonlop,Aij,cprjin,cprj_work)
2713 :
2714 : ! Xout = Xin + sum_ai pai cprj_work
2715 110182 : if (present(Xout)) then
2716 0 : call xgBlock_copy(Xin,Xout)
2717 0 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work,Xout,work_mpi)
2718 : else ! in place version
2719 110182 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work,Xin,work_mpi)
2720 : end if
2721 :
2722 110182 : end subroutine xg_nonlop_getAX
2723 :
2724 103665 : subroutine xg_nonlop_getDX(xg_nonlop,diag,Xin,cprjin,cprj_work,work_mpi,Xout)
2725 :
2726 : use iso_c_binding
2727 :
2728 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2729 : type(xgBlock_t), intent(in) :: cprjin,diag
2730 : type(xgBlock_t), intent(inout) :: Xin,cprj_work,work_mpi
2731 : type(xgBlock_t), optional, intent(inout) :: Xout
2732 :
2733 : integer :: nblocks,nspinor
2734 : integer :: nrows,nrows_cprj
2735 : integer :: ncols,ncols_cprj
2736 :
2737 34555 : nblocks = xmpi_comm_size(comm(cprjin))
2738 :
2739 34555 : nspinor = xg_nonlop%nspinor
2740 :
2741 : ! check sizes
2742 34555 : if (present(Xout)) then
2743 0 : call xgBlock_check(Xin,Xout)
2744 : end if
2745 34555 : call xgBlock_check(cprjin,cprj_work)
2746 34555 : call xgBlock_getsize(Xin,nrows,ncols)
2747 34555 : call xgBlock_getsize(cprjin,nrows_cprj,ncols_cprj)
2748 :
2749 : ! cprj_work = sum_j Saij cprjin
2750 34555 : call xgBlock_zero(cprj_work)
2751 34555 : call xg_nonlop_apply_diag(xg_nonlop,diag,cprjin,cprj_work)
2752 :
2753 : ! Xout = Xin + sum_ai pai cprj_work
2754 34555 : if (present(Xout)) then
2755 0 : call xgBlock_copy(Xin,Xout)
2756 0 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work,Xout,work_mpi)
2757 : else ! in place version
2758 34555 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work,Xin,work_mpi)
2759 : end if
2760 :
2761 34555 : end subroutine xg_nonlop_getDX
2762 :
2763 142425 : subroutine xg_nonlop_getHX(xg_nonlop,Xin,cprjin,cprj_work,work_mpi,Xout)
2764 :
2765 : use iso_c_binding
2766 :
2767 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2768 : type(xgBlock_t), intent(in) :: cprjin
2769 : type(xgBlock_t), intent(inout) :: Xin,cprj_work,work_mpi
2770 : type(xgBlock_t), optional, intent(inout) :: Xout
2771 :
2772 142425 : if (present(Xout)) then
2773 0 : if (xg_nonlop%paw) then
2774 0 : call xg_nonlop_getAX(xg_nonlop,xg_nonlop%Dij_spin,Xin,cprjin,cprj_work,work_mpi,Xout=Xout)
2775 : else
2776 0 : call xg_nonlop_getDX(xg_nonlop,xg_nonlop%ekb%self,Xin,cprjin,cprj_work,work_mpi,Xout=Xout)
2777 : end if
2778 : else
2779 142425 : if (xg_nonlop%paw) then
2780 107870 : call xg_nonlop_getAX(xg_nonlop,xg_nonlop%Dij_spin,Xin,cprjin,cprj_work,work_mpi)
2781 : else
2782 34555 : call xg_nonlop_getDX(xg_nonlop,xg_nonlop%ekb%self,Xin,cprjin,cprj_work,work_mpi)
2783 : end if
2784 : end if
2785 :
2786 142425 : end subroutine xg_nonlop_getHX
2787 :
2788 2312 : subroutine xg_nonlop_getSX(xg_nonlop,Xin,cprjin,cprj_work,work_mpi,Xout)
2789 :
2790 : use iso_c_binding
2791 :
2792 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2793 : type(xgBlock_t), intent(in) :: cprjin
2794 : type(xgBlock_t), intent(inout) :: Xin,cprj_work,work_mpi
2795 : type(xgBlock_t), optional, intent(inout) :: Xout
2796 :
2797 2312 : if (.not.xg_nonlop%paw) then
2798 0 : ABI_ERROR('Not implemented with paw=False.')
2799 : end if
2800 :
2801 2312 : call xg_nonlop_getAX(xg_nonlop,xg_nonlop%Sij%self,Xin,cprjin,cprj_work,work_mpi,Xout)
2802 :
2803 2312 : end subroutine xg_nonlop_getSX
2804 :
2805 284850 : subroutine xg_nonlop_getSm1X(xg_nonlop,Xin,cprjin,cprj_work1,cprj_work2,work_mpi,Xout)
2806 :
2807 : use iso_c_binding
2808 :
2809 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2810 : type(xgBlock_t), intent(in) :: cprjin
2811 : type(xgBlock_t), intent(inout) :: Xin,cprj_work1,cprj_work2,work_mpi
2812 : type(xgBlock_t), optional, intent(inout) :: Xout
2813 :
2814 : integer :: nblocks,nspinor
2815 : integer :: nrows,nrows_cprj
2816 : integer :: ncols,ncols_cprj
2817 :
2818 94950 : if (.not.xg_nonlop%paw) then
2819 0 : ABI_ERROR('Not implemented with paw=False.')
2820 : end if
2821 :
2822 94950 : nblocks = xmpi_comm_size(comm(cprjin))
2823 :
2824 94950 : nspinor = xg_nonlop%nspinor
2825 :
2826 94950 : if (present(Xout)) then
2827 0 : call xgBlock_check(Xin,Xout)
2828 : end if
2829 94950 : call xgBlock_check(cprjin,cprj_work1)
2830 94950 : call xgBlock_check(cprjin,cprj_work2)
2831 94950 : call xgBlock_getsize(Xin,nrows,ncols)
2832 94950 : call xgBlock_getsize(cprjin,nrows_cprj,ncols_cprj)
2833 94950 : if (.not.associated(xg_nonlop%invSij_approx_k)) then
2834 0 : ABI_ERROR('invSij_approx not associated')
2835 : end if
2836 94950 : if (.not.associated(xg_nonlop%gram_proj_k)) then
2837 0 : ABI_ERROR('gram_proj_k should be associated')
2838 : end if
2839 :
2840 : call xg_nonlop_precond_iterative_refinement(xg_nonlop,xg_nonlop%gram_proj_k%self,xg_nonlop%invSij_approx_k%self,&
2841 94950 : & cprjin,cprj_work1,cprj_work2)
2842 94950 : call xgBlock_scale(cprj_work1,-1.0d0,1)
2843 :
2844 : ! Xout = Xin + sum_ai pai cprj_work
2845 94950 : if (present(Xout)) then
2846 0 : call xgBlock_copy(Xin,Xout)
2847 0 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work1,Xout,work_mpi)
2848 : else ! in place version
2849 94950 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work1,Xin,work_mpi)
2850 : end if
2851 :
2852 94950 : end subroutine xg_nonlop_getSm1X
2853 :
2854 450870 : subroutine xg_nonlop_getHmeSX(xg_nonlop,Xin,cprjin,Xout,eigen,cprj_work,work_mpi,no_H)
2855 :
2856 : use iso_c_binding
2857 :
2858 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2859 : type(xgBlock_t), intent(in) :: Xin,eigen,cprjin
2860 : type(xgBlock_t), intent(inout) :: Xout,cprj_work,work_mpi
2861 : logical,optional,intent(in) :: no_H
2862 :
2863 : real(dp) :: tsec(2)
2864 : integer :: nblocks,shift
2865 : integer :: nrows,ncols,nrows_out,ncols_out
2866 : integer :: nrows_diag,ncols_diag
2867 : integer :: nrows_cprj,ncols_cprj
2868 : integer :: nrows_cprj_work,ncols_cprj_work
2869 : integer :: nspinor
2870 : logical :: no_H_
2871 : type(xgBlock_t) :: cprj_work_spinor
2872 :
2873 75145 : call timab(tim_getHmeSX,1,tsec)
2874 :
2875 75145 : if (.not.xg_nonlop%paw) then
2876 0 : ABI_ERROR('Not implemented with paw=False.')
2877 : end if
2878 :
2879 75145 : nblocks = xmpi_comm_size(comm(cprjin))
2880 :
2881 75145 : nspinor = xg_nonlop%nspinor
2882 :
2883 75145 : call xgBlock_getsize(Xin,nrows,ncols)
2884 75145 : call xgBlock_getsize(Xout,nrows_out,ncols_out)
2885 75145 : call xgBlock_getsize(eigen,nrows_diag,ncols_diag)
2886 75145 : call xgBlock_getsize(cprjin,nrows_cprj,ncols_cprj)
2887 75145 : call xgBlock_getsize(cprj_work,nrows_cprj_work,ncols_cprj_work)
2888 75145 : if (ncols/=nrows_diag.or.ncols*nspinor/=nblocks*ncols_cprj.or.ncols_cprj/=ncols_cprj_work.or.ncols/=ncols_out) then
2889 0 : ABI_ERROR('wrong ncols')
2890 : end if
2891 75145 : if (nrows/=nrows_out) then
2892 0 : ABI_ERROR('nrows/=nrows_out')
2893 : end if
2894 75145 : if (nrows_cprj/=nrows_cprj_work) then
2895 0 : ABI_ERROR('nrows_cprj/=nrows_cprj_work')
2896 : end if
2897 75145 : if (ncols_diag/=1) then
2898 0 : ABI_ERROR('ncols_diag should be one')
2899 : end if
2900 :
2901 : ! cprj_work = sum_j Saij cprjin
2902 75145 : call xgBlock_zero(cprj_work)
2903 75145 : call xg_nonlop_apply_Aij(xg_nonlop,xg_nonlop%Sij%self,cprjin,cprj_work)
2904 :
2905 : ! cprj_work = - e cprj_work = -e sum_j Saij cprjin
2906 75145 : shift = xg_nonlop%me_band*ncols_cprj/xg_nonlop%nspinor
2907 75145 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
2908 75145 : call xgBlock_ymax(cprj_work_spinor,eigen,shift,nblocks)
2909 :
2910 75145 : no_H_=.False.
2911 75145 : if (present(no_H)) then
2912 12920 : no_H_ = no_H
2913 : end if
2914 12920 : if (.not.no_H_) then
2915 : ! cprj_work = sum_j Daij cprjin + cprj_work
2916 62225 : call xg_nonlop_apply_Aij(xg_nonlop,xg_nonlop%Dij_spin,cprjin,cprj_work)
2917 : end if
2918 :
2919 : ! Xout = Xout + sum_ai pai cprj_work
2920 75145 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work,Xout,work_mpi)
2921 :
2922 : ! Xout = Xout - e Xin
2923 75145 : call xgBlock_yxmax(Xout,eigen,Xin)
2924 :
2925 75145 : call timab(tim_getHmeSX,2,tsec)
2926 :
2927 75145 : end subroutine xg_nonlop_getHmeSX
2928 : !!***
2929 :
2930 : !!****f* m_xg_nonlop/xg_nonlop_forces
2931 : !! NAME
2932 : !!
2933 : !! FUNCTION
2934 : !!
2935 : !! INPUTS
2936 : !!
2937 19100 : subroutine xg_nonlop_forces_stress(xg_nonlop,Xin,cprjin,cprj_work,eigen,forces,stress,gprimd)
2938 :
2939 : use m_geometry, only : strconv
2940 :
2941 : type(xg_nonlop_t), intent(in) :: xg_nonlop
2942 : type(xgBlock_t), intent(in) :: Xin,eigen
2943 : type(xgBlock_t), intent(in) :: cprjin
2944 : type(xgBlock_t), intent(inout) :: cprj_work
2945 : type(xgBlock_t), optional, intent(inout) :: forces
2946 : type(xgBlock_t), optional, intent(inout) :: stress
2947 : real(dp), optional, intent(in) :: gprimd(:,:)
2948 :
2949 : real(dp) :: tsec(2)
2950 : integer :: cplex,iband
2951 : integer :: nmpi,shift
2952 : integer :: nrows,ncols,ncols_cprj,ncols_cprj_nospin
2953 : integer :: nrows_diag,ncols_diag
2954 : integer :: nspinor,space_cprj
2955 :
2956 : type(xg_t) :: cprj_deriv,work_mpi_deriv,dot,dot_all
2957 19100 : real(dp), pointer :: dot_(:,:),dot_all_(:,:)
2958 : logical :: do_forces,do_stress
2959 : real(dp) :: work(6)
2960 19100 : real(dp), pointer :: stress_(:,:)
2961 : type(xgBlock_t) :: cprjin_spinor,cprj_work_spinor
2962 :
2963 19100 : call timab(tim_forces_stress,1,tsec)
2964 :
2965 19100 : call timab(tim_fst_start,1,tsec)
2966 :
2967 19100 : if (comm(Xin)/=xg_nonlop%comm_band) then
2968 0 : ABI_ERROR('wrong communicator')
2969 : end if
2970 :
2971 19100 : nmpi = xmpi_comm_size(comm(Xin))
2972 19100 : nspinor = xg_nonlop%nspinor
2973 :
2974 19100 : call xgBlock_getsize(Xin,nrows,ncols)
2975 19100 : call xgBlock_getsize(eigen,nrows_diag,ncols_diag)
2976 19100 : if (ncols/=nrows_diag) then
2977 0 : ABI_ERROR('wrong ncols')
2978 : end if
2979 19100 : if (ncols_diag/=1) then
2980 0 : ABI_ERROR('ncols_diag should be one')
2981 : end if
2982 :
2983 19100 : ncols_cprj = ncols*nspinor/nmpi
2984 19100 : ncols_cprj_nospin = ncols/nmpi
2985 :
2986 19100 : space_cprj = xg_nonlop%space_cprj
2987 :
2988 19100 : do_forces = .false.
2989 19100 : if (present(forces)) then
2990 18300 : do_forces = .true.
2991 18300 : if (rows(forces)/=3*xg_nonlop%natom) then
2992 0 : ABI_ERROR('rows(forces)/=3*natom')
2993 : end if
2994 18300 : if (cols(forces)/=ncols) then
2995 0 : ABI_ERROR('cols(forces)/=ncols')
2996 : end if
2997 : end if
2998 :
2999 19100 : do_stress = .false.
3000 19100 : if (present(stress)) then
3001 1140 : do_stress = .true.
3002 1140 : if (.not.present(gprimd)) then
3003 0 : ABI_ERROR('If stress is present, gprimd must be present too')
3004 : end if
3005 1140 : if (rows(stress)/=6) then
3006 0 : ABI_ERROR('rows(stress)/=6')
3007 : end if
3008 1140 : if (cols(stress)/=ncols) then
3009 0 : ABI_ERROR('cols(stress)/=ncols')
3010 : end if
3011 : end if
3012 :
3013 19100 : call xgBlock_zero(cprj_work)
3014 19100 : if (xg_nonlop%paw) then
3015 : ! cprj_work = sum_j Saij cprjin
3016 17212 : call xg_nonlop_apply_Aij(xg_nonlop,xg_nonlop%Sij%self,cprjin,cprj_work)
3017 : ! cprj_work = - e cprj_work = -e sum_j Saij cprjin
3018 17212 : shift = xg_nonlop%me_band*ncols_cprj_nospin
3019 17212 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
3020 17212 : call xgBlock_ymax(cprj_work_spinor,eigen,shift,nmpi)
3021 : end if
3022 :
3023 : ! cprj_work = sum_j Daij cprjin + cprj_work
3024 19100 : if (xg_nonlop%paw) then
3025 17212 : call xg_nonlop_apply_Aij(xg_nonlop,xg_nonlop%Dij_spin,cprjin,cprj_work)
3026 : else
3027 1888 : call xg_nonlop_apply_diag(xg_nonlop,xg_nonlop%ekb%self,cprjin,cprj_work)
3028 : end if
3029 :
3030 19100 : call timab(tim_fst_start,2,tsec)
3031 :
3032 19100 : if (do_forces) then
3033 :
3034 18300 : call xg_init(cprj_deriv ,space_cprj,3*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3035 18300 : call xg_init(work_mpi_deriv,space_cprj,3*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3036 :
3037 18300 : call timab(tim_fst_cprj_deriv_f,1,tsec)
3038 18300 : call xg_nonlop_getcprj_deriv(xg_nonlop,Xin,cprj_deriv%self,work_mpi_deriv%self,DERIV_ATOM)
3039 18300 : call timab(tim_fst_cprj_deriv_f,2,tsec)
3040 :
3041 18300 : call xg_free(work_mpi_deriv)
3042 :
3043 18300 : call xgBlock_zero(forces)
3044 :
3045 18300 : call timab(tim_fst_mult_cprj_f,1,tsec)
3046 18300 : call xg_nonlop_mult_cprj_forces(xg_nonlop,cprj_work,cprj_deriv%self,forces)
3047 18300 : call timab(tim_fst_mult_cprj_f,2,tsec)
3048 :
3049 18300 : call xg_free(cprj_deriv)
3050 :
3051 : end if
3052 :
3053 19100 : if (do_stress) then
3054 :
3055 1140 : call xg_init(cprj_deriv ,space_cprj,6*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3056 1140 : call xg_init(work_mpi_deriv,space_cprj,6*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3057 :
3058 1140 : call timab(tim_fst_cprj_deriv_str,1,tsec)
3059 1140 : call xg_nonlop_getcprj_deriv(xg_nonlop,Xin,cprj_deriv%self,work_mpi_deriv%self,DERIV_STRESS)
3060 1140 : call timab(tim_fst_cprj_deriv_str,2,tsec)
3061 :
3062 1140 : call xg_free(work_mpi_deriv)
3063 :
3064 1140 : call xgBlock_zero(stress)
3065 :
3066 1140 : call timab(tim_fst_mult_cprj_str,1,tsec)
3067 1140 : call xg_nonlop_mult_cprj_stress(xg_nonlop,cprj_work,cprj_deriv%self,stress)
3068 1140 : call timab(tim_fst_mult_cprj_str,2,tsec)
3069 :
3070 1140 : call xg_free(cprj_deriv)
3071 :
3072 1140 : call timab(tim_fst_work_str,1,tsec)
3073 1140 : call xg_init(dot,space_cprj,ncols_cprj_nospin,1)
3074 :
3075 1140 : call xgBlock_reshape_spinor(cprjin,cprjin_spinor,nspinor,COLS2ROWS)
3076 1140 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
3077 1140 : call xgBlock_colwiseDotProduct(cprjin_spinor,cprj_work_spinor,dot%self,comm_loc=xmpi_comm_null)
3078 :
3079 1140 : cplex=1
3080 1140 : if (space_cprj==SPACE_C) cplex=2
3081 :
3082 1140 : call xg_init(dot_all,space_cprj,ncols,1)
3083 1140 : if (xmpi_comm_size(xg_nonlop%comm_band)>1) then
3084 768 : call xgBlock_zero(dot_all%self)
3085 768 : call xgBlock_reverseMap(dot%self,dot_)
3086 768 : call xgBlock_reverseMap(dot_all%self,dot_all_)
3087 768 : shift = cplex*xg_nonlop%me_band*ncols_cprj_nospin
3088 2752 : do iband=1,ncols_cprj_nospin
3089 2752 : dot_all_(1+cplex*(iband-1)+shift,1) = dot_(1+cplex*(iband-1),1)
3090 : end do
3091 768 : call xgBlock_mpi_sum(dot_all%self,comm=xg_nonlop%comm_band)
3092 : else
3093 372 : call xgBlock_reverseMap(dot%self,dot_all_)
3094 : end if
3095 :
3096 1140 : call xgBlock_reverseMap(stress,stress_)
3097 :
3098 11208 : do iband=1,ncols
3099 70476 : work = stress_(:,iband)
3100 10068 : call strconv(work,gprimd,work)
3101 40272 : stress_(1:3,iband) = work(1:3) - dot_all_(1+cplex*(iband-1),1)
3102 41412 : stress_(4:6,iband) = work(4:6)
3103 : end do
3104 :
3105 1140 : call xg_free(dot)
3106 1140 : call xg_free(dot_all)
3107 :
3108 1140 : call timab(tim_fst_work_str,2,tsec)
3109 :
3110 : end if
3111 :
3112 19100 : call timab(tim_forces_stress,2,tsec)
3113 :
3114 19100 : end subroutine xg_nonlop_forces_stress
3115 : !!***
3116 :
3117 18300 : subroutine xg_nonlop_mult_cprj_forces(xg_nonlop,cprj,cprj_deriv,forces)
3118 :
3119 : type(xg_nonlop_t), intent(in) :: xg_nonlop
3120 : type(xgBlock_t), intent(in) :: cprj,cprj_deriv
3121 : type(xgBlock_t), intent(inout) :: forces
3122 :
3123 : !real(dp) :: tsec(2)
3124 : integer :: ia,idir,ilmn,iband,iband_spinor,my_iband,itypat,nlmn,nattyp_i
3125 : integer :: ispinor,iforces,icprj,icprj_deriv
3126 : integer :: ncols_cprj,ncols_cprj_nospin
3127 : integer :: nspinor
3128 : integer :: shift_itypat,shift_itypat_nlmn,shift_itypat_3nlmn
3129 :
3130 18300 : real(dp), pointer :: forces_(:,:)
3131 18300 : complex(dp), pointer :: cprj_(:,:),cprj_deriv_(:,:)
3132 18300 : real(dp), pointer :: cprj_real(:,:),cprj_deriv_real(:,:)
3133 : real(dp) :: forces_tmp
3134 :
3135 18300 : nspinor = xg_nonlop%nspinor
3136 :
3137 18300 : ncols_cprj = cols(cprj)
3138 18300 : ncols_cprj_nospin = ncols_cprj/nspinor
3139 :
3140 18300 : call xgBlock_reverseMap(forces,forces_)
3141 :
3142 28903 : select case(xg_nonlop%space_cprj)
3143 :
3144 : case (SPACE_C)
3145 :
3146 10603 : call xgBlock_reverseMap(cprj,cprj_)
3147 10603 : call xgBlock_reverseMap(cprj_deriv,cprj_deriv_)
3148 :
3149 10603 : shift_itypat=0
3150 10603 : shift_itypat_nlmn=0
3151 10603 : shift_itypat_3nlmn=0
3152 : !$omp parallel default (none) &
3153 : !$omp& shared(xg_nonlop,forces_,cprj_deriv_,cprj_), &
3154 : !$omp& firstprivate(shift_itypat,shift_itypat_nlmn,shift_itypat_3nlmn,ncols_cprj_nospin,nspinor), &
3155 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,idir,iforces,my_iband), &
3156 : !$omp& private(iband_spinor,icprj,icprj_deriv,forces_tmp)
3157 32376 : do itypat = 1, xg_nonlop%ntypat
3158 21773 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
3159 21773 : nattyp_i = xg_nonlop%nattyp(itypat)
3160 32376 : if (nattyp_i>0) then
3161 : !$omp do collapse(3)
3162 91565 : do iband=1,ncols_cprj_nospin
3163 163021 : do ia = 1, nattyp_i
3164 356232 : do idir=1,3
3165 : forces_tmp = zero
3166 512880 : do ispinor=1,nspinor
3167 3146064 : do ilmn=1,nlmn
3168 2633184 : iband_spinor = ispinor + nspinor*(iband-1)
3169 2633184 : icprj = ilmn + nlmn*(ia-1) + shift_itypat_nlmn
3170 2633184 : icprj_deriv = ilmn + nlmn*(idir-1) + 3*nlmn*(ia-1) + shift_itypat_3nlmn
3171 : forces_tmp = forces_tmp &
3172 2931696 : & + 2 * dble(conjg(cprj_deriv_(icprj_deriv,iband_spinor))*cprj_(icprj,iband_spinor))
3173 : end do
3174 : end do
3175 214368 : iforces = idir + 3*(ia-1) + shift_itypat
3176 214368 : my_iband = iband + xg_nonlop%me_band*ncols_cprj_nospin
3177 285824 : forces_(iforces,my_iband) = forces_(iforces,my_iband) + forces_tmp
3178 : end do
3179 : end do
3180 : end do
3181 : !$omp end do
3182 21157 : shift_itypat = shift_itypat + 3*nattyp_i
3183 21157 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
3184 21157 : shift_itypat_3nlmn = shift_itypat_3nlmn + nattyp_i*3*nlmn
3185 : end if
3186 : end do
3187 : !$omp end parallel
3188 :
3189 : case (SPACE_R)
3190 :
3191 7697 : call xgBlock_reverseMap(cprj,cprj_real)
3192 7697 : call xgBlock_reverseMap(cprj_deriv,cprj_deriv_real)
3193 :
3194 7697 : shift_itypat=0
3195 7697 : shift_itypat_nlmn=0
3196 7697 : shift_itypat_3nlmn=0
3197 : !$omp parallel default (none) &
3198 : !$omp& shared(xg_nonlop,forces_,cprj_deriv_real,cprj_real), &
3199 : !$omp& firstprivate(shift_itypat,shift_itypat_nlmn,shift_itypat_3nlmn,ncols_cprj_nospin,nspinor), &
3200 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,idir,iforces,my_iband), &
3201 : !$omp& private(iband_spinor,icprj,icprj_deriv,forces_tmp)
3202 23074 : do itypat = 1, xg_nonlop%ntypat
3203 15377 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
3204 15377 : nattyp_i = xg_nonlop%nattyp(itypat)
3205 23074 : if (nattyp_i>0) then
3206 : !$omp do collapse(3)
3207 50357 : do iband=1,ncols_cprj_nospin
3208 85373 : do ia = 1, nattyp_i
3209 175044 : do idir=1,3
3210 : forces_tmp = zero
3211 210096 : do ispinor=1,nspinor
3212 1648224 : do ilmn=1,nlmn
3213 1438128 : iband_spinor = ispinor + nspinor*(iband-1)
3214 1438128 : icprj = ilmn + nlmn*(ia-1) + shift_itypat_nlmn
3215 1438128 : icprj_deriv = ilmn + nlmn*(idir-1) + 3*nlmn*(ia-1) + shift_itypat_3nlmn
3216 : forces_tmp = forces_tmp &
3217 1543176 : & + 2 * cprj_deriv_real(icprj_deriv,iband_spinor)*cprj_real(icprj,iband_spinor)
3218 : end do
3219 : end do
3220 105048 : iforces = idir + 3*(ia-1) + shift_itypat
3221 105048 : my_iband = iband + xg_nonlop%me_band*ncols_cprj_nospin
3222 140064 : forces_(iforces,my_iband) = forces_(iforces,my_iband) + forces_tmp
3223 : end do
3224 : end do
3225 : end do
3226 : !$omp end do
3227 15377 : shift_itypat = shift_itypat + 3*nattyp_i
3228 15377 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
3229 15377 : shift_itypat_3nlmn = shift_itypat_3nlmn + nattyp_i*3*nlmn
3230 : end if
3231 : end do
3232 : !$omp end parallel
3233 :
3234 : case default
3235 18300 : ABI_ERROR("Wrong space")
3236 :
3237 : end select
3238 :
3239 18300 : call xgBlock_mpi_sum(forces,comm=xg_nonlop%comm_band)
3240 :
3241 18300 : end subroutine xg_nonlop_mult_cprj_forces
3242 :
3243 : !!****f* m_xg_nonlop/xg_nonlop_stress
3244 : !! NAME
3245 : !!
3246 : !! FUNCTION
3247 : !!
3248 : !! INPUTS
3249 : !!
3250 1140 : subroutine xg_nonlop_mult_cprj_stress(xg_nonlop,cprj,cprj_deriv,stress)
3251 :
3252 : type(xg_nonlop_t), intent(in) :: xg_nonlop
3253 : type(xgBlock_t), intent(in) :: cprj
3254 : type(xgBlock_t), intent(in) :: cprj_deriv
3255 : type(xgBlock_t), intent(inout) :: stress
3256 :
3257 : !real(dp) :: tsec(2)
3258 : integer :: ia,idir,ilmn,iband,iband_spinor,my_iband,itypat,nlmn,nattyp_i
3259 : integer :: ispinor,icprj,icprj_deriv
3260 : integer :: ncols_cprj,ncols_cprj_nospin,nspinor
3261 : integer :: shift_itypat,shift_itypat_nlmn,shift_itypat_6nlmn
3262 :
3263 1140 : complex(dp), pointer :: cprj_(:,:),cprj_deriv_(:,:)
3264 1140 : real(dp), pointer :: cprj_real(:,:),cprj_deriv_real(:,:)
3265 1140 : real(dp), pointer :: stress_(:,:)
3266 : real(dp) :: stress_tmp
3267 :
3268 : ! call timab(tim_getHmeSX,1,tsec)
3269 :
3270 1140 : nspinor = xg_nonlop%nspinor
3271 :
3272 1140 : ncols_cprj = cols(cprj)
3273 1140 : ncols_cprj_nospin = ncols_cprj/nspinor
3274 :
3275 1140 : call xgBlock_reverseMap(stress,stress_)
3276 :
3277 1871 : select case(xg_nonlop%space_cprj)
3278 :
3279 : case (SPACE_C)
3280 :
3281 731 : call xgBlock_reverseMap(cprj,cprj_)
3282 731 : call xgBlock_reverseMap(cprj_deriv,cprj_deriv_)
3283 :
3284 731 : shift_itypat=0
3285 731 : shift_itypat_nlmn=0
3286 731 : shift_itypat_6nlmn=0
3287 : !$omp parallel default (none) &
3288 : !$omp& shared(xg_nonlop,stress_,cprj_deriv_,cprj_), &
3289 : !$omp& firstprivate(shift_itypat,shift_itypat_nlmn,shift_itypat_6nlmn,ncols_cprj_nospin,nspinor), &
3290 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,idir,my_iband), &
3291 : !$omp& private(iband_spinor,icprj,icprj_deriv,stress_tmp)
3292 2176 : do itypat = 1, xg_nonlop%ntypat
3293 1445 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
3294 1445 : nattyp_i = xg_nonlop%nattyp(itypat)
3295 2176 : if (nattyp_i>0) then
3296 : !$omp do collapse(2)
3297 7309 : do iband=1,ncols_cprj_nospin
3298 42685 : do idir=1,6
3299 : stress_tmp = zero
3300 83472 : do ispinor=1,nspinor
3301 138288 : do ia = 1, nattyp_i
3302 634560 : do ilmn=1,nlmn
3303 531648 : iband_spinor = ispinor + nspinor*(iband-1)
3304 531648 : icprj = ilmn + nlmn*(ia-1) + shift_itypat_nlmn
3305 531648 : icprj_deriv = ilmn + nlmn*(idir-1) + 6*nlmn*(ia-1) + shift_itypat_6nlmn
3306 : stress_tmp = stress_tmp &
3307 586464 : & + 2 * dble(conjg(cprj_deriv_(icprj_deriv,iband_spinor))*cprj_(icprj,iband_spinor))
3308 : end do
3309 : end do
3310 : end do
3311 35376 : my_iband = iband + xg_nonlop%me_band*ncols_cprj_nospin
3312 41272 : stress_(idir,my_iband) = stress_(idir,my_iband) + stress_tmp
3313 : end do
3314 : end do
3315 : !$omp end do
3316 1413 : shift_itypat = shift_itypat + 6*nattyp_i
3317 1413 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
3318 1413 : shift_itypat_6nlmn = shift_itypat_6nlmn + nattyp_i*6*nlmn
3319 : end if
3320 : end do
3321 : !$omp end parallel
3322 :
3323 : case (SPACE_R)
3324 :
3325 409 : call xgBlock_reverseMap(cprj,cprj_real)
3326 409 : call xgBlock_reverseMap(cprj_deriv,cprj_deriv_real)
3327 :
3328 409 : shift_itypat=0
3329 409 : shift_itypat_nlmn=0
3330 409 : shift_itypat_6nlmn=0
3331 : !$omp parallel default (none) &
3332 : !$omp& shared(xg_nonlop,stress_,cprj_deriv_real,cprj_real), &
3333 : !$omp& firstprivate(shift_itypat,shift_itypat_nlmn,shift_itypat_6nlmn,ncols_cprj_nospin,nspinor), &
3334 : !$omp& private(itypat,nattyp_i,nlmn,ia,ilmn,idir,my_iband), &
3335 : !$omp& private(iband_spinor,icprj,icprj_deriv,stress_tmp)
3336 1210 : do itypat = 1, xg_nonlop%ntypat
3337 801 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
3338 801 : nattyp_i = xg_nonlop%nattyp(itypat)
3339 1210 : if (nattyp_i>0) then
3340 : !$omp do collapse(2)
3341 3653 : do iband=1,ncols_cprj_nospin
3342 20765 : do idir=1,6
3343 : stress_tmp = zero
3344 34224 : do ispinor=1,nspinor
3345 51552 : do ia = 1, nattyp_i
3346 303720 : do ilmn=1,nlmn
3347 269280 : iband_spinor = ispinor + nspinor*(iband-1)
3348 269280 : icprj = ilmn + nlmn*(ia-1) + shift_itypat_nlmn
3349 269280 : icprj_deriv = ilmn + nlmn*(idir-1) + 6*nlmn*(ia-1) + shift_itypat_6nlmn
3350 : stress_tmp = stress_tmp &
3351 286608 : & + 2 * cprj_deriv_real(icprj_deriv,iband_spinor)*cprj_real(icprj,iband_spinor)
3352 : end do
3353 : end do
3354 : end do
3355 17112 : my_iband = iband + xg_nonlop%me_band*ncols_cprj_nospin
3356 19964 : stress_(idir,my_iband) = stress_(idir,my_iband) + stress_tmp
3357 : end do
3358 : end do
3359 : !$omp end do
3360 801 : shift_itypat = shift_itypat + 6*nattyp_i
3361 801 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
3362 801 : shift_itypat_6nlmn = shift_itypat_6nlmn + nattyp_i*6*nlmn
3363 : end if
3364 : end do
3365 : !$omp end parallel
3366 :
3367 : case default
3368 1140 : ABI_ERROR("Wrong space")
3369 :
3370 : end select
3371 :
3372 1140 : call xgBlock_mpi_sum(stress,comm=xg_nonlop%comm_band)
3373 :
3374 :
3375 : ! call timab(tim_getHmeSX,2,tsec)
3376 :
3377 1140 : end subroutine xg_nonlop_mult_cprj_stress
3378 : !!***
3379 :
3380 0 : end module m_xg_nonlop
3381 : !!***
|
