Branch data Line data Source code
1 : : !!****m* ABINIT/m_nonlop_ylm
2 : : !! NAME
3 : : !! m_nonlop_ylm
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !! COPYRIGHT
8 : : !! Copyright (C) 1998-2026 ABINIT group (MT)
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 : : module m_nonlop_ylm
22 : :
23 : : use defs_basis
24 : : use m_xmpi
25 : : use m_abicore
26 : : use m_errors
27 : :
28 : : use defs_abitypes, only : MPI_type
29 : : use m_geometry, only : strconv
30 : : use m_kg, only : ph1d3d, mkkpg
31 : : use m_pawcprj, only : pawcprj_type
32 : : use m_opernla_ylm, only : opernla_ylm,opernla_counter
33 : : use m_opernla_ylm_mv, only : opernla_ylm_mv,opernla_mv_counter,opernla_mv_dgemv_counter
34 : : use m_opernlb_ylm, only : opernlb_ylm,opernlb_counter
35 : : use m_opernlb_ylm_mv, only : opernlb_ylm_mv,opernlb_mv_counter,opernlb_mv_dgemv_counter
36 : : use m_opernlc_ylm, only : opernlc_ylm
37 : : use m_opernld_ylm, only : opernld_ylm
38 : : use m_kg, only : mkkpgcart
39 : : ! use m_time, only : timab
40 : :
41 : : implicit none
42 : :
43 : : private
44 : : !!***
45 : :
46 : : public :: nonlop_ylm
47 : : public :: nonlop_ylm_init_counters
48 : : public :: nonlop_ylm_stop_counters
49 : : public :: nonlop_ylm_output_counters
50 : : !!***
51 : :
52 : : contains
53 : : !!***
54 : :
55 : : !!****f* ABINIT/nonlop_ylm
56 : : !! NAME
57 : : !! nonlop_ylm
58 : : !!
59 : : !! FUNCTION
60 : : !! * Compute application of a nonlocal operator Vnl in order to get:
61 : : !! - contracted elements (energy, forces, stresses, ...), if signs=1
62 : : !! - a function in reciprocal space (|out> = Vnl|in>), if signs=2
63 : : !! Operator Vnl, as the following general form:
64 : : !! $Vnl=sum_{R,lmn,l''m''n''} {|P_{Rlmn}> Enl^{R}_{lmn,l''m''n''} <P_{Rl''m''n''}|}$
65 : : !! Operator Vnl is -- in the typical case -- the nonlocal potential.
66 : : !! - With norm-conserving pseudopots, $Enl^{R}_{lmn,l''m''n''}$ is the
67 : : !! Kleinmann-Bylander energy $Ekb^{R}_{ln}$.
68 : : !! - In a PAW calculation, $Enl^{R}_{lmn,l''m''n''}$ are the nonlocal
69 : : !! coefficients to connect projectors $D_{ij}$.
70 : : !! - The |P_{Rlmn}> are the projector functions.
71 : : !! * Optionnaly, in case of PAW calculation, compute:
72 : : !! - Application of the overlap matrix in reciprocal space
73 : : !! (<in|S|in> or (I+S)|in>).
74 : : !! - Application of (Vnl-lambda.S) in reciprocal space
75 : : !! (<in|Vnl-lambda.S|in> and derivatives or (Vnl-lambda.S)|in>).
76 : : !! * This routine uses spherical harmonics Ylm to express Vnl.
77 : : !!
78 : : !! INPUTS
79 : : !! atindx1(natom)=index table for atoms, inverse of atindx
80 : : !! choice: chooses possible output:
81 : : !! choice=0 => do nothing (only compute WF projected with NL projectors)
82 : : !! =1 => non-local energy contribution
83 : : !! =2 => 1st derivative(s) with respect to atomic position(s)
84 : : !! =3 => 1st derivative(s) with respect to strain(s)
85 : : !! =22=> mixed 2nd derivative(s) with respect to atomic pos. and q vector (at q=0)
86 : : !! =25=> mixed 3rd derivative(s) with respect to atomic pos. and two q vectors (at q=0)
87 : : !! =23=> 1st derivative(s) with respect to atomic pos. and
88 : : !! 1st derivative(s) with respect to atomic pos. and strains
89 : : !! =4 => 2nd derivative(s) with respect to 2 atomic pos.
90 : : !! =24=> 1st derivative(s) with respect to atm. pos. and
91 : : !! =33=> mixed 2nd derivative(s) with respect to strain and q vector (at q=0)
92 : : !! 2nd derivative(s) with respect to 2 atomic pos.
93 : : !! =5 => 1st derivative(s) with respect to k wavevector, typically
94 : : !! sum_ij [ |p_i> D_ij <dp_j/dk| + |dp_i/dk> D_ij < p_j| ]
95 : : !! =6 => 2nd derivative(s) with respect to 2 strains and
96 : : !! mixed 2nd derivative(s) with respect to strains & atomic pos.
97 : : !! =51 =>right 1st derivative(s) with respect to k wavevector, typically
98 : : !! sum_ij [ |p_i> D_ij <dp_j/dk| ]
99 : : !! =52 =>left 1st derivative(s) with respect to k wavevector, typically
100 : : !! sum_ij [ |dp_i/dk> D_ij < p_j| ]
101 : : !! =53 =>twist 1st derivative(s) with respect to k, typically
102 : : !! sum_ij [ |dp_i/dk_(idir+1)> D_ij <dp_j//dk_(idir+2)|
103 : : !! =54=> mixed 2nd derivative(s) with respect to atomic pos. and left k wavevector
104 : : !! =55=> mixed 2nd derivative(s) with respect to strain and right k wavevector
105 : : !! =7 => apply operator $\sum_i [ |p_i> <p_i| ],
106 : : !! same as overlap operator with s_ij=identity (paw_opt==3 only)
107 : : !! =8 => 2nd derivatives with respect to 2 k wavevectors
108 : : !! =81=> partial 2nd derivatives with respect to 2 k wavevectors,
109 : : !! full derivative with respect to k1, right derivative with respect to k2,
110 : : !! (derivative with respect to k of choice 51), typically
111 : : !! sum_ij [ |dp_i/dk1> D_ij <dp_j/dk2| + |p_i> D_ij < d2p_j/dk1dk2| ]
112 : : !! Only choices 1,2,3,23,4,5,6 are compatible with useylm=0.
113 : : !! Only choices 1,2,22,25,3,5,33,51,52,53,7,8,81 are compatible with signs=2
114 : : !! cpopt=flag defining the status of cprjin%cp(:)=<Proj_i|Cnk> scalars (see below, side effects)
115 : : !! dimenl1,dimenl2=dimensions of enl (see enl)
116 : : !! dimekbq=1 if enl factors do not contain a exp(-iqR) phase, 2 is they do
117 : : !! dimffnlin=second dimension of ffnlin (1+number of derivatives)
118 : : !! dimffnlout=second dimension of ffnlout (1+number of derivatives)
119 : : !! enl(cplex_enl*dimenl1,dimenl2,nspinortot**2,dimekbq)=
120 : : !! ->Norm conserving : ==== when paw_opt=0 ====
121 : : !! (Real) Kleinman-Bylander energies (hartree)
122 : : !! dimenl1=lmnmax - dimenl2=ntypat
123 : : !! dimekbq is 2 if Enl contains a exp(-iqR) phase, 1 otherwise
124 : : !! ->PAW : ==== when paw_opt=1, 2 or 4 ====
125 : : !! (Real or complex, hermitian) Dij coefs to connect projectors
126 : : !! dimenl1=cplex_enl*lmnmax*(lmnmax+1)/2 - dimenl2=natom
127 : : !! These are complex numbers if cplex_enl=2
128 : : !! enl(:,:,1) contains Dij^up-up
129 : : !! enl(:,:,2) contains Dij^dn-dn
130 : : !! enl(:,:,3) contains Dij^up-dn (only if nspinor=2)
131 : : !! enl(:,:,4) contains Dij^dn-up (only if nspinor=2)
132 : : !! dimekbq is 2 if Dij contains a exp(-iqR) phase, 1 otherwise
133 : : !! ffnlin(npwin,dimffnlin,lmnmax,ntypat)=nonlocal form factors to be used
134 : : !! for the application of the nonlocal operator to the |in> vector
135 : : !! ffnlout(npwout,dimffnlout,lmnmax,ntypat)=nonlocal form factors to be used
136 : : !! for the application of the nonlocal operator to the |out> vector
137 : : !! -----------------------------------------------------------
138 : : !! gprimd(3,3)=dimensional reciprocal space primitive translations
139 : : !! idir=direction of the - atom to be moved in the case (choice=2,signs=2) or (choice=22,signs=2)
140 : : !! - k point direction in the case (choice=5,signs=2S)
141 : : !! for choice 53, twisted derivative involves idir+1 and idir+2 (mod 3)
142 : : !! - strain component (1:6) in the case (choice=3,signs=2) or (choice=6,signs=1)
143 : : !! - strain component (1:9) in the case (choice=33,signs=2)
144 : : !! - (1:9) components to specify the atom to be moved and the second q-gradient
145 : : !! direction in the case (choice=25,signs=2)
146 : : !! indlmn(6,i,ntypat)= array giving l,m,n,lm,ln,s for i=lmn
147 : : !! istwf_k=option parameter that describes the storage of wfs
148 : : !! kgin(3,npwin)=integer coords of planewaves in basis sphere, for the |in> vector
149 : : !! kgout(3,npwout)=integer coords of planewaves in basis sphere, for the |out> vector
150 : : !! kpgin(npw,npkgin)= (k+G) components and related data, for the |in> vector
151 : : !! kpgout(npw,nkpgout)=(k+G) components and related data, for the |out> vector
152 : : !! kptin(3)=k point in terms of recip. translations, for the |in> vector
153 : : !! kptout(3)=k point in terms of recip. translations, for the |out> vector
154 : : !! lambda=factor to be used when computing (Vln-lambda.S) - only for paw_opt=2
155 : : !! Typically lambda is the eigenvalue (or its guess)
156 : : !! lmnmax=max. number of (l,m,n) components over all types of atoms
157 : : !! matblk=dimension of the arrays ph3din and ph3dout
158 : : !! mgfft=maximum size of 1D FFTs
159 : : !! mpi_enreg=information about MPI parallelization
160 : : !! natom=number of atoms in cell
161 : : !! nattyp(ntypat)=number of atoms of each type
162 : : !! ngfft(18)=contain all needed information about 3D FFT, see ~abinit/doc/variables/vargs.htm#ngfft
163 : : !! nkpgin,nkpgout=second sizes of arrays kpgin/kpgout
164 : : !! nloalg(3)=governs the choice of the algorithm for nonlocal operator
165 : : !! nnlout=dimension of enlout (when signs=1 and choice>0):
166 : : !! ==== if paw_opt=0, 1 or 2 ====
167 : : !! choice nnlout | choice nnlout
168 : : !! 1 1 | 51 6 (complex)
169 : : !! 2 3*natom | 52 6 (complex)
170 : : !! 3 6 | 53 6 (complex)
171 : : !! 4 6*natom | 54 9*natom
172 : : !! 23 6+3*natom | 55 36 (complex)
173 : : !! 24 9*natom | 6 36+18*natom
174 : : !! 5 3 | 8 6
175 : : !! | 81 18 (complex)
176 : : !! ==== if paw_opt=3 ====
177 : : !! choice nnlout
178 : : !! 1 1
179 : : !! 2 3*natom
180 : : !! 5 3
181 : : !! 51 3
182 : : !! 52 3
183 : : !! 54 9*natom
184 : : !! 55 36
185 : : !! 7 1
186 : : !! 8 6
187 : : !! 81 9
188 : : !! ==== if paw_opt=4 ====
189 : : !! not available
190 : : !! npwin=number of planewaves for given k point, for the |in> vector
191 : : !! npwout=number of planewaves for given k point, for the |out> vector
192 : : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
193 : : !! nspinortot=total number of spinorial components of the wavefunctions
194 : : !! ntypat=number of types of atoms in cell
195 : : !! paw_opt= define the nonlocal operator concerned with:
196 : : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
197 : : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
198 : : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
199 : : !! paw_opt=3 : PAW overlap matrix (Sij)
200 : : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
201 : : !! phkxredin(2,natom)=phase factors exp(2 pi kptin.xred)
202 : : !! phkxredout(2,natom)=phase factors exp(2 pi kptout.xred)
203 : : !! ph1d(2,3*(2*mgfft+1)*natom)=1D structure factors phase information
204 : : !! ph3din(2,npwin,matblk)=3D structure factors, for each atom and plane wave (in)
205 : : !! ph3dout(2,npwout,matblk)=3-dim structure factors, for each atom and plane wave (out)
206 : : !! [qdir]= optional,direction of the q-gradient (only for choice=22 choice=25 and choice=33)
207 : : !! signs= if 1, get contracted elements (energy, forces, stress, ...)
208 : : !! if 2, applies the non-local operator to a function in reciprocal space
209 : : !! sij(dimenl1,ntypat*(paw_opt/3))=overlap matrix components (only if paw_opt=2, 3 or 4)
210 : : !! ucvol=unit cell volume (bohr^3)
211 : : !! vectin(2,npwin*nspinor)=input cmplx wavefunction coefficients <G|in>
212 : : !! [cprjin_left(natom,nspinor)]=The projected input wave function <p_nlm|in_left>
213 : : !! for the left wavefunction. Data are assumed to be in memory, they are NOT recalculated here.
214 : : !!
215 : : !! OUTPUT
216 : : !! ==== if (signs==1) ====
217 : : !! --If (paw_opt==0, 1 or 2)
218 : : !! enlout(nnlout)= contribution to the non-local part of the following properties:
219 : : !! if choice=1 : enlout(1) -> the energy
220 : : !! if choice=2 : enlout(3*natom) -> 1st deriv. of energy wrt atm. pos (forces)
221 : : !! if choice=3 : enlout(6) -> 1st deriv. of energy wrt strain (stresses)
222 : : !! if choice=4 : enlout(6*natom) -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
223 : : !! if choice=23: enlout(6+3*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
224 : : !! 1st deriv. of energy wrt strain (stresses)
225 : : !! if choice=24: enlout(9*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
226 : : !! 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
227 : : !! if choice=5 : enlout(3) -> 1st deriv. of energy wrt k
228 : : !! if choice=51: enlout(3) -> 1st deriv. (right) of energy wrt k
229 : : !! if choice=52: enlout(3) -> 1st deriv. (left) of energy wrt k
230 : : !! if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
231 : : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
232 : : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
233 : : !! if choice=6 : enlout(36+18*natom) -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
234 : : !! 2nd deriv. of energy wrt to atm. pos and strain (internal strain)
235 : : !! if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
236 : : !! if choice=81: enlout(9) -> 2nd deriv.of E: full derivative w.r.t. k1, right derivative w.r.t k2
237 : : !! --If (paw_opt==3)
238 : : !! if choice=1 : enlout(1) -> contribution to <c|S|c> (note: not including <c|c>)
239 : : !! if choice=2 : enlout(3*natom) -> contribution to <c|dS/d_atm.pos|c>
240 : : !! if choice=51: enlout(3) -> contribution to <c|d(right)S/d_k|c>
241 : : !! if choice=52: enlout(3) -> contribution to <c|d(left)S/d_k|c>
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=7 : enlout(1) -> contribution to <c|sum_i[p_i><p_i]|c>
245 : : !! if choice=8 : enlout(6) -> contribution to <c|d2S/d_k1d_k2|c>
246 : : !! if choice=81: enlout(9) -> contribution to <c|dS/d_k1[d(right)d_k2]|c>
247 : : !! --If (paw_opt==4)
248 : : !! not available
249 : : !! ==== if (signs==2) ====
250 : : !! --if (paw_opt=0)
251 : : !! vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
252 : : !! or one of its derivatives to the input vect.
253 : : !! if (choice=22) <G|d2V_nonlocal/d(atm. pos)dq|vect_in> (at q=0)
254 : : !! if (choice=25) <G|d3V_nonlocal/d(atm. pos)dqdq|vect_in> (at q=0)
255 : : !! if (choice=33) <G|d2V_nonlocal/d(strain)dq|vect_in> (at q=0)
256 : : !! --if (paw_opt=0, 1 or 4)
257 : : !! vectout(2,npwout*my_nspinor*ndat)=result of the aplication of the concerned operator
258 : : !! or one of its derivatives to the input vect.:
259 : : !! if (choice=1) <G|V_nonlocal|vect_in>
260 : : !! if (choice=2) <G|dV_nonlocal/d(atm. pos)|vect_in>
261 : : !! if (choice=3) <G|dV_nonlocal/d(strain)|vect_in>
262 : : !! if (choice=5) <G|dV_nonlocal/d(k)|vect_in>
263 : : !! if (choice=51) <G|d(right)V_nonlocal/d(k)|vect_in>
264 : : !! if (choice=52) <G|d(left)V_nonlocal/d(k)|vect_in>
265 : : !! if (choice=53) <G|d(twist)V_nonlocal/d(k)|vect_in>
266 : : !! if (choice=8) <G|d2V_nonlocal/d(k)d(k)|vect_in>
267 : : !! if (choice=81) <G|d[d(right)V_nonlocal/d(k)]/d(k)|vect_in>
268 : : !! --if (paw_opt=2)
269 : : !! vectout(2,npwout*my_nspinor*ndat)=final vector in reciprocal space:
270 : : !! if (choice=1) <G|V_nonlocal-lambda.(I+S)|vect_in>
271 : : !! if (choice=2) <G|d[V_nonlocal-lambda.(I+S)]/d(atm. pos)|vect_in>
272 : : !! if (choice=3) <G|d[V_nonlocal-lambda.(I+S)]/d(strain)|vect_in>
273 : : !! if (choice=5) <G|d[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
274 : : !! if (choice=51) <G|d(right)[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
275 : : !! if (choice=52) <G|d(left)[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
276 : : !! if (choice=53) <G|d(twist)[V_nonlocal-lambda.(I+S)]/d(k)|vect_in>
277 : : !! if (choice=8) <G|d2[V_nonlocal-lambda.(I+S)]/d(k)d(k)|vect_in>
278 : : !! if (choice=81) <G|d[d(right[V_nonlocal-lambda.(I+S)]/d(k)]/d(k)|vect_in>
279 : : !! --if (paw_opt=3 or 4)
280 : : !! svectout(2,npwout*my_nspinor*ndat)=result of the aplication of Sij (overlap matrix)
281 : : !! or one of its derivatives to the input vect.:
282 : : !! if (choice=1) <G|I+S|vect_in>
283 : : !! if (choice=2) <G|dS/d(atm. pos)|vect_in>
284 : : !! if (choice=3) <G|dS/d(strain)|vect_in>
285 : : !! if (choice=5) <G|dS/d(k)|vect_in>
286 : : !! if (choice=51) <G|d(right)S/d(k)|vect_in>
287 : : !! if (choice=52) <G|d(left)S/d(k)|vect_in>
288 : : !! if (choice=53) <G|d(twist)S/d(k)|vect_in>
289 : : !! if (choice=3) <G|d[V_nonlocal-lambda.(I+S)]/d(strain)|vect_in>
290 : : !! if (choice=7) <G|sum_i[p_i><p_i]|vect_in>
291 : : !! if (choice=8) <G|d2S/d(k)d(k)|vect_in>
292 : : !! if (choice=81) <G|d[d(right)S/d(k)]/d(k)|vect_in>
293 : : !!
294 : : !! SIDE EFFECTS
295 : : !! cprjin(natom,nspinor) <type(pawcprj_type)>=projected input wave function |in> on non-local projectors
296 : : !! =<p_lmn|in> and derivatives
297 : : !! Treatment depends on cpopt parameter:
298 : : !! if cpopt=-1, <p_lmn|in> (and derivatives)
299 : : !! are computed here (and not saved)
300 : : !! if cpopt= 0, <p_lmn|in> are computed here and saved
301 : : !! derivatives are eventually computed but not saved
302 : : !! if cpopt= 1, <p_lmn|in> and first derivatives are computed here and saved
303 : : !! other derivatives are eventually computed but not saved
304 : : !! if cpopt= 2 <p_lmn|in> are already in memory;
305 : : !! first (and 2nd) derivatives are computed here and not saved
306 : : !! if cpopt= 3 <p_lmn|in> are already in memory;
307 : : !! first derivatives are computed here and saved
308 : : !! other derivatives are eventually computed but not saved
309 : : !! if cpopt= 4 <p_lmn|in> and first derivatives are already in memory;
310 : : !! other derivatives are not computed (except when choice=8 or 81)
311 : : !! This option is not compatible with choice=4,24 or 6
312 : : !! Warning: for cpopt= 1 or 3, derivatives wrt strains do not contain
313 : : !! the contribution due to the volume change;
314 : : !! i.e. <dp_lmn/dEps|in> are incomplete.
315 : : !!
316 : : !! NOTES
317 : : !! This application of the nonlocal operator is programmed using a direct
318 : : !! implementation of spherical harmonics (Ylm). Abinit used historically
319 : : !! Legendre polynomials for the application of nonlocal operator; but the
320 : : !! implementation of PAW algorithm enforced the use of Ylm.
321 : : !!
322 : : !! In the case signs=1, the array vectout is not used, nor modified
323 : : !! so that the same array as vectin can be used as a dummy argument;
324 : : !! the same is true for the pairs npwin-npwout, ffnlin-ffnlout,
325 : : !! kgin-kgout, ph3din-ph3dout, phkredin-phkxredout).
326 : : !!
327 : : !! Notes about choice==33:
328 : : !! **Since the 2nd derivative w.r.t q-vector is calculated along cartesian
329 : : !! directions, the 1/twopi**2 factor (that in the rest of the code is applied
330 : : !! in the reduced to cartesian derivative conversion process) is here
331 : : !! explicictly included in the formulas.
332 : : !!
333 : : !! **Notice that idir=1-9, in contrast to the strain perturbation (idir=1-6),
334 : : !! because this term is not symmetric w.r.t permutations of the two strain
335 : : !! indices.(Also applies for choice=25)
336 : : !!
337 : : !! **A -i factor has been factorized out in all the contributions of the second
338 : : !! q-gradient of the metric Hamiltonian and in the first and second q-gradients
339 : : !! of the atomic displacement Hamiltonian. This is lately included in the
340 : : !! matrix element calculation.
341 : : !!
342 : : !! TODO
343 : : !! * Complete implementation of spin-orbit
344 : : !!
345 : : !! SOURCE
346 : :
347 : 103206688 : subroutine nonlop_ylm(atindx1,choice,cpopt,cprjin,dimenl1,dimenl2,dimekbq,dimffnlin,dimffnlout,&
348 [ - + ]: 25801672 : & enl,enlout,ffnlin,ffnlout,gprimd,idir,indlmn,istwf_k,&
349 : 25801672 : & kgin,kgout,kpgin,kpgout,kptin,kptout,lambda,lmnmax,matblk,mgfft,&
350 [ - + ]: 25801672 : & mpi_enreg,natom,nattyp,ngfft,nkpgin,nkpgout,nloalg,nnlout,&
351 : 25801672 : & npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,phkxredin,phkxredout,ph1d,&
352 [ + + + - : 51603344 : & ph3din,ph3dout,signs,sij,svectout,ucvol,vectin,vectout,cprjin_left,&
- + - + -
+ - + ]
353 [ + + + - : 25801672 : & enlout_im,ndat_left,qdir)
- + ]
354 : :
355 : : !Arguments ------------------------------------
356 : : !scalars
357 : : integer,intent(in) :: choice,cpopt,dimenl1,dimenl2,dimekbq,dimffnlin,dimffnlout,idir
358 : : integer,intent(in) :: istwf_k,lmnmax,matblk,mgfft,natom,nkpgin,nkpgout,nnlout
359 : : integer,intent(in) :: npwin,npwout,nspinor,nspinortot,ntypat,paw_opt,signs
360 : : integer,intent(in),optional :: qdir,ndat_left
361 : : real(dp),intent(in) :: lambda,ucvol
362 : : type(MPI_type),intent(in) :: mpi_enreg
363 : : !arrays
364 : : integer,intent(in) :: atindx1(natom),kgin(3,npwin)
365 : : integer,intent(in),target :: indlmn(6,lmnmax,ntypat)
366 : : integer,intent(in) :: kgout(3,npwout),nattyp(ntypat),ngfft(18),nloalg(3)
367 : : real(dp),intent(in) :: enl(dimenl1,dimenl2,nspinortot**2,dimekbq)
368 : : real(dp),intent(in),target :: ffnlin(npwin,dimffnlin,lmnmax,ntypat)
369 : : real(dp),intent(in),target :: ffnlout(npwout,dimffnlout,lmnmax,ntypat)
370 : : real(dp),intent(in) :: gprimd(3,3)
371 : : real(dp),intent(in),target :: kpgin(npwin,nkpgin),kpgout(npwout,nkpgout)
372 : : real(dp),intent(in) :: kptin(3),kptout(3),ph1d(2,3*(2*mgfft+1)*natom)
373 : : real(dp),intent(in) :: phkxredin(2,natom),phkxredout(2,natom)
374 : : real(dp),intent(in) :: sij(dimenl1,ntypat*((paw_opt+1)/3))
375 : : real(dp),intent(inout) :: ph3din(2,npwin,matblk),ph3dout(2,npwout,matblk)
376 : : real(dp),intent(inout) :: vectin(:,:)
377 : : real(dp),intent(out) :: enlout(:)
378 : : real(dp),intent(out),optional :: enlout_im(:)
379 : : real(dp),intent(out) :: svectout(:,:)
380 : : real(dp),intent(inout) :: vectout (:,:)
381 : : type(pawcprj_type),intent(inout) :: cprjin(:,:)
382 : : type(pawcprj_type),optional,intent(in) :: cprjin_left(:,:)
383 : :
384 : : !Local variables-------------------------------
385 : : !scalars
386 : : integer :: choice_a,choice_b,cplex,cplex_enl,cplex_fac,ia,ia1,ia2,ia3,ia4,ia5
387 : : integer :: iatm,ic,idir1,idir2,ii,ierr,ilmn,ishift,ispinor,itypat,jc,mincat,mu,mua,mub,mu0
388 : : integer :: n1,n2,n3,nd2gxdt,ndat_left_,ndgxdt,ndgxdt_stored,nd2gxdtfac,ndgxdtfac
389 : : integer :: nincat,nkpgin_,nkpgout_,nlmn,nu,nua1,nua2,nub1,nub2,optder
390 : : real(dp) :: enlk!, tsec(2)
391 : : logical :: check,testnl,no_opernla_mv,no_opernlb_mv
392 : : character(len=500) :: message
393 : : !arrays
394 : : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
395 : : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
396 : 25801672 : integer,allocatable :: cplex_dgxdt(:),cplex_d2gxdt(:)
397 : 25801672 : integer,contiguous, pointer :: indlmn_typ(:,:)
398 : 25801672 : real(dp),allocatable :: d2gxdt(:,:,:,:,:),d2gxdtfac(:,:,:,:,:),d2gxdtfac_sij(:,:,:,:,:)
399 : 25801672 : real(dp),allocatable :: dgxdt(:,:,:,:,:),dgxdtfac(:,:,:,:,:),dgxdtfac_sij(:,:,:,:,:)
400 : 25801672 : real(dp),allocatable :: ddkk(:),fnlk(:),gmet(:,:)
401 : 25801672 : real(dp),allocatable :: gx(:,:,:,:),gxfac(:,:,:,:),gxfac_sij(:,:,:,:),gx_left(:,:,:,:)
402 : 25801672 : real(dp),allocatable :: sij_typ(:),strnlk(:)
403 : 25801672 : real(dp),allocatable :: work1(:),work2(:),work3(:,:),work4(:,:),work5(:,:,:),work6(:,:,:),work7(:,:,:)
404 : 25801672 : real(dp),contiguous, pointer :: ffnlin_typ(:,:,:),ffnlout_typ(:,:,:),kpgin_(:,:),kpgout_(:,:)
405 : :
406 : : ! **********************************************************************
407 : :
408 : : DBG_ENTER("COLL")
409 : :
410 : : ! call timab(1100,1,tsec)
411 : :
412 : : !Check consistency of arguments
413 : : !==============================================================
414 : :
415 : : !signs=1, almost all choices
416 : 25801672 : if (signs==1) then
417 [ + + ]: 1410081 : if(paw_opt<3) then
418 : : check=(choice==0 .or.choice==1 .or.choice==2 .or.choice==3 .or.choice==4 .or.&
419 : : & choice==23.or.choice==24.or.choice==5 .or.choice==51.or.choice==52.or.&
420 : : & choice==53.or.choice==54.or.choice==55.or.&
421 [ - + ]: 1167089 : & choice==6 .or.choice==8 .or.choice==81)
422 [ + - ]: 242992 : else if (paw_opt==3) then
423 : : check=(choice== 0.or.choice== 1.or.choice== 2.or.choice==3.or.choice==5.or.&
424 : : & choice==23.or.choice==51.or.choice==52.or.choice==53.or.choice==54.or.choice==55.or.&
425 [ - + ]: 242992 : & choice== 8.or.choice==81)
426 : : else
427 : : check = .false.
428 : : end if
429 : 0 : ABI_CHECK(check,'BUG: choice not compatible (for signs=1)')
430 : : end if
431 : :
432 : : !signs=2, less choices
433 [ + + ]: 25801672 : if (signs==2) then
434 : : check=(choice==0.or.choice==1.or.choice==2.or.choice==22.or.choice==25.or.choice==3 .or.&
435 : : & choice==5.or.choice==33.or.choice==51.or.choice==52.or.choice==53.or.choice==54.or.&
436 [ - + ]: 24391591 : & choice==7.or.choice==8.or.choice==81)
437 : 0 : ABI_CHECK(check,'BUG: choice not compatible (for signs=2)')
438 : : end if
439 : : !1<=idir<=6 is required when choice=3 and signs=2
440 [ + + + + ]: 25801672 : if (choice==3.and.signs==2) then
441 : 261762 : check=(idir>=1.and.idir<=6)
442 [ - + ]: 261762 : ABI_CHECK(check,'BUG: choice=3 and signs=2 requires 1<=idir<=6')
443 : : !1<=idir<=9 is required when choice= 25 or 33 and signs=2
444 [ + + + - ]: 25539910 : else if ((choice==25.or.choice==33).and.signs==2) then
445 : 622080 : check=(idir>=1.and.idir<=9)
446 [ - + ]: 622080 : ABI_CHECK(check,'BUG: choice= 25 or 33 and signs=2 requires 1<=idir<=9')
447 : : !1<=idir<=9 is required when choice==8/81 and signs=2
448 [ + + + + : 24917830 : else if ((choice==8.or.choice==81.or.choice==54).and.signs==2) then
+ + ]
449 : 72807 : check=(idir>=1.and.idir<=9)
450 [ - + ]: 72807 : ABI_CHECK(check,'BUG: choice=8/81 and signs=2 requires 1<=idir<=9')
451 : : else
452 : : ! signs=2 requires 1<=idir<=3 when choice>1
453 [ + + + + : 24845023 : check=(signs/=2.or.choice<=1.or.choice==7.or.(idir>=1.and.idir<=3))
+ + - + ]
454 : 0 : ABI_CHECK(check,'BUG: signs=2 requires 1<=idir<=3')
455 : : end if
456 : : !1<=qdir<=3 is required when choice==22 or choice==25 or choice==33 and signs=2
457 [ + + + - ]: 25801672 : if ((choice==22.or.choice==25.or.choice==33).and.signs==2) then
458 : 1002208 : check=(qdir>=1.and.qdir<=3)
459 [ - + ]: 1002208 : ABI_CHECK(check,'BUG: choice=22,25 or 33 and signs=2 requires 1<=qdir<=3')
460 : : end if
461 : :
462 : : !check allowed values for cpopt
463 : 25801672 : check=(cpopt>=-1.and.cpopt<=4)
464 [ - + ]: 25801672 : ABI_CHECK(check,'bad value for cpopt')
465 [ + + + - : 25801672 : check=(cpopt/=4.or.(choice/=4.and.choice/=24.and.choice/=6))
- + ]
466 : 0 : ABI_CHECK(check,'BUG: cpopt=4 not allowed for 2nd derivatives')
467 [ + + - + ]: 25801672 : check=(cpopt/=2.or.(choice/=8.and.choice/=81))
468 : 0 : ABI_CHECK(check,'BUG: cpopt=2 not allowed for choice=8,81, use cpopt=4 instead')
469 : : !check conditions for optional arguments
470 [ + + + - : 25801672 : check=((.not.present(cprjin_left)).or.(signs==1.and.choice==1))
- + ]
471 : 0 : ABI_CHECK(check,'BUG: when cprjin_left is present, must have choice=1,signs=1')
472 : : !protect special case choice==7
473 [ + + - + ]: 25801672 : check=(choice/=7.or.paw_opt==3)
474 : 0 : ABI_CHECK(check,'BUG: when choice=7, paw_opt must be 3')
475 : : !
476 : :
477 : : !Test: size of blocks of atoms
478 [ + + ]: 65888006 : mincat=min(NLO_MINCAT,maxval(nattyp))
479 [ + + - + ]: 25801672 : if (nloalg(2)<=0.and.mincat>matblk) then
480 : : write(message, '(a,a,a,i4,a,i4,a)' ) &
481 : 0 : & 'With nloc_mem<=0, mincat must be less than matblk.',ch10,&
482 : 0 : & 'Their value is ',mincat,' and ',matblk,'.'
483 : 0 : ABI_BUG(message)
484 : : end if
485 : 25801672 : ndat_left_=1
486 [ + + ]: 25801672 : if (present(ndat_left)) then
487 : 220960 : ndat_left_=ndat_left
488 : : end if
489 [ - + ]: 25801672 : if (nloalg(1)<2.or.nloalg(1)>10) then
490 : 0 : ABI_ERROR('nloalg(1) should be between 2 and 10.')
491 : : end if
492 : : ! Determine which implementation to use : matrix-vector (mv), matrix-vector with dgmev (mv-dgemv), or native
493 : : !nloalg(1)| opernla | opernlb
494 : : !------------------------------
495 : : ! 2 | mv-dgemv | mv-dgemv
496 : : ! 3 | mv | mv
497 : : ! 4(def) | native | native
498 : : ! 5 | mv-dgemv | mv
499 : : ! 6 | mv | mv-dgemv
500 : : ! 7 | mv-dgemv | native
501 : : ! 8 | native | mv
502 : : ! 9 | mv | native
503 : : ! 10 | native | mv-dgemv
504 [ + - ]: 25801672 : no_opernla_mv = nloalg(1)==4.or.nloalg(1)==8.or.nloalg(1)==10 ! have to be consistent with getcprj
505 [ + - ]: 25801672 : no_opernlb_mv = nloalg(1)==4.or.nloalg(1)==7.or.nloalg(1)==9
506 : :
507 : : !Define dimensions of projected scalars
508 : : !==============================================================
509 : :
510 : : !Define some useful variables
511 : 25801672 : n1=ngfft(1);n2=ngfft(2);n3=ngfft(3)
512 [ + + ]: 25801672 : choice_a=merge(choice,1,choice/=7);choice_b=choice_a
513 [ + + ]: 25801672 : if (cpopt>=2) choice_a=-choice_a
514 [ + + ]: 25801672 : cplex=2;if (istwf_k>1) cplex=1 !Take into account TR-symmetry
515 [ + + ]: 25801672 : cplex_enl=1;if (paw_opt>0) cplex_enl=2*dimenl1/(lmnmax*(lmnmax+1))
516 : 25801672 : cplex_fac=max(cplex,dimekbq)
517 [ + + + + : 25801672 : if ((nspinortot==2.or.cplex_enl==2).and.paw_opt>0.and.choice/=7) cplex_fac=2
+ + + + ]
518 : :
519 : : !Define dimensions of projected scalars
520 : 25801672 : ndgxdt=0;ndgxdtfac=0;nd2gxdt=0;nd2gxdtfac=0
521 [ + + ]: 25801672 : if (choice==2) then
522 [ + + ]: 1766744 : if (signs==1) ndgxdt=3
523 [ + + ]: 1766744 : if (signs==2) ndgxdt=1
524 [ + + ]: 1766744 : if (signs==2) ndgxdtfac=1
525 : : end if
526 [ + + ]: 25801672 : if (choice==22) then
527 [ + - ]: 380128 : if (signs==2) ndgxdt=1
528 [ + - ]: 380128 : if (signs==2) ndgxdtfac=1
529 : : end if
530 [ + + ]: 25801672 : if (choice==25) then
531 [ + - ]: 248832 : if (signs==2) ndgxdt=1
532 [ + - ]: 248832 : if (signs==2) ndgxdtfac=1
533 : : end if
534 [ + + ]: 25801672 : if (choice==3) then
535 [ + + ]: 288373 : if (signs==1) ndgxdt=6
536 [ + + ]: 288373 : if (signs==2) ndgxdt=1
537 [ + + ]: 288373 : if (signs==2) ndgxdtfac=1
538 : : end if
539 [ + + ]: 25801672 : if (choice==23) then
540 [ + - ]: 43756 : if (signs==1) ndgxdt=9
541 : : end if
542 [ + + ]: 25801672 : if (choice==4) then
543 [ + - ]: 24336 : if(signs==1) ndgxdt=3
544 [ + - ]: 24336 : if(signs==1) ndgxdtfac=3
545 [ + - ]: 24336 : if(signs==1) nd2gxdt=6
546 : : end if
547 [ - + ]: 25801672 : if (choice==24) then
548 [ # # ]: 0 : if(signs==1) ndgxdt=3
549 [ # # ]: 0 : if(signs==1) ndgxdtfac=3
550 [ # # ]: 0 : if(signs==1) nd2gxdt=6
551 : : end if
552 [ + + ]: 25801672 : if (choice==33) then
553 [ + - ]: 373248 : if(signs==2) ndgxdt=2
554 [ + - ]: 373248 : if(signs==2) ndgxdtfac=2
555 [ + - ]: 373248 : if(signs==2) nd2gxdt=3
556 [ + - ]: 373248 : if(signs==2) nd2gxdtfac=3
557 : : end if
558 [ + + ]: 25801672 : if (choice==5) then
559 [ + + ]: 852101 : if(signs==1) ndgxdt=3
560 [ + + ]: 852101 : if(signs==2) ndgxdt=1
561 [ + + ]: 852101 : if(signs==2) ndgxdtfac=1
562 : : end if
563 [ + + ]: 25801672 : if (choice==51) then
564 [ - + ]: 106864 : if(signs==1) ndgxdt=3
565 [ + - ]: 106864 : if(signs==2) ndgxdt=1
566 [ + - ]: 106864 : if(signs==2) ndgxdtfac=1
567 : : end if
568 [ - + ]: 25801672 : if (choice==52) then
569 [ # # ]: 0 : if(signs==1) ndgxdt=3
570 [ # # ]: 0 : if(signs==2) ndgxdt=1
571 [ # # ]: 0 : if(signs==2) ndgxdtfac=1
572 : : end if
573 [ - + ]: 25801672 : if (choice==53) then
574 [ # # ]: 0 : if(signs==1) ndgxdt=3
575 [ # # ]: 0 : if(signs==1) ndgxdtfac=3
576 [ # # ]: 0 : if(signs==2) ndgxdt=1
577 [ # # ]: 0 : if(signs==2) ndgxdtfac=1
578 : : end if
579 [ + + ]: 25801672 : if (choice==54) then
580 [ + + ]: 34432 : if(signs==1) ndgxdt=6
581 [ + + ]: 34432 : if(signs==1) ndgxdtfac=6
582 [ + + ]: 34432 : if(signs==1) nd2gxdt=9
583 [ + + ]: 34432 : if(signs==2) ndgxdt=1
584 [ + + ]: 34432 : if(signs==2) nd2gxdt=1
585 [ + + ]: 34432 : if(signs==2) ndgxdtfac=1
586 [ + + ]: 34432 : if(signs==2) nd2gxdtfac=1
587 : : end if
588 [ + + ]: 25801672 : if (choice==55) then
589 [ + - ]: 2576 : if(signs==1) ndgxdt=9
590 [ + - ]: 2576 : if(signs==1) ndgxdtfac=9
591 [ + - ]: 2576 : if(signs==1) nd2gxdt=18
592 : : end if
593 [ + + ]: 25801672 : if (choice==6) then
594 [ + - ]: 2400 : if(signs==1) ndgxdt=9
595 [ + - ]: 2400 : if(signs==1) ndgxdtfac=9
596 [ + - ]: 2400 : if(signs==1) nd2gxdt=54
597 : : end if
598 [ + + ]: 25801672 : if (choice==8) then
599 [ + + ]: 47271 : if(signs==1) ndgxdt=3
600 [ + + ]: 47271 : if(signs==1) ndgxdtfac=3
601 [ + + ]: 47271 : if(signs==1) nd2gxdt=6
602 [ + + ]: 47271 : if(signs==2) ndgxdt=2
603 [ + + ]: 47271 : if(signs==2) ndgxdtfac=2
604 [ + + ]: 47271 : if(signs==2) nd2gxdt=1
605 [ + + ]: 47271 : if(signs==2) nd2gxdtfac=1
606 : : end if
607 [ + + ]: 25801672 : if (choice==81) then
608 [ + + ]: 3456 : if(signs==1) ndgxdt=3
609 [ + + ]: 3456 : if(signs==1) ndgxdtfac=3
610 [ + + ]: 3456 : if(signs==1) nd2gxdt=6
611 [ + + ]: 3456 : if(signs==2) ndgxdt=1
612 [ + + ]: 3456 : if(signs==2) ndgxdtfac=1
613 [ + + ]: 3456 : if(signs==2) nd2gxdt=1
614 [ + + ]: 3456 : if(signs==2) nd2gxdtfac=1
615 : : end if
616 [ - + ]: 25801672 : ABI_CHECK(ndgxdtfac<=ndgxdt,"BUG: ndgxdtfac>ndgxdt!")
617 [ + + ]: 25801672 : optder=0;if (ndgxdtfac>0) optder=1
618 [ + + ]: 25801672 : if (nd2gxdtfac>0) optder=2
619 : :
620 : : !Consistency tests
621 [ + + ]: 25801672 : if (cpopt==4) then
622 [ + + - + ]: 2170126 : if (ndgxdt>0.and.cprjin(1,1)%ncpgr<=0) then
623 : 0 : message='cprjin%ncpgr=0 not allowed with cpopt=4 and these (choice,signs) !'
624 : 0 : ABI_BUG(message)
625 : : end if
626 : : end if
627 [ + + ]: 25801672 : if (cpopt==1.or.cpopt==3) then
628 [ - + ]: 97076 : if (cprjin(1,1)%ncpgr<ndgxdt) then
629 : 0 : message='should have cprjin%ncpgr>=ndgxdt with cpopt=1 or 3 !'
630 : 0 : ABI_BUG(message)
631 : : end if
632 : : end if
633 : :
634 : :
635 : : !Additional steps before calculation
636 : : !==============================================================
637 : :
638 : : !Initialize output arrays
639 [ + + ]: 25801672 : if (signs==1) then
640 [ + - - + : 4230243 : ABI_MALLOC(fnlk,(3*natom))
+ - + - +
- - + -
+ ]
641 [ - + ]: 1410081 : ABI_MALLOC(ddkk,(6))
642 [ - + ]: 1410081 : ABI_MALLOC(strnlk,(6))
643 [ + + + + : 32032749 : enlk=zero;fnlk=zero;ddkk=zero;strnlk=zero
+ + ]
644 [ + + ]: 7679992 : enlout(:)=zero
645 [ + + ]: 1410081 : if (present(enlout_im)) then
646 [ + + ]: 2217120 : enlout_im(:)=zero
647 : : end if
648 : : end if
649 [ + + ]: 25801672 : if (signs==2) then
650 [ + + + + : 15081112032 : if (paw_opt==0.or.paw_opt==1.or.paw_opt==4) vectout(:,:)=zero
+ + ]
651 [ + + - + : 24391591 : if (paw_opt==2.and.choice==1) vectout(:,:)=-lambda*vectin(:,:)
- - - - ]
652 [ + + + - : 24415373 : if (paw_opt==2.and.choice> 1) vectout(:,:)=zero
+ + + + ]
653 [ + + ]: 24391591 : if (paw_opt==3.or.paw_opt==4) then
654 [ + + + + : 9841767948 : if (choice==1) svectout(:,:)=vectin(:,:)
+ + ]
655 [ + + + + : 1369560298 : if (choice> 1) svectout(:,:)=zero
+ + ]
656 : : end if
657 : : end if
658 : :
659 : : !Eventually re-compute (k+G) vectors (and related data)
660 : 25801672 : nkpgin_=0
661 [ + + ]: 25801672 : if (choice==2.or.choice==22.or.choice==25.or.choice==33.or.choice==54) nkpgin_=3
662 [ + + ]: 25801672 : if (signs==1) then
663 [ + + ]: 1410081 : if (choice==4.or.choice==24) nkpgin_=9
664 [ + + ]: 1410081 : if (choice==3.or.choice==23.or.choice==6) nkpgin_=3
665 [ + + ]: 1410081 : if (choice==55) nkpgin_=3
666 : : end if
667 [ + + ]: 25801672 : if (nkpgin<nkpgin_) then
668 [ - + + - : 7568828 : ABI_MALLOC(kpgin_,(npwin,nkpgin_))
- + + - +
- + - - +
- + ]
669 : :
670 : : !For the metric derivatives we need kpg in Cartesian coordinates
671 [ + + ]: 1892207 : if (choice==33) then
672 : 373248 : call mkkpgcart(gprimd,kgin,kpgin_,kptin,nkpgin_,npwin)
673 : : else
674 : 1518959 : call mkkpg(kgin,kpgin_,kptin,nkpgin_,npwin)
675 : : end if
676 : :
677 : : else
678 : 23909465 : nkpgin_ = nkpgin
679 : 23909465 : kpgin_ => kpgin
680 : : end if
681 : :
682 [ + + ]: 25801672 : nkpgout_=0
683 [ + + ]: 3091757 : if ((choice==2.or.choice==22.or.choice==25.or.choice==3.or.choice==33.or.choice==54).and.signs==2) nkpgout_=3
684 [ + + ]: 25801672 : if (nkpgout<nkpgout_) then
685 [ - + + - : 7875440 : ABI_MALLOC(kpgout_,(npwout,nkpgout_))
- + + - +
- + - - +
- + ]
686 : :
687 : : !For the metric derivatives we need kpg in Cartesian coordinates
688 [ + + ]: 1968860 : if (choice==33) then
689 : 373248 : call mkkpgcart(gprimd,kgout,kpgout_,kptout,nkpgout_,npwout)
690 : : else
691 : 1595612 : call mkkpg(kgout,kpgout_,kptout,nkpgout_,npwout)
692 : : end if
693 : :
694 : : else
695 : 23832812 : nkpgout_ = nkpgout
696 : 23832812 : kpgout_ => kpgout
697 : : end if
698 : :
699 : : !Big loop on atom types
700 : : !==============================================================
701 : :
702 : 25801672 : ia1=1;iatm=0
703 [ + + ]: 65888006 : do itypat=1,ntypat
704 : :
705 : : ! Get atom loop indices for different types:
706 : 40086334 : ia2=ia1+nattyp(itypat)-1;ia5=1
707 : :
708 : : ! Select quantities specific to the current type of atom
709 [ + + + + ]: 501003762 : nlmn=count(indlmn(3,:,itypat)>0)
710 : :
711 : : ! Test on local part
712 : 40086334 : testnl=(paw_opt/=0)
713 [ + + + + : 49699147 : if (paw_opt==0) testnl=any(abs(enl(:,:,:,:))>tol10)
+ + + + +
+ + + ]
714 : :
715 : : ! Some non-local part is to be applied for that type of atom
716 [ + + ]: 40086334 : if (testnl) then
717 : :
718 : : ! Store some quantities depending only of the atom type
719 : 40084990 : ffnlin_typ => ffnlin(:,:,:,itypat)
720 : 40084990 : indlmn_typ => indlmn(:,:,itypat)
721 [ + + ]: 40084990 : if (signs==2) then
722 : 38055821 : ffnlout_typ => ffnlout(:,:,:,itypat)
723 : : end if
724 [ + + ]: 40084990 : if (paw_opt>=2) then
725 [ + - - + : 68151843 : ABI_MALLOC(sij_typ,(nlmn*(nlmn+1)/2))
+ - + - +
- - + -
+ ]
726 [ + + ]: 22717281 : if (cplex_enl==1) then
727 [ + + ]: 1500580140 : do ilmn=1,nlmn*(nlmn+1)/2
728 : 1500580140 : sij_typ(ilmn)=sij(ilmn,itypat)
729 : : end do
730 : : else
731 [ + + ]: 348226337 : do ilmn=1,nlmn*(nlmn+1)/2
732 : 348226337 : sij_typ(ilmn)=sij(2*ilmn-1,itypat)
733 : : end do
734 : : end if
735 : : else
736 [ - + ]: 17367709 : ABI_MALLOC(sij_typ,(0))
737 : : end if
738 : :
739 : : ! Loop over atoms of the same type
740 : : ! ==============================================================
741 : :
742 : : ! Cut the sum on different atoms in blocks, to allow memory saving.
743 : : ! Inner summations on atoms will be done from ia3 to ia4.
744 : : ! Note: the maximum range from ia3 to ia4 is mincat (max. increment of atoms).
745 : :
746 [ - + - - : 40084990 : do ia3=ia1,ia2,mincat
+ + ]
747 : 40217298 : ia4=min(ia2,ia3+mincat-1)
748 : : ! Give the increment of number of atoms in this subset.
749 : 40217298 : nincat=ia4-ia3+1
750 : :
751 : : ! Prepare the phase factors if they were not already computed
752 [ + + ]: 40217298 : if (nloalg(2)<=0) then
753 : 362256 : call ph1d3d(ia3,ia4,kgin,matblk,natom,npwin,n1,n2,n3,phkxredin,ph1d,ph3din)
754 : : end if
755 : :
756 : : ! Allocate memory for projected scalars
757 [ - + + - : 241303788 : ABI_MALLOC(gx,(cplex,nlmn,nincat,nspinor))
+ + + - -
+ + - - +
+ - + - +
+ - + -
+ ]
758 [ - + + - : 241303788 : ABI_MALLOC(gxfac,(cplex_fac,nlmn,nincat,nspinor))
+ + + - -
+ + - - +
+ - + - +
+ - + -
+ ]
759 [ - + + - : 281521086 : ABI_MALLOC(dgxdt,(cplex,ndgxdt,nlmn,nincat,nspinor))
+ + + - +
+ + - - +
+ - - + +
- + - + +
- + - + ]
760 [ - + + - : 281521086 : ABI_MALLOC(d2gxdt,(cplex,nd2gxdt,nlmn,nincat,nspinor))
+ + + - +
+ + - - +
+ - - + +
- + - + +
- + - + ]
761 [ - + + - : 281521086 : ABI_MALLOC(d2gxdtfac,(cplex_fac,nd2gxdtfac,nlmn,nincat,nspinor))
+ + + - +
+ + - - +
+ - - + +
- + - + +
- + - + ]
762 [ - + + - : 281521086 : ABI_MALLOC(dgxdtfac,(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor))
+ + + - +
+ + - - +
+ - - + +
- + - + +
- + - + ]
763 [ + + + + : 3957159764 : gx(:,:,:,:)=zero;gxfac(:,:,:,:)=zero
+ + + + +
+ + + + +
+ + ]
764 [ + + + + : 438619774 : if (ndgxdt>0) dgxdt(:,:,:,:,:)=zero
+ + + + +
+ + + ]
765 [ + + + + : 296514634 : if (ndgxdtfac>0) dgxdtfac(:,:,:,:,:)=zero
+ + + + +
+ + + ]
766 [ + + + + : 136998391 : if (nd2gxdt>0) d2gxdt(:,:,:,:,:)=zero
+ + + + +
+ + + ]
767 [ + + + + : 111542559 : if (nd2gxdtfac>0) d2gxdtfac(:,:,:,:,:)=zero
+ + + + +
+ + + ]
768 [ + + ]: 40217298 : if (paw_opt>=3) then
769 [ - + + - : 105730395 : ABI_MALLOC(gxfac_sij,(cplex,nlmn,nincat,nspinor))
- + + - -
+ + - - +
+ - - + -
+ ]
770 [ - + + - : 148022553 : ABI_MALLOC(dgxdtfac_sij,(cplex,ndgxdtfac,nlmn,nincat,nspinor))
+ + + - -
+ + - - +
+ - - + +
- + - + +
- + - + ]
771 [ - + + - : 148022553 : ABI_MALLOC(d2gxdtfac_sij,(cplex,nd2gxdtfac,nlmn,nincat,nspinor))
+ + + - -
+ + - - +
+ - - + +
- + - + +
- + + - ]
772 [ + + + + : 1168000089 : gxfac_sij(:,:,:,:)=zero
+ + + + ]
773 [ + + + + : 101255467 : if (ndgxdtfac>0) dgxdtfac_sij(:,:,:,:,:)=zero
+ + + + +
+ + + ]
774 [ + + + + : 22353333 : if (nd2gxdtfac>0) d2gxdtfac_sij(:,:,:,:,:) = zero
+ + + + +
+ + + ]
775 : : else
776 [ - + ]: 19071219 : ABI_MALLOC(gxfac_sij,(0,0,0,0))
777 [ - + ]: 19071219 : ABI_MALLOC(dgxdtfac_sij,(0,0,0,0,0))
778 [ - + ]: 19071219 : ABI_MALLOC(d2gxdtfac_sij,(0,0,0,0,0))
779 : : end if
780 : :
781 : : ! When istwf_k > 1, gx derivatives can be real or pure imaginary
782 : : ! cplex_dgxdt(i) = 1 if dgxdt(1,i,:,:) is real, 2 if it is pure imaginary
783 : : ! cplex_d2gxdt(i) = 1 if d2gxdt(1,i,:,:) is real, 2 if it is pure imaginary
784 [ + + + - : 120651894 : ABI_MALLOC(cplex_dgxdt,(ndgxdt))
+ - + + -
+ - + ]
785 [ + + + - : 120651894 : ABI_MALLOC(cplex_d2gxdt,(nd2gxdt))
+ - + + -
+ + - ]
786 [ + + + + ]: 89643256 : cplex_dgxdt(:) = 1 ; cplex_d2gxdt(:) = 1
787 [ + + ]: 40217298 : if(ndgxdt > 0) then
788 : : if (choice==5.or.choice==51.or.choice==52.or.choice==53.or. &
789 [ + + + + ]: 8258442 : & choice==8.or.choice==81) cplex_dgxdt(:) = 2
790 [ + + + + : 5144153 : if (choice==54.and.signs==1) cplex_dgxdt(4:6) = 2
+ + ]
791 [ + + + + : 5099737 : if (choice==54.and.signs==2) cplex_dgxdt(:) = 2
+ + ]
792 [ + + + - : 5074265 : if (choice==55.and.signs==1) cplex_dgxdt(7:9) = 2
+ + ]
793 : : end if
794 [ + + ]: 40217298 : if(nd2gxdt > 0) then
795 [ + + + + ]: 830429 : if (choice==54) cplex_d2gxdt(:) = 2
796 [ + + + - : 655997 : if (choice==55.and.signs==1) cplex_d2gxdt(1:18)= 2
+ + ]
797 : : end if
798 : :
799 : : ! Compute projection of current wave function |c> on each
800 : : ! non-local projector: <p_lmn|c>
801 : : ! ==============================================================
802 : :
803 : : ! Retrieve eventually <p_lmn|c> coeffs (and derivatives)
804 [ + + ]: 40217298 : if (cpopt>=2) then
805 [ + + ]: 21904115 : do ispinor=1,nspinor
806 [ + + ]: 34996768 : do ia=1,nincat
807 [ + + + + ]: 353684593 : gx(1:cplex,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)
808 : : end do
809 : : end do
810 : : end if
811 [ + + + + ]: 40217298 : if (cpopt==4.and.ndgxdt>0) then
812 : 1321288 : ndgxdt_stored = cprjin(1,1)%ncpgr
813 : 1321288 : ishift=0
814 [ + + + + : 1321288 : if (((choice==2).or.(choice==3)).and.(ndgxdt_stored>ndgxdt).and.(signs==2)) ishift=idir-ndgxdt
+ - ]
815 [ + + + + : 1321288 : if ((choice==2).and.(ndgxdt_stored==9).and.(signs==2)) ishift=ishift+6
+ - ]
816 [ + + + + : 1321288 : if (choice==2.and.(ndgxdt_stored>ndgxdt).and.(signs==1)) ishift=ndgxdt_stored-ndgxdt
- + ]
817 [ - + ]: 1321288 : if(cplex == 2) then
818 [ + + ]: 2661008 : do ispinor=1,nspinor
819 [ + + ]: 4184436 : do ia=1,nincat
820 [ + + + - : 2863148 : if (ndgxdt_stored==ndgxdt.or.(ndgxdt_stored>ndgxdt.and.((choice==2).or.(choice==3)))) then
+ + ]
821 [ + + + + : 52044884 : dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)
+ + ]
822 [ + - + - ]: 92304 : else if (signs==2.and.ndgxdt_stored==3) then
823 [ + - + + : 92304 : if (choice==5.or.choice==51.or.choice==52) then ! ndgxdt=1
- ]
824 [ + + + + ]: 2259600 : dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir,1:nlmn)
825 : : else if (choice==53) then ! ndgxdt=1
826 : 0 : idir1 = mod(idir+1,3)+1
827 [ # # # # ]: 0 : dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir1,1:nlmn)
828 : : else if (choice==8) then ! ndgxdt=2
829 : 768 : idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
830 [ + + + + ]: 19200 : dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir1,1:nlmn)
831 [ + + + + ]: 19200 : dgxdt(1:2,2,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir2,1:nlmn)
832 : : else if (choice==81) then ! ndgxdt=1
833 : 1152 : idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
834 [ + + + + ]: 28800 : dgxdt(1:2,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(1:2,idir2,1:nlmn)
835 : : end if
836 : : end if
837 : : end do
838 : : end do
839 : : else ! cplex != 2
840 [ # # ]: 0 : do ispinor=1,nspinor
841 [ # # ]: 0 : do ia=1,nincat
842 [ # # ]: 0 : do ilmn=1,nlmn
843 [ # # # # : 0 : if (ndgxdt_stored==ndgxdt.or.(ndgxdt_stored>ndgxdt.and.((choice==2).or.(choice==3)))) then
# # ]
844 [ # # ]: 0 : do ii=1,ndgxdt
845 : 0 : ic = cplex_dgxdt(ii)
846 : 0 : dgxdt(1,ii,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)
847 : : end do
848 [ # # # # ]: 0 : else if (signs==2.and.ndgxdt_stored==3) then
849 [ # # # # : 0 : if (choice==5.or.choice==51.or.choice==52) then ! ndgxdt=1
# ]
850 : 0 : dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir,ilmn)
851 : : else if (choice==53) then ! ndgxdt=1
852 : 0 : idir1 = mod(idir+1,3)+1
853 [ # # ]: 0 : dgxdt(1,1,1:nlmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir1,1:nlmn)
854 : : else if (choice==8) then ! ndgxdt=2
855 : 0 : idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
856 : 0 : dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir1,ilmn)
857 : 0 : dgxdt(1,2,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(2),idir2,ilmn)
858 : : else if (choice==81) then ! ndgxdt=1
859 : 0 : idir1=(idir-1)/3+1; idir2=mod((idir-1),3)+1
860 : 0 : dgxdt(1,1,ilmn,ia,ispinor)=cprjin(iatm+ia,ispinor)%dcp(cplex_dgxdt(1),idir2,ilmn)
861 : : end if
862 : : end if
863 : : end do
864 : : end do
865 : : end do
866 : : end if ! cplex == 2
867 : : end if ! cpopt==4 and ndgxdt>0
868 : :
869 : : ! Computation or <p_lmn|c> (and derivatives) for this block of atoms if :
870 : : ! <p_lmn|c> are not in memory : cpopt<=1
871 : : ! OR <p_lmn|c> are in memory, but we need derivatives : cpopt<=3 and abs(choice_a)>1
872 : : ! OR <p_lmn|c> and first derivatives are in memory, but we need second derivatives : choice=8 or 81
873 [ + + + + : 40217298 : if (cpopt<=1.or.(cpopt<=3.and.abs(choice_a)>1).or.choice==8.or.choice==81) then
+ + + + +
+ ]
874 : : ! if ((cpopt<4.and.choice_a/=-1).or.choice==8.or.choice==81) then
875 [ + + + - ]: 29411818 : if (abs(choice_a)>1.or.no_opernla_mv) then
876 : : ! call timab(1101,1,tsec)
877 : : call opernla_ylm(choice_a,cplex,cplex_dgxdt,cplex_d2gxdt,dimffnlin,d2gxdt,dgxdt,ffnlin_typ,gx,&
878 : : & ia3,idir,indlmn_typ,istwf_k,kpgin_,matblk,mpi_enreg,nd2gxdt,ndgxdt,nincat,nkpgin_,nlmn,&
879 : 29411818 : & nloalg,npwin,nspinor,ph3din,signs,ucvol,vectin,qdir=qdir)
880 : : ! call timab(1101,2,tsec)
881 : : else
882 : : ! call timab(1102,1,tsec)
883 : : call opernla_ylm_mv(choice_a,cplex,dimffnlin,ffnlin_typ,gx,&
884 : : & ia3,indlmn_typ,istwf_k,matblk,mpi_enreg,nincat,nlmn,&
885 : 0 : & nloalg,npwin,nspinor,ph3din,ucvol,vectin)
886 : : ! call timab(1102,2,tsec)
887 : : end if
888 : : end if
889 : :
890 : : ! Transfer result to output variable cprj (if requested)
891 : : ! cprj(:)%cp receive the <p_i|Psi> factors (p_i: non-local projector)
892 : : ! Be careful: cprj(:)%dcp does not exactly contain the derivative of cprj(:)%cp.
893 : : ! - Volume contributions (in case of strain derivative) are not included,
894 : : ! - Global coordinate transformation are not applied,
895 : : ! cprj(:)%dcp is meant to be a restart argument of the present nonlop routine.
896 [ + + ]: 40217298 : if (cpopt==0.or.cpopt==1) then
897 [ + + ]: 18661806 : do ispinor=1,nspinor
898 [ + + ]: 32110120 : do ia=1,nincat
899 : 13448314 : cprjin(iatm+ia,ispinor)%nlmn=nlmn
900 [ + + + + ]: 388048460 : cprjin(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)=gx(1:cplex,1:nlmn,ia,ispinor)
901 [ + + + + ]: 30788371 : if (cplex==1) cprjin(iatm+ia,ispinor)%cp(2,1:nlmn)=zero
902 : : end do
903 : : end do
904 : : end if
905 [ + + + + ]: 40217298 : if ((cpopt==1.or.cpopt==3).and.ndgxdt>0) then
906 : 173888 : ishift=0
907 [ + + + + ]: 173888 : if ((choice==2).and.(cprjin(1,1)%ncpgr>ndgxdt)) ishift=cprjin(1,1)%ncpgr-ndgxdt
908 [ - + ]: 173888 : if(cplex==2)then
909 [ + + ]: 347776 : do ispinor=1,nspinor
910 [ + + ]: 542416 : do ia=1,nincat
911 [ + + + + : 15554976 : cprjin(iatm+ia,ispinor)%dcp(1:2,1+ishift:ndgxdt+ishift,1:nlmn)=dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)
+ + ]
912 : : ! cprjin(iatm+ia,ispinor)%dcp(1:2,1:ndgxdt,1:nlmn)=dgxdt(1:2,1:ndgxdt,1:nlmn,ia,ispinor)
913 : : end do
914 : : end do
915 : : else
916 [ # # ]: 0 : do ispinor=1,nspinor
917 [ # # ]: 0 : do ia=1,nincat
918 [ # # ]: 0 : do ilmn=1,nlmn
919 [ # # ]: 0 : do ii=1,ndgxdt
920 : 0 : ic = cplex_dgxdt(ii) ; jc = 3 - ic
921 : 0 : cprjin(iatm+ia,ispinor)%dcp(ic,ii+ishift,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor)
922 : 0 : cprjin(iatm+ia,ispinor)%dcp(jc,ii+ishift,ilmn)=zero
923 : : ! cprjin(iatm+ia,ispinor)%dcp(ic,ii,ilmn)=dgxdt(1,ii,ilmn,ia,ispinor)
924 : : ! cprjin(iatm+ia,ispinor)%dcp(jc,ii,ilmn)=zero
925 : : end do
926 : : end do
927 : : end do
928 : : end do
929 : : end if
930 : : end if
931 : :
932 : : ! If choice==0, that's all for these atoms !
933 [ + + ]: 40217298 : if (choice>0) then
934 [ + + ]: 37303032 : if(choice/=7) then
935 : : ! call timab(1105,1,tsec)
936 : : ! Contraction from <p_i|c> to Sum_j[Dij.<p_j|c>] (and derivatives)
937 : : call opernlc_ylm(atindx1,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_enl,cplex_fac,dgxdt,dgxdtfac,dgxdtfac_sij,&
938 : : & d2gxdt,d2gxdtfac,d2gxdtfac_sij,dimenl1,dimenl2,dimekbq,enl,gx,gxfac,gxfac_sij,&
939 : : & iatm,indlmn_typ,itypat,lambda,mpi_enreg,natom,ndgxdt,ndgxdtfac,nd2gxdt,nd2gxdtfac,&
940 : 35673080 : & nincat,nlmn,nspinor,nspinortot,optder,paw_opt,sij_typ)
941 : : ! call timab(1105,2,tsec)
942 : : else
943 [ + - + - : 79693664 : gxfac_sij=gx
+ - + - -
+ - - + +
+ + + + +
+ ]
944 : : end if
945 : :
946 : : ! Operate with the non-local potential on the projected scalars,
947 : : ! in order to get contributions to energy/forces/stress/dyn.mat
948 : : ! ==============================================================
949 [ + + ]: 37303032 : if (signs==1) then
950 [ + + ]: 768799 : if (.not.present(cprjin_left)) then
951 : : ! call timab(1106,1,tsec)
952 : : call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
953 : : & enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,ndat_left_,nd2gxdt,ndgxdt,ndgxdtfac,&
954 [ - + - - : 547839 : & nincat,nlmn,nnlout,nspinor,paw_opt,strnlk)
- - - + -
- ]
955 : : ! call timab(1106,2,tsec)
956 : : else
957 [ - + + - : 1325760 : ABI_MALLOC(gx_left,(cplex,nlmn,nincat,nspinor*ndat_left_))
- + + - -
+ + - - +
+ - + - +
- - + -
+ ]
958 : : ! Retrieve <p_lmn|c> coeffs
959 [ + + ]: 2217120 : do ispinor=1,nspinor*ndat_left_
960 [ + + ]: 6773920 : do ia=1,nincat
961 [ + + + + ]: 149784960 : gx_left(1:cplex,1:nlmn,ia,ispinor)=cprjin_left(iatm+ia,ispinor)%cp(1:cplex,1:nlmn)
962 : : end do
963 : : end do
964 : : ! TODO
965 : : ! if (cpopt==4.and.ndgxdt>0) then
966 : : ! do ispinor=1,nspinor
967 : : ! do ia=1,nincat
968 : : ! dgxdt_left(1:cplex,1:ndgxdt,1:nlmn,ia,ispinor)=cprjin_left(iatm+ia,ispinor)%dcp(1:cplex,1:ndgxdt,1:nlmn)
969 : : ! end do
970 : : ! end do
971 : : ! end if
972 [ + - ]: 220960 : if (present(enlout_im)) then
973 : : ! call timab(1108,1,tsec)
974 : : call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
975 : : & enlk,enlout,fnlk,gx_left,gxfac,gxfac_sij,ia3,natom,ndat_left_,nd2gxdt,ndgxdt,ndgxdtfac,&
976 [ - + - - : 220960 : & nincat,nlmn,nnlout,nspinor,paw_opt,strnlk,enlout_im=enlout_im)
- - - + -
- - - - +
- - - + -
- ]
977 : : ! call timab(1108,2,tsec)
978 : : else
979 : : ! call timab(1107,1,tsec)
980 : : call opernld_ylm(choice_b,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
981 : : & enlk,enlout,fnlk,gx_left,gxfac,gxfac_sij,ia3,natom,ndat_left_,nd2gxdt,ndgxdt,ndgxdtfac,&
982 [ # # # # : 0 : & nincat,nlmn,nnlout,nspinor,paw_opt,strnlk)
# # # # #
# ]
983 : : ! call timab(1107,2,tsec)
984 : : end if
985 [ - + ]: 220960 : ABI_FREE(gx_left)
986 : : end if
987 : : end if ! signs == 1
988 : :
989 : : ! Operate with the non-local potential on the projected scalars,
990 : : ! in order to get matrix element
991 : : ! ==============================================================
992 [ + + ]: 37303032 : if (signs==2) then
993 : : ! Prepare the phase factors if they were not already computed
994 [ + + ]: 36534233 : if(nloalg(2)<=0) then
995 : 362256 : call ph1d3d(ia3,ia4,kgout,matblk,natom,npwout,n1,n2,n3,phkxredout,ph1d,ph3dout)
996 : : end if
997 [ + + + - ]: 36534233 : if (abs(choice_b)>1.or.no_opernlb_mv) then
998 : : ! call timab(1103,1,tsec)
999 : : call opernlb_ylm(choice_b,cplex,cplex_dgxdt,cplex_d2gxdt,cplex_fac,&
1000 : : & d2gxdtfac,d2gxdtfac_sij,dgxdtfac,dgxdtfac_sij,dimffnlout,ffnlout_typ,gxfac,gxfac_sij,ia3,&
1001 : : & idir,indlmn_typ,kpgout_,matblk,ndgxdtfac,nd2gxdtfac,nincat,nkpgout_,nlmn,&
1002 : 36534233 : & nloalg,npwout,nspinor,paw_opt,ph3dout,svectout,ucvol,vectout,qdir=qdir)
1003 : : ! call timab(1103,2,tsec)
1004 : : else
1005 : : ! call timab(1104,1,tsec)
1006 : : call opernlb_ylm_mv(choice_b,cplex,cplex_fac,&
1007 : : & dimffnlout,ffnlout_typ,gxfac,gxfac_sij,ia3,indlmn_typ,matblk,nincat,nlmn,&
1008 : 0 : & nloalg,npwout,nspinor,paw_opt,ph3dout,svectout,ucvol,vectout)
1009 : : ! call timab(1104,2,tsec)
1010 : : end if
1011 : : end if ! signs == 2
1012 : :
1013 : : end if ! choice>=0
1014 : :
1015 : : ! Deallocate temporary projected scalars
1016 [ - + ]: 40217298 : ABI_FREE(gx)
1017 [ - + ]: 40217298 : ABI_FREE(gxfac)
1018 [ - + ]: 40217298 : ABI_FREE(dgxdt)
1019 [ - + ]: 40217298 : ABI_FREE(dgxdtfac)
1020 [ - + ]: 40217298 : ABI_FREE(d2gxdt)
1021 [ - + ]: 40217298 : ABI_FREE(d2gxdtfac)
1022 [ - + ]: 40217298 : ABI_FREE(dgxdtfac_sij)
1023 [ - + ]: 40217298 : ABI_FREE(d2gxdtfac_sij)
1024 [ - + ]: 40217298 : ABI_FREE(gxfac_sij)
1025 [ - + ]: 40217298 : ABI_FREE(cplex_dgxdt)
1026 [ - + ]: 40217298 : ABI_FREE(cplex_d2gxdt)
1027 : :
1028 : : ! End sum on atom subset loop
1029 [ + + ]: 40355411 : iatm=iatm+nincat;ia5=ia5+nincat
1030 : : end do
1031 : : !if (paw_opt>=2) then
1032 [ - + ]: 40084990 : ABI_FREE(sij_typ)
1033 : : !end if
1034 : :
1035 : : ! End condition of existence of a non-local part
1036 : : else
1037 [ - + ]: 1344 : if (cpopt==0.or.cpopt==1) then
1038 [ # # ]: 0 : do ispinor=1,nspinor
1039 [ # # ]: 0 : do ia=1,nattyp(itypat)
1040 [ # # # # ]: 0 : cprjin(iatm+ia,ispinor)%cp(:,1:nlmn)=zero
1041 : : end do
1042 : : end do
1043 : : end if
1044 [ - + - - ]: 1344 : if ((cpopt==1.or.cpopt==3).and.ndgxdt>0) then
1045 : 0 : ishift=0
1046 [ # # # # ]: 0 : if ((choice==2).and.(cprjin(1,1)%ncpgr>ndgxdt)) ishift=cprjin(1,1)%ncpgr-ndgxdt
1047 [ # # ]: 0 : do ispinor=1,nspinor
1048 [ # # ]: 0 : do ia=1,nattyp(itypat)
1049 [ # # # # : 0 : cprjin(iatm+ia,ispinor)%dcp(:,1+ishift:ndgxdt+ishift,1:nlmn)=zero
# # ]
1050 : : ! cprjin(iatm+ia,ispinor)%dcp(:,1:ndgxdt,1:nlmn)=zero
1051 : : end do
1052 : : end do
1053 : : end if
1054 : 1344 : iatm=iatm+nattyp(itypat)
1055 : : end if
1056 : :
1057 : : ! End atom type loop
1058 : 65888006 : ia1=ia2+1
1059 : : end do
1060 : :
1061 : :
1062 : : !Reduction in case of parallelism
1063 : : !==============================================================
1064 : :
1065 [ + + + + ]: 25801672 : if (signs==1.and.mpi_enreg%paral_spinor==1) then
1066 [ + + ]: 10872 : if (nnlout/=0) then
1067 : 1304 : call xmpi_sum(enlout,mpi_enreg%comm_spinor,ierr)
1068 : : end if
1069 [ + + ]: 10872 : if (choice==3.or.choice==6.or.choice==23) then
1070 : 1304 : call xmpi_sum(enlk,mpi_enreg%comm_spinor,ierr)
1071 : : end if
1072 [ - + ]: 10872 : if (choice==6) then
1073 : 0 : call xmpi_sum(fnlk,mpi_enreg%comm_spinor,ierr)
1074 : 0 : call xmpi_sum(strnlk,mpi_enreg%comm_spinor,ierr)
1075 : : end if
1076 [ - + ]: 10872 : if (choice==55) then
1077 : 0 : call xmpi_sum(ddkk,mpi_enreg%comm_spinor,ierr)
1078 : : end if
1079 : : end if
1080 : :
1081 : : !Coordinate transformations
1082 : : !==============================================================
1083 : :
1084 : : !Need sometimes gmet
1085 [ + + - + : 25801672 : if ((signs==1.and.paw_opt<=3).and. &
+ + + - +
- + + +
+ ]
1086 : : & (choice==5 .or.choice==51.or.choice==52.or.choice==53.or.&
1087 : : & choice==54.or.choice==55)) then
1088 [ - + ]: 15120 : ABI_MALLOC(gmet,(3,3))
1089 [ - + - - : 26587912 : gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
- - + + +
+ + + + +
+ + ]
1090 : : end if
1091 : :
1092 : : !1st derivative wrt to strain (stress tensor):
1093 : : ! - convert from reduced to cartesian coordinates
1094 : : ! - substract volume contribution
1095 [ + + + + : 25801672 : if ((choice==3.or.choice==23).and.signs==1.and.paw_opt<=3) then
+ - ]
1096 : 70367 : mu0=0 ! Shift to be applied in enlout array
1097 [ + - ]: 70367 : ABI_MALLOC(work1,(6))
1098 [ + + ]: 492569 : work1(1:6)=enlout(mu0+1:mu0+6)
1099 : 70367 : call strconv(work1,gprimd,work1)
1100 [ + + ]: 281468 : enlout(mu0+1:mu0+3)=(work1(1:3)-enlk)
1101 [ + + ]: 281468 : enlout(mu0+4:mu0+6)= work1(4:6)
1102 [ - + ]: 70367 : ABI_FREE(work1)
1103 : : end if
1104 : :
1105 : : !1st derivative wrt to k wave vector (ddk):
1106 : : ! - convert from cartesian to reduced coordinates
1107 [ + + + + : 25801672 : if ((choice==5.or.choice==53).and.signs==1.and.paw_opt<=3) then
+ - ]
1108 : 1152 : mu0=0 ! Shift to be applied in enlout array
1109 [ + - ]: 1152 : ABI_MALLOC(work1,(3))
1110 [ + + ]: 4608 : work1(:)=enlout(mu0+1:mu0+3)
1111 [ + + ]: 4608 : enlout(mu0+1:mu0+3)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3)
1112 : 1152 : ABI_FREE(work1)
1113 : : end if
1114 [ + + - + : 25801672 : if ((choice==51.or.choice==52).and.signs==1.and.paw_opt<=3) then
- - ]
1115 : 0 : mu0=0 ! Shift to be applied in enlout array
1116 [ # # ]: 0 : ABI_MALLOC(work1,(3))
1117 [ # # ]: 0 : do mu=1,2 ! Loop for Re,Im
1118 [ # # ]: 0 : work1(1:3)=(/enlout(mu0+1),enlout(mu0+3),enlout(mu0+5)/)
1119 : 0 : enlout(mu0+1)=gmet(1,1)*work1(1)+gmet(1,2)*work1(2)+gmet(1,3)*work1(3)
1120 : 0 : enlout(mu0+3)=gmet(2,1)*work1(1)+gmet(2,2)*work1(2)+gmet(2,3)*work1(3)
1121 : 0 : enlout(mu0+5)=gmet(3,1)*work1(1)+gmet(3,2)*work1(2)+gmet(3,3)*work1(3)
1122 : 0 : mu0=mu0+1
1123 : : end do
1124 [ # # ]: 0 : ABI_FREE(work1)
1125 : : end if
1126 : :
1127 : : !2nd derivative wrt to k wave vector and atomic position (effective charges):
1128 : : ! - convert from cartesian to reduced coordinates
1129 [ + + + + : 25801672 : if (choice==54.and.signs==1.and.paw_opt<=3) then
+ - ]
1130 : 11392 : mu0=0 ! Shift to be applied in enlout array
1131 [ - + ]: 11392 : ABI_MALLOC(work1,(3))
1132 [ - + ]: 11392 : ABI_MALLOC(work2,(3))
1133 [ + + ]: 81664 : do mu=1,3*natom
1134 : : ! First, real part
1135 : 70272 : work1(1)=enlout(mu0+1);work1(2)=enlout(mu0+3);work1(3)=enlout(mu0+5)
1136 [ + + ]: 281088 : work2(:)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3)
1137 : 70272 : enlout(mu0+1)=work2(1);enlout(mu0+3)=work2(2);enlout(mu0+5)=work2(3)
1138 : : ! Then imaginary part
1139 : 70272 : work1(1)=enlout(mu0+2);work1(2)=enlout(mu0+4);work1(3)=enlout(mu0+6)
1140 [ + + ]: 281088 : work2(:)=gmet(:,1)*work1(1)+gmet(:,2)*work1(2)+gmet(:,3)*work1(3)
1141 : 70272 : enlout(mu0+2)=work2(1);enlout(mu0+4)=work2(2);enlout(mu0+6)=work2(3)
1142 : 81664 : mu0=mu0+6
1143 : : end do
1144 [ - + ]: 11392 : ABI_FREE(work1)
1145 [ - + ]: 11392 : ABI_FREE(work2)
1146 : : end if
1147 : :
1148 : : !2nd derivative wrt to k wave vector and strain (piezoelectric tensor):
1149 : : ! - convert from cartesian to reduced coordinates (k point)
1150 : : ! - convert from reduced to cartesian coordinates (strain)
1151 : : ! - substract volume contribution
1152 : : ! - symetrize strain components
1153 [ + + + - : 25801672 : if (choice==55.and.signs==1.and.paw_opt<=3) then
+ - ]
1154 [ - + ]: 2576 : ABI_MALLOC(work3,(2,3))
1155 [ - + ]: 2576 : ABI_MALLOC(work4,(2,3))
1156 [ - + ]: 2576 : ABI_MALLOC(work5,(2,3,6))
1157 [ - + ]: 2576 : ABI_MALLOC(work7,(2,3,6))
1158 [ + - ]: 2576 : ABI_MALLOC(work6,(2,3,3))
1159 [ + + ]: 10304 : do ic=1,3 ! gamma
1160 [ + + + + : 471408 : work5=zero
+ + ]
1161 [ + + ]: 30912 : do jc=1,3 ! nu
1162 [ + + ]: 100464 : do ii=1,3 ! lambda
1163 : 69552 : mu=(gamma(jc,ii)-1)*3+1
1164 : : work5(1,jc,ii)=gmet(ic,1)*enlout(2*mu-1)+gmet(ic,2)*enlout(2*mu+1) &
1165 : 69552 : & +gmet(ic,3)*enlout(2*mu+3)
1166 : : work5(2,jc,ii)=gmet(ic,1)*enlout(2*mu )+gmet(ic,2)*enlout(2*mu+2) &
1167 : 92736 : & +gmet(ic,3)*enlout(2*mu+4)
1168 : : end do
1169 : : end do
1170 [ + + + + : 239568 : work6=zero
+ + ]
1171 [ + + ]: 30912 : do jc=1,3 ! nu
1172 [ + + ]: 100464 : do ii=1,3 ! beta
1173 : : work6(1:cplex,ii,jc)=gprimd(ii,1)*work5(1:cplex,jc,1)+gprimd(ii,2)*work5(1:cplex,jc,2) &
1174 [ + + ]: 231840 : & +gprimd(ii,3)*work5(1:cplex,jc,3)
1175 : : end do
1176 : : end do
1177 [ + + ]: 33488 : do jc=1,3 ! alpha
1178 [ + + ]: 100464 : do ii=1,3 ! beta
1179 : 69552 : mu=gamma(jc,ii)
1180 : : work7(1:cplex,ic,mu)=gprimd(jc,1)*work6(1:cplex,ii,1)+gprimd(jc,2)*work6(1:cplex,ii,2) &
1181 [ + + ]: 231840 : & +gprimd(jc,3)*work6(1:cplex,ii,3)
1182 : : end do
1183 : : end do
1184 : : end do ! gamma
1185 : :
1186 [ + + ]: 10304 : do ii=1,3 ! alpha
1187 : : work3(1,ii)=gprimd(ii,1)*ddkk(2*1-1)+gprimd(ii,2)*ddkk(2*2-1) &
1188 : 7728 : & +gprimd(ii,3)*ddkk(2*3-1)
1189 : : work3(2,ii)=gprimd(ii,1)*ddkk(2*1 )+gprimd(ii,2)*ddkk(2*2 ) &
1190 : 10304 : & +gprimd(ii,3)*ddkk(2*3 )
1191 : : end do
1192 [ + + ]: 10304 : do ii=1,3 ! gamma
1193 : : work4(1,ii)=gmet(ii,1)*ddkk(2*1-1)+gmet(ii,2)*ddkk(2*2-1) &
1194 : 7728 : & +gmet(ii,3)*ddkk(2*3-1)
1195 : : work4(2,ii)=gmet(ii,1)*ddkk(2*1 )+gmet(ii,2)*ddkk(2*2 ) &
1196 : 10304 : & +gmet(ii,3)*ddkk(2*3 )
1197 : : end do
1198 : :
1199 [ + + ]: 18032 : do mu=1,6
1200 : 15456 : ii=alpha(mu) ! alpha
1201 : 15456 : ic=beta(mu) ! beta
1202 [ + + ]: 64400 : do jc=1,3 ! gamma
1203 : : work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-half &
1204 [ + + ]: 139104 : & *(gprimd(ic,jc)*work3(1:cplex,ii)+gprimd(ii,jc)*work3(1:cplex,ic))
1205 [ + + + + ]: 108192 : if (ii==ic) work7(1:cplex,jc,mu)=work7(1:cplex,jc,mu)-work4(1:cplex,jc)
1206 : : end do
1207 : : end do
1208 [ + + ]: 18032 : do mu=1,6 ! alpha,beta
1209 [ + + ]: 64400 : do nu=1,3 ! gamma
1210 : 46368 : mu0=3*(mu-1)+nu
1211 : 46368 : enlout(2*mu0-1)=work7(1,nu,mu)
1212 : 61824 : enlout(2*mu0 )=work7(2,nu,mu)
1213 : : end do
1214 : : end do
1215 [ - + ]: 2576 : ABI_FREE(gmet)
1216 [ - + ]: 2576 : ABI_FREE(work3)
1217 [ - + ]: 2576 : ABI_FREE(work4)
1218 [ - + ]: 2576 : ABI_FREE(work5)
1219 [ - + ]: 2576 : ABI_FREE(work6)
1220 [ - + ]: 2576 : ABI_FREE(work7)
1221 : : end if
1222 : :
1223 : : !2nd derivative wrt to 2 k wave vectors (effective mass):
1224 : : ! - convert from cartesian to reduced coordinates
1225 [ + + + + : 25801672 : if ((choice==8.or.choice==81).and.signs==1.and.paw_opt<=3) then
+ - ]
1226 : 960 : mu0=0 ! Shift to be applied in enlout array
1227 [ - + ]: 960 : ABI_MALLOC(work3,(3,3))
1228 [ - + ]: 960 : ABI_MALLOC(work4,(3,3))
1229 [ + + ]: 960 : mua=1;if (choice==81) mua=2
1230 [ + + ]: 2496 : do ii=1,mua ! Loop Re,Im
1231 [ + + ]: 1536 : if (choice==8) then ! enlout is real in Voigt notation
1232 : 384 : work3(1,1)=enlout(mu0+1) ; work3(1,2)=enlout(mu0+6) ; work3(1,3)=enlout(mu0+5)
1233 : 384 : work3(2,1)=enlout(mu0+6) ; work3(2,2)=enlout(mu0+2) ; work3(2,3)=enlout(mu0+4)
1234 : 384 : work3(3,1)=enlout(mu0+5) ; work3(3,2)=enlout(mu0+4) ; work3(3,3)=enlout(mu0+3)
1235 : : else ! enlout is complex in matrix notation
1236 : 1152 : work3(1,1)=enlout(mu0+1 ) ; work3(1,2)=enlout(mu0+3 ) ; work3(1,3)=enlout(mu0+5 )
1237 : 1152 : work3(2,1)=enlout(mu0+7 ) ; work3(2,2)=enlout(mu0+9 ) ; work3(2,3)=enlout(mu0+11)
1238 : 1152 : work3(3,1)=enlout(mu0+13) ; work3(3,2)=enlout(mu0+15) ; work3(3,3)=enlout(mu0+17)
1239 : : end if
1240 [ + + ]: 6144 : do mu=1,3
1241 [ + + ]: 19968 : work4(:,mu)=gprimd(:,1)*work3(mu,1)+gprimd(:,2)*work3(mu,2)+gprimd(:,3)*work3(mu,3)
1242 : : end do
1243 [ + + ]: 6144 : do mu=1,3
1244 [ + + ]: 19968 : work3(:,mu)=gprimd(:,1)*work4(mu,1)+gprimd(:,2)*work4(mu,2)+gprimd(:,3)*work4(mu,3)
1245 : : end do
1246 [ + + ]: 6144 : do mu=1,3
1247 [ + + ]: 19968 : work4(:,mu)=gprimd(1,:)*work3(mu,1)+gprimd(2,:)*work3(mu,2)+gprimd(3,:)*work3(mu,3)
1248 : : end do
1249 [ + + ]: 6144 : do mu=1,3
1250 [ + + ]: 19968 : work3(:,mu)=gprimd(1,:)*work4(mu,1)+gprimd(2,:)*work4(mu,2)+gprimd(3,:)*work4(mu,3)
1251 : : end do
1252 [ + + ]: 1536 : if (choice==8) then ! enlout is real in Voigt notation
1253 : 384 : enlout(mu0+1) = work3(1,1) ; enlout(mu0+2) = work3(2,2) ; enlout(mu0+3) = work3(3,3)
1254 : 384 : enlout(mu0+4) = work3(3,2) ; enlout(mu0+5) = work3(1,3) ; enlout(mu0+6) = work3(2,1)
1255 : : else ! enlout is complex in matrix notation
1256 : 1152 : enlout(mu0+1 )=work3(1,1) ; enlout(mu0+3 )=work3(1,2) ; enlout(mu0+5 )=work3(1,3)
1257 : 1152 : enlout(mu0+7 )=work3(2,1) ; enlout(mu0+9 )=work3(2,2) ; enlout(mu0+11)=work3(2,3)
1258 : 1152 : enlout(mu0+13)=work3(3,1) ; enlout(mu0+15)=work3(3,2) ; enlout(mu0+17)=work3(3,3)
1259 : : end if
1260 : 2496 : mu0=mu0+1
1261 : : end do
1262 [ - + ]: 960 : ABI_FREE(work3)
1263 [ - + ]: 960 : ABI_FREE(work4)
1264 : : end if
1265 : :
1266 : : !2nd derivative wrt to 2 strains (elastic tensor):
1267 : : ! - convert from reduced to cartesian coordinates
1268 : : ! - substract volume contribution
1269 [ + + + - : 25801672 : if (choice==6.and.signs==1.and.paw_opt<=3) then
+ - ]
1270 : 2400 : mu0=0 ! Shift to be applied in enlout array
1271 [ - + ]: 2400 : ABI_MALLOC(work1,(6))
1272 [ - + ]: 2400 : ABI_MALLOC(work2,(6))
1273 [ + - - + : 9600 : ABI_MALLOC(work3,(6+3*natom,6))
+ - + - +
- + - - +
- + ]
1274 [ + + ]: 7200 : work3(:,:)=reshape(enlout(mu0+1:mu0+6*(6+3*natom)),(/6+3*natom,6/))
1275 [ + + ]: 16800 : do mu=1,6
1276 : 16800 : call strconv(work3(1:6,mu),gprimd,work3(1:6,mu))
1277 : : end do
1278 [ + + ]: 32928 : do mu=1,6+3*natom
1279 [ + + ]: 213696 : work1(1:6)=work3(mu,1:6)
1280 : 30528 : call strconv(work1,gprimd,work2)
1281 [ + + ]: 216096 : work3(mu,1:6)=work2(1:6)
1282 : : end do
1283 [ + + ]: 4800 : enlout(mu0+1:mu0+6*(6+3*natom))=reshape(work3(:,:),(/6*(6+3*natom)/))
1284 [ - + ]: 2400 : ABI_FREE(work1)
1285 [ - + ]: 2400 : ABI_FREE(work2)
1286 [ - + ]: 2400 : ABI_FREE(work3)
1287 : 2400 : call strconv(strnlk,gprimd,strnlk)
1288 [ + + ]: 16800 : do mub=1,6
1289 : 14400 : nub1=alpha(mub);nub2=beta(mub)
1290 [ + + ]: 100800 : do mua=1,6
1291 : 86400 : mu=mu0+mua+(3*natom+6)*(mub-1)
1292 : 86400 : nua1=alpha(mua);nua2=beta(mua)
1293 [ + + ]: 86400 : if (mua<=3.and.mub<=3) enlout(mu)=enlout(mu)+enlk
1294 [ + + ]: 86400 : if (mua<=3) enlout(mu)=enlout(mu)-strnlk(mub)
1295 [ + + ]: 86400 : if (mub<=3) enlout(mu)=enlout(mu)-strnlk(mua)
1296 [ + + ]: 86400 : if (nub1==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub2))
1297 [ + + ]: 86400 : if (nub2==nua2) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua1,nub1))
1298 [ + + ]: 86400 : if (nub1==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub2))
1299 [ + + ]: 100800 : if (nub2==nua1) enlout(mu)=enlout(mu)-0.25d0*strnlk(gamma(nua2,nub1))
1300 : : end do
1301 [ + + ]: 16800 : if (mub<=3) then
1302 [ + + ]: 23328 : do nua1=1,natom
1303 : 16128 : nua2=3*(nua1-1);mu=mu0+nua2+6+(3*natom+6)*(mub-1)
1304 [ + + ]: 71712 : enlout(mu+1:mu+3)=enlout(mu+1:mu+3)-fnlk(nua2+1:nua2+3)
1305 : : end do
1306 : : end if
1307 : : end do
1308 : : end if
1309 : :
1310 [ + + ]: 25801672 : if (allocated(gmet)) then
1311 : 12544 : ABI_FREE(gmet)
1312 : : end if
1313 : :
1314 : : !Final deallocations
1315 : : !==============================================================
1316 : :
1317 [ + + ]: 25801672 : if (signs==1) then
1318 [ - + ]: 1410081 : ABI_FREE(fnlk)
1319 [ - + ]: 1410081 : ABI_FREE(ddkk)
1320 [ - + ]: 1410081 : ABI_FREE(strnlk)
1321 : : end if
1322 : :
1323 [ + + ]: 25801672 : if (nkpgin<nkpgin_) then
1324 [ - + ]: 1892207 : ABI_FREE(kpgin_)
1325 : : end if
1326 [ + + ]: 25801672 : if (nkpgout<nkpgout_) then
1327 [ - + ]: 1968860 : ABI_FREE(kpgout_)
1328 : : end if
1329 : :
1330 : : ! call timab(1100,2,tsec)
1331 : :
1332 : : DBG_EXIT("COLL")
1333 : :
1334 [ + + - + : 51603344 : end subroutine nonlop_ylm
- + - + ]
1335 : : !!***
1336 : :
1337 : : !!****f* ABINIT/nonlop_ylm_init_counters
1338 : : !! NAME
1339 : : !! nonlop_ylm_init_counters
1340 : : !!
1341 : : !! FUNCTION
1342 : : !!
1343 : : !! SOURCE
1344 : :
1345 : 0 : subroutine nonlop_ylm_init_counters()
1346 : :
1347 : 0 : opernla_counter = 0
1348 : 0 : opernlb_counter = 0
1349 : 0 : opernla_mv_counter = 0
1350 : 0 : opernlb_mv_counter = 0
1351 : 0 : opernla_mv_dgemv_counter = 0
1352 : 0 : opernlb_mv_dgemv_counter = 0
1353 : :
1354 : 0 : end subroutine nonlop_ylm_init_counters
1355 : : !!***
1356 : :
1357 : : !!****f* ABINIT/nonlop_ylm_stop_counters
1358 : : !! NAME
1359 : : !! nonlop_ylm_stop_counters
1360 : : !!
1361 : : !! FUNCTION
1362 : : !!
1363 : : !! SOURCE
1364 : :
1365 : 0 : subroutine nonlop_ylm_stop_counters()
1366 : :
1367 : 0 : opernla_counter = -1
1368 : 0 : opernlb_counter = -1
1369 : 0 : opernla_mv_counter = -1
1370 : 0 : opernlb_mv_counter = -1
1371 : 0 : opernla_mv_dgemv_counter = -1
1372 : 0 : opernlb_mv_dgemv_counter = -1
1373 : :
1374 : 0 : end subroutine nonlop_ylm_stop_counters
1375 : : !!***
1376 : :
1377 : : !!****f* ABINIT/nonlop_ylm_output_counters
1378 : : !! NAME
1379 : : !! nonlop_ylm_output_counters
1380 : : !!
1381 : : !! FUNCTION
1382 : : !!
1383 : : !! SOURCE
1384 : :
1385 [ # # ]: 0 : subroutine nonlop_ylm_output_counters(natom,nbandtot,ntypat,typat,mpi_enreg)
1386 : :
1387 : : !Arguments ------------------------------------
1388 : : !scalars
1389 : : integer,intent(in) :: natom,nbandtot,ntypat
1390 : : integer,intent(in) :: typat(:)
1391 : : type(MPI_type),intent(in) :: mpi_enreg
1392 : : !arrays
1393 : :
1394 : : !Local variables-------------------------------
1395 : : !scalars
1396 : : character(len=500) :: msg
1397 : : integer :: cnt,ia1,ia2,ia3,ia4,ia5,iatm,ierr,itypat,mincat,nincat,opernl_calls
1398 : : !arrays
1399 : 0 : integer :: nattyp(ntypat)
1400 : :
1401 [ # # ]: 0 : do itypat=1,ntypat
1402 : 0 : nattyp(itypat)=0
1403 [ # # ]: 0 : do iatm=1,natom
1404 [ # # ]: 0 : if(typat(iatm)==itypat)then
1405 : : ! atindx(iatom)=indx
1406 : : ! atindx1(indx)=iatom
1407 : : ! indx=indx+1
1408 : 0 : nattyp(itypat)=nattyp(itypat)+1
1409 : : end if
1410 : : end do
1411 : : end do
1412 [ # # ]: 0 : call wrtout([std_out,ab_out],'','COLL')
1413 : 0 : write(msg,'(a)') ' --- NONLOP YLM COUNTERS -----------------------------------------------------'
1414 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1415 [ # # ]: 0 : mincat=min(NLO_MINCAT,maxval(nattyp))
1416 : 0 : ia1=1;iatm=0;opernl_calls=0
1417 [ # # ]: 0 : do itypat=1,ntypat
1418 : : ! Get atom loop indices for different types:
1419 : 0 : ia2=ia1+nattyp(itypat)-1;ia5=1
1420 [ # # # # : 0 : do ia3=ia1,ia2,mincat
# # ]
1421 : 0 : ia4=min(ia2,ia3+mincat-1)
1422 : : ! Give the increment of number of atoms in this subset.
1423 : 0 : nincat=ia4-ia3+1
1424 : 0 : opernl_calls=opernl_calls+1
1425 : : ! End sum on atom subset loop
1426 [ # # ]: 0 : iatm=iatm+nincat;ia5=ia5+nincat
1427 : : end do
1428 : : ! End atom type loop
1429 : 0 : ia1=ia2+1
1430 : : end do
1431 [ # # ]: 0 : if (iatm/=natom) then
1432 : 0 : ABI_ERROR('iatm should be equal to natom!')
1433 : : end if
1434 : 0 : write(msg,'(a,i6)') ' Number of Calls in nonlop_ylm : NC = ',opernl_calls
1435 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1436 : 0 : write(msg,'(a,i6)') ' total Number of Bands : NB = ',nbandtot
1437 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1438 : 0 : write(msg,'(a)') ' | total count (TC) | TC/NC | TC/NC/NB'
1439 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1440 : 0 : write(msg,'(a)') ' -----------------------------------------------------------------------------'
1441 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1442 : 0 : call xmpi_sum(opernla_counter,mpi_enreg%comm_kpt,ierr)
1443 : 0 : call xmpi_sum(opernlb_counter,mpi_enreg%comm_kpt,ierr)
1444 : 0 : call xmpi_sum(opernla_mv_counter,mpi_enreg%comm_kpt,ierr)
1445 : 0 : call xmpi_sum(opernlb_mv_counter,mpi_enreg%comm_kpt,ierr)
1446 : 0 : call xmpi_sum(opernla_mv_dgemv_counter,mpi_enreg%comm_kpt,ierr)
1447 : 0 : call xmpi_sum(opernlb_mv_dgemv_counter,mpi_enreg%comm_kpt,ierr)
1448 : 0 : cnt=opernla_counter
1449 [ # # ]: 0 : if (cnt>0) then
1450 : 0 : write(msg,'(2(a,i16),a,f16.1)') '_opernla_ylm | ',&
1451 : 0 : & cnt,' | ',cnt/opernl_calls,' | ',dble(cnt)/opernl_calls/nbandtot
1452 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1453 : : end if
1454 : 0 : cnt=opernla_mv_counter
1455 [ # # ]: 0 : if (cnt>0) then
1456 : 0 : write(msg,'(2(a,i16),a,f16.1)') '_opernla_ylm_mv | ',&
1457 : 0 : & cnt,' | ',cnt/opernl_calls,' | ',dble(cnt)/opernl_calls/nbandtot
1458 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1459 : : end if
1460 : 0 : cnt=opernla_mv_dgemv_counter
1461 [ # # ]: 0 : if (cnt>0) then
1462 : 0 : write(msg,'(2(a,i16),a,f16.1)') '_opernla_ylm_mv(dgemv)| ',&
1463 : 0 : & cnt,' | ',cnt/opernl_calls,' | ',dble(cnt)/opernl_calls/nbandtot
1464 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1465 : : end if
1466 : 0 : cnt=opernlb_counter
1467 [ # # ]: 0 : if (cnt>0) then
1468 : 0 : write(msg,'(2(a,i16),a,f16.1)') '_opernlb_ylm | ',&
1469 : 0 : & cnt,' | ',cnt/opernl_calls,' | ',dble(cnt)/opernl_calls/nbandtot
1470 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1471 : : end if
1472 : 0 : cnt=opernlb_mv_counter
1473 [ # # ]: 0 : if (cnt>0) then
1474 : 0 : write(msg,'(2(a,i16),a,f16.1)') '_opernlb_ylm_mv | ',&
1475 : 0 : & cnt,' | ',cnt/opernl_calls,' | ',dble(cnt)/opernl_calls/nbandtot
1476 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1477 : : end if
1478 : 0 : cnt=opernlb_mv_dgemv_counter
1479 [ # # ]: 0 : if (cnt>0) then
1480 : 0 : write(msg,'(2(a,i16),a,f16.1)') '_opernlb_ylm_mv(dgemv)| ',&
1481 : 0 : & cnt,' | ',cnt/opernl_calls,' | ',dble(cnt)/opernl_calls/nbandtot
1482 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1483 : : end if
1484 : 0 : write(msg,'(a)') ' -----------------------------------------------------------------------------'
1485 [ # # ]: 0 : call wrtout([std_out,ab_out],msg,'COLL')
1486 : :
1487 : 0 : end subroutine nonlop_ylm_output_counters
1488 : : !!***
1489 : :
1490 : : end module m_nonlop_ylm
1491 : : !!***
|
