Branch data Line data Source code
1 : : !!****m* ABINIT/m_nonlop
2 : : !! NAME
3 : : !! m_nonlop
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (MT, FDahm)
9 : : !! This file is distributed under the terms of the
10 : : !! GNU General Public License, see ~abinit/COPYING
11 : : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : : !!
13 : : !! SOURCE
14 : :
15 : : #if defined HAVE_CONFIG_H
16 : : #include "config.h"
17 : : #endif
18 : :
19 : : #include "abi_common.h"
20 : :
21 : : ! nvtx related macro definition
22 : : #include "nvtx_macros.h"
23 : :
24 : : module m_nonlop
25 : :
26 : : use, intrinsic :: iso_c_binding, only: c_loc, c_associated
27 : : use defs_basis
28 : : use m_errors
29 : : use m_abicore
30 : : use m_xmpi
31 : : use m_xomp
32 : : use m_cgtools
33 : : use m_gemm_nonlop
34 : : use m_gemm_nonlop_gpu
35 : : use m_gemm_nonlop_projectors
36 : :
37 : : use defs_abitypes, only : MPI_type
38 : : use m_time, only : timab
39 : : use m_fstrings, only : sjoin, itoa, ftoa
40 : : use m_hamiltonian, only : gs_hamiltonian_type, KPRIME_H_K, K_H_KPRIME, K_H_K, KPRIME_H_KPRIME
41 : : use m_pawcprj, only : pawcprj_type, pawcprj_alloc, pawcprj_free, pawcprj_copy
42 : : use m_nonlop_pl, only : nonlop_pl
43 : : use m_nonlop_ylm, only : nonlop_ylm
44 : :
45 : : #if defined HAVE_GPU_CUDA
46 : : use m_manage_cuda
47 : : #endif
48 : : #if defined(HAVE_GPU_MARKERS)
49 : : use m_nvtx_data
50 : : #endif
51 : :
52 : : implicit none
53 : :
54 : : private
55 : : !!***
56 : :
57 : : public :: nonlop
58 : : integer,public,save :: nonlop_counter = 0
59 : : !!***
60 : :
61 : : contains
62 : : !!***
63 : :
64 : : !!****f* ABINIT/nonlop
65 : : !! NAME
66 : : !! nonlop
67 : : !!
68 : : !! FUNCTION
69 : : !! This routine is a driver to compute:
70 : : !! * Application of a nonlocal operator Vnl_k_k^prime in order to get:
71 : : !! - contracted elements (energy, forces, stresses, ...), if signs=1
72 : : !! - a function in reciprocal space (|out> = Vnl|in>), if signs=2
73 : : !! * Optionally, in case of PAW calculation:
74 : : !! - Application of the overlap matrix in reciprocal space (<in|S|in> or (I+S)|in>).
75 : : !! - Application of (Vnl-lambda.S) in reciprocal space
76 : : !! According to user's choice, the routine calls a subroutine, computing all quantities:
77 : : !! - using Legendre Polynomials Pl (Norm-conserving psps only)
78 : : !! - using Spherical Harmonics Ylm (N-conserving or PAW; compulsory for PAW)
79 : : !! - using GPUs (N-conserving or PAW)
80 : : !!
81 : : !! INPUTS
82 : : !! choice: chooses possible output:
83 : : !! choice=0 => do nothing (only compute WF projected with NL projectors)
84 : : !! =1 => non-local energy contribution
85 : : !! =2 => 1st derivative(s) with respect to atomic position(s)
86 : : !! =3 => 1st derivative(s) with respect to strain(s)
87 : : !! =22=> mixed 2nd derivative(s) with respect to atomic pos. and q vector (at q=0)
88 : : !! =25=> mixed 3rd derivative(s) with respect to atomic pos. and two q vectors (at q=0)
89 : : !! =23=> 1st derivative(s) with respect to atomic pos. and
90 : : !! 1st derivative(s) with respect to atomic pos. and strains
91 : : !! =4 => 2nd derivative(s) with respect to 2 atomic pos.
92 : : !! =24=> 1st derivative(s) with respect to atm. pos. and
93 : : !! 2nd derivative(s) with respect to 2 atomic pos.
94 : : !! =33=> mixed 2nd derivative(s) with respect to strain and q vector (at q=0)
95 : : !! =5 => 1st derivative(s) with respect to k wavevector, typically
96 : : !! sum_ij [ |p_i> D_ij <dp_j/dk| + |dp_i/dk> D_ij < p_j| ]
97 : : !! =6 => 2nd derivative(s) with respect to 2 strains and
98 : : !! mixed 2nd derivative(s) with respect to strains & atomic pos.
99 : : !! =51 =>right 1st derivative(s) with respect to k wavevector, typically
100 : : !! sum_ij [ |p_i> D_ij <dp_j/dk| ]
101 : : !! =52 =>left 1st derivative(s) with respect to k wavevector, typically
102 : : !! sum_ij [ |dp_i/dk> D_ij < p_j| ]
103 : : !! =53 =>twist 1st derivative(s) with respect to k, typically
104 : : !! sum_ij [ |dp_i/dk_(idir+1)> D_ij <dp_j//dk_(idir+2)|
105 : : !! where idir, idir+1, idir+2 taken mod 3
106 : : !! =54=> mixed 2nd derivative(s) with respect to atomic pos. and left k wavevector
107 : : !! =55=> mixed 2nd derivative(s) with respect to strain and right k wavevector
108 : : !! =7 => apply operator $\sum_i [ |p_i> <p_i| ],
109 : : !! same as overlap operator with s_ij=identity (paw_opt==3 only)
110 : : !! =8 => 2nd derivatives with respect to 2 k wavevectors
111 : : !! =81=> partial 2nd derivatives with respect to 2 k wavevectors,
112 : : !! full derivative with respect to k1, right derivative with respect to k2,
113 : : !! (derivative with respect to k of choice 51), typically
114 : : !! sum_ij [ |dp_i/dk1> D_ij <dp_j/dk2| + |p_i> D_ij < d2p_j/dk1dk2| ]
115 : : !! Only choices 1,2,3,23,4,5,6 are compatible with useylm=0.
116 : : !! Only choices 1,2,22,25,3,5,33,51,52,53,7,8,81 are compatible with signs=2
117 : : !! cpopt=flag defining the status of cprjin%cp(:)=<Proj_i|Cnk> scalars (see below, side effects)
118 : : !! [enl]=optional (if not present, use hamk%ekb); non-local coeffs connecting projectors
119 : : !! see hamk%ekb description
120 : : !! hamk <type(gs_hamiltonian_type)>=data defining the Hamiltonian at a given k (NL part involved here)
121 : : !! | atindx1(natom)=index table for atoms, inverse of atindx
122 : : !! | dimekb1,dimekb2=dimensions of ekb (see ham%ekb)
123 : : !! | dimekbq=1 if enl factors do not contain a exp(-iqR) phase, 2 is they do
124 : : !! | ekb(dimekb1,dimekb2,nspinor**2,dimekbq)=
125 : : !! | ->NC psps (paw_opt=0): Kleinman-Bylander energies (hartree)
126 : : !! | dimekb1=lmnmax, dimekb2=ntypat
127 : : !! | ->PAW (paw_opt=1 or 4): Dij coeffs connecting projectors (ij symmetric)
128 : : !! | dimekb1=cplex_ekb*lmnmax*(lmnmax+1)/2, dimekb2=natom
129 : : !! | Complex numbers if cplex_ekb=2
130 : : !! | ekb(:,:,1)= Dij^up-up, ekb(:,:,2)= Dij^dn-dn
131 : : !! | ekb(:,:,3)= Dij^up-dn, ekb(:,:,4)= Dij^dn-up (only if nspinor=2)
132 : : !! | ffnl_k(npw_k,dimffnl_k,lmnmax,ntypat)=nonlocal form factors at k
133 : : !! | ffnl_kp(npw_kp,dimffnl_kp,lmnmax,ntypat)=nonlocal form factors at k^prime
134 : : !! | gmet(3,3)=metric tensor for G vecs (in bohr**-2)
135 : : !! | gprimd(3,3)=dimensional reciprocal space primitive translations
136 : : !! | indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=ln (useylm=0) or i=lmn (useylm=1)
137 : : !! | istwf_k=option parameter that describes the storage of wfs at k
138 : : !! | istwf_kp=option parameter that describes the storage of wfs at k^prime
139 : : !! | lmnmax=max. number of (l,m,n) components over all types of atoms
140 : : !! | matblk=dimension of the arrays ph3d_k and ph3d_kp
141 : : !! | mgfft=maximum size of 1D FFTs
142 : : !! | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
143 : : !! | mpssoang= 1+max(spin*angular momentum) for nonlocal pseudopotentials
144 : : !! | natom=number of atoms in cell
145 : : !! | nattyp(ntypat)=number of atoms of each type
146 : : !! | ngfft(18)=contain all needed information about 3D FFT
147 : : !! ! kg_k(3,npw_k)=integer coords of planewaves in basis sphere, for k
148 : : !! ! kg_kp(3,npw_kp)=integer coords of planewaves in basis sphere, for k^prime
149 : : !! ! kpg_k(npw_k,:)= (k+G) components and related data
150 : : !! ! kpg_kp(npw_kp,:)=(k^prime+G) components and related data,
151 : : !! ! kpt(3)=k point in terms of recip. translations
152 : : !! ! kptp(3)=k^prime point in terms of recip. translations
153 : : !! | nloalg(3)=governs the choice of the algorithm for nonlocal operator
154 : : !! ! npw_k=number of (k+G) planewaves
155 : : !! ! npw_kp=number of (k^prime+G) planewaves
156 : : !! | ntypat=number of types of atoms in cell
157 : : !! | nspinor=total number of spinorial components of the wavefunctions
158 : : !! | ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
159 : : !! | ph3d_k(2,npw_k,matblk)=3D structure factors, for each atom and (k+g) plane wave
160 : : !! | ph3d_kp(2,npw_kp,matblk)=3-dim structure factors, for each atom and (k^prime+g) plane wave
161 : : !! | phkxred(2,natom)=phase factors exp(2 pi k.xred)
162 : : !! | phkpxred(2,natom)=phase factors exp(2 pi k^prime.xred)
163 : : !! | sij(dimekb1,ntypat)=overlap matrix components (only if paw_opt=2, 3 or 4)
164 : : !! | ucvol=unit cell volume (bohr^3)
165 : : !! | gpu_option= GPU implementation to use, i.e. cuda, openMP, ... (0=not using GPU)
166 : : !! | useylm=how the NL operator is to be applied: 1=using Ylm, 0=using Legendre polynomials
167 : : !! [iatom_only]=optional. If present (and >0), only projectors related to atom of index iatom_only
168 : : !! will be applied. (used fi to apply derivative of NL operator wrt an atomic displacement)
169 : : !! idir=direction of the - atom to be moved in the case (choice=2,signs=2) or (choice=22,signs=2)
170 : : !! - k point direction in the case (choice=5,51,or 52)
171 : : !! for choice 53 signs=2, cross derivatives are in idir+1 and idir+2 directions (mod 3)
172 : : !! - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
173 : : !! - strain component (1:9) in the case (choice=33,signs=2)
174 : : !! - (1:9) components to specify the atom to be moved and the second q-gradient
175 : : !! direction in the case (choice=25,signs=2)
176 : : !! lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
177 : : !! Typically lambda is the eigenvalue (or its guess)
178 : : !! mpi_enreg=information about MPI parallelization
179 : : !! ndat=number of wavefunctions on which to apply nonlop
180 : : !! nnlout=dimension of enlout (when signs=1 and choice>0):
181 : : !! ==== if paw_opt=0, 1 or 2 ====
182 : : !! choice nnlout | choice nnlout
183 : : !! 1 1 | 51 6 (complex)
184 : : !! 2 3*natom | 52 6 (complex)
185 : : !! 3 6 | 53 6 (complex)
186 : : !! 4 6*natom | 54 9*natom
187 : : !! 23 6+3*natom | 55 36 (complex)
188 : : !! 24 9*natom | 6 36+18*natom
189 : : !! 5 3 | 8 6
190 : : !! | 81 18 (complex)
191 : : !! ==== if paw_opt=3 ====
192 : : !! choice nnlout
193 : : !! 1 1
194 : : !! 2 3*natom
195 : : !! 5 3
196 : : !! 51 3
197 : : !! 52 3
198 : : !! 54 9*natom
199 : : !! 55 36
200 : : !! 7 1
201 : : !! 8 6
202 : : !! 81 9
203 : : !! ==== if paw_opt=4 ====
204 : : !! not available
205 : : !! [only_SO]=optional, flag to calculate only the SO part in nonlop
206 : : !! paw_opt= define the nonlocal operator concerned with:
207 : : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
208 : : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
209 : : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
210 : : !! paw_opt=3 : PAW overlap matrix (Sij)
211 : : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
212 : : !! [qdir]= optional, direction of the q-gradient (only for choice=22, choice=25 and choice=33)
213 : : !! [select_k]=optional, option governing the choice of k points to be used.
214 : : !! hamk datastructure contains quantities needed to apply NL operator
215 : : !! in reciprocal space between 2 kpoints, k and k^prime (equal in most cases);
216 : : !! if select_k=1, <k^prime|Vnl|k> is applied [default]
217 : : !! if select_k=2, <k|Vnl|k^prime> is applied
218 : : !! if select_k=3, <k|Vnl|k> is applied
219 : : !! if select_k=4, <k^prime|Vnl|k^prime> is applied
220 : : !! signs= if 1, get contracted elements (energy, forces, stress, ...)
221 : : !! if 2, applies the non-local operator to a function in reciprocal space
222 : : !! tim_nonlop=timing code of the calling routine (can be set to 0 if not attributed)
223 : : !! vectin(2,npwin*my_nspinor*ndat)=input cmplx wavefunction coefficients <G|Cnk>
224 : : !! vectproj(2,nprojs,my_nspinor*ndat)=Optional, vector to be used instead of cprjin%cp when provided
225 : : !!
226 : : !! OUTPUT
227 : : !! ==== if (signs==1) ====
228 : : !! --If (paw_opt==0, 1 or 2)
229 : : !! enlout(nnlout)= contribution to the non-local part of the following properties:
230 : : !! if choice=1 : enlout(1) -> the energy
231 : : !! if choice=2 : enlout(3*natom) -> 1st deriv. of energy wrt atm. pos (forces)
232 : : !! if choice=3 : enlout(6) -> 1st deriv. of energy wrt strain (stresses)
233 : : !! if choice=4 : enlout(6*natom) -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
234 : : !! if choice=23: enlout(6+3*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
235 : : !! 1st deriv. of energy wrt strain (stresses)
236 : : !! if choice=24: enlout(9*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
237 : : !! 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
238 : : !! if choice=5 : enlout(3) -> 1st deriv. of energy wrt k
239 : : !! if choice=51: enlout(3) -> 1st deriv. (right) of energy wrt k
240 : : !! if choice=52: enlout(3) -> 1st deriv. (left) of energy wrt k
241 : : !! if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
242 : : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
243 : : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
244 : : !! if choice=6 : enlout(36+18*natom) -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
245 : : !! 2nd deriv. of energy wrt to atm. pos and strain (internal strain)
246 : : !! if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
247 : : !! if choice=81: enlout(9) -> 2nd deriv.of E: full derivative w.r.t. k1, right derivative w.r.t k2
248 : : !! --If (paw_opt==3)
249 : : !! if choice=1 : enlout(1) -> contribution to <c|S|c> (note: not including <c|c>)
250 : : !! if choice=2 : enlout(3*natom) -> contribution to <c|dS/d_atm.pos|c>
251 : : !! if choice=51: enlout(3) -> contribution to <c|d(right)S/d_k|c>
252 : : !! if choice=52: enlout(3) -> contribution to <c|d(left)S/d_k|c>
253 : : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
254 : : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
255 : : !! if choice=7 : enlout(1) -> contribution to <c|sum_i[p_i><p_i]|c>
256 : : !! if choice=8 : enlout(6) -> contribution to <c|d2S/d_k1d_k2|c>
257 : : !! if choice=81: enlout(9) -> contribution to <c|dS/d_k1[d(right)d_k2]|c>
258 : : !! --If (paw_opt==4)
259 : : !! not available
260 : : !! ==== if (signs==2) ====
261 : : !! --if (paw_opt=0)
262 : : !! vectout(2,npwout*my_nspinor*ndat)=result of the application of the concerned operator
263 : : !! or one of its derivatives to the input vect.
264 : : !! if (choice=22) <G|d2V_nonlocal/d(atm. pos)dq|vect_in> (at q=0)
265 : : !! if (choice=25) <G|d3V_nonlocal/d(atm. pos)dqdq|vect_in> (at q=0)
266 : : !! if (choice=33) <G|d2V_nonlocal/d(strain)dq|vect_in> (at q=0)
267 : : !! --if (paw_opt=0, 1 or 4)
268 : : !! vectout(2,npwout*my_nspinor*ndat)=result of the application of the concerned operator
269 : : !! or one of its derivatives to the input vect.:
270 : : !! if (choice=1) <G|V_nonlocal|vect_in>
271 : : !! if (choice=2) <G|dV_nonlocal/d(atm. pos)|vect_in>
272 : : !! if (choice=3) <G|dV_nonlocal/d(strain)|vect_in>
273 : : !! if (choice=5) <G|dV_nonlocal/d(k)|vect_in>
274 : : !! if (choice=51) <G|d(right)V_nonlocal/d(k)|vect_in>
275 : : !! if (choice=52) <G|d(left)V_nonlocal/d(k)|vect_in>
276 : : !! if (choice=53) <G|d(twist)V_nonlocal/d(k)|vect_in>
277 : : !! if (choice=8) <G|d2V_nonlocal/d(k)d(k)|vect_in>
278 : : !! if (choice=81) <G|d[d(right)V_nonlocal/d(k)]/d(k)|vect_in>
279 : : !! --if (paw_opt=2)
280 : : !! vectout(2,npwout*my_nspinor*ndat)=final vector in reciprocal space:
281 : : !! if (choice=1) <G|V_nonlocal-lambda.(I+S)|vect_in>
282 : : !! if (choice=2) <G|d[V_nonlocal-lambda.(I+S)]/d(atm. pos)|vect_in>
283 : : !! if (choice=3) <G|d[V_nonlocal-lambda.(I+S)]/d(strain)|vect_in>
284 : : !! if (choice=5) <G|d[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
285 : : !! if (choice=51) <G|d(right)[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
286 : : !! if (choice=52) <G|d(left)[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
287 : : !! if (choice=53) <G|d(twist)[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
288 : : !! if (choice=8) <G|d2[V_nonlocal-lambda.(I+S)]/d(k)d(k)|vect_in>
289 : : !! if (choice=81) <G|d[d(right[V_nonlocal-lambda.(I+S)]/d(k)]/d(k)|vect_in>
290 : : !! --if (paw_opt=3 or 4)
291 : : !! svectout(2,npwout*my_nspinor*ndat)=result of the application of Sij (overlap matrix)
292 : : !! or one of its derivatives to the input vect.:
293 : : !! if (choice=1) <G|I+S|vect_in>
294 : : !! if (choice=2) <G|dS/d(atm. pos)|vect_in>
295 : : !! if (choice=3) <G|dS/d(strain)|vect_in>
296 : : !! if (choice=5) <G|dS/d(k)|vect_in>
297 : : !! if (choice=51) <G|d(right)S/d(k)|vect_in>
298 : : !! if (choice=52) <G|d(left)S/d(k)|vect_in>
299 : : !! if (choice=53) <G|d(twist)S/d(k)|vect_in>
300 : : !! if (choice=3) <G|d[V_nonlocal-lambda.(I+S)]/d(strain)|vect_in>
301 : : !! if (choice=7) <G|sum_i[p_i><p_i]|vect_in>
302 : : !! if (choice=8) <G|d2S/d(k)d(k)|vect_in>
303 : : !! if (choice=81) <G|d[d(right)S/d(k)]/d(k)|vect_in>
304 : : !!
305 : : !! SIDE EFFECTS
306 : : !! ==== ONLY IF useylm=1
307 : : !! cprjin(natom,my_nspinor*ndat) <type(pawcprj_type)>=projected input wave function |in> on non-local projectors
308 : : !! =<p_lmn|in> and derivatives
309 : : !! Treatment depends on cpopt parameter:
310 : : !! if cpopt=-1, <p_lmn|in> (and derivatives)
311 : : !! are computed here (and not saved)
312 : : !! if cpopt= 0, <p_lmn|in> are computed here and saved
313 : : !! derivatives are eventually computed but not saved
314 : : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
315 : : !! other derivatives are eventually computed but not saved
316 : : !! if cpopt= 2 <p_lmn|in> are already in memory;
317 : : !! first (and 2nd) derivatives are computed here and not saved
318 : : !! if cpopt= 3 <p_lmn|in> are already in memory;
319 : : !! first derivatives are computed here and saved
320 : : !! other derivatives are eventually computed but not saved
321 : : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
322 : : !! other derivatives are not computed
323 : : !! This option is not compatible with choice=4,24 or 6
324 : : !! If useylm=0, must have cpopt=-1!
325 : : !! Warning: for cpopt= 1 or 3, derivatives wrt strains do not contain
326 : : !! the contribution due to the volume change;
327 : : !! i.e. <dp_lmn/dEps|in> are incomplete.
328 : : !!
329 : : !! NOTES
330 : : !! * See nonlop_pl and nonlop_ylm to have more comments...
331 : : !! * In the case signs=1, the array vectout is not used.
332 : : !!
333 : : !! SOURCE
334 : :
335 [ - + - + ]: 102392042 : subroutine nonlop(choice,cpopt,cprjin,enlout,hamk,idir,lambda,mpi_enreg,ndat,nnlout,&
336 [ - + - + : 51196021 : paw_opt,signs,svectout,tim_nonlop,vectin,vectout,&
- + ]
337 [ + + + - : 204784084 : cprjin_left,enl,enl_ndat,enlout_im,iatom_only,ndat_left,only_SO,qdir,select_k,vectproj) !optional arguments
+ + + - +
+ + - + +
+ - - + +
+ + - -
+ ]
338 : :
339 : : !Arguments ------------------------------------
340 : : !scalars
341 : : integer,intent(in) :: choice,cpopt,idir,ndat,nnlout,paw_opt,signs,tim_nonlop
342 : : integer,intent(in),optional :: iatom_only,only_SO,qdir,ndat_left,select_k
343 : : type(MPI_type),intent(in) :: mpi_enreg
344 : : type(gs_hamiltonian_type),intent(in),target :: hamk
345 : : !arrays
346 : : real(dp),intent(in) :: lambda(ndat)
347 : : real(dp),contiguous, intent(in),target,optional :: enl(:,:,:,:),enl_ndat(:,:,:,:,:)
348 : : real(dp),intent(inout),target :: vectin(:,:)
349 : : real(dp),intent(out),target :: enlout(:),svectout(:,:)
350 : : real(dp),intent(out),optional :: enlout_im(:)
351 : : real(dp),intent(inout),target :: vectout(:,:)
352 : : type(pawcprj_type),intent(inout),target :: cprjin(:,:)
353 : : type(pawcprj_type),intent(inout),target,optional :: cprjin_left(:,:)
354 : : real(dp),intent(inout), contiguous, optional :: vectproj(:,:,:)
355 : :
356 : : !Local variables-------------------------------
357 : : !scalars
358 : : integer :: dimenl1,dimenl2,dimenl2_,dimekbq,dimffnlin,dimffnlout,dimsij,iatm,iatom_only_,idat
359 : : integer :: ii,ispden,ispinor,istwf_k,itypat,jspinor,matblk_,my_nspinor,n1,n2,n3,natom_,ncpgr_atm,ndat_left_
360 : : integer :: nkpgin,nkpgout,npwin,npwout,ntypat_,only_SO_,select_k_,shift1,shift2,shift3
361 : : logical :: atom_pert,force_recompute_ph3d,kpgin_allocated,kpgout_allocated, use_gemm_nonlop
362 : : !character(len=500) :: msg
363 : : !arrays
364 : : integer :: nlmn_atm(1),nloalg_(3)
365 : 51196021 : integer,pointer :: kgin(:,:),kgout(:,:)
366 : 51196021 : integer, contiguous, pointer :: atindx1_(:),indlmn_(:,:,:),nattyp_(:)
367 : : real(dp) :: tsec(2)
368 : 51196021 : real(dp), contiguous, pointer :: enl_ptr(:,:,:,:),enl_ndat_ptr(:,:,:,:,:)
369 : 51196021 : real(dp),pointer :: ffnlin(:,:,:,:),ffnlin_(:,:,:,:),ffnlout(:,:,:,:),ffnlout_(:,:,:,:)
370 : 51196021 : real(dp),pointer :: kpgin(:,:),kpgout(:,:)
371 : : real(dp) :: kptin(3),kptout(3)
372 : 51196021 : real(dp),pointer :: ph3din(:,:,:),ph3din_(:,:,:),ph3dout(:,:,:),ph3dout_(:,:,:)
373 : 51196021 : real(dp),pointer :: phkxredin(:,:),phkxredin_(:,:),phkxredout(:,:),phkxredout_(:,:)
374 : 51196021 : real(dp), contiguous, pointer :: ph1d_(:,:),sij_(:,:)
375 : 51196021 : real(dp), contiguous, pointer :: enl__(:,:,:,:),enl_ndat_(:,:,:,:,:)
376 : 51196021 : type(pawcprj_type),pointer :: cprjin_(:,:)
377 : : integer :: b0,b1,b2,b3,b4,e0,e1,e2,e3,e4
378 : : integer :: proj_shift,ia,nlmn
379 : : integer :: shift,shift_forces,shift_stress
380 : : integer :: nnlout_forces,nnlout_stress
381 : 51196021 : real(dp), allocatable :: enlout_forces(:),enlout_stress(:)
382 : : ! **********************************************************************
383 : :
384 : : DBG_ENTER("COLL")
385 : :
386 : : !Keep track of time spent in this routine (selection of different slots for different choices)
387 : 51196021 : call timab(220+tim_nonlop,1,tsec)
388 : :
389 : : ! Increment global counter
390 : : !$OMP MASTER
391 : 51196021 : nonlop_counter = nonlop_counter + ndat
392 : : !$OMP END MASTER
393 : :
394 [ + + ]: 51196021 : only_SO_=0; if (present(only_SO)) only_SO_=only_SO
395 : 51196021 : my_nspinor=max(1,hamk%nspinor/mpi_enreg%nproc_spinor)
396 : :
397 : 51196021 : force_recompute_ph3d=.false.
398 : :
399 : : ! Error(s) on incorrect input
400 [ + + ]: 51196021 : if (hamk%useylm==0) then
401 [ - + ]: 29075308 : if (paw_opt>0) then
402 : 0 : ABI_BUG('When paw_opt>0 you must use ylm version of nonlop! Set useylm 1.')
403 : : end if
404 [ - + ]: 29075308 : if (cpopt/=-1) then
405 : 0 : ABI_BUG('If useylm=0, ie no PAW, then cpopt/=-1 is not allowed !')
406 : : end if
407 [ - + ]: 29075308 : if (hamk%dimekbq/=1) then
408 : 0 : ABI_BUG('If useylm=0, ie no PAW, then dimekbq/=-1 is not allowed !')
409 : : end if
410 [ - + ]: 29075308 : if (hamk%gpu_option/=ABI_GPU_DISABLED) then
411 : 0 : ABI_BUG('When gpu_option/=0 you must use ylm version of nonlop! Set useylm to 1.')
412 : : end if
413 : : end if
414 [ - + - - ]: 51196021 : if (hamk%gpu_option/=ABI_GPU_DISABLED.and.hamk%dimekbq/=1) then
415 : 0 : ABI_BUG('GPU version of nonlop not compatible with a exp(-iqR) phase!')
416 : : end if
417 [ + - - + ]: 51196021 : if ((.not.associated(hamk%kg_k)).or.(.not.associated(hamk%kg_kp))) then
418 : 0 : ABI_BUG('kg_k/kg_kp should be associated!')
419 : : end if
420 [ + - - + ]: 51196021 : if ((.not.associated(hamk%ffnl_k)).or.(.not.associated(hamk%ffnl_kp))) then
421 : 0 : ABI_BUG('ffnl_k/ffnl_kp should be associated!')
422 : : end if
423 : : !if (hamk%istwf_k/=hamk%istwf_kp) then
424 : : ! ABI_BUG('istwf has to be the same for both k-points.')
425 : : !end if
426 : :
427 [ - + ]: 51196021 : if (present(enl) .and. present(enl_ndat)) then
428 : 0 : ABI_BUG("enl and enl_ndat cannot be specified concurrently !")
429 : : end if
430 : :
431 : : !Select k-dependent objects according to select_k input parameter
432 [ + + ]: 51196021 : select_k_=KPRIME_H_K;if (present(select_k)) select_k_=select_k
433 : : ! If both K-Kprime variant of each attribute of hamiltonian share the same
434 : : ! address, we can assume select_k==K_H_K.
435 : : if ( c_associated(c_loc(hamk%ffnl_k), c_loc(hamk%ffnl_kp)) &
436 [ + - - - : 82583513 : .and. c_associated(c_loc(hamk%kg_k), c_loc(hamk%kg_kp))) then
- - - - -
- - - + -
- - - - -
- - - - -
+ - - - -
- - - + -
- - - - -
- + + + +
- + - - -
- - + - -
- - - + -
- - - - -
- - - + -
- - - - -
- - + + ]
437 [ + + + - ]: 19808529 : if (associated(hamk%ph3d_k).and.associated(hamk%ph3d_kp)) then
438 [ - + - - : 19788137 : if (c_associated(c_loc(hamk%ph3d_k), c_loc(hamk%ph3d_kp))) then
- - - - -
- - + - -
- - - - -
- - + - -
- - - - -
+ - - - -
- - + - ]
439 : 19788137 : select_k_=K_H_K
440 : : end if
441 : : else
442 : 20392 : select_k_=K_H_K
443 : : end if
444 : : end if
445 : 51196021 : nkpgin=0;nkpgout=0;nullify(kpgin);nullify(kpgout)
446 : 51196021 : nullify(ph3din);nullify(ph3dout)
447 : : !print *, "in nonlop with select_k:", select_k_
448 : :
449 [ + + ]: 51196021 : if (select_k_==KPRIME_H_K) then
450 : : ! ===== <k^prime|Vnl|k> =====
451 [ + + + + ]: 69041420 : kptin = hamk%kpt_k ; kptout = hamk%kpt_kp
452 : 9863060 : npwin=hamk%npw_fft_k ; npwout=hamk%npw_fft_kp
453 : 9863060 : kgin => hamk%kg_k ; kgout => hamk%kg_kp
454 [ + + ]: 9863060 : if (associated(hamk%kpg_k)) then
455 : 9862676 : kpgin => hamk%kpg_k ; nkpgin=size(kpgin,2)
456 : : end if
457 [ + + ]: 9863060 : if (associated(hamk%kpg_kp)) then
458 : 9862676 : kpgout => hamk%kpg_kp ; nkpgout=size(kpgout,2)
459 : : end if
460 : 9863060 : phkxredin => hamk%phkxred ; phkxredout => hamk%phkpxred
461 : 9863060 : ffnlin => hamk%ffnl_k ; ffnlout => hamk%ffnl_kp
462 [ + + ]: 9863060 : if (associated(hamk%ph3d_k )) ph3din => hamk%ph3d_k
463 [ + - ]: 9863060 : if (associated(hamk%ph3d_kp)) ph3dout => hamk%ph3d_kp
464 [ + + - + ]: 9863060 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_k).and.associated(hamk%ph3d_kp)))
465 : 9863060 : istwf_k=hamk%istwf_k
466 [ + + ]: 41332961 : else if (select_k_==K_H_KPRIME) then
467 : : ! ===== <k|Vnl|k^prime> =====
468 [ + + + + ]: 250292 : kptin = hamk%kpt_kp ; kptout = hamk%kpt_k
469 : 35756 : npwin=hamk%npw_fft_kp ; npwout=hamk%npw_fft_k
470 : 35756 : kgin => hamk%kg_kp ; kgout => hamk%kg_k
471 [ + - ]: 35756 : if (associated(hamk%kpg_kp)) then
472 : 35756 : kpgin => hamk%kpg_kp ; nkpgin=size(kpgin,2)
473 : : end if
474 [ + - ]: 35756 : if (associated(hamk%kpg_k)) then
475 : 35756 : kpgout => hamk%kpg_k ; nkpgout=size(kpgout,2)
476 : : end if
477 : 35756 : phkxredin => hamk%phkpxred ; phkxredout => hamk%phkxred
478 : 35756 : ffnlin => hamk%ffnl_kp ; ffnlout => hamk%ffnl_k
479 [ + - ]: 35756 : if (associated(hamk%ph3d_kp)) ph3din => hamk%ph3d_kp
480 [ + - ]: 35756 : if (associated(hamk%ph3d_k )) ph3dout => hamk%ph3d_k
481 [ + - - + ]: 35756 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_kp).and.associated(hamk%ph3d_k)))
482 : 35756 : istwf_k=hamk%istwf_kp
483 [ + + ]: 41297205 : else if (select_k_==K_H_K) then
484 : : ! ===== <k|Vnl|k> =====
485 [ + + + + ]: 140109844 : kptin = hamk%kpt_k ; kptout = hamk%kpt_k
486 : 20015692 : npwin=hamk%npw_fft_k ; npwout=hamk%npw_fft_k
487 : 20015692 : kgin => hamk%kg_k ; kgout => hamk%kg_k
488 [ + + ]: 20015692 : if (associated(hamk%kpg_k)) then
489 : 19765130 : kpgin => hamk%kpg_k ; nkpgin=size(kpgin,2)
490 : : end if
491 [ + + ]: 20015692 : if (associated(hamk%kpg_k)) then
492 : 19765130 : kpgout => hamk%kpg_k ; nkpgout=size(kpgout,2)
493 : : end if
494 : 20015692 : phkxredin => hamk%phkxred ; phkxredout => hamk%phkxred
495 : 20015692 : ffnlin => hamk%ffnl_k ; ffnlout => hamk%ffnl_k
496 [ + + ]: 20015692 : if (associated(hamk%ph3d_k)) ph3din => hamk%ph3d_k
497 [ + + ]: 20015692 : if (associated(hamk%ph3d_k)) ph3dout => hamk%ph3d_k
498 : 20015692 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_k)))
499 : 20015692 : istwf_k=hamk%istwf_k
500 [ + - ]: 21281513 : else if (select_k_==KPRIME_H_KPRIME) then
501 : : ! ===== <k^prime|Vnl|k^prime> =====
502 [ + + + + ]: 148970591 : kptin = hamk%kpt_kp ; kptout = hamk%kpt_kp
503 : 21281513 : npwin=hamk%npw_fft_kp ; npwout=hamk%npw_fft_kp
504 : 21281513 : kgin => hamk%kg_kp ; kgout => hamk%kg_kp
505 [ + + ]: 21281513 : if (associated(hamk%kpg_kp)) then
506 : 21280865 : kpgin => hamk%kpg_kp ; nkpgin=size(kpgin,2)
507 : : end if
508 [ + + ]: 21281513 : if (associated(hamk%kpg_kp)) then
509 : 21280865 : kpgout => hamk%kpg_kp ; nkpgout=size(kpgout,2)
510 : : end if
511 : 21281513 : phkxredin => hamk%phkpxred ; phkxredout => hamk%phkpxred
512 : 21281513 : ffnlin => hamk%ffnl_kp ; ffnlout => hamk%ffnl_kp
513 [ + - ]: 21281513 : if (associated(hamk%ph3d_kp)) ph3din => hamk%ph3d_kp
514 [ + - ]: 21281513 : if (associated(hamk%ph3d_kp)) ph3dout => hamk%ph3d_kp
515 : 21281513 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_kp)))
516 : 21281513 : istwf_k=hamk%istwf_kp
517 : : end if
518 : :
519 [ - + ]: 51196021 : if (npwin==0.or.npwout==0) return
520 : 51196021 : dimffnlin=size(ffnlin,2);dimffnlout=size(ffnlout,2)
521 : 51196021 : kpgin_allocated=(.not.associated(kpgin))
522 [ + + ]: 51196021 : if (kpgin_allocated) then
523 [ - + + - : 503188 : ABI_MALLOC(kpgin,(npwin,0))
- + - + ]
524 : : end if
525 : 51196021 : kpgout_allocated=(.not.associated(kpgout))
526 [ + + ]: 51196021 : if (kpgout_allocated) then
527 [ - + + - : 503188 : ABI_MALLOC(kpgout,(npwout,0))
- + - + ]
528 : : end if
529 : :
530 : : !Check some sizes for safety
531 : : !if (paw_opt==0.or.cpopt<2.or.((cpopt==2.or.cpopt==3).and.choice>1)) then
532 : : !if (size(ffnlin,1)/=npwin.or.size(ffnlin,3)/=hamk%lmnmax) then
533 : : ! ABI_BUG('Incorrect size for ffnlin!')
534 : : !end if
535 [ + + ]: 51196021 : if(signs==2) then
536 [ + - - + ]: 49353960 : if (size(ffnlout,1)/=npwout.or.size(ffnlout,3)/=hamk%lmnmax) then
537 : 0 : ABI_BUG('Incorrect size for ffnlout!')
538 : : end if
539 : : end if
540 [ + - ]: 51196021 : if (associated(kpgin)) then
541 [ + + + + : 153588063 : if (size(kpgin) > 0 .and. size(kpgin,1)/=npwin) then
- + ]
542 : 0 : ABI_BUG('Incorrect size for kpgin')
543 : : end if
544 : : end if
545 [ + - ]: 51196021 : if (associated(kpgout)) then
546 [ + + + + : 153588063 : if (size(kpgout) > 0 .and. size(kpgout,1)/=npwout) then
- + ]
547 : 0 : ABI_BUG('Incorrect size for kpgout')
548 : : end if
549 : : end if
550 : : !This test is OK only because explicit sizes are passed to nonlop_* routines
551 [ + + + + ]: 153588063 : if (size(vectin)<2*npwin*my_nspinor*ndat) then
552 : : !FB: Allow the usage of nonlop from the "linalg" representation where
553 : : !FB: the cg are distributed over the plane waves with npband > 1
554 : : !FB: in case signs=1 & choice=1
555 [ + - + - ]: 51200 : if (signs==1 .and. choice==1) then
556 : 51200 : npwin = size(vectin,2)/ndat/my_nspinor
557 : : else
558 : 0 : ABI_BUG('Incorrect size for vectin!')
559 : : end if
560 : : end if
561 [ + + + + ]: 51196021 : if(choice/=0.and.signs==2) then
562 [ + + ]: 49175298 : if(paw_opt/=3) then
563 : : ! This test is OK only because explicit sizes are passed to nonlop_* routines
564 [ + + - + : 143545212 : ABI_CHECK_IGEQ(size(vectout), 2*npwout*my_nspinor*ndat, 'Incorrect size for vectout!')
- - ]
565 : : end if
566 [ + + ]: 49175298 : if(paw_opt>=3) then
567 [ + + - + : 31678446 : ABI_CHECK_IGEQ(size(svectout), 2*npwout*my_nspinor*ndat, 'Incorrect size for svectout!')
- - ]
568 : : end if
569 : : end if
570 [ + + + + ]: 51196021 : if(cpopt>=0 .and. .not. present(vectproj)) then
571 [ + + - + : 26369340 : ABI_CHECK_IGEQ(size(cprjin), hamk%natom*my_nspinor*ndat, 'Incorrect size for cprjin!')
- - ]
572 : : end if
573 : 51196021 : ndat_left_ = 1
574 [ + + ]: 51196021 : if (present(ndat_left)) then
575 : 220960 : ndat_left_ = ndat_left
576 : : end if
577 [ + + ]: 51196021 : if(present(cprjin_left)) then
578 [ + + - + ]: 662880 : if (size(cprjin_left)/=hamk%natom*my_nspinor*ndat*ndat_left_) then
579 : 0 : ABI_BUG('Incorrect size for cprjin_left!')
580 : : end if
581 : : end if
582 : :
583 : : !Non-local coefficients connecting projectors:
584 : : !If enl is present in the arg list, use it; instead use hamk%ebk
585 [ + + ]: 51196021 : if (present(enl)) then
586 : 3276696 : enl_ptr => enl
587 : 3276696 : dimenl1=size(enl,1);dimenl2=size(enl,2);dimekbq=size(enl,4)
588 : : else
589 : 47919325 : enl_ptr => hamk%ekb
590 : 47919325 : dimenl1=hamk%dimekb1;dimenl2=hamk%dimekb2;dimekbq=1
591 : : end if
592 : :
593 : : ! If enl_ndat is present (meaning enl is absent), use it
594 [ + + ]: 51196021 : if (present(enl_ndat)) then
595 : 35756 : enl_ndat_ptr => enl_ndat
596 : 35756 : dimenl1=size(enl_ndat,1);dimenl2=size(enl_ndat,2);dimekbq=size(enl_ndat,5)
597 : : else
598 [ - + ]: 51160265 : ABI_MALLOC(enl_ndat_ptr, (0,0,0,0,0))
599 : : end if
600 : :
601 : :
602 : : !A specific version of nonlop based on BLAS3 can be used
603 : : !But there are several restrictions
604 : :
605 : 51196021 : use_gemm_nonlop=.false.
606 [ + + ]: 51196021 : if (gemm_nonlop_use_gemm) then
607 : 13942 : use_gemm_nonlop=.true.
608 [ + + ]: 13942 : if(signs==2) then
609 : : use_gemm_nonlop= ( use_gemm_nonlop .and. &
610 : : & ( paw_opt /= 2 .and. &
611 : : & hamk%useylm /= 0 .and.&
612 : : & ((cpopt < 3 .and. (choice < 1 .or. choice == 7)) .or.&
613 [ + - + - : 11478 : & (choice==1 .or. choice==2 .or. choice==3 .or. choice==5 .or. choice==51))))
+ - + + +
- ]
614 : : !FIXME Derivatives of any kind not handled in CUDA GEMM nonlop
615 [ - + - - ]: 11478 : if(choice > 1 .and. choice/=7 .and. (hamk%gpu_option==ABI_GPU_LEGACY .or. hamk%gpu_option==ABI_GPU_KOKKOS)) use_gemm_nonlop=.false.
616 : : end if
617 [ + + ]: 13942 : if(signs==1) then
618 : : use_gemm_nonlop= ( use_gemm_nonlop .and. hamk%useylm/=0 .and. &
619 : : ! Forces and stress (forstr)
620 : : & ( ((choice >= 1 .and. choice <= 3) .or. choice == 23) ) .or. &
621 : : ! Rho ij
622 : : & choice == 0 .or.&
623 [ + - + - : 2464 : ( (choice == 54 .or. choice == 55 .or. choice == 4 .or. choice==6) ) )
+ + ]
624 : : !FIXME forces and constraints computation not handled in CUDA GEMM nonlop
625 [ + + + - ]: 2464 : if(choice > 0 .and. (hamk%gpu_option==ABI_GPU_LEGACY .or. hamk%gpu_option==ABI_GPU_KOKKOS)) use_gemm_nonlop=.false.
626 : : end if
627 : : end if
628 [ + - + + : 51196021 : if(gemm_nonlop_gpu_option/=hamk%gpu_option .or. force_recompute_ph3d .or. my_nspinor/=hamk%nspinor) use_gemm_nonlop=.false.
+ + ]
629 : :
630 : :
631 : : !In the case of a derivative with respect to an atomic displacement,
632 : : !and if <g|dVnl/dR|c> is required (signs=2), we only need to compute the
633 : : !derivatives of the projectors associated with the displaced atom.
634 [ + + ]: 51196021 : iatom_only_=-1;if (present(iatom_only)) iatom_only_=iatom_only
635 [ + + + + ]: 51196021 : atom_pert=((signs==2).and.(choice==2.or.choice==4.or.choice==22.or.choice==24.or.choice==25.or.choice==54))
636 : 51196021 : proj_shift=0
637 : :
638 [ + + ]: 51196021 : if (iatom_only_>0 .and. atom_pert) then
639 : : ! Handling atomic displacement with GEMM variant.
640 : : ! Arrays are fully passed as argument as when treating all atoms.
641 : : ! An atom offset computed below is passed to gemm_nonlop instead.
642 [ - + ]: 7110385 : if (use_gemm_nonlop) then
643 : 0 : iatm=1; proj_shift=0
644 [ # # ]: 0 : do itypat=1, hamk%ntypat
645 [ # # # # ]: 0 : nlmn=count(hamk%indlmn(3,:,itypat)>0)
646 [ # # ]: 0 : do ia=1,hamk%nattyp(itypat)
647 [ # # ]: 0 : if(iatm/=iatom_only_) then
648 : 0 : proj_shift = proj_shift + nlmn
649 : 0 : iatm = iatm + 1
650 : : end if
651 : : end do
652 [ # # ]: 0 : if(iatm==iatom_only) exit
653 : : end do
654 : : end if
655 : : ! We consider only atom with index iatom_only
656 : 7110385 : iatm=hamk%atindx(iatom_only_);itypat=hamk%typat(iatom_only_)
657 : 7110385 : natom_=1 ; ntypat_=1 ; dimenl2_=1 ; matblk_=1
658 [ + + ]: 28441540 : nloalg_(:)=hamk%nloalg(:)
659 [ - + ]: 7110385 : ABI_MALLOC(atindx1_,(1))
660 [ - + ]: 7110385 : ABI_MALLOC(nattyp_,(1))
661 : 7110385 : atindx1_(1)=1 ; nattyp_(1)=1
662 : : ! Store at the right place the 1d phases
663 : 7110385 : n1=hamk%ngfft(1);n2=hamk%ngfft(2);n3=hamk%ngfft(3)
664 [ - + + - : 21331155 : ABI_MALLOC(ph1d_,(2,(2*n1+1)+(2*n2+1)+(2*n3+1)))
+ - + - -
+ - + ]
665 : 7110385 : shift1=(iatm-1)*(2*n1+1)
666 [ + + + + ]: 700264348 : ph1d_(:,1:2*n1+1)=hamk%ph1d(:,1+shift1:2*n1+1+shift1)
667 : 7110385 : shift2=(iatm-1)*(2*n2+1)+hamk%natom*(2*n1+1)
668 [ + + + + ]: 698814772 : ph1d_(:,1+2*n1+1:2*n2+1+2*n1+1)=hamk%ph1d(:,1+shift2:2*n2+1+shift2)
669 : 7110385 : shift3=(iatm-1)*(2*n3+1)+hamk%natom*(2*n1+1+2*n2+1)
670 [ + + + + ]: 731364052 : ph1d_(:,1+2*n1+1+2*n2+1:2*n3+1+2*n2+1+2*n1+1)=hamk%ph1d(:,1+shift3:2*n3+1+shift3)
671 [ - + ]: 7110385 : ABI_MALLOC(phkxredin_,(2,1))
672 [ - + ]: 7110385 : ABI_MALLOC(phkxredout_,(2,1))
673 [ + - + + : 35551925 : phkxredin_(:,1)=phkxredin(:,iatm)
+ + ]
674 [ + - + + : 35551925 : phkxredout_(:,1)=phkxredout(:,iatm)
+ + ]
675 [ - + + - : 21331155 : ABI_MALLOC(ph3din_,(2,npwin,1))
+ - + - -
+ - + ]
676 [ - + + - : 21331155 : ABI_MALLOC(ph3dout_,(2,npwout,1))
+ - + - -
+ - + ]
677 [ + + - + ]: 7110385 : if (force_recompute_ph3d.or.hamk%matblk<hamk%natom) then
678 : 2606200 : nloalg_(2)=-abs(nloalg_(2)) !Will compute the 3D phase factors inside nonlop
679 : : else
680 [ + - + + : 8086866981 : ph3din_(:,1:npwin,1)=ph3din(:,1:npwin,iatm)
+ + + + +
+ ]
681 [ + - + + : 8066921277 : ph3dout_(:,1:npwout,1)=ph3dout(:,1:npwout,iatm)
+ + + + +
+ ]
682 : : end if
683 [ - + + - : 35551925 : ABI_MALLOC(ffnlin_,(npwin,dimffnlin,hamk%lmnmax,1))
- + + - -
+ + - + -
+ - - + -
+ ]
684 [ - + + - : 35551925 : ABI_MALLOC(ffnlout_,(npwout,dimffnlout,hamk%lmnmax,1))
- + + - -
+ + - + -
+ - - + -
+ ]
685 [ + - + + : 24739423013 : ffnlin_(:,:,:,1)=ffnlin(:,:,:,itypat)
+ + + + +
+ + + +
+ ]
686 [ + - + + : 31131811487 : ffnlout_(:,:,:,1)=ffnlout(:,:,:,itypat)
+ + + + +
+ + + +
+ ]
687 [ + + + - : 22754651 : ABI_MALLOC(cprjin_,(1,my_nspinor*ndat*((cpopt+5)/5)))
+ - + + -
+ - + + +
+ + ]
688 [ + + ]: 7110385 : if (cpopt>=0) then
689 : 663340 : nlmn_atm(1)=cprjin(iatm,1)%nlmn
690 : 663340 : ncpgr_atm=cprjin(iatm,1)%ncpgr
691 : 663340 : call pawcprj_alloc(cprjin_,ncpgr_atm,nlmn_atm)
692 [ + + ]: 1375088 : do idat=1,ndat
693 [ + + ]: 2086836 : do ispinor=1,my_nspinor
694 : 711748 : jspinor=ispinor+(idat-1)*my_nspinor
695 : 1423496 : call pawcprj_copy(cprjin(iatm:iatm,jspinor:jspinor),cprjin_(1:1,jspinor:jspinor))
696 : : end do
697 : : end do
698 : : end if
699 [ + + ]: 7110385 : if (present(enl_ndat)) then
700 [ - + + - : 100944 : ABI_MALLOC(enl_ndat_,(size(enl_ndat_ptr,1),1,hamk%nspinor**2,size(enl_ndat_ptr,5),ndat))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
701 [ + + ]: 82056 : do idat=1,ndat
702 [ + + ]: 212520 : do ii=1,size(enl_ndat_ptr,5)
703 [ + + ]: 326160 : do ispden=1,hamk%nspinor**2
704 [ + - + - ]: 260928 : if (dimenl2==hamk%natom .and. hamk%usepaw==1) then
705 [ + - + + : 9523872 : enl_ndat_(:,1,ispden,ii,idat)=enl_ndat_ptr(:,iatom_only_,ispden,idat,ii)
+ + ]
706 [ # # ]: 0 : else if (dimenl2==hamk%ntypat) then
707 [ # # # # : 0 : enl_ndat_(:,1,ispden,ii,idat)=enl_ndat_ptr(:,itypat,ispden,idat,ii)
# # ]
708 : : else
709 [ # # # # : 0 : enl_ndat_(:,1,ispden,ii,idat)=enl_ndat_ptr(:,1,ispden,idat,ii)
# # ]
710 : : end if
711 : : end do
712 : : end do
713 : : end do
714 [ + + + - ]: 35467805 : else if (size(enl_ptr)>0) then
715 [ - + + - : 35467805 : ABI_MALLOC(enl__,(size(enl_ptr,1),1,hamk%nspinor**2,size(enl_ptr,4)))
- + + - -
+ + - + -
+ - - + -
+ ]
716 [ + + ]: 14187122 : do ii=1,size(enl_ptr,4)
717 [ + + ]: 22017117 : do ispden=1,hamk%nspinor**2
718 [ + + + + ]: 14923556 : if (dimenl2==hamk%natom .and. hamk%usepaw==1) then
719 [ + - + + : 76442300 : enl__(:,1,ispden,ii)=enl_ptr(:,iatom_only_,ispden,ii)
+ + ]
720 [ + - ]: 7011063 : else if (dimenl2==hamk%ntypat) then
721 [ + - + + : 79000291 : enl__(:,1,ispden,ii)=enl_ptr(:,itypat,ispden,ii)
+ + ]
722 : : else
723 [ # # # # : 0 : enl__(:,1,ispden,ii)=enl_ptr(:,1,ispden,ii)
# # ]
724 : : end if
725 : : end do
726 : : end do
727 : : else
728 [ # # ]: 0 : ABI_MALLOC(enl__,(0,0,0,0))
729 : : end if
730 [ + - ]: 7110385 : if (allocated(hamk%sij)) then
731 : 7110385 : dimsij=size(hamk%sij,1)
732 [ + + + - : 21331155 : ABI_MALLOC(sij_,(dimsij,1))
+ - + + -
+ - + ]
733 [ + + ]: 7110385 : if (size(hamk%sij,2)==hamk%ntypat) then
734 [ + + ]: 39253104 : sij_(:,1)=hamk%sij(:,itypat)
735 [ + + - + ]: 18823887 : else if (size(hamk%sij)>0) then
736 [ # # ]: 0 : sij_(:,1)=hamk%sij(:,1)
737 : : end if
738 : : end if
739 [ - + + - : 21331155 : ABI_MALLOC(indlmn_,(6,hamk%lmnmax,1))
+ - + - -
+ - + ]
740 [ + + + + ]: 329346257 : indlmn_(:,:,1)=hamk%indlmn(:,:,itypat)
741 : :
742 : : else
743 : : ! Usual case: all atoms are processed
744 : 44085636 : natom_ =hamk%natom; ntypat_=hamk%ntypat
745 : 44085636 : dimenl2_=dimenl2 ; matblk_=hamk%matblk
746 [ + + ]: 176342544 : nloalg_(:) = hamk%nloalg(:)
747 : 44085636 : atindx1_ => hamk%atindx1
748 : 44085636 : nattyp_ => hamk%nattyp
749 : 44085636 : ph1d_ => hamk%ph1d
750 : 44085636 : phkxredin_ => phkxredin
751 : 44085636 : phkxredout_ => phkxredout
752 : 44085636 : ffnlin_ => ffnlin
753 : 44085636 : ffnlout_ => ffnlout
754 : 44085636 : cprjin_ => cprjin
755 : :
756 : 44085636 : enl__ => enl_ptr
757 : :
758 [ + + ]: 44085636 : if (present(enl_ndat)) then
759 [ + - ]: 18932 : if (.not. use_gemm_nonlop) then
760 : : ! An issue with Intel 16 forces to do this conversion
761 [ - + + - : 132524 : ABI_MALLOC(enl_ndat_,(size(enl_ndat_ptr,1),natom_,hamk%nspinor**2,size(enl_ndat_ptr,5),ndat))
- + + - -
+ + - - +
+ - - + +
- + - + -
- + - + ]
762 [ + + ]: 94346 : do idat=1,ndat
763 [ + - + + : 23020202 : enl_ndat_(:,:,:,:,idat)=enl_ndat_ptr(:,:,:,idat,:)
+ + + + +
+ + + + +
+ + + + ]
764 : : end do
765 : : else
766 : 0 : enl_ndat_ => enl_ndat_ptr
767 : : end if
768 : : end if
769 : 44085636 : sij_ => hamk%sij
770 : 44085636 : indlmn_ => hamk%indlmn
771 [ - + ]: 44085636 : if (force_recompute_ph3d) then
772 : 0 : nloalg_(2)=-abs(nloalg_(2)) !Will compute the 3D phase factors inside nonlop
773 [ # # # # : 0 : ABI_MALLOC(ph3din_,(2,npwin,hamk%matblk))
# # # # #
# # # # #
# # ]
774 [ # # # # : 0 : ABI_MALLOC(ph3dout_,(2,npwout,hamk%matblk))
# # # # #
# # # # #
# # ]
775 : : else
776 : 44085636 : ph3din_ => ph3din
777 : 44085636 : ph3dout_ => ph3dout
778 : : end if
779 : :
780 : : end if
781 : :
782 [ + + ]: 51196021 : if(use_gemm_nonlop) then
783 : :
784 [ + - ]: 13942 : if(hamk%gpu_option==ABI_GPU_DISABLED .or. hamk%gpu_option==ABI_GPU_OPENMP) then
785 : :
786 : : ! If forces and stresses are both asked, compute them separately if set to (choice=={2,3})
787 [ + + + - : 13942 : if(choice==23 .and. signs==1 .and. (gemm_nonlop_split_choice23 .or. hamk%gpu_option==ABI_GPU_DISABLED)) then
+ - + - ]
788 : 212 : nnlout_forces = 3*hamk%natom
789 : 212 : nnlout_stress = 6
790 [ - + + - : 636 : ABI_MALLOC(enlout_forces,(nnlout_forces*ndat))
+ - + - -
+ - + ]
791 [ - + + - : 636 : ABI_MALLOC(enlout_stress,(nnlout_stress*ndat))
+ - + - -
+ - + ]
792 : : call gemm_nonlop(hamk%atindx1, 2,cpopt,cprjin,dimenl1,dimenl2,dimekbq,&
793 : : dimffnlin,dimffnlout,enl_ptr,enl_ndat_ptr,enlout_forces,ffnlin,ffnlout,&
794 : : hamk%gmet,hamk%gprimd,&
795 : : idir,hamk%indlmn,istwf_k,kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,&
796 : : hamk%lmnmax,hamk%matblk,hamk%mgfft,mpi_enreg,&
797 : : hamk%natom,hamk%nattyp,ndat,hamk%ngfft,nkpgin,nkpgout,nloalg_,&
798 : : nnlout_forces,npwin,npwout,my_nspinor,hamk%nspinor,hamk%ntypat,only_SO_,paw_opt,&
799 : : ph3din,ph3dout,signs,hamk%sij,svectout,&
800 : : tim_nonlop,hamk%ucvol,hamk%useylm,vectin,vectout,proj_shift,select_k_,&
801 : : iatom_only_,hamk%typat,hamk%usepaw,&
802 [ - + - - : 424 : vectproj=vectproj,gpu_option=hamk%gpu_option)
- - - - -
+ + - - +
- + - + -
+ - + - +
- + - + -
+ - - -
- ]
803 : : call gemm_nonlop(hamk%atindx1, 3,cpopt,cprjin,dimenl1,dimenl2,dimekbq,&
804 : : dimffnlin,dimffnlout,enl_ptr,enl_ndat_ptr,enlout_stress,ffnlin,ffnlout,&
805 : : hamk%gmet,hamk%gprimd,&
806 : : idir,hamk%indlmn,istwf_k,kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,&
807 : : hamk%lmnmax,hamk%matblk,hamk%mgfft,mpi_enreg,&
808 : : hamk%natom,hamk%nattyp,ndat,hamk%ngfft,nkpgin,nkpgout,nloalg_,&
809 : : nnlout_stress,npwin,npwout,my_nspinor,hamk%nspinor,hamk%ntypat,only_SO_,paw_opt,&
810 : : ph3din,ph3dout,signs,hamk%sij,svectout,&
811 : : tim_nonlop,hamk%ucvol,hamk%useylm,vectin,vectout,proj_shift,select_k_,&
812 : : iatom_only_,hamk%typat,hamk%usepaw,&
813 [ - + - - : 424 : vectproj=vectproj,gpu_option=hamk%gpu_option)
- - - - +
- - + - +
- + - + -
+ - + - +
- + - + -
- - - ]
814 [ + + ]: 980 : do idat=1,ndat
815 : 768 : shift = (idat-1)*nnlout
816 : 768 : shift_forces = (idat-1)*nnlout_forces
817 : 768 : shift_stress = (idat-1)*nnlout_stress
818 [ + + ]: 5376 : enlout(shift+1:shift+nnlout_stress) = enlout_stress(shift_stress+1:shift_stress+nnlout_stress)
819 [ + + ]: 13652 : enlout(shift+nnlout_stress+1:shift+nnlout_stress+nnlout_forces) = enlout_forces(shift_forces+1:shift_forces+nnlout_forces)
820 : : end do
821 [ - + ]: 212 : ABI_FREE(enlout_forces)
822 [ - + ]: 212 : ABI_FREE(enlout_stress)
823 : : else
824 : : call gemm_nonlop(hamk%atindx1,choice,cpopt,cprjin,dimenl1,dimenl2,dimekbq,&
825 : : dimffnlin,dimffnlout,enl_ptr,enl_ndat_ptr,enlout,ffnlin,ffnlout,hamk%gmet,hamk%gprimd,&
826 : : idir,hamk%indlmn,istwf_k,kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,&
827 : : hamk%lmnmax,hamk%matblk,hamk%mgfft,mpi_enreg,&
828 : : hamk%natom,hamk%nattyp,ndat,hamk%ngfft,nkpgin,nkpgout,nloalg_,&
829 : : nnlout,npwin,npwout,my_nspinor,hamk%nspinor,hamk%ntypat,only_SO_,paw_opt,&
830 : : ph3din,ph3dout,signs,hamk%sij,svectout,&
831 : : tim_nonlop,hamk%ucvol,hamk%useylm,vectin,vectout,proj_shift,select_k_,&
832 : : iatom_only_,hamk%typat,hamk%usepaw,&
833 [ - + - - : 26308 : vectproj=vectproj,gpu_option=hamk%gpu_option)
- - - + -
- - - - -
+ + + + -
+ - - - +
- + - + -
+ - + - +
- + - + -
+ - - -
- ]
834 : : end if
835 : :
836 : : else if (hamk%gpu_option==ABI_GPU_LEGACY .or. hamk%gpu_option==ABI_GPU_KOKKOS) then
837 : :
838 : : #if defined HAVE_GPU_CUDA
839 : : call gemm_nonlop_gpu(atindx1_, choice, cpopt, cprjin_, dimenl1, dimenl2_, dimekbq, &
840 : : dimffnlin, dimffnlout, &
841 : : enl_ptr, ffnlin, ffnlout, indlmn_, istwf_k, &
842 : : lambda, hamk%lmnmax, matblk_, &
843 : : mpi_enreg, natom_, nattyp_, ndat, nkpgin, nkpgout, &
844 : : nnlout, npwin, npwout, my_nspinor, hamk%nspinor, ntypat_, paw_opt, &
845 : : ph3din, ph3dout, sij_, svectout, &
846 : : hamk%ucvol, hamk%useylm, vectin, vectout, select_k_, &
847 : : hamk%gpu_option,vectproj=vectproj)
848 : : #endif
849 : :
850 : : end if
851 : :
852 : : else
853 : :
854 : : #ifdef HAVE_OPENMP_OFFLOAD
855 : : if(hamk%gpu_option==ABI_GPU_OPENMP) then
856 : : if(xomp_target_is_present(c_loc(vectin))) then
857 : : !$OMP TARGET UPDATE FROM(vectin)
858 : : end if
859 : : end if
860 : : #endif
861 : :
862 : : !$omp parallel do default(shared), &
863 : : !$omp& firstprivate(ndat,npwin,my_nspinor,choice,signs,paw_opt,npwout,cpopt,nnlout,enl__), &
864 : : !$omp& private(b0,b1,b2,b3,b4,e0,e1,e2,e3,e4)
865 : : !!$omp& schedule(static), if(hamk%gpu_option==ABI_GPU_DISABLED)
866 [ + + ]: 106492360 : do idat=1, ndat
867 [ + + ]: 55310281 : if(present(enl_ndat)) enl__ => enl_ndat_(:,:,:,:,idat)
868 : : !vectin_idat => vectin(:,1+npwin*my_nspinor*(idat-1):npwin*my_nspinor*idat)
869 : 55310281 : b0 = 1+npwin*my_nspinor*(idat-1)
870 : 55310281 : e0 = npwin*my_nspinor*idat
871 [ + + + + : 55310281 : if (choice/=0.and.signs==2.and.paw_opt/=3) then
+ + ]
872 : : !vectout_idat => vectout(:,1+npwout*my_nspinor*(idat-1):npwout*my_nspinor*idat)
873 : 49766902 : b1 = 1+npwout*my_nspinor*(idat-1)
874 : 49766902 : e1 = npwout*my_nspinor*idat
875 : : else
876 : : !vectout_idat => vectout
877 : 5543379 : b1 = lbound(vectout,dim=2)
878 [ + + ]: 10796582 : e1 = ubound(vectout,dim=2)
879 : : end if
880 [ + + + + : 55310281 : if (choice/=0.and.signs==2.and.paw_opt>=3) then
+ + ]
881 : : !svectout_idat => svectout(:,1+npwout*my_nspinor*(idat-1):npwout*my_nspinor*idat)
882 : 12246567 : b2 = 1+npwout*my_nspinor*(idat-1)
883 : 12246567 : e2 = npwout*my_nspinor*idat
884 : : else
885 : : !svectout_idat => svectout
886 : 43063714 : b2 = lbound(svectout,dim=2)
887 [ + + ]: 82443476 : e2 = ubound(svectout,dim=2)
888 : : end if
889 : :
890 [ + + ]: 55310281 : if (cpopt>=0) then
891 : : !cprjin_idat => cprjin_(:,my_nspinor*(idat-1)+1:my_nspinor*(idat))
892 : 12082595 : b3 = my_nspinor*(idat-1)+1
893 : 12082595 : e3 = my_nspinor*(idat)
894 : : else
895 : : !cprjin_idat => cprjin_
896 [ + + ]: 43227686 : b3 = lbound(cprjin_,dim=2)
897 [ + + ]: 43227686 : e3 = ubound(cprjin_,dim=2)
898 : : end if
899 [ + + ]: 55310281 : if (nnlout>0) then
900 : : !enlout_idat => enlout((idat-1)*nnlout+1:(idat*nnlout))
901 : 50026612 : b4 = (idat-1)*nnlout*ndat_left_+1
902 : 50026612 : e4 = (idat*nnlout*ndat_left_)
903 : : else
904 : : !enlout_idat => enlout
905 : 5283669 : b4 = lbound(enlout,dim=1)
906 [ + + ]: 7298140 : e4 = ubound(enlout,dim=1)
907 : : end if
908 : :
909 : : ! Legendre Polynomials version
910 [ + + ]: 106492360 : if (hamk%useylm==0) then
911 : : call nonlop_pl(choice,dimenl1,dimenl2_,dimffnlin,dimffnlout,enl__,&
912 : : & enlout(b4:e4),ffnlin_,ffnlout_,hamk%gmet,hamk%gprimd,idir,indlmn_,hamk%ispin_gbt,istwf_k,&
913 : : & kgin,kgout,kpgin,kpgout,kptin,kptout,hamk%lmnmax,matblk_,hamk%mgfft,&
914 : : & mpi_enreg,hamk%mpsang,hamk%mpssoang,natom_,nattyp_,hamk%ngfft,&
915 : : & nkpgin,nkpgout,nloalg_,npwin,npwout,my_nspinor,hamk%nspinor,&
916 : : & ntypat_,only_SO_,phkxredin_,phkxredout_,ph1d_,ph3din_,ph3dout_,signs,hamk%spinaxis,hamk%ucvol,hamk%use_gbt,&
917 [ - + - + : 29595225 : & vectin(:,b0:e0),vectout(:,b1:e1))
- + - + -
+ - + - +
- + - + -
+ ]
918 : : ! Spherical Harmonics version
919 [ + - ]: 25715056 : else if (hamk%gpu_option==ABI_GPU_DISABLED .or. hamk%gpu_option==ABI_GPU_OPENMP) then
920 [ + + ]: 25715056 : if (present(cprjin_left).and.present(enlout_im)) then
921 : : call nonlop_ylm(atindx1_,choice,cpopt,cprjin_(:,b3:e3),dimenl1,dimenl2_,dimekbq,&
922 : : & dimffnlin,dimffnlout,enl__,enlout(b4:e4),ffnlin_,ffnlout_,hamk%gprimd,idir,&
923 : : & indlmn_,istwf_k,kgin,kgout,kpgin,kpgout,kptin,kptout,lambda(idat),&
924 : : & hamk%lmnmax,matblk_,hamk%mgfft,mpi_enreg,natom_,nattyp_,hamk%ngfft,&
925 : : & nkpgin,nkpgout,nloalg_,nnlout,npwin,npwout,my_nspinor,hamk%nspinor,&
926 : : & ntypat_,paw_opt,phkxredin_,phkxredout_,ph1d_,ph3din_,ph3dout_,signs,sij_,&
927 : : & svectout(:,b2:e2),hamk%ucvol,vectin(:,b0:e0),vectout(:,b1:e1),qdir=qdir,&
928 [ - + - + : 220960 : cprjin_left=cprjin_left,enlout_im=enlout_im,ndat_left=ndat_left_)
- + - + -
+ - + - +
- + - + -
+ ]
929 : : else
930 : : call nonlop_ylm(atindx1_,choice,cpopt,cprjin_(:,b3:e3),dimenl1,dimenl2_,dimekbq,&
931 : : & dimffnlin,dimffnlout,enl__,enlout(b4:e4),ffnlin_,ffnlout_,hamk%gprimd,idir,&
932 : : & indlmn_,istwf_k,kgin,kgout,kpgin,kpgout,kptin,kptout,lambda(idat),&
933 : : & hamk%lmnmax,matblk_,hamk%mgfft,mpi_enreg,natom_,nattyp_,hamk%ngfft,&
934 : : & nkpgin,nkpgout,nloalg_,nnlout,npwin,npwout,my_nspinor,hamk%nspinor,&
935 : : & ntypat_,paw_opt,phkxredin_,phkxredout_,ph1d_,ph3din_,ph3dout_,signs,sij_,&
936 [ - + - + : 25494096 : & svectout(:,b2:e2),hamk%ucvol,vectin(:,b0:e0),vectout(:,b1:e1),qdir=qdir)
- + - + -
+ - + - +
- + - + -
+ ]
937 : : end if
938 : : ! GPU version
939 : : else
940 : : call nonlop_gpu(atindx1_,choice,cpopt,cprjin(:,b3:e3),dimenl1,dimenl2_,&
941 : : & dimffnlin,dimffnlout,enl__,enlout(b4:e4),ffnlin_,ffnlout_,hamk%gprimd,idir,&
942 : : & indlmn_,istwf_k,kgin,kgout,kpgin,kpgout,kptin,kptout,lambda(idat),&
943 : : & hamk%lmnmax,matblk_,hamk%mgfft,mpi_enreg,natom_,nattyp_,hamk%ngfft,&
944 : : & nkpgin,nkpgout,nloalg_,nnlout,npwin,npwout,my_nspinor,hamk%nspinor,&
945 : : & ntypat_,paw_opt,phkxredin_,phkxredout_,ph1d_,ph3din_,ph3dout_,signs,sij_,&
946 [ # # # # : 0 : & svectout(:,b2:e2),hamk%ucvol,vectin(:,b0:e0),vectout(:,b1:e1))
# # # # #
# # # # #
# # # # #
# ]
947 : : end if
948 : :
949 : : end do
950 : : !$omp end parallel do
951 : :
952 : : #ifdef HAVE_OPENMP_OFFLOAD
953 : : if(hamk%gpu_option==ABI_GPU_OPENMP) then
954 : : if(signs==2 .and. (paw_opt <= 2 .or. paw_opt == 4) .and. xomp_target_is_present(c_loc(vectout))) then
955 : : !$OMP TARGET UPDATE TO(vectout)
956 : : end if
957 : : if(signs==2 .and. paw_opt>=3 .and. xomp_target_is_present(c_loc(svectout))) then
958 : : !$OMP TARGET UPDATE TO(svectout)
959 : : end if
960 : : end if
961 : : #endif
962 : :
963 : : end if
964 : :
965 : : !Release temporary storage
966 [ + + + - ]: 51196021 : if (iatom_only_>0.and.atom_pert) then
967 [ + + ]: 7110385 : if (cpopt>=0) then
968 : 663340 : call pawcprj_free(cprjin_)
969 : : end if
970 [ - + ]: 7110385 : ABI_FREE(atindx1_)
971 [ - + ]: 7110385 : ABI_FREE(nattyp_)
972 [ - + ]: 7110385 : ABI_FREE(ph1d_)
973 [ - + ]: 7110385 : ABI_FREE(ph3din_)
974 [ - + ]: 7110385 : ABI_FREE(ph3dout_)
975 [ - + ]: 7110385 : ABI_FREE(phkxredin_)
976 [ - + ]: 7110385 : ABI_FREE(phkxredout_)
977 [ - + ]: 7110385 : ABI_FREE(ffnlin_)
978 [ - + ]: 7110385 : ABI_FREE(ffnlout_)
979 [ + + ]: 7110385 : if (present(enl_ndat)) then
980 [ - + ]: 16824 : ABI_FREE(enl_ndat_)
981 : : else
982 [ - + ]: 7093561 : ABI_FREE(enl__)
983 : : end if
984 [ - + ]: 7110385 : ABI_FREE(indlmn_)
985 [ - + ]: 7110385 : ABI_FREE(cprjin_)
986 [ + - ]: 7110385 : if (allocated(hamk%sij)) then
987 [ - + ]: 7110385 : ABI_FREE(sij_)
988 : : end if
989 : : else
990 [ - + ]: 44085636 : if (force_recompute_ph3d) then
991 [ # # ]: 0 : ABI_FREE(ph3din_)
992 [ # # ]: 0 : ABI_FREE(ph3dout_)
993 : : end if
994 : :
995 [ + + ]: 44085636 : if (present(enl_ndat) .and. .not. use_gemm_nonlop) then
996 [ - + ]: 18932 : ABI_FREE(enl_ndat_)
997 : : end if
998 : :
999 : : end if
1000 : :
1001 [ + + ]: 51196021 : if (kpgin_allocated) then
1002 [ - + ]: 251594 : ABI_FREE(kpgin)
1003 : : end if
1004 [ + + ]: 51196021 : if (kpgout_allocated) then
1005 [ - + ]: 251594 : ABI_FREE(kpgout)
1006 : : end if
1007 [ + + ]: 51196021 : if (.not. present(enl_ndat)) then
1008 [ - + ]: 51160265 : ABI_FREE(enl_ndat_ptr)
1009 : : end if
1010 : :
1011 : 51196021 : call timab(220+tim_nonlop,2,tsec)
1012 : :
1013 : : DBG_EXIT("COLL")
1014 : :
1015 : 51196021 : end subroutine nonlop
1016 : : !!***
1017 : :
1018 : : !!****f* ABINIT/nonlop_gpu
1019 : : !! NAME
1020 : : !! nonlop_gpu
1021 : : !!
1022 : : !! FUNCTION
1023 : : !! Compute application of a nonlocal operator, using GPU (NVidia Cuda)
1024 : : !! This routine is an interface to Cuda Kernel gpu_nonlop.cu
1025 : : !!
1026 : : !! INPUTS
1027 : : !! atindx1(natom)=index table for atoms, inverse of atindx
1028 : : !! choice: chooses possible output:
1029 : : !! choice=0 => do nothing (only compute WF projected with NL projectors)
1030 : : !! =1 => a non-local energy contribution
1031 : : !! =2 => a gradient with respect to atomic position(s)
1032 : : !! =3 => a gradient with respect to strain(s)
1033 : : !! =23=> a gradient with respect to atm. pos. and strain(s)
1034 : : !! cpopt=flag defining the status of cprjin%cp(:)=<Proj_i|Cnk> scalars (see below, side effects)
1035 : : !! dimenl1,dimenl2=dimensions of enl (see enl)
1036 : : !! dimffnlin=second dimension of ffnlin (1+number of derivatives)
1037 : : !! dimffnlout=second dimension of ffnlout (1+number of derivatives)
1038 : : !! enl(dimenl1,dimenl2,nspinortot**2)=
1039 : : !! ->Norm conserving : ==== when paw_opt=0 ====
1040 : : !! (Real) Kleinman-Bylander energies (hartree)
1041 : : !! dimenl1=lmnmax - dimenl2=ntypat
1042 : : !! ->PAW : ==== when paw_opt=1, 2 or 4 ====
1043 : : !! (Real or complex, hermitian) Dij coefs to connect projectors
1044 : : !! dimenl1=cplex_enl*lmnmax*(lmnmax+1)/2 - dimenl2=natom
1045 : : !! ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used
1046 : : !! for the application of the nonlocal operator to the |in> vector
1047 : : !! ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used
1048 : : !! for the application of the nonlocal operator to the |out> vector
1049 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
1050 : : !! idir=direction of the - atom to be moved in the case (choice=2,signs=2),
1051 : : !! - k point direction in the case (choice=5,signs=2)
1052 : : !! for choice 53, twisted derivative involves idir+1 and idir+2 (mod 3)
1053 : : !! - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
1054 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
1055 : : !! istwf_k=option parameter that describes the storage of wfs
1056 : : !! kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector
1057 : : !! kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector
1058 : : !! kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector
1059 : : !! kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector
1060 : : !! kptin(3)=k point in terms of recip. translations, for the |in> vector
1061 : : !! kptout(3)=k point in terms of recip. translations, for the |out> vector
1062 : : !! lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
1063 : : !! Typically lambda is the eigenvalue (or its guess)
1064 : : !! lmnmax=max. number of (l,m,n) components over all types of atoms
1065 : : !! matblk=dimension of the arrays ph3din and ph3dout
1066 : : !! mgfft=maximum size of 1D FFTs
1067 : : !! mpi_enreg=information about MPI parallelization
1068 : : !! natom=number of atoms in cell
1069 : : !! nattyp(ntypat)=number of atoms of each type
1070 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
1071 : : !! nkpgin,nkpgout=second sizes of arrays kpgin/kpgout
1072 : : !! nloalg(3)=governs the choice of the algorithm for nonlocal operator
1073 : : !! nnlout=dimension of enlout (when signs=1 and choice>0):
1074 : : !! ==== if paw_opt=0, 1 or 2 ====
1075 : : !! choice=1=>nnlout=1 choice=2=>nnlout=3*natom choice=3=>nnlout=6
1076 : : !! ==== if paw_opt=3 ====
1077 : : !! choice=1 =>nnlout=1
1078 : : !! ==== if paw_opt=4 ====
1079 : : !! not available
1080 : : !! npwin=number of planewaves for given k point, for the |in> vector
1081 : : !! npwout=number of planewaves for given k point, for the |out> vector
1082 : : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
1083 : : !! nspinortot=number of spinorial components of the wavefunctions on current proc
1084 : : !! ntypat=number of types of atoms in cell
1085 : : !! paw_opt= define the nonlocal operator concerned with:
1086 : : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
1087 : : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
1088 : : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
1089 : : !! paw_opt=3 : PAW overlap matrix (Sij)
1090 : : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
1091 : : !! phkxredin(2,natom)=phase factors exp(2 pi kptin.xred)
1092 : : !! phkxredout(2,natom)=phase factors exp(2 pi kptout.xred)
1093 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
1094 : : !! ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in)
1095 : : !! ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out)
1096 : : !! signs= if 1, get contracted elements (energy, forces, stress, ...)
1097 : : !! if 2, applies the non-local operator to a function in reciprocal space
1098 : : !! sij(dimenl1,ntypat*(paw_opt/3))=overlap matrix components (only if paw_opt=2, 3 or 4)
1099 : : !! ucvol=unit cell volume (bohr^3)
1100 : : !! vectin(2,npwin*nspinor)=input cmplx wavefunction coefficients <G|Cnk>
1101 : : !! [cprjin_left(natom,nspinor)]=The projected input wave function <p_nlm|in_left>
1102 : : !! for the left wavefunction. Data are assumed to be in memory, they are NOT recalculated here.
1103 : : !! Only signs==1 and choice==1 are supported.
1104 : : !!
1105 : : !! OUTPUT
1106 : : !! ==== if (signs==1) ====
1107 : : !! --If (paw_opt==0, 1 or 2)
1108 : : !! enlout(nnlout)= contribution to the non-local part of the following properties:
1109 : : !! if choice=1 : enlout(1) -> the energy
1110 : : !! if choice=2 : enlout(1:3*natom) -> the forces
1111 : : !! if choice=3 : enlout(1:6) -> the stresses
1112 : : !! if choice=23: enlout(1:6+3*natom) -> the forces and the stresses
1113 : : !! --If (paw_opt==3)
1114 : : !! if choice=1 : enlout(nnlout)= contribution to <c|S|c> (nnlout=1)
1115 : : !! --If (paw_opt==4)
1116 : : !! not available
1117 : : !! ==== if (signs==2) ====
1118 : : !! --if (paw_opt=0, 1 or 4)
1119 : : !! vectout(2,npwout*nspinor)=result of the application of the concerned operator
1120 : : !! or one of its derivatives to the input vect.:
1121 : : !! if (choice=1) <G|V_nonlocal|vect_start>
1122 : : !! if (choice=2) <G|dV_nonlocal/d(atm coord)|vect_start>
1123 : : !! if (choice=3) <G|dV_nonlocal/d(strain)|vect_start>
1124 : : !! if (paw_opt=2)
1125 : : !! vectout(2,npwout*nspinor)=final vector in reciprocal space:
1126 : : !! if (choice=1) <G|V_nonlocal-lambda.(I+S)|vect_start>
1127 : : !! if (choice=2) <G|d[V_nonlocal-lambda.(I+S)]/d(atm coord)|vect_start>
1128 : : !! if (choice=3) <G|d[V_nonlocal-lambda.(I+S)]/d(strain)|vect_start>
1129 : : !! --if (paw_opt=3 or 4)
1130 : : !! svectout(2,npwout*nspinor)=result of the application of Sij (overlap matrix)
1131 : : !! or one of its derivatives to the input vect.:
1132 : : !! if (choice=1) <G|I+S|vect_start>
1133 : : !! if (choice=2) <G|dS/d(atm coord)|vect_start>
1134 : : !! if (choice=3) <G|dS/d(strain)|vect_start>
1135 : : !!
1136 : : !! SIDE EFFECTS
1137 : : !! cprjin(natom,nspinor) <type(pawcprj_type)>=projected input wave function |in> on non-local projectors
1138 : : !! =<p_lmn|in> and derivatives
1139 : : !! Treatment depends on cpopt parameter:
1140 : : !! if cpopt=-1, <p_lmn|in> (and derivatives)
1141 : : !! are computed here (and not saved)
1142 : : !! if cpopt= 0, <p_lmn|in> are computed here and saved
1143 : : !! derivatives are eventually computed but not saved
1144 : : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
1145 : : !! other derivatives are eventually computed but not saved
1146 : : !!
1147 : : !! TODO
1148 : : !! * Implementation for spinorial wave functions (nspinor=2)
1149 : : !! * Implementation for response function (phonons, ddk, elastic tensor, ...)
1150 : : !!
1151 : : !! SOURCE
1152 : :
1153 : :
1154 : 0 : subroutine nonlop_gpu(atindx1,choice,cpopt,cprjin,dimenl1,dimenl2,dimffnlin,dimffnlout,&
1155 : 0 : & enl,enlout,ffnlin,ffnlout,gprimd,idir,indlmn,istwf_k,&
1156 : 0 : & kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,lmnmax,matblk,mgfft,&
1157 : 0 : & mpi_enreg,natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,nnlout,&
1158 : 0 : & npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,phkxredin,phkxredout,ph1d,&
1159 [ # # # # : 0 : & ph3din,ph3dout,signs,sij,svectout,ucvol,vectin,vectout)
# # # # ]
1160 : :
1161 : : !Arguments ------------------------------------
1162 : : !scalars
1163 : : integer,intent(in) :: choice,cpopt,dimenl1,dimenl2,dimffnlin,dimffnlout,idir
1164 : : integer,intent(in) :: istwf_k,lmnmax,matblk,mgfft,natom,nkpgin,nkpgout,nnlout
1165 : : integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,signs
1166 : : real(dp),intent(in) :: lambda,ucvol
1167 : : type(MPI_type),intent(in) :: mpi_enreg
1168 : : !arrays
1169 : : integer,intent(in) :: atindx1(natom),indlmn(6,lmnmax,ntypat),kgin(3,npwin)
1170 : : integer,intent(in) :: kgout(3,npwout),nattyp(ntypat),ngfft(18),nloalg(3)
1171 : : real(dp),intent(in) :: enl(dimenl1,dimenl2,nspinortot**2)
1172 : : real(dp),intent(in) :: ffnlin(npwin,dimffnlin,lmnmax,ntypat)
1173 : : real(dp),intent(in) :: ffnlout(npwout,dimffnlout,lmnmax,ntypat) !,gmet(3,3)
1174 : : real(dp),intent(in) :: gprimd(3,3),kpgin(npwin,nkpgin),kpgout(npwout,nkpgout)
1175 : : real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom)
1176 : : real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom)
1177 : : real(dp),intent(in) :: sij(dimenl1,ntypat*((paw_opt+1)/3))
1178 : : real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk)
1179 : : real(dp),intent(inout) :: vectin(:,:)
1180 : : real(dp),intent(out) :: enlout(:)
1181 : : real(dp),intent(out),target :: svectout(:,:)
1182 : : real(dp),intent(out),target :: vectout (:,:)
1183 : : type(pawcprj_type),intent(inout) :: cprjin(:,:)
1184 : :
1185 : : !Local variables-------------------------------
1186 : : !scalars
1187 : : integer :: ia,iatom,ilmn,iproj,ispinor,itypat,signs_
1188 : : real(dp) :: doti
1189 : : character(len=500) :: msg
1190 : : !arrays
1191 : 0 : real(dp),allocatable :: proj(:,:)
1192 : 0 : real(dp),pointer :: svectout_(:,:),vectout_(:,:)
1193 : : ! **********************************************************************
1194 : :
1195 : : DBG_ENTER("COLL")
1196 : :
1197 : : !Error on bad choice
1198 [ # # # # : 0 : if ((choice<0 .or. (choice>3.and.choice/=7)).and. choice/=23 .and. choice/=24) then
# # ]
1199 : 0 : write(msg,'(a,i0,a)')'Does not presently support this choice=',choice,'.'
1200 : 0 : ABI_BUG(msg)
1201 : : end if
1202 [ # # ]: 0 : if (cpopt<-1.or.cpopt>2) then
1203 : 0 : ABI_BUG('Bad value for cpopt !')
1204 : : end if
1205 [ # # ]: 0 : if (nspinor==2) then
1206 : 0 : ABI_ERROR('nspinor=2 (spinorial WF) not yet allowed !')
1207 : : end if
1208 : :
1209 [ # # ]: 0 : if ((cpopt==0).or.(cpopt==1).or.(cpopt==2)) then
1210 [ # # # # : 0 : ABI_MALLOC(proj,(2,lmnmax*natom))
# # # # #
# # # ]
1211 [ # # # # ]: 0 : proj=zero;
1212 : : end if
1213 : :
1214 : : !Workaround to get choice=1/signs=1 working
1215 [ # # # # ]: 0 : if (choice==1.and.signs==1) then
1216 : 0 : signs_=2
1217 [ # # # # : 0 : ABI_MALLOC(vectout_,(2,npwin*nspinor))
# # # # #
# # # ]
1218 [ # # # # : 0 : ABI_MALLOC(svectout_,(2,npwin*nspinor*(paw_opt/3)))
# # # # #
# # # ]
1219 : : else
1220 : 0 : signs_=signs;vectout_=>vectout;svectout_=>svectout
1221 : : end if
1222 : :
1223 : : !if cpot==2, the projections are already in memory
1224 [ # # ]: 0 : if (cpopt>=2) then
1225 : : iproj=0
1226 [ # # ]: 0 : do ispinor=1,nspinor
1227 : : iatom=0
1228 [ # # ]: 0 : do itypat=1,ntypat
1229 [ # # ]: 0 : do ia=1,nattyp(itypat)
1230 : 0 : iatom=iatom+1
1231 [ # # ]: 0 : do ilmn=1,cprjin(iatom,1)%nlmn
1232 : 0 : iproj=iproj+1
1233 [ # # ]: 0 : proj(:,iproj)=cprjin(iatom,1)%cp(:,ilmn)
1234 : : end do
1235 : : end do
1236 : : end do
1237 : : end do
1238 : : end if
1239 : :
1240 : : #if defined HAVE_GPU_CUDA
1241 : : call gpu_nonlop(atindx1,choice,cpopt,proj,dimenl1,dimenl2,dimffnlin,dimffnlout,&
1242 : : & enl,enlout,ffnlin,ffnlout,gprimd,idir,indlmn,istwf_k,&
1243 : : & kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,lmnmax,matblk,mgfft,&
1244 : : & mpi_enreg%me_g0_fft,natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,nnlout,&
1245 : : & npwin,npwout,nspinor,ntypat,paw_opt,phkxredin,phkxredout,ph1d,&
1246 : : & ph3din,ph3dout,signs_,sij,svectout_,pi,ucvol,vectin,vectout_)
1247 : : #else
1248 : : ABI_UNUSED(nnlout)
1249 : : #endif
1250 : :
1251 [ # # # # ]: 0 : if (choice==1.and.signs==1) then
1252 [ # # ]: 0 : if (paw_opt/=3) then
1253 [ # # # # : 0 : call dotprod_g(enlout(1),doti,istwf_k,npwin*nspinor,1,vectin,vectout_,mpi_enreg%me_g0_fft,mpi_enreg%comm_spinorfft)
# # # # #
# # # ]
1254 : : else
1255 [ # # # # : 0 : call dotprod_g(enlout(1),doti,istwf_k,npwin*nspinor,1,vectin,svectout_,mpi_enreg%me_g0_fft,mpi_enreg%comm_spinorfft)
# # # # #
# # # ]
1256 : : end if
1257 [ # # ]: 0 : ABI_FREE(vectout_)
1258 [ # # ]: 0 : ABI_FREE(svectout_)
1259 : : else
1260 : 0 : nullify(vectout_,svectout_)
1261 : : end if
1262 : :
1263 [ # # ]: 0 : if ((cpopt==0).or.(cpopt==1)) then
1264 : : iproj=0
1265 [ # # ]: 0 : do ispinor=1,nspinor
1266 : : iatom=0
1267 [ # # ]: 0 : do itypat=1,ntypat
1268 [ # # ]: 0 : do ia=1,nattyp(itypat)
1269 [ # # # # ]: 0 : iatom=iatom+1;cprjin(iatom,1)%nlmn=count(indlmn(3,:,itypat)>0)
1270 [ # # ]: 0 : do ilmn=1,cprjin(iatom,1)%nlmn
1271 [ # # ]: 0 : iproj=iproj+1;cprjin(iatom,1)%cp(:,ilmn)=proj(:,iproj)
1272 : : end do
1273 : : end do
1274 : : end do
1275 : : end do
1276 : : end if
1277 : :
1278 [ # # ]: 0 : if (allocated(proj)) then
1279 : 0 : ABI_FREE(proj)
1280 : : end if
1281 : :
1282 : : DBG_EXIT("COLL")
1283 : :
1284 : : !Fake statements to satisfy ABI rules
1285 : : #if ! defined HAVE_GPU_CUDA
1286 : : if (.false.) then
1287 : : write(std_out,*) atindx1,enl,ffnlin,ffnlout,gprimd
1288 : : write(std_out,*) idir,istwf_k,kgin,kgout,kpgin,kpgout
1289 : : write(std_out,*) kptin,kptout,lambda,mpi_enreg%me
1290 : : write(std_out,*) ngfft,nloalg,ph1d,ph3din,ph3dout
1291 : : write(std_out,*) phkxredin,phkxredout,signs,sij
1292 : : write(std_out,*) ucvol,vectin
1293 : : end if
1294 : : #endif
1295 : :
1296 [ # # ]: 0 : end subroutine nonlop_gpu
1297 : : !!***
1298 : :
1299 : : end module m_nonlop
1300 : : !!***
|
