Branch data 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 [ - + - + ]: 163 : subroutine xg_nonlop_init(xg_nonlop,indlmn,my_natom,nattyp,mkmem,ntypat,nspinor,ucvol,usepaw,&
212 [ + - + + : 163 : 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 : 163 : call timab(tim_init,1,tsec)
232 : :
233 : 163 : xg_nonlop%mkmem=mkmem
234 : 163 : xg_nonlop%nspinor=nspinor
235 : 163 : xg_nonlop%comm_atom=comm_atom
236 [ + + ]: 163 : 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 : 55 : xg_nonlop%mpi_atmtab=>null()
244 : : end if
245 : 163 : xg_nonlop%my_natom=my_natom
246 : 163 : xg_nonlop%me_band=me_band
247 : 163 : xg_nonlop%comm_band=comm_band
248 : :
249 [ + - ]: 163 : if (xg_nonlop_option==0.or.xg_nonlop_option==1) then
250 : 163 : xg_nonlop%option = xg_nonlop_option
251 : : else
252 : 0 : ABI_ERROR('Wrong value of xg_nonlop_option')
253 : : end if
254 : :
255 : 163 : xg_nonlop%paw=usepaw==1
256 : :
257 : 163 : xg_nonlop%space_pw=0
258 : 163 : xg_nonlop%space_cprj=0
259 : 163 : xg_nonlop%space_Dij=0
260 : :
261 : 163 : xg_nonlop%weight=four_pi/sqrt(ucvol)
262 : 163 : xg_nonlop%ntypat=ntypat
263 : :
264 : 163 : xg_nonlop%nattyp=>nattyp
265 : 163 : xg_nonlop%indlmn=>indlmn
266 : :
267 : 163 : natom=0
268 [ + + ]: 479 : do itypat=1,ntypat
269 : 479 : natom = natom + nattyp(itypat)
270 : : end do
271 : :
272 [ - + + - : 489 : ABI_MALLOC(xg_nonlop%nlmn_ntypat,(ntypat))
+ - + - -
+ - + -
+ ]
273 [ - + + - : 489 : ABI_MALLOC(xg_nonlop%nlmn_natom,(natom))
+ - + - -
+ - + -
+ ]
274 : :
275 : 163 : cprjdim=0
276 : 163 : shift=0
277 : 163 : nlmn_max=0
278 [ + + ]: 479 : do itypat=1,ntypat
279 [ + + + + ]: 4860 : nlmn = count(indlmn(3,:,itypat)>0)
280 : 316 : if (nlmn>nlmn_max) nlmn_max=nlmn
281 : 316 : xg_nonlop%nlmn_ntypat(itypat) = nlmn
282 : 316 : nattyp_i = nattyp(itypat)
283 [ + + ]: 479 : if (nattyp_i>0) then
284 : 312 : cprjdim = cprjdim + nattyp_i*nlmn
285 [ + + ]: 678 : xg_nonlop%nlmn_natom(1+shift:nattyp_i+shift) = nlmn
286 : : shift = shift + nattyp_i
287 : : end if
288 : : end do
289 : 163 : xg_nonlop%nlmn_max = nlmn_max
290 : 163 : xg_nonlop%natom = natom
291 : 163 : xg_nonlop%cprjdim = cprjdim
292 : :
293 [ - + + - : 1589 : ABI_MALLOC(xg_nonlop%projectors,(mkmem))
+ - + - -
+ - + +
+ ]
294 [ - + + - : 1426 : ABI_MALLOC(xg_nonlop%ffnl_gather,(mkmem))
- + - + +
+ ]
295 [ - + + - : 1426 : ABI_MALLOC(xg_nonlop%ph3d_gather,(mkmem))
- + - + +
+ ]
296 [ + + ]: 163 : if (xg_nonlop%paw) then
297 [ - + + - : 1186 : ABI_MALLOC(xg_nonlop%gram_proj,(mkmem))
- + - + +
+ ]
298 [ - + + - : 1186 : ABI_MALLOC(xg_nonlop%invSij_approx,(mkmem))
- + - + +
+ ]
299 : : end if
300 : :
301 : 163 : nmpi = xmpi_comm_size(xg_nonlop%comm_band)
302 [ - + + - : 489 : ABI_MALLOC(xg_nonlop%l_npw_k,(nmpi))
+ - + - -
+ - + -
+ ]
303 [ - + + - : 326 : ABI_MALLOC(xg_nonlop%l_shift_npw_k,(nmpi))
- + - + -
+ ]
304 : :
305 : 163 : call timab(tim_init,2,tsec)
306 : :
307 : 163 : 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 : 2588 : 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 : 2588 : xg_nonlop%weight=four_pi/sqrt(ucvol)
326 : :
327 : 2588 : 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 [ - + ]: 143 : 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 : 143 : call timab(tim_init,1,tsec)
350 : :
351 : 143 : paral_atom=(xmpi_comm_size(xg_nonlop%comm_atom)>1)
352 : :
353 [ - + + - : 429 : ABI_MALLOC(l_cplex,(xg_nonlop%natom))
+ - + - -
+ - + ]
354 [ + + ]: 469 : l_cplex=0
355 [ + + ]: 283 : do iatom=1,xg_nonlop%my_natom ! loop over atoms treated by this proc
356 [ + + ]: 140 : iatom_tot=iatom;if (paral_atom) iatom_tot=xg_nonlop%mpi_atmtab(iatom)
357 : 283 : l_cplex(iatom_tot)=paw_ij(iatom)%cplex_dij
358 : : end do
359 : 143 : call xmpi_sum(l_cplex,xg_nonlop%comm_atom,ierr)
360 : 143 : xg_nonlop%cplex_alldij = 0
361 [ + + ]: 469 : do iatom=1,xg_nonlop%natom ! loop over all atoms
362 [ + + ]: 469 : if (xg_nonlop%cplex_alldij<l_cplex(iatom)) xg_nonlop%cplex_alldij=l_cplex(iatom)
363 : : end do
364 [ - + ]: 143 : ABI_FREE(l_cplex)
365 : :
366 : 143 : call timab(tim_init,2,tsec)
367 : :
368 : 143 : end subroutine xg_nonlop_init_cplex_alldij
369 : : !!***
370 : :
371 : 163 : subroutine xg_nonlop_destroy(xg_nonlop)
372 : :
373 : : type(xg_nonlop_t),intent(inout) :: xg_nonlop
374 : :
375 : : integer :: ikpt
376 : :
377 [ + + ]: 163 : if (xg_nonlop%paw) then
378 : 139 : call xg_nonlop_destroy_Sij(xg_nonlop) ! Can be destroyed before
379 : 139 : 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 [ - + ]: 163 : ABI_FREE(xg_nonlop%nlmn_ntypat)
385 [ - + ]: 163 : ABI_FREE(xg_nonlop%nlmn_natom)
386 : :
387 [ - + ]: 163 : ABI_FREE(xg_nonlop%l_npw_k)
388 [ - + ]: 163 : ABI_FREE(xg_nonlop%l_shift_npw_k)
389 : :
390 [ + + ]: 1263 : do ikpt=1,xg_nonlop%mkmem
391 : 1100 : call xg_free(xg_nonlop%projectors(ikpt))
392 : 1100 : call xg_free(xg_nonlop%ffnl_gather(ikpt))
393 : 1100 : call xg_free(xg_nonlop%ph3d_gather(ikpt))
394 [ + + ]: 1263 : if (xg_nonlop%paw) then
395 : 908 : call xg_free(xg_nonlop%gram_proj(ikpt))
396 : 908 : call xg_free(xg_nonlop%invSij_approx(ikpt))
397 : : end if
398 : : end do
399 [ - + ]: 163 : ABI_FREE(xg_nonlop%projectors)
400 [ - + ]: 163 : ABI_FREE(xg_nonlop%ffnl_gather)
401 [ - + ]: 163 : ABI_FREE(xg_nonlop%ph3d_gather)
402 [ + + ]: 163 : if (xg_nonlop%paw) then
403 [ - + ]: 139 : ABI_FREE(xg_nonlop%gram_proj)
404 [ - + ]: 139 : ABI_FREE(xg_nonlop%invSij_approx)
405 : : end if
406 : :
407 : 163 : 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 [ - + - + ]: 2204 : 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 : 2204 : real(dp),pointer :: Dij_iatom_(:,:)
462 : : real(dp) :: tsec(2)
463 : : type(xgBlock_t) :: Dij_iatom
464 : :
465 : : ! *************************************************************************
466 : :
467 : 2204 : call timab(tim_make_Dij,1,tsec)
468 : :
469 [ - + ]: 2204 : if (.not.xg_nonlop%paw) then
470 : 0 : ABI_ERROR('Not implemented with paw=False.')
471 : : end if
472 : :
473 : 2204 : nspinor = xg_nonlop%nspinor
474 : 2204 : natom = xg_nonlop%natom
475 : 2204 : nlmn_max = xg_nonlop%nlmn_max
476 : :
477 : 2204 : paral_atom=(xmpi_comm_size(xg_nonlop%comm_atom)>1)
478 : :
479 : 2204 : cplex_alldij = xg_nonlop%cplex_alldij
480 [ + + ]: 2204 : if (cplex_alldij==1) then
481 : 1661 : 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 : 2204 : call xg_init(xg_nonlop%Dij,xg_nonlop%space_Dij,nspinor*nlmn_max,nspinor*nlmn_max*natom*nsppol,xmpi_comm_null)
489 : :
490 [ + + ]: 4472 : do isppol=1,nsppol
491 : :
492 [ + + ]: 6582 : do iatom=1,xg_nonlop%my_natom ! loop over atoms treated by this proc
493 : 2110 : iatom_input = iatom
494 [ + + ]: 2110 : if (paral_atom) iatom_input=xg_nonlop%mpi_atmtab(iatom) ! mpi_atmtab(iatom) has the ordering of the input file
495 : 2110 : iatom_type = atindx(iatom_input) ! convert iatom from input file to iatom ordered by type
496 : 2110 : nlmn=xg_nonlop%nlmn_natom(iatom_type)
497 : 2110 : shift=1+(iatom_type-1)*nspinor*nlmn_max+(isppol-1)*natom*nspinor*nlmn_max
498 : 2110 : call xg_setBlock(xg_nonlop%Dij,Dij_iatom,nspinor*nlmn,nspinor*nlmn,fcol=shift)
499 : 2110 : cplex_dij=paw_ij(iatom)%cplex_dij
500 : 2110 : call xgBlock_reverseMap(Dij_iatom,Dij_iatom_)
501 [ + + ]: 6902 : do jsp=1,nspinor
502 : 2524 : jsps = (jsp-1)*nlmn
503 [ + + ]: 27566 : do jlmn=1,nlmn
504 : 22932 : j0lmn=jlmn*(jlmn-1)/2
505 : 22932 : jjlmn=j0lmn+jlmn
506 [ + + ]: 56252 : do isp=1,nspinor
507 : 30796 : isps = (isp-1)*nlmn
508 [ + + ]: 30796 : 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 [ + + ]: 30796 : if (cplex_dij==1) then
521 : 15068 : Dij_iatom_(cplex_alldij*(jlmn+isps-1)+1,jlmn+jsps) = paw_ij(iatom)%dij(jjlmn,ijsp)
522 [ - + ]: 15068 : 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 [ + + ]: 197334 : do ilmn=1,jlmn-1
528 : : ! see m_hamiltonian:pawdij2ekb and opernlc_ylm
529 : 143606 : ijlmn=j0lmn+ilmn
530 [ + + ]: 174402 : if (cplex_dij==1) then
531 : 66238 : Dij_iatom_(cplex_alldij*(ilmn+isps-1)+1,jlmn+jsps) = paw_ij(iatom)%dij(ijlmn,ijsp)
532 [ - + ]: 66238 : if (cplex_alldij==2) Dij_iatom_(2*(ilmn+isps),jlmn+jsps) = zero
533 : 66238 : Dij_iatom_(cplex_alldij*(jlmn+jsps-1)+1,ilmn+isps) = paw_ij(iatom)%dij(ijlmn,ijsp)
534 [ - + ]: 66238 : 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 [ + + ]: 2204 : if (paral_atom) then
551 : 1692 : call xgBlock_mpi_sum(xg_nonlop%Dij%self,comm=xg_nonlop%comm_atom)
552 : : end if
553 : :
554 : 2204 : call timab(tim_make_Dij,2,tsec)
555 : :
556 : 2204 : end subroutine xg_nonlop_make_Dij
557 : : !!***
558 : :
559 : 2268 : 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 : 2268 : nspinor = xg_nonlop%nspinor
567 : 2268 : natom = xg_nonlop%natom
568 : 2268 : nlmn_max = xg_nonlop%nlmn_max
569 : :
570 [ + + ]: 2268 : 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 : 2268 : shift=1+(isppol-1)*natom*nspinor*nlmn_max
577 : 2268 : call xg_setBlock(xg_nonlop%Dij,xg_nonlop%Dij_spin,nspinor*nlmn_max,nspinor*nlmn_max*natom,fcol=shift)
578 : :
579 : 2268 : end subroutine xg_nonlop_set_Dij_spin
580 : : !!***
581 : :
582 [ - + ]: 139 : 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 : 139 : 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 : 139 : call timab(tim_make_Sij,1,tsec)
599 : :
600 [ - + ]: 139 : if (.not.xg_nonlop%paw) then
601 : 0 : ABI_ERROR('Not implemented with paw=False.')
602 : : end if
603 : :
604 : 139 : ntypat = xg_nonlop%ntypat
605 : 139 : nlmn_max = xg_nonlop%nlmn_max
606 : :
607 : 139 : call xg_init(xg_nonlop%Sij,SPACE_R,nlmn_max,nlmn_max*ntypat,xmpi_comm_self)
608 : :
609 : 139 : inv_sij_ = .false.
610 [ + - ]: 139 : if (present(inv_sij)) inv_sij_ = inv_sij
611 : :
612 [ + + ]: 139 : 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 [ + + ]: 407 : do itypat=1,ntypat
617 : :
618 : 268 : nlmn=xg_nonlop%nlmn_ntypat(itypat)
619 : :
620 : 268 : shift=1+(itypat-1)*nlmn_max
621 : 268 : call xg_setBlock(xg_nonlop%Sij,Sij_itypat,nlmn,nlmn,fcol=shift)
622 : 268 : call xgBlock_reverseMap(Sij_itypat,Sij_itypat_)
623 [ + + ]: 2602 : do jlmn=1,nlmn
624 : 2334 : j0lmn=jlmn*(jlmn-1)/2
625 : 2334 : jjlmn=j0lmn+jlmn
626 : 2334 : Sij_itypat_(jlmn,jlmn) = pawtab(itypat)%sij(jjlmn)
627 [ + + ]: 12481 : do ilmn=1,jlmn-1
628 : 9879 : ijlmn=j0lmn+ilmn
629 : 9879 : Sij_itypat_(ilmn,jlmn) = pawtab(itypat)%sij(ijlmn)
630 : 12213 : Sij_itypat_(jlmn,ilmn) = pawtab(itypat)%sij(ijlmn)
631 : : end do
632 : : end do
633 [ + + ]: 407 : 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 : 139 : call timab(tim_make_Sij,2,tsec)
643 : :
644 : 139 : 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 : 2343 : subroutine xg_nonlop_destroy_Dij(xg_nonlop)
663 : :
664 : : type(xg_nonlop_t),intent(inout) :: xg_nonlop
665 : :
666 : : ! *************************************************************************
667 : :
668 [ - + ]: 2343 : if (.not.xg_nonlop%paw) then
669 : 0 : ABI_ERROR('Not implemented with paw=False.')
670 : : end if
671 : :
672 : 2343 : call xg_free(xg_nonlop%Dij)
673 : :
674 : 2343 : end subroutine xg_nonlop_destroy_Dij
675 : : !!***
676 : :
677 : 139 : subroutine xg_nonlop_destroy_Sij(xg_nonlop)
678 : :
679 : : type(xg_nonlop_t),intent(inout) :: xg_nonlop
680 : :
681 : : ! *************************************************************************
682 : :
683 [ - + ]: 139 : if (.not.xg_nonlop%paw) then
684 : 0 : ABI_ERROR('Not implemented with paw=False.')
685 : : end if
686 : :
687 : 139 : call xg_free(xg_nonlop%Sij)
688 : 139 : call xg_free(xg_nonlop%Sijm1)
689 : :
690 : 139 : end subroutine xg_nonlop_destroy_Sij
691 : : !!***
692 : :
693 : 2322 : subroutine xg_nonlop_compute_projs(xg_nonlop)
694 : :
695 : : type(xg_nonlop_t),intent(inout) :: xg_nonlop
696 : :
697 : 2322 : complex(dp),pointer :: projectors_k_(:,:)
698 : 2322 : real(dp),pointer :: projectors_k_real(:,:)
699 : 2322 : real(dp),pointer :: ffnl_gather_k_(:,:)
700 : 2322 : complex(dp),pointer :: ph3d_gather_k_(:,:)
701 : 2322 : 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 : 2322 : ntypat = xg_nonlop%ntypat
710 : :
711 : : ! 4pi/sqrt(ucvol) * (-i)^l
712 : 2322 : cil(1) = ( 1.0_DP, 0.0_DP) * xg_nonlop%weight
713 : 2322 : cil(2) = ( 0.0_DP,-1.0_DP) * xg_nonlop%weight
714 : 2322 : cil(3) = (-1.0_DP, 0.0_DP) * xg_nonlop%weight
715 : 2322 : cil(4) = ( 0.0_DP, 1.0_DP) * xg_nonlop%weight
716 : :
717 : 2322 : compute_gather = xg_nonlop%option==0
718 : :
719 : 2322 : if (compute_gather) then
720 : 2194 : call xgBlock_zero(xg_nonlop%ffnl_gather_k%self)
721 : 2194 : call xgBlock_zero(xg_nonlop%ph3d_gather_k%self)
722 : 2194 : call xgBlock_reverseMap(xg_nonlop%ffnl_gather_k%self,ffnl_gather_k_)
723 : : end if
724 : :
725 : 2322 : shift_ipw = xg_nonlop%l_shift_npw_k(xg_nonlop%me_band+1)
726 : :
727 : 3658 : select case(xg_nonlop%space_pw)
728 : :
729 : : case (SPACE_C)
730 : :
731 : 1336 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_)
732 [ + + ]: 1336 : if (compute_gather) then
733 : 1208 : 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 [ + + ]: 4054 : do itypat = 1, ntypat
742 : 2718 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
743 : 2718 : nattyp_i = xg_nonlop%nattyp(itypat)
744 [ + + ]: 4054 : if (nattyp_i>0) then
745 : : !! projectors = 4pi/sqrt(ucvol) * (-i)^l * conj(ph3d) * ffnl
746 : : !$omp do collapse(3)
747 [ + + ]: 5358 : do ia = 1, nattyp_i
748 [ + + ]: 30518 : do ilmn=1,nlmn
749 [ + + ]: 2193560 : do ipw=1,xg_nonlop%npw_k
750 : 2165696 : iatom = ia + shift_itypat
751 : 2165696 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
752 : 2165696 : ph3d_ipw = cmplx( xg_nonlop%ph3d_k(1,ipw,iatom), xg_nonlop%ph3d_k(2,ipw,iatom), kind=DP)
753 : 2165696 : ffnl_ipw = xg_nonlop%ffnl_k(ipw, 1, ilmn, itypat)
754 : : !
755 [ + + ]: 2165696 : if (compute_gather) then
756 : 2097408 : ph3d_gather_k_(ipw+shift_ipw,iatom) = ph3d_ipw
757 : 2097408 : ffnl_gather_k_(ipw+shift_ipw,ilmn+(itypat-1)*xg_nonlop%nlmn_max) = ffnl_ipw
758 : : end if
759 : : !
760 : 2165696 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
761 : 2190856 : 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 : 2654 : shift_itypat = shift_itypat + nattyp_i
767 : 2654 : 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 [ + + - ]: 2322 : ABI_ERROR("Wrong space")
820 : :
821 : : end select
822 : :
823 [ + + ]: 2322 : if (compute_gather) then
824 : 2194 : call xgBlock_mpi_sum(xg_nonlop%ffnl_gather_k%self,comm=xg_nonlop%comm_band)
825 : 2194 : call xgBlock_mpi_sum(xg_nonlop%ph3d_gather_k%self,comm=xg_nonlop%comm_band)
826 : : end if
827 : :
828 [ + + ]: 2322 : end subroutine xg_nonlop_compute_projs
829 : : !!***
830 : :
831 : 18301 : 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 : 18301 : complex(dp),pointer :: projectors_k_(:,:)
837 : 18301 : complex(dp),pointer :: projectors_deriv_atom_k_(:,:)
838 : 18301 : real(dp),pointer :: projectors_k_real(:,:)
839 : 18301 : 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 : 18301 : if (.not.associated(xg_nonlop%projectors_k)) then
847 : 0 : ABI_ERROR('projectors_k should be associated')
848 : : end if
849 : :
850 [ - + ]: 18301 : if (rows(projs_deriv_atom)/=xg_nonlop%npw_k) then
851 : 0 : ABI_ERROR('rows(projs_deriv_atom)/=npw_k')
852 : : end if
853 [ - + ]: 18301 : 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 : 18301 : ntypat = xg_nonlop%ntypat
858 : :
859 : 28905 : select case(xg_nonlop%space_pw)
860 : :
861 : : case (SPACE_C)
862 : :
863 : 10604 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_)
864 : 10604 : call xgBlock_reverseMap(projs_deriv_atom,projectors_deriv_atom_k_)
865 : 10604 : shift_itypat_nlmn=0
866 : 10604 : 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 [ + + ]: 32379 : do itypat = 1, ntypat
872 : 21775 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
873 : 21775 : nattyp_i = xg_nonlop%nattyp(itypat)
874 [ + + ]: 32379 : 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 [ + + ]: 42383 : do ia = 1, nattyp_i
878 [ + + ]: 230907 : do ilmn=1,nlmn
879 [ + + ]: 16229123 : do ipw=1,xg_nonlop%npw_k
880 [ + + ]: 64266024 : do idir=1,3
881 : 48058125 : icol = ilmn + (ia-1)*nlmn + shift_itypat_nlmn
882 : 48058125 : ctmp = ( 0.0_DP, -1.0_DP) * two_pi * xg_nonlop%kpg_k(ipw,idir)
883 : 48058125 : icol_deriv = ilmn + (idir-1)*nlmn + (ia-1)*3*nlmn + shift_itypat_3nlmn
884 : 64077500 : 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 : 21159 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
891 : 21159 : 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 [ + + - ]: 18301 : ABI_ERROR("Wrong space")
939 : :
940 : : end select
941 : :
942 [ - + ]: 18301 : end subroutine xg_nonlop_compute_projs_deriv_atom
943 : : !!***
944 : :
945 : 1141 : 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 : 1141 : complex(dp),pointer :: projectors_k_(:,:)
951 : 1141 : complex(dp),pointer :: projectors_deriv_stress_k_(:,:)
952 : 1141 : real(dp),pointer :: projectors_k_real(:,:)
953 : 1141 : 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 : 1141 : if (rows(projs_deriv_stress)/=xg_nonlop%npw_k) then
961 : 0 : ABI_ERROR('rows(projs_deriv_atom)/=npw_k')
962 : : end if
963 [ - + ]: 1141 : 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 : 1141 : ntypat = xg_nonlop%ntypat
968 : :
969 : : ! 4pi/sqrt(ucvol) * (-i)^l
970 : 1141 : cil(1) = ( 1.0_DP, 0.0_DP) * xg_nonlop%weight
971 : 1141 : cil(2) = ( 0.0_DP,-1.0_DP) * xg_nonlop%weight
972 : 1141 : cil(3) = (-1.0_DP, 0.0_DP) * xg_nonlop%weight
973 : 1141 : cil(4) = ( 0.0_DP, 1.0_DP) * xg_nonlop%weight
974 : :
975 : 1873 : select case(xg_nonlop%space_pw)
976 : :
977 : : case (SPACE_C)
978 : :
979 : 732 : call xgBlock_reverseMap(xg_nonlop%projectors_k%self,projectors_k_)
980 : 732 : call xgBlock_reverseMap(projs_deriv_stress,projectors_deriv_stress_k_)
981 : 732 : shift_itypat=0
982 : 732 : shift_itypat_nlmn=0
983 : 732 : 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 [ + + ]: 2179 : do itypat = 1, ntypat
990 : 1447 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
991 : 1447 : nattyp_i = xg_nonlop%nattyp(itypat)
992 [ + + ]: 2179 : 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 [ + + ]: 2895 : do ia = 1, nattyp_i
996 [ + + ]: 17299 : do ilmn=1,nlmn
997 [ + + ]: 1382148 : do ipw=1,xg_nonlop%npw_k
998 : 1366264 : iatom = ia + shift_itypat
999 : 1366264 : il=mod(xg_nonlop%indlmn(1,ilmn,itypat),4)+1
1000 : 1366264 : ph3d_ipw = cmplx( xg_nonlop%ph3d_k(1,ipw,iatom), xg_nonlop%ph3d_k(2,ipw,iatom), kind=DP)
1001 [ + + ]: 5465056 : ffnl_ipw(:) = xg_nonlop%ffnl_k(ipw, 2:4, ilmn, itypat)
1002 [ + + ]: 5465056 : ctmp(:) = - cil(il) * conjg(ph3d_ipw) * xg_nonlop%kpg_k(ipw,:)
1003 : 1366264 : icol_shift = ilmn + (ia-1)*6*nlmn + shift_itypat_6nlmn
1004 : : ! diagonal part
1005 [ + + ]: 5465056 : do idir=1,3
1006 : 4098792 : icol_deriv = icol_shift + (idir-1)*nlmn
1007 : 5465056 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(idir) * ffnl_ipw(idir)
1008 : : end do
1009 : : ! off-diagonal part (which is symmetric)
1010 [ + + ]: 5465056 : ctmp(:) = half*ctmp(:)
1011 : 1366264 : icol_deriv = icol_shift + (4-1)*nlmn
1012 : 1366264 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(2) * ffnl_ipw(3) + ctmp(3) * ffnl_ipw(2)
1013 : 1366264 : icol_deriv = icol_shift + (5-1)*nlmn
1014 : 1366264 : projectors_deriv_stress_k_(ipw,icol_deriv) = ctmp(1) * ffnl_ipw(3) + ctmp(3) * ffnl_ipw(1)
1015 : 1366264 : icol_deriv = icol_shift + (6-1)*nlmn
1016 : 1380668 : 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 : 1415 : shift_itypat = shift_itypat + nattyp_i
1022 : 1415 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
1023 : 1415 : 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 [ + + - ]: 1141 : ABI_ERROR("Wrong space")
1092 : :
1093 : : end select
1094 : :
1095 [ - + ]: 1141 : 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 [ - + - + : 18260 : 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 : 18260 : 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 : 18260 : call timab(tim_make_k,1,tsec)
1260 : :
1261 : 18260 : me_g0_loc = -1
1262 : 18260 : me_g0_fft_loc = -1
1263 [ + + ]: 18260 : if (istwf_k==1) then
1264 : 10373 : xg_nonlop%cplex=2
1265 : 10373 : xg_nonlop%space_pw=SPACE_C
1266 : 10373 : 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 : 18260 : cplex = xg_nonlop%cplex
1280 : 18260 : space_cprj = xg_nonlop%space_cprj
1281 : :
1282 : 18260 : xg_nonlop%npw_k = npw_k
1283 : :
1284 [ + + ]: 73448 : xg_nonlop%l_npw_k(:) = 0
1285 : 18260 : xg_nonlop%l_npw_k(xg_nonlop%me_band+1) = npw_k
1286 : 18260 : call xmpi_sum(xg_nonlop%l_npw_k,xg_nonlop%comm_band,ierr)
1287 : :
1288 [ + + ]: 73448 : xg_nonlop%total_npw_k = sum(xg_nonlop%l_npw_k)
1289 [ + + ]: 73448 : xg_nonlop%max_npw_k = maxval(xg_nonlop%l_npw_k)
1290 : :
1291 : 18260 : xg_nonlop%l_shift_npw_k(1) = 0
1292 : 18260 : nmpi = xmpi_comm_size(xg_nonlop%comm_band)
1293 [ + + ]: 55188 : do iblock=2,nmpi
1294 : 55188 : 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 : 18260 : ntypat = xg_nonlop%ntypat
1298 : 18260 : nlmn_max = xg_nonlop%nlmn_max
1299 : :
1300 : 18260 : xg_nonlop%ph3d_k => ph3d_k
1301 : 18260 : xg_nonlop%ffnl_k => ffnl_k
1302 : :
1303 : 18260 : xg_nonlop%kpg_k => kpg_k
1304 : :
1305 : 18260 : xg_nonlop%projectors_k => xg_nonlop%projectors(ikpt)
1306 [ + + ]: 18260 : if (xg_nonlop%option==0) then
1307 : 17236 : xg_nonlop%ffnl_gather_k => xg_nonlop%ffnl_gather(ikpt)
1308 : 17236 : xg_nonlop%ph3d_gather_k => xg_nonlop%ph3d_gather(ikpt)
1309 : : end if
1310 : :
1311 [ + + ]: 18260 : if (xg_nonlop%paw) then
1312 : 15188 : xg_nonlop%gram_proj_k => xg_nonlop%gram_proj(ikpt)
1313 : 15188 : xg_nonlop%invSij_approx_k => xg_nonlop%invSij_approx(ikpt)
1314 : : end if
1315 : :
1316 [ + + ]: 18260 : if (compute_proj) then
1317 : :
1318 : : call xg_init(xg_nonlop%projectors_k,xg_nonlop%space_pw,npw_k,xg_nonlop%cprjdim,&
1319 : 2322 : xg_nonlop%comm_band,me_g0=me_g0_loc)
1320 : :
1321 [ + + ]: 2322 : if (xg_nonlop%option==0) then
1322 : 2194 : 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 : 2194 : xmpi_comm_null,me_g0=me_g0_fft_loc)
1325 : : end if
1326 : :
1327 : 2322 : call xg_nonlop_compute_projs(xg_nonlop)
1328 : :
1329 : 2322 : compute_invS_approx_=.false.
1330 [ + + ]: 2322 : if (present(compute_invS_approx)) compute_invS_approx_ = compute_invS_approx
1331 [ + + ]: 1161 : 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 : 2322 : compute_gram_=.false.
1366 [ + + ]: 2322 : if (present(compute_gram)) compute_gram_ = compute_gram
1367 [ + + ]: 1161 : 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 : 18260 : call timab(tim_make_k,2,tsec)
1412 : :
1413 : 18260 : end subroutine xg_nonlop_make_k
1414 : : !!***
1415 : :
1416 : 679302 : 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 [ - + ]: 679302 : if (.not.associated(xg_nonlop%projectors_k)) then
1428 : 0 : ABI_ERROR('projectors_k should be associated')
1429 : : end if
1430 : :
1431 : 679302 : projs = xg_nonlop%projectors_k%self
1432 : :
1433 : 679302 : npw = xg_nonlop%npw_k
1434 : 679302 : nspinor = xg_nonlop%nspinor
1435 : 679302 : blocksize = cols(cprjX)
1436 : :
1437 : : ! Check projector sizes
1438 [ - + ]: 679302 : if (npw/=rows(projs)) then
1439 : 0 : ABI_ERROR('npw/=rows(projs)')
1440 : : end if
1441 [ - + ]: 679302 : if (xg_nonlop%cprjdim/=cols(projs)) then
1442 : 0 : ABI_ERROR('cols(projs)/=cprjdim')
1443 : : end if
1444 : : ! Check cprj sizes
1445 [ - + ]: 679302 : 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 : 679302 : nmpi = xmpi_comm_size(comm(cprjX))
1451 [ + + ]: 679302 : if (nmpi==1) then ! sequential
1452 [ - + ]: 197026 : if (xmpi_comm_size(comm(X))/=1) then
1453 : 0 : ABI_ERROR('size(comm(X))/=1')
1454 : : end if
1455 [ - + ]: 197026 : if (rows(X)/=npw*nspinor) then
1456 : 0 : ABI_ERROR('rows(X)/=npw*nspinor')
1457 : : end if
1458 [ - + ]: 197026 : 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 : 679302 : 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 : 728794 : 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 : 364397 : call timab(tim_getcprj,1,tsec)
1509 : :
1510 : 364397 : call xg_nonlop_set_nmpi(xg_nonlop,X,cprjX,nmpi,blocksize,fft_representation)
1511 : :
1512 : 364397 : nspinor = xg_nonlop%nspinor
1513 : 364397 : projs = xg_nonlop%projectors_k%self
1514 : :
1515 : 364397 : call xgBlock_reshape_spinor(X,X_spinor,nspinor,ROWS2COLS)
1516 : :
1517 [ + + ]: 364397 : if (nmpi==1) then
1518 : :
1519 : 106483 : call timab(tim_getcprj_gemm,1,tsec)
1520 : 106483 : call xgBlock_gemm('t','n',1.0d0,projs,X_spinor,0.d0,cprjX)
1521 : 106483 : 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 : 364397 : call timab(tim_getcprj,2,tsec)
1641 : :
1642 : 364397 : 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 : 19442 : 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 : 19442 : npw = xg_nonlop%npw_k
1668 : 19442 : nspinor = xg_nonlop%nspinor
1669 : 19442 : blocksize = cols(cprjX)
1670 : :
1671 : 19442 : if (option==DERIV_ATOM) then
1672 : 18301 : proj_size = 3*xg_nonlop%cprjdim
1673 [ + - ]: 1141 : else if (option==DERIV_STRESS) then
1674 : 1141 : proj_size = 6*xg_nonlop%cprjdim
1675 : : else
1676 : 0 : ABI_ERROR('Bad option')
1677 : : end if
1678 : : ! Check cprj sizes
1679 [ - + ]: 19442 : 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 : 19442 : nmpi = xmpi_comm_size(comm(cprjX))
1684 [ + + ]: 19442 : if (nmpi==1) then ! sequential
1685 [ - + ]: 6498 : if (xmpi_comm_size(comm(X))/=1) then
1686 : 0 : ABI_ERROR('size(comm(X))/=1')
1687 : : end if
1688 [ - + ]: 6498 : if (rows(X)/=npw*nspinor) then
1689 : 0 : ABI_ERROR('rows(X)/=npw*nspinor')
1690 : : end if
1691 [ - + ]: 6498 : 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 : 19442 : call xg_init(projs_deriv,xg_nonlop%space_pw,npw,proj_size,comm=xg_nonlop%comm_band,me_g0=me_g0(X))
1708 : :
1709 : 18301 : select case (option)
1710 : : case (DERIV_ATOM)
1711 : 18301 : call xg_nonlop_compute_projs_deriv_atom(xg_nonlop,projs_deriv%self)
1712 : : case (DERIV_STRESS)
1713 [ + + - ]: 19442 : call xg_nonlop_compute_projs_deriv_stress(xg_nonlop,projs_deriv%self)
1714 : : end select
1715 : :
1716 : 19442 : call xgBlock_reshape_spinor(X,X_spinor,nspinor,ROWS2COLS)
1717 : :
1718 [ + + ]: 19442 : if (nmpi==1) then
1719 : :
1720 : : ! call timab(tim_getcprj_gemm,1,tsec)
1721 : 6498 : 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 : 19442 : call xg_free(projs_deriv)
1753 : : ! call timab(tim_getcprj,2,tsec)
1754 : :
1755 [ + + ]: 19442 : end subroutine xg_nonlop_getcprj_deriv
1756 : : !!***
1757 : :
1758 : 629810 : 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 : 314905 : call timab(tim_apply_prj,1,tsec)
1774 : :
1775 : 314905 : call xg_nonlop_set_nmpi(xg_nonlop,X,cprjX,nmpi,blocksize,fft_representation)
1776 : :
1777 : 314905 : projs = xg_nonlop%projectors_k%self
1778 : :
1779 : 314905 : nspinor = xg_nonlop%nspinor
1780 : 314905 : call xgBlock_reshape_spinor(X,X_spinor,nspinor,ROWS2COLS)
1781 : :
1782 [ + + ]: 314905 : if (nmpi==1) then
1783 : :
1784 : 90543 : call timab(tim_apply_prj_gemm,1,tsec)
1785 : 90543 : call xgBlock_gemm('n','n',1.0d0,projs,cprjX,1.d0,X_spinor)
1786 : 90543 : 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 : 314905 : call timab(tim_apply_prj,2,tsec)
1922 : :
1923 : 314905 : 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 : 2828322 : 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 : 1414161 : real(dp),pointer :: cprjin_nlmn_(:,:),cprjout_nlmn_(:,:)
2039 : 1414161 : real(dp),pointer :: cprjin_(:,:),cprjout_(:,:)
2040 : : real(dp) :: tsec(2)
2041 : : logical :: aij_r2c,A_with_spin_
2042 : :
2043 : 1414161 : call timab(tim_apply_Aij,1,tsec)
2044 : :
2045 : 1414161 : call xgBlock_getsize(cprjin,nrows,ncols)
2046 : :
2047 : 1414161 : cprjdim = xg_nonlop%cprjdim
2048 [ - + ]: 1414161 : if (nrows/=cprjdim) then
2049 : 0 : ABI_ERROR('nrows/=cprjdim')
2050 : : end if
2051 : 1414161 : call xgBlock_check(cprjin,cprjout)
2052 : :
2053 : 1414161 : space_cprj = space(cprjin)
2054 [ + + ]: 1414161 : if (space_cprj==SPACE_C) then
2055 : : cplex=2
2056 : : else
2057 : 522161 : cplex=1
2058 : : end if
2059 : 1414161 : space_aij = space(Aij)
2060 : 1414161 : aij_r2c = .false.
2061 [ + + ]: 1414161 : if (space_cprj/=space_aij) then
2062 [ - + ]: 345960 : 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 : 1414161 : nspinor=xg_nonlop%nspinor
2070 : 1414161 : nlmn_max=xg_nonlop%nlmn_max
2071 : :
2072 : 1414161 : A_with_spin_=.true.
2073 [ + + ]: 1414161 : 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 : 706136 : nlmn_max_1atom = nlmn_max * nspinor
2078 : 706136 : ncols_1atom = ncols / nspinor
2079 : : else
2080 : 708025 : nlmn_max_1atom = nlmn_max
2081 : 708025 : ncols_1atom = ncols
2082 : : end if
2083 : :
2084 [ + + ]: 1414161 : if (aij_r2c) then
2085 [ - + ]: 345960 : if (rows(Aij) /= nlmn_max) then
2086 : 0 : ABI_ERROR('wrong rows for Aij')
2087 : : end if
2088 [ + + ]: 345960 : 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 : 1414161 : call xg_init(cprjin_nlmn_max ,space_cprj,nlmn_max_1atom,ncols_1atom)
2110 : 1414161 : call xg_init(cprjout_nlmn_max,space_cprj,nlmn_max_1atom,ncols_1atom)
2111 : :
2112 : 1414161 : call xgBlock_reverseMap(cprjin ,cprjin_ )
2113 : 1414161 : call xgBlock_reverseMap(cprjout,cprjout_)
2114 : :
2115 : 1414161 : shift = 0
2116 : 1414161 : shift_itypat = 0
2117 [ + + ]: 4233366 : do itypat=1,xg_nonlop%ntypat
2118 : 2819205 : nlmn=xg_nonlop%nlmn_ntypat(itypat)
2119 [ + + ]: 2819205 : if (A_with_spin_) then
2120 : 1378591 : nlmn_1atom = nlmn*nspinor
2121 : : else
2122 : 1440614 : nlmn_1atom = nlmn
2123 : : end if
2124 : 2819205 : nattyp_i=xg_nonlop%nattyp(itypat)
2125 [ + + ]: 4233366 : if (nattyp_i>0) then
2126 : 2769106 : call xg_setBlock(cprjin_nlmn_max ,cprjin_nlmn ,nlmn_1atom,ncols_1atom)
2127 : 2769106 : call xg_setBlock(cprjout_nlmn_max,cprjout_nlmn,nlmn_1atom,ncols_1atom)
2128 : 2769106 : call xgBlock_reverseMap(cprjin_nlmn ,cprjin_nlmn_ )
2129 : 2769106 : call xgBlock_reverseMap(cprjout_nlmn,cprjout_nlmn_)
2130 [ + + ]: 2769106 : if (aij_r2c) call xg_init(Aij_complex,SPACE_C,nlmn_1atom,nlmn_1atom)
2131 [ + + ]: 2769106 : if (.not.loop_over_atoms) then
2132 [ + + ]: 2720058 : if (aij_r2c) then
2133 : 628717 : 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 [ + + ]: 5611940 : do ia=1,nattyp_i
2139 [ + + ]: 2842834 : 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 : 2842834 : call xgBlock_getsize(Aij_iatom,nrows_A,ncols_A)
2148 [ + + ]: 2842834 : 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 [ - + ]: 693232 : if (nrows_A/=nlmn) then
2157 : 0 : ABI_ERROR('nrows_A/=nlmn')
2158 : : end if
2159 [ - + ]: 693232 : 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 [ + + ]: 2842834 : if (aij_r2c) then
2165 : 693232 : call xgBlock_r2c(Aij_iatom,Aij_complex%self,nspinor)
2166 : 693232 : 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 [ + + ]: 2842834 : if (A_with_spin_) then
2172 [ + + ]: 5343620 : do iband=1,ncols_1atom
2173 [ + - + + : 121614500 : cprjin_nlmn_(1:cplex*nlmn,iband) = cprjin_(1+shift:cplex*nlmn+shift,1+nspinor*(iband-1))
+ + ]
2174 : : end do
2175 [ + + ]: 1413584 : 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 : 2842834 : call xgBlock_gemm('n','n',1.0d0,Aij_iatom_,cprjin_nlmn,0.d0,cprjout_nlmn,timing=.false.)
2187 : :
2188 [ + + ]: 2842834 : if (A_with_spin_) then
2189 [ + + ]: 5343620 : 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 [ + - + + : 121614500 : & + cprjout_nlmn_(1:cplex*nlmn,iband)
+ + ]
2192 : : end do
2193 [ + + ]: 1413584 : 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 : 8454774 : shift=shift+cplex*nlmn
2206 : : end do
2207 : :
2208 [ + + ]: 2769106 : if (aij_r2c) call xg_free(Aij_complex)
2209 : 2769106 : shift_itypat = shift_itypat + nattyp_i
2210 : : end if
2211 : :
2212 : : end do
2213 : :
2214 : 1414161 : call xg_free(cprjin_nlmn_max)
2215 : 1414161 : call xg_free(cprjout_nlmn_max)
2216 : :
2217 : 1414161 : call timab(tim_apply_Aij,2,tsec)
2218 : :
2219 : 1414161 : 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 : 442645 : 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 : 442645 : call timab(tim_mult_cprj,1,tsec)
2319 : :
2320 : 442645 : comm_cprj = comm(cprj_right)
2321 : 442645 : nblocks_mpi = xmpi_comm_size(comm_cprj)
2322 : :
2323 : 442645 : cprjdim = xg_nonlop%cprjdim
2324 : 442645 : nspinor = xg_nonlop%nspinor
2325 : :
2326 : 442645 : call xgBlock_reshape_spinor(cprj_right,cprj_right_spinor,nspinor,COLS2ROWS)
2327 : 442645 : call xgBlock_reshape_spinor(cprj_left ,cprj_left_spinor ,nspinor,COLS2ROWS)
2328 : :
2329 : 442645 : nrows_r = rows(cprj_right_spinor)
2330 : 442645 : nrows_l = rows(cprj_left_spinor)
2331 : :
2332 : 442645 : ncols_r = cols(cprj_right_spinor)
2333 : 442645 : ncols_l = cols(cprj_left_spinor)
2334 : :
2335 : 442645 : me_band = xg_nonlop%me_band
2336 : :
2337 [ - + ]: 442645 : if (nrows_r/=nspinor*cprjdim) then
2338 : 0 : ABI_ERROR("rows(cprj_right)/=nspinor*cprjdim")
2339 : : end if
2340 [ - + ]: 442645 : if (nrows_l/=nspinor*cprjdim) then
2341 : 0 : ABI_ERROR("rows(cprj_left)/=nspinor*cprjdim")
2342 : : end if
2343 [ - + ]: 442645 : if (rows(res)/=nblocks_mpi*ncols_l) then
2344 : 0 : ABI_ERROR("rows(res)/=nblocks_mpi*cols(cprj_left)")
2345 : : end if
2346 [ - + ]: 442645 : if (cols(res)/=nblocks_mpi*ncols_r) then
2347 : 0 : ABI_ERROR("cols(res)/=nblocks_mpi*cols(cprj_right)")
2348 : : end if
2349 : :
2350 : 442645 : blocksize_ = ncols_r
2351 [ + - ]: 442645 : if (present(blocksize)) then
2352 [ - + ]: 442645 : if (mod(blocksize,nspinor)/=0) then
2353 : 0 : ABI_ERROR("wrong blocksize (nspinor)")
2354 : : end if
2355 [ - + ]: 442645 : 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 [ + + ]: 442645 : if (nblocks_mpi==1) then
2362 : :
2363 : 159861 : call timab(tim_mult_cprj_gemm,1,tsec)
2364 : 159861 : call xgBlock_gemm('t','n',1.0d0,cprj_left_spinor,cprj_right_spinor,1.d0,res)
2365 : 159861 : 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 : 442645 : call timab(tim_mult_cprj,2,tsec)
2473 : :
2474 : 442645 : 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 : 411092 : 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 : 411092 : call xgBlock_zero(cprj_work)
2599 : 411092 : call xg_nonlop_apply_Aij(xg_nonlop,Aij,cprj_right,cprj_work)
2600 : :
2601 : 411092 : blocksize_ = cols(cprj_right)
2602 [ + - ]: 411092 : if (present(blocksize)) then
2603 : 411092 : blocksize_ = blocksize
2604 : : end if
2605 : 411092 : call xg_nonlop_mult_cprj(xg_nonlop,cprj_left,cprj_work,res,blocksize=blocksize_)
2606 : :
2607 : 411092 : 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 : 255320 : 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 : 255320 : call timab(tim_getXSY,1,tsec)
2642 : :
2643 [ - + ]: 255320 : if (.not.xg_nonlop%paw) then
2644 : 0 : ABI_ERROR('Not implemented with paw=False.')
2645 : : end if
2646 : :
2647 : 255320 : blocksize_ = cols(cprj_right)
2648 [ + + ]: 255320 : if (present(blocksize)) then
2649 : 217666 : blocksize_ = blocksize
2650 : : end if
2651 : 255320 : call xg_nonlop_getXAY(xg_nonlop,xg_nonlop%Sij%self,cprj_left,cprj_right,cprj_work,res,blocksize=blocksize_)
2652 : :
2653 : 255320 : call timab(tim_getXSY,2,tsec)
2654 : :
2655 : 255320 : end subroutine xg_nonlop_getXSY
2656 : : !!***
2657 : :
2658 : 187325 : 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 : 187325 : call timab(tim_getXHY,1,tsec)
2669 : :
2670 : 187325 : blocksize_ = cols(cprj_right)
2671 [ + - ]: 187325 : if (present(blocksize)) then
2672 : 187325 : blocksize_ = blocksize
2673 : : end if
2674 [ + + ]: 187325 : if (xg_nonlop%paw) then
2675 : 155772 : 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 : 187325 : call timab(tim_getXHY,2,tsec)
2681 : :
2682 : 187325 : 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 : 451308 : 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 : 75218 : call timab(tim_getHmeSX,1,tsec)
2874 : :
2875 [ - + ]: 75218 : if (.not.xg_nonlop%paw) then
2876 : 0 : ABI_ERROR('Not implemented with paw=False.')
2877 : : end if
2878 : :
2879 : 75218 : nblocks = xmpi_comm_size(comm(cprjin))
2880 : :
2881 : 75218 : nspinor = xg_nonlop%nspinor
2882 : :
2883 : 75218 : call xgBlock_getsize(Xin,nrows,ncols)
2884 : 75218 : call xgBlock_getsize(Xout,nrows_out,ncols_out)
2885 : 75218 : call xgBlock_getsize(eigen,nrows_diag,ncols_diag)
2886 : 75218 : call xgBlock_getsize(cprjin,nrows_cprj,ncols_cprj)
2887 : 75218 : call xgBlock_getsize(cprj_work,nrows_cprj_work,ncols_cprj_work)
2888 [ + - + - : 75218 : 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 [ - + ]: 75218 : if (nrows/=nrows_out) then
2892 : 0 : ABI_ERROR('nrows/=nrows_out')
2893 : : end if
2894 [ - + ]: 75218 : if (nrows_cprj/=nrows_cprj_work) then
2895 : 0 : ABI_ERROR('nrows_cprj/=nrows_cprj_work')
2896 : : end if
2897 [ - + ]: 75218 : 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 : 75218 : call xgBlock_zero(cprj_work)
2903 : 75218 : 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 : 75218 : shift = xg_nonlop%me_band*ncols_cprj/xg_nonlop%nspinor
2907 : 75218 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
2908 : 75218 : call xgBlock_ymax(cprj_work_spinor,eigen,shift,nblocks)
2909 : :
2910 : 75218 : no_H_=.False.
2911 [ + + ]: 75218 : 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 : 62298 : 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 : 75218 : call xg_nonlop_apply_prj(xg_nonlop,cprj_work,Xout,work_mpi)
2921 : :
2922 : : ! Xout = Xout - e Xin
2923 : 75218 : call xgBlock_yxmax(Xout,eigen,Xin)
2924 : :
2925 : 75218 : call timab(tim_getHmeSX,2,tsec)
2926 : :
2927 : 75218 : 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 [ + + + - : 19101 : 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 : 19101 : real(dp), pointer :: dot_(:,:),dot_all_(:,:)
2958 : : logical :: do_forces,do_stress
2959 : : real(dp) :: work(6)
2960 : 19101 : real(dp), pointer :: stress_(:,:)
2961 : : type(xgBlock_t) :: cprjin_spinor,cprj_work_spinor
2962 : :
2963 : 19101 : call timab(tim_forces_stress,1,tsec)
2964 : :
2965 : 19101 : call timab(tim_fst_start,1,tsec)
2966 : :
2967 [ - + ]: 19101 : if (comm(Xin)/=xg_nonlop%comm_band) then
2968 : 0 : ABI_ERROR('wrong communicator')
2969 : : end if
2970 : :
2971 : 19101 : nmpi = xmpi_comm_size(comm(Xin))
2972 : 19101 : nspinor = xg_nonlop%nspinor
2973 : :
2974 : 19101 : call xgBlock_getsize(Xin,nrows,ncols)
2975 : 19101 : call xgBlock_getsize(eigen,nrows_diag,ncols_diag)
2976 [ - + ]: 19101 : if (ncols/=nrows_diag) then
2977 : 0 : ABI_ERROR('wrong ncols')
2978 : : end if
2979 [ - + ]: 19101 : if (ncols_diag/=1) then
2980 : 0 : ABI_ERROR('ncols_diag should be one')
2981 : : end if
2982 : :
2983 : 19101 : ncols_cprj = ncols*nspinor/nmpi
2984 : 19101 : ncols_cprj_nospin = ncols/nmpi
2985 : :
2986 : 19101 : space_cprj = xg_nonlop%space_cprj
2987 : :
2988 : 19101 : do_forces = .false.
2989 [ + + ]: 19101 : if (present(forces)) then
2990 : 18301 : do_forces = .true.
2991 [ - + ]: 18301 : if (rows(forces)/=3*xg_nonlop%natom) then
2992 : 0 : ABI_ERROR('rows(forces)/=3*natom')
2993 : : end if
2994 [ - + ]: 18301 : if (cols(forces)/=ncols) then
2995 : 0 : ABI_ERROR('cols(forces)/=ncols')
2996 : : end if
2997 : : end if
2998 : :
2999 : 19101 : do_stress = .false.
3000 [ + + ]: 19101 : if (present(stress)) then
3001 : 1141 : do_stress = .true.
3002 [ - + ]: 1141 : if (.not.present(gprimd)) then
3003 : 0 : ABI_ERROR('If stress is present, gprimd must be present too')
3004 : : end if
3005 [ - + ]: 1141 : if (rows(stress)/=6) then
3006 : 0 : ABI_ERROR('rows(stress)/=6')
3007 : : end if
3008 [ - + ]: 1141 : if (cols(stress)/=ncols) then
3009 : 0 : ABI_ERROR('cols(stress)/=ncols')
3010 : : end if
3011 : : end if
3012 : :
3013 : 19101 : call xgBlock_zero(cprj_work)
3014 [ + + ]: 19101 : if (xg_nonlop%paw) then
3015 : : ! cprj_work = sum_j Saij cprjin
3016 : 17213 : 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 : 17213 : shift = xg_nonlop%me_band*ncols_cprj_nospin
3019 : 17213 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
3020 : 17213 : call xgBlock_ymax(cprj_work_spinor,eigen,shift,nmpi)
3021 : : end if
3022 : :
3023 : : ! cprj_work = sum_j Daij cprjin + cprj_work
3024 [ + + ]: 19101 : if (xg_nonlop%paw) then
3025 : 17213 : 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 : 19101 : call timab(tim_fst_start,2,tsec)
3031 : :
3032 [ + + ]: 19101 : if (do_forces) then
3033 : :
3034 : 18301 : call xg_init(cprj_deriv ,space_cprj,3*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3035 : 18301 : call xg_init(work_mpi_deriv,space_cprj,3*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3036 : :
3037 : 18301 : call timab(tim_fst_cprj_deriv_f,1,tsec)
3038 : 18301 : call xg_nonlop_getcprj_deriv(xg_nonlop,Xin,cprj_deriv%self,work_mpi_deriv%self,DERIV_ATOM)
3039 : 18301 : call timab(tim_fst_cprj_deriv_f,2,tsec)
3040 : :
3041 : 18301 : call xg_free(work_mpi_deriv)
3042 : :
3043 : 18301 : call xgBlock_zero(forces)
3044 : :
3045 : 18301 : call timab(tim_fst_mult_cprj_f,1,tsec)
3046 : 18301 : call xg_nonlop_mult_cprj_forces(xg_nonlop,cprj_work,cprj_deriv%self,forces)
3047 : 18301 : call timab(tim_fst_mult_cprj_f,2,tsec)
3048 : :
3049 : 18301 : call xg_free(cprj_deriv)
3050 : :
3051 : : end if
3052 : :
3053 [ + + ]: 19101 : if (do_stress) then
3054 : :
3055 : 1141 : call xg_init(cprj_deriv ,space_cprj,6*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3056 : 1141 : call xg_init(work_mpi_deriv,space_cprj,6*xg_nonlop%cprjdim,ncols_cprj,comm=xg_nonlop%comm_band)
3057 : :
3058 : 1141 : call timab(tim_fst_cprj_deriv_str,1,tsec)
3059 : 1141 : call xg_nonlop_getcprj_deriv(xg_nonlop,Xin,cprj_deriv%self,work_mpi_deriv%self,DERIV_STRESS)
3060 : 1141 : call timab(tim_fst_cprj_deriv_str,2,tsec)
3061 : :
3062 : 1141 : call xg_free(work_mpi_deriv)
3063 : :
3064 : 1141 : call xgBlock_zero(stress)
3065 : :
3066 : 1141 : call timab(tim_fst_mult_cprj_str,1,tsec)
3067 : 1141 : call xg_nonlop_mult_cprj_stress(xg_nonlop,cprj_work,cprj_deriv%self,stress)
3068 : 1141 : call timab(tim_fst_mult_cprj_str,2,tsec)
3069 : :
3070 : 1141 : call xg_free(cprj_deriv)
3071 : :
3072 : 1141 : call timab(tim_fst_work_str,1,tsec)
3073 : 1141 : call xg_init(dot,space_cprj,ncols_cprj_nospin,1)
3074 : :
3075 : 1141 : call xgBlock_reshape_spinor(cprjin,cprjin_spinor,nspinor,COLS2ROWS)
3076 : 1141 : call xgBlock_reshape_spinor(cprj_work,cprj_work_spinor,nspinor,COLS2ROWS)
3077 : 1141 : call xgBlock_colwiseDotProduct(cprjin_spinor,cprj_work_spinor,dot%self,comm_loc=xmpi_comm_null)
3078 : :
3079 : 1141 : cplex=1
3080 [ + + ]: 1141 : if (space_cprj==SPACE_C) cplex=2
3081 : :
3082 : 1141 : call xg_init(dot_all,space_cprj,ncols,1)
3083 [ + + ]: 1141 : 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 : 373 : call xgBlock_reverseMap(dot%self,dot_all_)
3094 : : end if
3095 : :
3096 : 1141 : call xgBlock_reverseMap(stress,stress_)
3097 : :
3098 [ + + ]: 11217 : do iband=1,ncols
3099 [ + + ]: 70532 : work = stress_(:,iband)
3100 [ - + - - : 10076 : call strconv(work,gprimd,work)
- - - - -
+ ]
3101 [ + + ]: 40304 : stress_(1:3,iband) = work(1:3) - dot_all_(1+cplex*(iband-1),1)
3102 [ + + ]: 41445 : stress_(4:6,iband) = work(4:6)
3103 : : end do
3104 : :
3105 : 1141 : call xg_free(dot)
3106 : 1141 : call xg_free(dot_all)
3107 : :
3108 : 1141 : call timab(tim_fst_work_str,2,tsec)
3109 : :
3110 : : end if
3111 : :
3112 : 19101 : call timab(tim_forces_stress,2,tsec)
3113 : :
3114 : 19101 : end subroutine xg_nonlop_forces_stress
3115 : : !!***
3116 : :
3117 : 18301 : 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 : 18301 : real(dp), pointer :: forces_(:,:)
3131 : 18301 : complex(dp), pointer :: cprj_(:,:),cprj_deriv_(:,:)
3132 : 18301 : real(dp), pointer :: cprj_real(:,:),cprj_deriv_real(:,:)
3133 : : real(dp) :: forces_tmp
3134 : :
3135 : 18301 : nspinor = xg_nonlop%nspinor
3136 : :
3137 : 18301 : ncols_cprj = cols(cprj)
3138 : 18301 : ncols_cprj_nospin = ncols_cprj/nspinor
3139 : :
3140 : 18301 : call xgBlock_reverseMap(forces,forces_)
3141 : :
3142 : 28905 : select case(xg_nonlop%space_cprj)
3143 : :
3144 : : case (SPACE_C)
3145 : :
3146 : 10604 : call xgBlock_reverseMap(cprj,cprj_)
3147 : 10604 : call xgBlock_reverseMap(cprj_deriv,cprj_deriv_)
3148 : :
3149 : 10604 : shift_itypat=0
3150 : 10604 : shift_itypat_nlmn=0
3151 : 10604 : 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 [ + + ]: 32379 : do itypat = 1, xg_nonlop%ntypat
3158 : 21775 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
3159 : 21775 : nattyp_i = xg_nonlop%nattyp(itypat)
3160 [ + + ]: 32379 : if (nattyp_i>0) then
3161 : : !$omp do collapse(3)
3162 [ + + ]: 91583 : do iband=1,ncols_cprj_nospin
3163 [ + + ]: 163055 : do ia = 1, nattyp_i
3164 [ + + ]: 356312 : do idir=1,3
3165 : : forces_tmp = zero
3166 [ + + ]: 512976 : do ispinor=1,nspinor
3167 [ + + ]: 3146544 : do ilmn=1,nlmn
3168 : 2633568 : iband_spinor = ispinor + nspinor*(iband-1)
3169 : 2633568 : icprj = ilmn + nlmn*(ia-1) + shift_itypat_nlmn
3170 : 2633568 : icprj_deriv = ilmn + nlmn*(idir-1) + 3*nlmn*(ia-1) + shift_itypat_3nlmn
3171 : : forces_tmp = forces_tmp &
3172 : 2932128 : & + 2 * dble(conjg(cprj_deriv_(icprj_deriv,iband_spinor))*cprj_(icprj,iband_spinor))
3173 : : end do
3174 : : end do
3175 : 214416 : iforces = idir + 3*(ia-1) + shift_itypat
3176 : 214416 : my_iband = iband + xg_nonlop%me_band*ncols_cprj_nospin
3177 : 285888 : forces_(iforces,my_iband) = forces_(iforces,my_iband) + forces_tmp
3178 : : end do
3179 : : end do
3180 : : end do
3181 : : !$omp end do
3182 : 21159 : shift_itypat = shift_itypat + 3*nattyp_i
3183 : 21159 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
3184 : 21159 : 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 [ + + - ]: 18301 : ABI_ERROR("Wrong space")
3236 : :
3237 : : end select
3238 : :
3239 : 18301 : call xgBlock_mpi_sum(forces,comm=xg_nonlop%comm_band)
3240 : :
3241 : 18301 : 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 : 1141 : 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 : 1141 : complex(dp), pointer :: cprj_(:,:),cprj_deriv_(:,:)
3264 : 1141 : real(dp), pointer :: cprj_real(:,:),cprj_deriv_real(:,:)
3265 : 1141 : real(dp), pointer :: stress_(:,:)
3266 : : real(dp) :: stress_tmp
3267 : :
3268 : : ! call timab(tim_getHmeSX,1,tsec)
3269 : :
3270 : 1141 : nspinor = xg_nonlop%nspinor
3271 : :
3272 : 1141 : ncols_cprj = cols(cprj)
3273 : 1141 : ncols_cprj_nospin = ncols_cprj/nspinor
3274 : :
3275 : 1141 : call xgBlock_reverseMap(stress,stress_)
3276 : :
3277 : 1873 : select case(xg_nonlop%space_cprj)
3278 : :
3279 : : case (SPACE_C)
3280 : :
3281 : 732 : call xgBlock_reverseMap(cprj,cprj_)
3282 : 732 : call xgBlock_reverseMap(cprj_deriv,cprj_deriv_)
3283 : :
3284 : 732 : shift_itypat=0
3285 : 732 : shift_itypat_nlmn=0
3286 : 732 : 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 [ + + ]: 2179 : do itypat = 1, xg_nonlop%ntypat
3293 : 1447 : nlmn = xg_nonlop%nlmn_ntypat(itypat)
3294 : 1447 : nattyp_i = xg_nonlop%nattyp(itypat)
3295 [ + + ]: 2179 : if (nattyp_i>0) then
3296 : : !$omp do collapse(2)
3297 [ + + ]: 7327 : do iband=1,ncols_cprj_nospin
3298 [ + + ]: 42799 : do idir=1,6
3299 : : stress_tmp = zero
3300 [ + + ]: 83664 : do ispinor=1,nspinor
3301 [ + + ]: 138576 : do ia = 1, nattyp_i
3302 [ + + ]: 635520 : do ilmn=1,nlmn
3303 : 532416 : iband_spinor = ispinor + nspinor*(iband-1)
3304 : 532416 : icprj = ilmn + nlmn*(ia-1) + shift_itypat_nlmn
3305 : 532416 : icprj_deriv = ilmn + nlmn*(idir-1) + 6*nlmn*(ia-1) + shift_itypat_6nlmn
3306 : : stress_tmp = stress_tmp &
3307 : 587328 : & + 2 * dble(conjg(cprj_deriv_(icprj_deriv,iband_spinor))*cprj_(icprj,iband_spinor))
3308 : : end do
3309 : : end do
3310 : : end do
3311 : 35472 : my_iband = iband + xg_nonlop%me_band*ncols_cprj_nospin
3312 : 41384 : stress_(idir,my_iband) = stress_(idir,my_iband) + stress_tmp
3313 : : end do
3314 : : end do
3315 : : !$omp end do
3316 : 1415 : shift_itypat = shift_itypat + 6*nattyp_i
3317 : 1415 : shift_itypat_nlmn = shift_itypat_nlmn + nattyp_i*nlmn
3318 : 1415 : 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 [ + + - ]: 1141 : ABI_ERROR("Wrong space")
3369 : :
3370 : : end select
3371 : :
3372 : 1141 : call xgBlock_mpi_sum(stress,comm=xg_nonlop%comm_band)
3373 : :
3374 : :
3375 : : ! call timab(tim_getHmeSX,2,tsec)
3376 : :
3377 : 1141 : end subroutine xg_nonlop_mult_cprj_stress
3378 : : !!***
3379 : :
3380 [ # # # # : 0 : end module m_xg_nonlop
# # # # #
# # # # #
# # # # #
# ]
3381 : : !!***
|
