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 101619076 : subroutine nonlop(choice,cpopt,cprjin,enlout,hamk,idir,lambda,mpi_enreg,ndat,nnlout,&
336 50809538 : paw_opt,signs,svectout,tim_nonlop,vectin,vectout,&
337 203238152 : 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 50809538 : integer,pointer :: kgin(:,:),kgout(:,:)
366 50809538 : integer, contiguous, pointer :: atindx1_(:),indlmn_(:,:,:),nattyp_(:)
367 : real(dp) :: tsec(2)
368 50809538 : real(dp), contiguous, pointer :: enl_ptr(:,:,:,:),enl_ndat_ptr(:,:,:,:,:)
369 50809538 : real(dp),pointer :: ffnlin(:,:,:,:),ffnlin_(:,:,:,:),ffnlout(:,:,:,:),ffnlout_(:,:,:,:)
370 50809538 : real(dp),pointer :: kpgin(:,:),kpgout(:,:)
371 : real(dp) :: kptin(3),kptout(3)
372 50809538 : real(dp),pointer :: ph3din(:,:,:),ph3din_(:,:,:),ph3dout(:,:,:),ph3dout_(:,:,:)
373 50809538 : real(dp),pointer :: phkxredin(:,:),phkxredin_(:,:),phkxredout(:,:),phkxredout_(:,:)
374 50809538 : real(dp), contiguous, pointer :: ph1d_(:,:),sij_(:,:)
375 50809538 : real(dp), contiguous, pointer :: enl__(:,:,:,:),enl_ndat_(:,:,:,:,:)
376 50809538 : 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 50809538 : 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 50809538 : call timab(220+tim_nonlop,1,tsec)
388 :
389 : ! Increment global counter
390 : !$OMP MASTER
391 50809538 : nonlop_counter = nonlop_counter + ndat
392 : !$OMP END MASTER
393 :
394 50809538 : only_SO_=0; if (present(only_SO)) only_SO_=only_SO
395 50809538 : my_nspinor=max(1,hamk%nspinor/mpi_enreg%nproc_spinor)
396 :
397 50809538 : force_recompute_ph3d=.false.
398 :
399 : ! Error(s) on incorrect input
400 50809538 : if (hamk%useylm==0) then
401 28937731 : 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 28937731 : if (cpopt/=-1) then
405 0 : ABI_BUG('If useylm=0, ie no PAW, then cpopt/=-1 is not allowed !')
406 : end if
407 28937731 : 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 28937731 : 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 50809538 : 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 50809538 : 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 50809538 : 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 50809538 : 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 50809538 : 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 81931562 : .and. c_associated(c_loc(hamk%kg_k), c_loc(hamk%kg_kp))) then
437 19687514 : if (associated(hamk%ph3d_k).and.associated(hamk%ph3d_kp)) then
438 19667122 : if (c_associated(c_loc(hamk%ph3d_k), c_loc(hamk%ph3d_kp))) then
439 19667122 : select_k_=K_H_K
440 : end if
441 : else
442 20392 : select_k_=K_H_K
443 : end if
444 : end if
445 50809538 : nkpgin=0;nkpgout=0;nullify(kpgin);nullify(kpgout)
446 50809538 : nullify(ph3din);nullify(ph3dout)
447 : !print *, "in nonlop with select_k:", select_k_
448 :
449 50809538 : if (select_k_==KPRIME_H_K) then
450 : ! ===== <k^prime|Vnl|k> =====
451 68614616 : kptin = hamk%kpt_k ; kptout = hamk%kpt_kp
452 9802088 : npwin=hamk%npw_fft_k ; npwout=hamk%npw_fft_kp
453 9802088 : kgin => hamk%kg_k ; kgout => hamk%kg_kp
454 9802088 : if (associated(hamk%kpg_k)) then
455 9802088 : kpgin => hamk%kpg_k ; nkpgin=size(kpgin,2)
456 : end if
457 9802088 : if (associated(hamk%kpg_kp)) then
458 9802088 : kpgout => hamk%kpg_kp ; nkpgout=size(kpgout,2)
459 : end if
460 9802088 : phkxredin => hamk%phkxred ; phkxredout => hamk%phkpxred
461 9802088 : ffnlin => hamk%ffnl_k ; ffnlout => hamk%ffnl_kp
462 9802088 : if (associated(hamk%ph3d_k )) ph3din => hamk%ph3d_k
463 9802088 : if (associated(hamk%ph3d_kp)) ph3dout => hamk%ph3d_kp
464 9802088 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_k).and.associated(hamk%ph3d_kp)))
465 9802088 : istwf_k=hamk%istwf_k
466 41007450 : 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 40971694 : else if (select_k_==K_H_K) then
484 : ! ===== <k|Vnl|k> =====
485 138796616 : kptin = hamk%kpt_k ; kptout = hamk%kpt_k
486 19828088 : npwin=hamk%npw_fft_k ; npwout=hamk%npw_fft_k
487 19828088 : kgin => hamk%kg_k ; kgout => hamk%kg_k
488 19828088 : if (associated(hamk%kpg_k)) then
489 19578558 : kpgin => hamk%kpg_k ; nkpgin=size(kpgin,2)
490 : end if
491 19828088 : if (associated(hamk%kpg_k)) then
492 19578558 : kpgout => hamk%kpg_k ; nkpgout=size(kpgout,2)
493 : end if
494 19828088 : phkxredin => hamk%phkxred ; phkxredout => hamk%phkxred
495 19828088 : ffnlin => hamk%ffnl_k ; ffnlout => hamk%ffnl_k
496 19828088 : if (associated(hamk%ph3d_k)) ph3din => hamk%ph3d_k
497 19828088 : if (associated(hamk%ph3d_k)) ph3dout => hamk%ph3d_k
498 19828088 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_k)))
499 19828088 : istwf_k=hamk%istwf_k
500 21143606 : else if (select_k_==KPRIME_H_KPRIME) then
501 : ! ===== <k^prime|Vnl|k^prime> =====
502 148005242 : kptin = hamk%kpt_kp ; kptout = hamk%kpt_kp
503 21143606 : npwin=hamk%npw_fft_kp ; npwout=hamk%npw_fft_kp
504 21143606 : kgin => hamk%kg_kp ; kgout => hamk%kg_kp
505 21143606 : if (associated(hamk%kpg_kp)) then
506 21143606 : kpgin => hamk%kpg_kp ; nkpgin=size(kpgin,2)
507 : end if
508 21143606 : if (associated(hamk%kpg_kp)) then
509 21143606 : kpgout => hamk%kpg_kp ; nkpgout=size(kpgout,2)
510 : end if
511 21143606 : phkxredin => hamk%phkpxred ; phkxredout => hamk%phkpxred
512 21143606 : ffnlin => hamk%ffnl_kp ; ffnlout => hamk%ffnl_kp
513 21143606 : if (associated(hamk%ph3d_kp)) ph3din => hamk%ph3d_kp
514 21143606 : if (associated(hamk%ph3d_kp)) ph3dout => hamk%ph3d_kp
515 21143606 : force_recompute_ph3d=(.not.(associated(hamk%ph3d_kp)))
516 21143606 : istwf_k=hamk%istwf_kp
517 : end if
518 :
519 50809538 : if (npwin==0.or.npwout==0) return
520 50809538 : dimffnlin=size(ffnlin,2);dimffnlout=size(ffnlout,2)
521 50809538 : kpgin_allocated=(.not.associated(kpgin))
522 50809538 : if (kpgin_allocated) then
523 499060 : ABI_MALLOC(kpgin,(npwin,0))
524 : end if
525 50809538 : kpgout_allocated=(.not.associated(kpgout))
526 50809538 : if (kpgout_allocated) then
527 499060 : 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 50809538 : if(signs==2) then
536 48982968 : 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 50809538 : if (associated(kpgin)) then
541 152428614 : 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 50809538 : if (associated(kpgout)) then
546 152428614 : 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 152428614 : 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 50809538 : if(choice/=0.and.signs==2) then
562 48822162 : if(paw_opt/=3) then
563 : ! This test is OK only because explicit sizes are passed to nonlop_* routines
564 142582896 : ABI_CHECK_IGEQ(size(vectout), 2*npwout*my_nspinor*ndat, 'Incorrect size for vectout!')
565 : end if
566 48822162 : if(paw_opt>=3) then
567 31136925 : ABI_CHECK_IGEQ(size(svectout), 2*npwout*my_nspinor*ndat, 'Incorrect size for svectout!')
568 : end if
569 : end if
570 50809538 : if(cpopt>=0 .and. .not. present(vectproj)) then
571 26011137 : ABI_CHECK_IGEQ(size(cprjin), hamk%natom*my_nspinor*ndat, 'Incorrect size for cprjin!')
572 : end if
573 50809538 : ndat_left_ = 1
574 50809538 : if (present(ndat_left)) then
575 220960 : ndat_left_ = ndat_left
576 : end if
577 50809538 : 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 50809538 : if (present(enl)) then
586 3239256 : enl_ptr => enl
587 3239256 : dimenl1=size(enl,1);dimenl2=size(enl,2);dimekbq=size(enl,4)
588 : else
589 47570282 : enl_ptr => hamk%ekb
590 47570282 : 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 50809538 : 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 50773782 : 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 50809538 : use_gemm_nonlop=.false.
606 50809538 : 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 50809538 : 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 50809538 : iatom_only_=-1;if (present(iatom_only)) iatom_only_=iatom_only
635 50809538 : atom_pert=((signs==2).and.(choice==2.or.choice==4.or.choice==22.or.choice==24.or.choice==25.or.choice==54))
636 50809538 : proj_shift=0
637 :
638 50809538 : 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 7087345 : 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 7087345 : iatm=hamk%atindx(iatom_only_);itypat=hamk%typat(iatom_only_)
657 7087345 : natom_=1 ; ntypat_=1 ; dimenl2_=1 ; matblk_=1
658 28349380 : nloalg_(:)=hamk%nloalg(:)
659 7087345 : ABI_MALLOC(atindx1_,(1))
660 7087345 : ABI_MALLOC(nattyp_,(1))
661 7087345 : atindx1_(1)=1 ; nattyp_(1)=1
662 : ! Store at the right place the 1d phases
663 7087345 : n1=hamk%ngfft(1);n2=hamk%ngfft(2);n3=hamk%ngfft(3)
664 21262035 : ABI_MALLOC(ph1d_,(2,(2*n1+1)+(2*n2+1)+(2*n3+1)))
665 7087345 : shift1=(iatm-1)*(2*n1+1)
666 695748508 : ph1d_(:,1:2*n1+1)=hamk%ph1d(:,1+shift1:2*n1+1+shift1)
667 7087345 : shift2=(iatm-1)*(2*n2+1)+hamk%natom*(2*n1+1)
668 694298932 : ph1d_(:,1+2*n1+1:2*n2+1+2*n1+1)=hamk%ph1d(:,1+shift2:2*n2+1+shift2)
669 7087345 : shift3=(iatm-1)*(2*n3+1)+hamk%natom*(2*n1+1+2*n2+1)
670 726848212 : 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 7087345 : ABI_MALLOC(phkxredin_,(2,1))
672 7087345 : ABI_MALLOC(phkxredout_,(2,1))
673 35436725 : phkxredin_(:,1)=phkxredin(:,iatm)
674 35436725 : phkxredout_(:,1)=phkxredout(:,iatm)
675 21262035 : ABI_MALLOC(ph3din_,(2,npwin,1))
676 21262035 : ABI_MALLOC(ph3dout_,(2,npwout,1))
677 7087345 : 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 8019866661 : ph3din_(:,1:npwin,1)=ph3din(:,1:npwin,iatm)
681 7999920957 : ph3dout_(:,1:npwout,1)=ph3dout(:,1:npwout,iatm)
682 : end if
683 35436725 : ABI_MALLOC(ffnlin_,(npwin,dimffnlin,hamk%lmnmax,1))
684 35436725 : ABI_MALLOC(ffnlout_,(npwout,dimffnlout,hamk%lmnmax,1))
685 24447967013 : ffnlin_(:,:,:,1)=ffnlin(:,:,:,itypat)
686 30840355487 : ffnlout_(:,:,:,1)=ffnlout(:,:,:,itypat)
687 22639451 : ABI_MALLOC(cprjin_,(1,my_nspinor*ndat*((cpopt+5)/5)))
688 7087345 : if (cpopt>=0) then
689 640300 : nlmn_atm(1)=cprjin(iatm,1)%nlmn
690 640300 : ncpgr_atm=cprjin(iatm,1)%ncpgr
691 640300 : call pawcprj_alloc(cprjin_,ncpgr_atm,nlmn_atm)
692 1329008 : do idat=1,ndat
693 2017716 : do ispinor=1,my_nspinor
694 688708 : jspinor=ispinor+(idat-1)*my_nspinor
695 1377416 : call pawcprj_copy(cprjin(iatm:iatm,jspinor:jspinor),cprjin_(1:1,jspinor:jspinor))
696 : end do
697 : end do
698 : end if
699 7087345 : 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 35352605 : else if (size(enl_ptr)>0) then
715 35352605 : ABI_MALLOC(enl__,(size(enl_ptr,1),1,hamk%nspinor**2,size(enl_ptr,4)))
716 14141042 : do ii=1,size(enl_ptr,4)
717 21947997 : do ispden=1,hamk%nspinor**2
718 14877476 : if (dimenl2==hamk%natom .and. hamk%usepaw==1) then
719 72225980 : 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 7087345 : if (allocated(hamk%sij)) then
731 7087345 : dimsij=size(hamk%sij,1)
732 21262035 : ABI_MALLOC(sij_,(dimsij,1))
733 7087345 : if (size(hamk%sij,2)==hamk%ntypat) then
734 37133424 : 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 21262035 : ABI_MALLOC(indlmn_,(6,hamk%lmnmax,1))
740 327226577 : indlmn_(:,:,1)=hamk%indlmn(:,:,itypat)
741 :
742 : else
743 : ! Usual case: all atoms are processed
744 43722193 : natom_ =hamk%natom; ntypat_=hamk%ntypat
745 43722193 : dimenl2_=dimenl2 ; matblk_=hamk%matblk
746 174888772 : nloalg_(:) = hamk%nloalg(:)
747 43722193 : atindx1_ => hamk%atindx1
748 43722193 : nattyp_ => hamk%nattyp
749 43722193 : ph1d_ => hamk%ph1d
750 43722193 : phkxredin_ => phkxredin
751 43722193 : phkxredout_ => phkxredout
752 43722193 : ffnlin_ => ffnlin
753 43722193 : ffnlout_ => ffnlout
754 43722193 : cprjin_ => cprjin
755 :
756 43722193 : enl__ => enl_ptr
757 :
758 43722193 : 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 43722193 : sij_ => hamk%sij
770 43722193 : indlmn_ => hamk%indlmn
771 43722193 : 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 43722193 : ph3din_ => ph3din
777 43722193 : ph3dout_ => ph3dout
778 : end if
779 :
780 : end if
781 :
782 50809538 : 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 105284647 : do idat=1, ndat
867 54489051 : if(present(enl_ndat)) enl__ => enl_ndat_(:,:,:,:,idat)
868 : !vectin_idat => vectin(:,1+npwin*my_nspinor*(idat-1):npwin*my_nspinor*idat)
869 54489051 : b0 = 1+npwin*my_nspinor*(idat-1)
870 54489051 : e0 = npwin*my_nspinor*idat
871 54489051 : 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 49294683 : b1 = 1+npwout*my_nspinor*(idat-1)
874 49294683 : e1 = npwout*my_nspinor*idat
875 : else
876 : !vectout_idat => vectout
877 5194368 : b1 = lbound(vectout,dim=2)
878 10098560 : e1 = ubound(vectout,dim=2)
879 : end if
880 54489051 : 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 11780917 : b2 = 1+npwout*my_nspinor*(idat-1)
883 11780917 : e2 = npwout*my_nspinor*idat
884 : else
885 : !svectout_idat => svectout
886 42708134 : b2 = lbound(svectout,dim=2)
887 81733092 : e2 = ubound(svectout,dim=2)
888 : end if
889 :
890 54489051 : if (cpopt>=0) then
891 : !cprjin_idat => cprjin_(:,my_nspinor*(idat-1)+1:my_nspinor*(idat))
892 11645205 : b3 = my_nspinor*(idat-1)+1
893 11645205 : e3 = my_nspinor*(idat)
894 : else
895 : !cprjin_idat => cprjin_
896 42843846 : b3 = lbound(cprjin_,dim=2)
897 42843846 : e3 = ubound(cprjin_,dim=2)
898 : end if
899 54489051 : if (nnlout>0) then
900 : !enlout_idat => enlout((idat-1)*nnlout+1:(idat*nnlout))
901 49519234 : b4 = (idat-1)*nnlout*ndat_left_+1
902 49519234 : e4 = (idat*nnlout*ndat_left_)
903 : else
904 : !enlout_idat => enlout
905 4969817 : b4 = lbound(enlout,dim=1)
906 6681169 : e4 = ubound(enlout,dim=1)
907 : end if
908 :
909 : ! Legendre Polynomials version
910 105284647 : 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 29457648 : & vectin(:,b0:e0),vectout(:,b1:e1))
918 : ! Spherical Harmonics version
919 25031403 : else if (hamk%gpu_option==ABI_GPU_DISABLED .or. hamk%gpu_option==ABI_GPU_OPENMP) then
920 25031403 : 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 24810443 : & 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 50809538 : if (iatom_only_>0.and.atom_pert) then
967 7087345 : if (cpopt>=0) then
968 640300 : call pawcprj_free(cprjin_)
969 : end if
970 7087345 : ABI_FREE(atindx1_)
971 7087345 : ABI_FREE(nattyp_)
972 7087345 : ABI_FREE(ph1d_)
973 7087345 : ABI_FREE(ph3din_)
974 7087345 : ABI_FREE(ph3dout_)
975 7087345 : ABI_FREE(phkxredin_)
976 7087345 : ABI_FREE(phkxredout_)
977 7087345 : ABI_FREE(ffnlin_)
978 7087345 : ABI_FREE(ffnlout_)
979 7087345 : if (present(enl_ndat)) then
980 16824 : ABI_FREE(enl_ndat_)
981 : else
982 7070521 : ABI_FREE(enl__)
983 : end if
984 7087345 : ABI_FREE(indlmn_)
985 7087345 : ABI_FREE(cprjin_)
986 7087345 : if (allocated(hamk%sij)) then
987 7087345 : ABI_FREE(sij_)
988 : end if
989 : else
990 43722193 : if (force_recompute_ph3d) then
991 0 : ABI_FREE(ph3din_)
992 0 : ABI_FREE(ph3dout_)
993 : end if
994 :
995 43722193 : 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 50809538 : if (kpgin_allocated) then
1002 249530 : ABI_FREE(kpgin)
1003 : end if
1004 50809538 : if (kpgout_allocated) then
1005 249530 : ABI_FREE(kpgout)
1006 : end if
1007 50809538 : if (.not. present(enl_ndat)) then
1008 50773782 : ABI_FREE(enl_ndat_ptr)
1009 : end if
1010 :
1011 50809538 : call timab(220+tim_nonlop,2,tsec)
1012 :
1013 : DBG_EXIT("COLL")
1014 :
1015 50809538 : 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 : !!***
|
