Branch data Line data Source code
1 : : !!****m* ABINIT/m_mkffkg
2 : : !! NAME
3 : : !! m_mkffkg
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, MT, DRH)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! SOURCE
15 : :
16 : : #if defined HAVE_CONFIG_H
17 : : #include "config.h"
18 : : #endif
19 : :
20 : : #include "abi_common.h"
21 : :
22 : : module m_mkffkg
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : :
27 : : implicit none
28 : :
29 : : private
30 : : !!***
31 : :
32 : : public :: mkffkg
33 : : public :: dfpt_mkffkg
34 : : !!***
35 : :
36 : : contains
37 : : !!***
38 : :
39 : : !!****f* ABINIT/dfpt_mkffkg
40 : : !! NAME
41 : : !! dfpt_mkffkg
42 : : !!
43 : : !! FUNCTION
44 : : !! Prepare the application of the projectors to the shifted wavefunctions,
45 : : !! by precomputing the k+G factors and their product with the form factors
46 : : !! Do this on a block of plane waves.
47 : : !!
48 : : !! INPUTS
49 : : !! choice=governs the combination of k+G vectors to be computed
50 : : !! ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
51 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
52 : : !! nffnl=3rd dimension of ffnl(2, conventional, or 3 for 2nd derivatives)
53 : : !! idir=direction of the perturbation (needed if choice==2 and ndgxdt==1,
54 : : !! or if choice==5)
55 : : !! indlmn(6,i,ntypat)=array giving l,m,n,lm,ln,spin for i=ln
56 : : !! ipw1 = index of the first plane wave treated in this block
57 : : !! ispinor=1 or 2, gives the spinorial component of ffnl to be used
58 : : !! itypat = type of atom, needed for ffnl
59 : : !! kg_k(3,npw)=integer coords of planewaves in basis sphere
60 : : !! kpg_k(npw,npkg)= (k+G) components and related data
61 : : !! kpt(3)=real components of k point in terms of recip. translations
62 : : !! lmnmax=max. number of (l,n) components over all type of psps
63 : : !! mblkpw=first dimension of kpgx
64 : : !! ndgxdt=number of components of first order derivative
65 : : !! nffkg=number of products of ffnls with combinations of k+G
66 : : !! nincpw=number of plane waves in the block
67 : : !! nkpg=second size of array kpg_k
68 : : !! nlang = number of angular momenta to be treated = 1 + highest ang. mom.
69 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
70 : : !! npw = total number of plane waves in reciprocal space
71 : : !! ntens=second dimension of kpgx, number of distinct tensorial products
72 : : !! ntypat = number of type of atoms, dimension needed for ffnl
73 : : !!
74 : : !! OUTPUT
75 : : !! kpgx(mblkpw,ntens)=different tensorial products of k+G
76 : : !! ffkg(nffkg,mblkpw)=different products of ffnls with k+G
77 : : !! parity(nffkg)=parity of the tensorial product of k+G (2 if even, 1 of odd)
78 : : !!
79 : : !! NOTES
80 : : !! This routine must be thread-safe as it is called inside loops that are OpenMP parallelized.
81 : : !! Please, do not add variables with the save attribute or SIDE EFFECTS.
82 : : !!
83 : : !! SOURCE
84 : :
85 : 185127742 : subroutine dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
86 : 185127742 : & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
87 : 185127742 : & npw,ntens,ntypat,parity)
88 : :
89 : : !Arguments ------------------------------------
90 : : !scalars
91 : : integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt
92 : : integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat
93 : : !arrays
94 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw)
95 : : integer,intent(out) :: parity(nffkg)
96 : : real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
97 : : real(dp),intent(in) :: kpt(3)
98 : : real(dp),intent(out) :: ffkg(nffkg,mblkpw),kpgx(mblkpw,ntens)
99 : :
100 : : !Local variables-------------------------------
101 : : !scalars
102 : : integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj
103 : : integer :: nffkge
104 : : real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z
105 : :
106 : : ! *************************************************************************
107 : :
108 : 185127742 : jj=0;ilangx=0
109 : :
110 : : !This will be useless after all the modifications have been done
111 [ + + ]: 27293488195 : do ipw=1,nincpw
112 : 27293488195 : kpgx(ipw,1)=1.0d0
113 : : end do
114 : :
115 : : !Initialize kpgx array related to tensors defined below
116 : : if ( nlang>=2 .or. choice==2 .or. choice==3 .or. choice==4 .or. choice==5&
117 [ + + + + : 185127742 : & .or. choice==6 .or. choice==23) then
+ + + + +
+ + + +
+ ]
118 [ - + ]: 164914888 : if (nkpg>=3) then
119 [ # # ]: 0 : kpgx(1:nincpw,2)=kpg_k(ipw1+1:ipw1+nincpw,1)
120 [ # # ]: 0 : kpgx(1:nincpw,3)=kpg_k(ipw1+1:ipw1+nincpw,2)
121 [ # # ]: 0 : kpgx(1:nincpw,4)=kpg_k(ipw1+1:ipw1+nincpw,3)
122 : : else
123 : 164914888 : ig=ipw1
124 [ + + ]: 24623102307 : do ipw=1,nincpw
125 : 24458187419 : kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig))
126 : 24458187419 : kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig))
127 : 24458187419 : kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig))
128 : 24643315161 : ig=ig+1
129 : : end do
130 : : end if
131 : : end if
132 [ + + + + ]: 185127742 : if (nlang>=3 .or. choice==3 .or. choice==6 .or. choice==23) then
133 : : ! Define (k+G) part of rank 2 symmetric tensor (6 components), l=2
134 : : ! Compressed storage is 11 22 33 32 31 21
135 [ + - ]: 96961528 : if (nkpg>=9) then
136 [ # # ]: 0 : kpgx(1:nincpw,5) =kpg_k(ipw1+1:ipw1+nincpw,4)
137 [ # # ]: 0 : kpgx(1:nincpw,6) =kpg_k(ipw1+1:ipw1+nincpw,5)
138 [ # # ]: 0 : kpgx(1:nincpw,7) =kpg_k(ipw1+1:ipw1+nincpw,6)
139 [ # # ]: 0 : kpgx(1:nincpw,8) =kpg_k(ipw1+1:ipw1+nincpw,7)
140 [ # # ]: 0 : kpgx(1:nincpw,9) =kpg_k(ipw1+1:ipw1+nincpw,8)
141 [ # # ]: 0 : kpgx(1:nincpw,10)=kpg_k(ipw1+1:ipw1+nincpw,9)
142 : : else
143 [ + + ]: 14154800631 : do ipw=1,nincpw
144 : 14057839103 : kpgx(ipw, 5) = kpgx(ipw, 2)*kpgx(ipw, 2)
145 : 14057839103 : kpgx(ipw, 6) = kpgx(ipw, 3)*kpgx(ipw, 3)
146 : 14057839103 : kpgx(ipw, 7) = kpgx(ipw, 4)*kpgx(ipw, 4)
147 : 14057839103 : kpgx(ipw, 8) = kpgx(ipw, 4)*kpgx(ipw, 3)
148 : 14057839103 : kpgx(ipw, 9) = kpgx(ipw, 4)*kpgx(ipw, 2)
149 : 14154800631 : kpgx(ipw,10) = kpgx(ipw, 3)*kpgx(ipw, 2)
150 : : end do
151 : : end if
152 : : end if
153 [ + + + + : 185127742 : if (nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) .or. choice==6) then
+ + + + ]
154 : : ! Define (k+G) part of rank 3 symmetric tensor (10 components), l=3
155 : : ! Compressed storage is 111 221 331 321 311 211 222 332 322 333
156 [ + + ]: 488131196 : do ipw=1,nincpw
157 : 484214506 : kpgx(ipw,11) = kpgx(ipw, 5)*kpgx(ipw, 2)
158 : 484214506 : kpgx(ipw,12) = kpgx(ipw, 6)*kpgx(ipw, 2)
159 : 484214506 : kpgx(ipw,13) = kpgx(ipw, 7)*kpgx(ipw, 2)
160 : 484214506 : kpgx(ipw,14) = kpgx(ipw, 8)*kpgx(ipw, 2)
161 : 484214506 : kpgx(ipw,15) = kpgx(ipw, 9)*kpgx(ipw, 2)
162 : 484214506 : kpgx(ipw,16) = kpgx(ipw,10)*kpgx(ipw, 2)
163 : 484214506 : kpgx(ipw,17) = kpgx(ipw, 6)*kpgx(ipw, 3)
164 : 484214506 : kpgx(ipw,18) = kpgx(ipw, 7)*kpgx(ipw, 3)
165 : 484214506 : kpgx(ipw,19) = kpgx(ipw, 8)*kpgx(ipw, 3)
166 : 669342248 : kpgx(ipw,20) = kpgx(ipw, 7)*kpgx(ipw, 4)
167 : : end do
168 : : end if
169 [ + + + + : 185127742 : if (((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6) then
+ + ]
170 : : ! Add additional tensors for strain gradients
171 : : ! Define (k+G) part of rank 4 symmetric tensor (15 components), l=2
172 : : ! Compressed storage is 1111 2211 3311 3211 3111 2111 2221 3321 3221
173 : : ! 3331 2222 3322 3222 3332 3333
174 [ + + ]: 172388838 : do ipw=1,nincpw
175 : 170519538 : kpgx(ipw,21) = kpgx(ipw, 5)*kpgx(ipw, 5)
176 : 170519538 : kpgx(ipw,22) = kpgx(ipw, 6)*kpgx(ipw, 5)
177 : 170519538 : kpgx(ipw,23) = kpgx(ipw, 7)*kpgx(ipw, 5)
178 : 170519538 : kpgx(ipw,24) = kpgx(ipw, 8)*kpgx(ipw, 5)
179 : 170519538 : kpgx(ipw,25) = kpgx(ipw, 9)*kpgx(ipw, 5)
180 : 170519538 : kpgx(ipw,26) = kpgx(ipw,10)*kpgx(ipw, 5)
181 : 170519538 : kpgx(ipw,27) = kpgx(ipw, 6)*kpgx(ipw,10)
182 : 170519538 : kpgx(ipw,28) = kpgx(ipw, 7)*kpgx(ipw,10)
183 : 170519538 : kpgx(ipw,29) = kpgx(ipw, 8)*kpgx(ipw,10)
184 : 170519538 : kpgx(ipw,30) = kpgx(ipw, 7)*kpgx(ipw, 9)
185 : 170519538 : kpgx(ipw,31) = kpgx(ipw, 6)*kpgx(ipw, 6)
186 : 170519538 : kpgx(ipw,32) = kpgx(ipw, 7)*kpgx(ipw, 6)
187 : 170519538 : kpgx(ipw,33) = kpgx(ipw, 8)*kpgx(ipw, 6)
188 : 170519538 : kpgx(ipw,34) = kpgx(ipw, 7)*kpgx(ipw, 8)
189 : 355647280 : kpgx(ipw,35) = kpgx(ipw, 7)*kpgx(ipw, 7)
190 : : end do
191 : : end if
192 [ + + + + : 185127742 : if (((choice==3.or.choice==23) .and. nlang>=4) .or. (choice==6 .and. nlang>=2)) then
+ + + + ]
193 : : ! Define (k+G) part of rank 5 symmetric tensor (21 components), l=3
194 : : ! Compressed storage is 11111 22111 33111 32111 31111 21111
195 : : ! 22211 33211 32211 33311 22221 33221 32221 33321 33331
196 : : ! 22222 33222 32222 33322 33332 33333
197 [ + + ]: 10482041 : do ipw=1,nincpw
198 : 10403476 : kpgx(ipw,36) = kpgx(ipw,21)*kpgx(ipw, 2)
199 : 10403476 : kpgx(ipw,37) = kpgx(ipw,22)*kpgx(ipw, 2)
200 : 10403476 : kpgx(ipw,38) = kpgx(ipw,23)*kpgx(ipw, 2)
201 : 10403476 : kpgx(ipw,39) = kpgx(ipw,24)*kpgx(ipw, 2)
202 : 10403476 : kpgx(ipw,40) = kpgx(ipw,25)*kpgx(ipw, 2)
203 : 10403476 : kpgx(ipw,41) = kpgx(ipw,26)*kpgx(ipw, 2)
204 : 10403476 : kpgx(ipw,42) = kpgx(ipw,27)*kpgx(ipw, 2)
205 : 10403476 : kpgx(ipw,43) = kpgx(ipw,28)*kpgx(ipw, 2)
206 : 10403476 : kpgx(ipw,44) = kpgx(ipw,29)*kpgx(ipw, 2)
207 : 10403476 : kpgx(ipw,45) = kpgx(ipw,30)*kpgx(ipw, 2)
208 : 10403476 : kpgx(ipw,46) = kpgx(ipw,31)*kpgx(ipw, 2)
209 : 10403476 : kpgx(ipw,47) = kpgx(ipw,32)*kpgx(ipw, 2)
210 : 10403476 : kpgx(ipw,48) = kpgx(ipw,33)*kpgx(ipw, 2)
211 : 10403476 : kpgx(ipw,49) = kpgx(ipw,34)*kpgx(ipw, 2)
212 : 10403476 : kpgx(ipw,50) = kpgx(ipw,35)*kpgx(ipw, 2)
213 : 10403476 : kpgx(ipw,51) = kpgx(ipw,31)*kpgx(ipw, 3)
214 : 10403476 : kpgx(ipw,52) = kpgx(ipw,32)*kpgx(ipw, 3)
215 : 10403476 : kpgx(ipw,53) = kpgx(ipw,33)*kpgx(ipw, 3)
216 : 10403476 : kpgx(ipw,54) = kpgx(ipw,34)*kpgx(ipw, 3)
217 : 10403476 : kpgx(ipw,55) = kpgx(ipw,35)*kpgx(ipw, 3)
218 : 195531218 : kpgx(ipw,56) = kpgx(ipw,35)*kpgx(ipw, 4)
219 : : end do
220 : : end if
221 [ + + + + ]: 185127742 : if (choice==6 .and. nlang>=3) then
222 : : ! Define (k+G) part of rank 6 symmetric tensor (28 components)
223 : : ! Compressed storage is
224 : : ! 111111 221111 331111 321111 311111 211111 222111 332111 322111
225 : : ! 333111 222211 332211 322211 333211 333311 222221 332221 322221
226 : : ! 333221 333321 333331 222222 332222 322222 333222 333322 333332
227 : : ! 333333
228 [ + + ]: 3038944 : do ipw=1,nincpw
229 : 3012128 : kpgx(ipw,57) = kpgx(ipw,36)*kpgx(ipw, 2)
230 : 3012128 : kpgx(ipw,58) = kpgx(ipw,37)*kpgx(ipw, 2)
231 : 3012128 : kpgx(ipw,59) = kpgx(ipw,38)*kpgx(ipw, 2)
232 : 3012128 : kpgx(ipw,60) = kpgx(ipw,39)*kpgx(ipw, 2)
233 : 3012128 : kpgx(ipw,61) = kpgx(ipw,40)*kpgx(ipw, 2)
234 : 3012128 : kpgx(ipw,62) = kpgx(ipw,41)*kpgx(ipw, 2)
235 : 3012128 : kpgx(ipw,63) = kpgx(ipw,42)*kpgx(ipw, 2)
236 : 3012128 : kpgx(ipw,64) = kpgx(ipw,43)*kpgx(ipw, 2)
237 : 3012128 : kpgx(ipw,65) = kpgx(ipw,44)*kpgx(ipw, 2)
238 : 3012128 : kpgx(ipw,66) = kpgx(ipw,45)*kpgx(ipw, 2)
239 : 3012128 : kpgx(ipw,67) = kpgx(ipw,46)*kpgx(ipw, 2)
240 : 3012128 : kpgx(ipw,68) = kpgx(ipw,47)*kpgx(ipw, 2)
241 : 3012128 : kpgx(ipw,69) = kpgx(ipw,48)*kpgx(ipw, 2)
242 : 3012128 : kpgx(ipw,70) = kpgx(ipw,49)*kpgx(ipw, 2)
243 : 3012128 : kpgx(ipw,71) = kpgx(ipw,50)*kpgx(ipw, 2)
244 : 3012128 : kpgx(ipw,72) = kpgx(ipw,51)*kpgx(ipw, 2)
245 : 3012128 : kpgx(ipw,73) = kpgx(ipw,52)*kpgx(ipw, 2)
246 : 3012128 : kpgx(ipw,74) = kpgx(ipw,53)*kpgx(ipw, 2)
247 : 3012128 : kpgx(ipw,75) = kpgx(ipw,54)*kpgx(ipw, 2)
248 : 3012128 : kpgx(ipw,76) = kpgx(ipw,55)*kpgx(ipw, 2)
249 : 3012128 : kpgx(ipw,77) = kpgx(ipw,56)*kpgx(ipw, 2)
250 : 3012128 : kpgx(ipw,78) = kpgx(ipw,51)*kpgx(ipw, 3)
251 : 3012128 : kpgx(ipw,79) = kpgx(ipw,52)*kpgx(ipw, 3)
252 : 3012128 : kpgx(ipw,80) = kpgx(ipw,53)*kpgx(ipw, 3)
253 : 3012128 : kpgx(ipw,81) = kpgx(ipw,54)*kpgx(ipw, 3)
254 : 3012128 : kpgx(ipw,82) = kpgx(ipw,55)*kpgx(ipw, 3)
255 : 3012128 : kpgx(ipw,83) = kpgx(ipw,56)*kpgx(ipw, 3)
256 : 3038944 : kpgx(ipw,84) = kpgx(ipw,56)*kpgx(ipw, 4)
257 : : end do
258 : : end if
259 [ + + - + ]: 185127742 : if (choice==6 .and. nlang==4) then
260 : : ! Define (k+G) part of rank 7 symmetric tensor (36 components)
261 : : ! Compressed storage is
262 : : ! 1111111 2211111 3311111 3211111 3111111 2111111 2221111 3321111 3221111
263 : : ! 3331111 2222111 3322111 3222111 3332111 3333111 2222211 3322211 3222211
264 : : ! 3332211 3333211 3333311 2222221 3322221 3222221 3332221 3333221 3333321
265 : : ! 3333331 2222222 3322222 3222222 3332222 3333222 3333322 3333332 3333333
266 [ # # ]: 0 : do ipw=1,nincpw
267 : 0 : kpgx(ipw,85) = kpgx(ipw,57)*kpgx(ipw, 2)
268 : 0 : kpgx(ipw,86) = kpgx(ipw,58)*kpgx(ipw, 2)
269 : 0 : kpgx(ipw,87) = kpgx(ipw,59)*kpgx(ipw, 2)
270 : 0 : kpgx(ipw,88) = kpgx(ipw,60)*kpgx(ipw, 2)
271 : 0 : kpgx(ipw,89) = kpgx(ipw,61)*kpgx(ipw, 2)
272 : 0 : kpgx(ipw,90) = kpgx(ipw,62)*kpgx(ipw, 2)
273 : 0 : kpgx(ipw,91) = kpgx(ipw,63)*kpgx(ipw, 2)
274 : 0 : kpgx(ipw,92) = kpgx(ipw,64)*kpgx(ipw, 2)
275 : 0 : kpgx(ipw,93) = kpgx(ipw,65)*kpgx(ipw, 2)
276 : 0 : kpgx(ipw,94) = kpgx(ipw,66)*kpgx(ipw, 2)
277 : 0 : kpgx(ipw,95) = kpgx(ipw,67)*kpgx(ipw, 2)
278 : 0 : kpgx(ipw,96) = kpgx(ipw,68)*kpgx(ipw, 2)
279 : 0 : kpgx(ipw,97) = kpgx(ipw,69)*kpgx(ipw, 2)
280 : 0 : kpgx(ipw,98) = kpgx(ipw,70)*kpgx(ipw, 2)
281 : 0 : kpgx(ipw,99) = kpgx(ipw,71)*kpgx(ipw, 2)
282 : 0 : kpgx(ipw,100) = kpgx(ipw,72)*kpgx(ipw, 2)
283 : 0 : kpgx(ipw,101) = kpgx(ipw,73)*kpgx(ipw, 2)
284 : 0 : kpgx(ipw,102) = kpgx(ipw,74)*kpgx(ipw, 2)
285 : 0 : kpgx(ipw,103) = kpgx(ipw,75)*kpgx(ipw, 2)
286 : 0 : kpgx(ipw,104) = kpgx(ipw,76)*kpgx(ipw, 2)
287 : 0 : kpgx(ipw,105) = kpgx(ipw,77)*kpgx(ipw, 2)
288 : 0 : kpgx(ipw,106) = kpgx(ipw,78)*kpgx(ipw, 2)
289 : 0 : kpgx(ipw,107) = kpgx(ipw,79)*kpgx(ipw, 2)
290 : 0 : kpgx(ipw,108) = kpgx(ipw,80)*kpgx(ipw, 2)
291 : 0 : kpgx(ipw,109) = kpgx(ipw,81)*kpgx(ipw, 2)
292 : 0 : kpgx(ipw,110) = kpgx(ipw,82)*kpgx(ipw, 2)
293 : 0 : kpgx(ipw,111) = kpgx(ipw,83)*kpgx(ipw, 2)
294 : 0 : kpgx(ipw,112) = kpgx(ipw,84)*kpgx(ipw, 2)
295 : 0 : kpgx(ipw,113) = kpgx(ipw,78)*kpgx(ipw, 3)
296 : 0 : kpgx(ipw,114) = kpgx(ipw,79)*kpgx(ipw, 3)
297 : 0 : kpgx(ipw,115) = kpgx(ipw,80)*kpgx(ipw, 3)
298 : 0 : kpgx(ipw,116) = kpgx(ipw,81)*kpgx(ipw, 3)
299 : 0 : kpgx(ipw,117) = kpgx(ipw,82)*kpgx(ipw, 3)
300 : 0 : kpgx(ipw,118) = kpgx(ipw,83)*kpgx(ipw, 3)
301 : 0 : kpgx(ipw,119) = kpgx(ipw,84)*kpgx(ipw, 3)
302 : 0 : kpgx(ipw,120) = kpgx(ipw,84)*kpgx(ipw, 4)
303 : : end do
304 : : end if
305 : :
306 : : !*****************************************************************************
307 : : !
308 : : !Packing of composite projectors in ffkg
309 : :
310 : : iffkg=0
311 : :
312 : : !Treat composite projectors for the energy
313 : : iln0=0
314 [ + + ]: 1236587097 : do ilmn=1,lmnmax
315 : 1051459355 : iln=indlmn(5,ilmn,itypat)
316 [ + + ]: 1236587097 : if (iln>iln0) then
317 : 792397829 : iln0=iln
318 : 792397829 : ilang=1+indlmn(1,ilmn,itypat)
319 : 792397829 : iproj=indlmn(3,ilmn,itypat)
320 [ + - ]: 792397829 : if(iproj>0)then
321 : 792397829 : ilang2=(ilang*(ilang+1))/2
322 : :
323 [ + + ]: 792397829 : if(ilang==1)then
324 : : ! Treat s-component separately
325 : 294665150 : ig=ipw1
326 : 294665150 : iffkg=iffkg+1
327 [ + + ]: 44757284605 : do ipw=1,nincpw
328 : 44462619455 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)
329 : 44757284605 : ig=ig+1
330 : : end do
331 : 294665150 : parity(iffkg)=2
332 : : else
333 : : ! Treat other components (could be made faster by treating explicitely
334 : : ! each angular momentum)
335 [ + + ]: 2545085286 : do ii=1,ilang2
336 : : ! Get the starting address for the relevant tensor
337 : 2047352607 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
338 : 2047352607 : ig=ipw1
339 : 2047352607 : iffkg=iffkg+1
340 [ + + ]: >30317*10^7 : do ipw=1,nincpw
341 : >30112*10^7 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
342 : >30317*10^7 : ig=ig+1
343 : : end do
344 [ + + ]: 2047352607 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
345 [ + + ]: 2545085286 : if(ilang==3)parity(iffkg)=2
346 : : end do
347 : : end if
348 : :
349 : : ! End condition if(iproj>0)
350 : : end if
351 : :
352 : : ! End loop on ilang (ilmn)
353 : : end if
354 : : end do
355 : :
356 : : !This is the number of composite projectors for the energy
357 : 185127742 : nffkge=iffkg
358 : :
359 : : !Second, treat forces : actually, this part could be rationalized,
360 : : !since the outcome is a multiplication by three of the number
361 : : !of composite projectors for the energy, while less should be needed
362 [ + + + + ]: 185127742 : if((choice==2.or.choice==23) .and. ndgxdt/=1)then
363 [ + + ]: 12295439 : do ii=1,nffkge
364 [ + + ]: 1772449069 : do ipw=1,nincpw
365 : 1760843441 : ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,2)
366 : 1760843441 : ffkg(iffkg+2,ipw)=ffkg(ii,ipw)*kpgx(ipw,3)
367 : 1772449069 : ffkg(iffkg+3,ipw)=ffkg(ii,ipw)*kpgx(ipw,4)
368 : : end do
369 : 11605628 : parity(iffkg+1)=3-parity(ii)
370 : 11605628 : parity(iffkg+2)=parity(iffkg+1)
371 : 11605628 : parity(iffkg+3)=parity(iffkg+1)
372 : 12295439 : iffkg=iffkg+3
373 : : end do
374 : : end if
375 : : !Note that the additional number of projectors for forces is 3*nffkge
376 : :
377 : : !Third, treat first-derivative of the non-local operator
378 : : !with respect to an atomic displacement in one direction :
379 [ + + + + ]: 185127742 : if(choice==2 .and. ndgxdt==1)then
380 [ + + ]: 315096894 : do ii=1,nffkge
381 [ + + ]: 45437053629 : do ipw=1,nincpw
382 : 45437053629 : ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,idir+1)
383 : : end do
384 : 288735662 : parity(iffkg+1)=3-parity(ii)
385 : 315096894 : iffkg=iffkg+1
386 : : end do
387 : : end if
388 : : !Note that the additional number of projectors for this case is nffkge
389 : :
390 : :
391 : : !Fourth, treat dynamical matrices : like forces, this part could be rationalized.
392 [ + + ]: 185127742 : if(choice==4)then
393 [ + + ]: 3991862 : do ii=1,nffkge
394 [ + + ]: 526234992 : do ipw=1,nincpw
395 : 522560983 : kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4)
396 : 522560983 : ffkg_now=ffkg(ii,ipw)
397 : 522560983 : ffkg(iffkg+1,ipw)=ffkg_now*kpg_x
398 : 522560983 : ffkg(iffkg+2,ipw)=ffkg_now*kpg_y
399 : 522560983 : ffkg(iffkg+3,ipw)=ffkg_now*kpg_z
400 : 522560983 : ffkg(iffkg+4,ipw)=ffkg_now*kpg_x*kpg_x
401 : 522560983 : ffkg(iffkg+5,ipw)=ffkg_now*kpg_y*kpg_y
402 : 522560983 : ffkg(iffkg+6,ipw)=ffkg_now*kpg_z*kpg_z
403 : 522560983 : ffkg(iffkg+7,ipw)=ffkg_now*kpg_z*kpg_y
404 : 522560983 : ffkg(iffkg+8,ipw)=ffkg_now*kpg_z*kpg_x
405 : 526234992 : ffkg(iffkg+9,ipw)=ffkg_now*kpg_y*kpg_x
406 : : end do
407 [ + + ]: 14696036 : parity(iffkg+1:iffkg+3)=3-parity(ii)
408 [ + + ]: 25718063 : parity(iffkg+4:iffkg+9)=parity(ii)
409 : 3991862 : iffkg=iffkg+9
410 : : end do
411 : : end if
412 : : !Note that the additional number of projectors for dynamical matrices is 9*nffkge
413 : :
414 : : !Treat composite projectors for the stress or 1st derivative contribution
415 : : !to frozen-wavefunction part of elastic tensor
416 : : !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..)
417 [ + + ]: 185127742 : if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then
418 : :
419 : : iln0=0
420 [ + + ]: 34491119 : do ilmn=1,lmnmax
421 [ + + ]: 34491119 : if (ispinor==indlmn(6,ilmn,itypat)) then
422 : 24720203 : iln=indlmn(5,ilmn,itypat)
423 [ + - ]: 24720203 : if (iln>iln0) then
424 : 24720203 : iln0=iln
425 : 24720203 : ilang=1+indlmn(1,ilmn,itypat)
426 : 24720203 : iproj=indlmn(3,ilmn,itypat)
427 [ + - ]: 24720203 : if(iproj>0)then
428 : : ! number of unique tensor components
429 [ + + ]: 24720203 : if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2
430 [ + + ]: 24720203 : if(choice==5)ilangx=(ilang*(ilang+1))/2
431 : :
432 [ + + ]: 186359960 : do ii=1,ilangx
433 : : ! Get the starting address for the relevant tensor
434 : 123147811 : if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
435 [ + + ]: 161639757 : if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6
436 : 161639757 : ig=ipw1
437 [ + + ]: 161639757 : iffkg=iffkg+1
438 : : if(choice==3 .or. choice==6 .or. choice==23)then
439 [ + + ]: 11714975452 : do ipw=1,nincpw
440 : 11591827641 : ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
441 : 11714975452 : ig=ig+1
442 : : end do
443 : : else
444 [ + + ]: 3802795220 : do ipw=1,nincpw
445 : : ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*&
446 : : & (kpgx(ipw,2)*gmet(1,idir)+ &
447 : : & kpgx(ipw,3)*gmet(2,idir)+ &
448 : 3764303274 : & kpgx(ipw,4)*gmet(3,idir) )
449 : 3802795220 : ig=ig+1
450 : : end do
451 : : end if
452 [ + + ]: 161639757 : if(ilang==1 .or. ilang==3)parity(iffkg)=2
453 [ + + ]: 161639757 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
454 [ + + ]: 186359960 : if(choice==5)parity(iffkg)=3-parity(iffkg)
455 : : end do
456 : :
457 : : ! End condition if(iproj>0)
458 : : end if
459 : : ! End condition on iln
460 : : end if
461 : : ! End condition if(ispinor=indlmn(6,...))
462 : : end if
463 : : ! End loop on ilmn
464 : : end do
465 : :
466 : : ! End condition of stress
467 : : end if
468 : :
469 : : !Treat composite projectors for the 2nd derivative wrt 2 strains
470 : : !and wrt one strain and one atomic displacement (internal strain)
471 : : !contributions to frozen-wavefunction part of (generalized) elastic tensor.
472 : : !There are 3 sets on terms (in historical order):
473 : : !first, terms with ffnl(:,3,...) and rank+4 tensors.
474 : : !second, terms with ffnl(:,1,...) and rank+1 tensors.
475 : : !third, terms with ffnl(:,2,...) and rank+3 tensors.
476 : :
477 [ + + ]: 185127742 : if(choice==6)then
478 : :
479 : : iln0=0
480 [ + + ]: 198592 : do ilmn=1,lmnmax
481 [ + + ]: 198592 : if (ispinor==indlmn(6,ilmn,itypat)) then
482 : 165000 : iln=indlmn(5,ilmn,itypat)
483 [ + - ]: 165000 : if (iln>iln0) then
484 : 165000 : iln0=iln
485 : 165000 : ilang=1+indlmn(1,ilmn,itypat)
486 : 165000 : iproj=indlmn(3,ilmn,itypat)
487 : :
488 [ + - ]: 165000 : if(iproj>0)then
489 : : ! First set of terms
490 : : ! number of unique tensor components
491 : 165000 : ilangx=((ilang+4)*(ilang+5))/2
492 : :
493 [ + + ]: 3660496 : do ii=1,ilangx
494 : : ! Get the starting address for the relevant tensor
495 : 3495496 : jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6
496 : 3495496 : ig=ipw1
497 : 3495496 : iffkg=iffkg+1
498 [ + + ]: 394676100 : do ipw=1,nincpw
499 : 391180604 : ffkg(iffkg,ipw)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj)
500 : 394676100 : ig=ig+1
501 : : end do
502 [ + + ]: 3495496 : if(ilang==1 .or. ilang==3)parity(iffkg)=2
503 [ + + ]: 3660496 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
504 : : end do
505 : :
506 : : ! Second set of terms
507 : : ! number of unique tensor components
508 : 165000 : ilangx=((ilang+1)*(ilang+2))/2
509 : :
510 [ + + ]: 1196776 : do ii=1,ilangx
511 : : ! Get the starting address for the relevant tensor
512 : 1031776 : jj=ii+((ilang)*(ilang+1)*(ilang+2))/6
513 : 1031776 : ig=ipw1
514 : 1031776 : iffkg=iffkg+1
515 [ + + ]: 116507400 : do ipw=1,nincpw
516 : 115475624 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
517 : 116507400 : ig=ig+1
518 : : end do
519 [ + + ]: 1031776 : if(ilang==1 .or. ilang==3)parity(iffkg)=1
520 [ + + ]: 1196776 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
521 : : end do
522 : :
523 : : ! Third set of terms
524 : : ! number of unique tensor components
525 : 165000 : ilangx=((ilang+3)*(ilang+4))/2
526 : :
527 [ + + ]: 2674256 : do ii=1,ilangx
528 : : ! Get the starting address for the relevant tensor
529 : 2509256 : jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6
530 : 2509256 : ig=ipw1
531 : 2509256 : iffkg=iffkg+1
532 [ + + ]: 283325180 : do ipw=1,nincpw
533 : 280815924 : ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
534 : 283325180 : ig=ig+1
535 : : end do
536 [ + + ]: 2509256 : if(ilang==1 .or. ilang==3)parity(iffkg)=1
537 [ + + ]: 2674256 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
538 : : end do
539 : :
540 : : ! End condition if(iproj>0)
541 : : end if
542 : : ! End condition on iln
543 : : end if
544 : : ! End condition if(ispinor=indlmn(6,...))
545 : : end if
546 : : ! End loop on ilmn
547 : : end do
548 : :
549 : : ! End condition of 2nd strain derivatives
550 : : end if
551 : :
552 : : !For ddk perturbation, treat the part that depend on ffnl(:,1,..)
553 : : !no contribution from s state
554 [ + + + + ]: 185127742 : if(nlang>=2 .and. choice==5)then
555 : : iln0=0
556 [ + + ]: 17135036 : do ilmn=1,lmnmax
557 [ + + ]: 17135036 : if (ispinor==indlmn(6,ilmn,itypat)) then
558 : 11948832 : iln=indlmn(5,ilmn,itypat)
559 [ + - ]: 11948832 : if (iln>iln0) then
560 : 11948832 : iln0=iln
561 : 11948832 : ilang=1+indlmn(1,ilmn,itypat)
562 [ + + ]: 11948832 : if (ilang>=2) then
563 : 7858038 : iproj=indlmn(3,ilmn,itypat)
564 [ + - ]: 7858038 : if(iproj>0)then
565 : 7858038 : ilang2=(ilang*(ilang-1))/2
566 : :
567 [ + + ]: 22316068 : do ii=1,ilang2
568 : : ! Get the starting address for the relevant tensor
569 : 14458030 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
570 : 14458030 : ig=ipw1
571 : 14458030 : iffkg=iffkg+1
572 [ + + ]: 1394727738 : do ipw=1,nincpw
573 : 1380269708 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
574 : 1394727738 : ig=ig+1
575 : : end do
576 [ + + ]: 14458030 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
577 [ + + ]: 22316068 : if(ilang==3)parity(iffkg)=1
578 : : end do
579 : :
580 : : ! End condition if(iproj>0)
581 : : end if
582 : : ! End condition if(ilang>=2)
583 : : end if
584 : : ! End condition if(iln>iln0)
585 : : end if
586 : : ! End condition if(ispinor=indlmn(6,...))
587 : : end if
588 : : ! End loop on ilmn
589 : : end do
590 : : ! End condition of p,d or f state
591 : : end if
592 : :
593 : : !DEBUG
594 : : !write(std_out,*)' dfpt_mkffkg : exit '
595 : : !ENDDEBUG
596 : :
597 : 185127742 : end subroutine dfpt_mkffkg
598 : : !!***
599 : :
600 : : !!****f* ABINIT/mkffkg
601 : : !! NAME
602 : : !! mkffkg
603 : : !!
604 : : !! FUNCTION
605 : : !! Prepare the application of the projectors to the shifted wavefunctions,
606 : : !! by precomputing the k+G factors and their product with the form factors
607 : : !! Do this on a block of plane waves.
608 : : !!
609 : : !! INPUTS
610 : : !! choice=governs the combination of k+G vectors to be computed
611 : : !! ffnl(npw,nffnl,lmnmax,ntypat)=nonlocal form factors on basis sphere.
612 : : !! gmet(3,3)=metric tensor for G vecs (in bohr**-2)
613 : : !! nffnl=3rd dimension of ffnl(2, conventional, or 3 for 2nd derivatives)
614 : : !! idir=direction of the perturbation (needed if choice==2 and ndgxdt==1,
615 : : !! or if choice==5)
616 : : !! indlmn(6,i,ntypat)=array giving l,m,n,lm,ln,spin for i=ln
617 : : !! ipw1 = index of the first plane wave treated in this block
618 : : !! ispinor=1 or 2, gives the spinorial component of ffnl to be used
619 : : !! itypat = type of atom, needed for ffnl
620 : : !! kg_k(3,npw)=integer coords of planewaves in basis sphere
621 : : !! kpg_k(npw,npkg)= (k+G) components and related data
622 : : !! kpt(3)=real components of k point in terms of recip. translations
623 : : !! lmnmax=max. number of (l,n) components over all type of psps
624 : : !! mblkpw=first dimension of kpgx
625 : : !! ndgxdt=number of components of first order derivative
626 : : !! nffkg=number of products of ffnls with combinations of k+G
627 : : !! nincpw=number of plane waves in the block
628 : : !! nkpg=second size of array kpg_k
629 : : !! nlang = number of angular momenta to be treated = 1 + highest ang. mom.
630 : : !! nloalg(3)=governs the choice of the algorithm for non-local operator.
631 : : !! npw = total number of plane waves in reciprocal space
632 : : !! ntens=second dimension of kpgx, number of distinct tensorial products
633 : : !! ntypat = number of type of atoms, dimension needed for ffnl
634 : : !!
635 : : !! OUTPUT
636 : : !! kpgx(mblkpw,ntens)=different tensorial products of k+G
637 : : !! ffkg(mblkpw,nffkg)=different products of ffnls with k+G
638 : : !! parity(nffkg)=parity of the tensorial product of k+G (2 if even, 1 of odd)
639 : : !!
640 : : !! NOTES
641 : : !! This routine must be thread-safe as it is called inside loops that are OpenMP parallelized.
642 : : !! Please, do not add variables with the save attribute or SIDE EFFECTS.
643 : : !!
644 : : !! SOURCE
645 : :
646 : 91404 : subroutine mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
647 : 91404 : & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
648 : 91404 : & npw,ntens,ntypat,parity)
649 : :
650 : : !Arguments ------------------------------------
651 : : !scalars
652 : : integer,intent(in) :: choice,idir,ipw1,ispinor,itypat,lmnmax,mblkpw,ndgxdt
653 : : integer,intent(in) :: nffkg,nffnl,nincpw,nkpg,nlang,npw,ntens,ntypat
654 : : !arrays
655 : : integer,intent(in) :: indlmn(6,lmnmax,ntypat),kg_k(3,npw)
656 : : integer,intent(out) :: parity(nffkg)
657 : : real(dp),intent(in) :: ffnl(npw,nffnl,lmnmax,ntypat),gmet(3,3),kpg_k(npw,nkpg)
658 : : real(dp),intent(in) :: kpt(3)
659 : : real(dp),intent(out) :: ffkg(mblkpw,nffkg),kpgx(mblkpw,ntens)
660 : :
661 : : !Local variables-------------------------------
662 : : !scalars
663 : : integer :: iffkg,ig,ii,ilang,ilang2,ilangx,ilmn,iln,iln0,iproj,ipw,jj
664 : : integer :: nffkge
665 : : real(dp) :: ffkg_now,kpg_x,kpg_y,kpg_z
666 : : !arrays
667 : :
668 : : ! *************************************************************************
669 : :
670 : 91404 : jj=0;ilangx=0
671 : :
672 : : !This will be useless after all the modifications have been done
673 [ + + ]: 15487040 : do ipw=1,nincpw
674 : 15487040 : kpgx(ipw,1)=1.0d0
675 : : end do
676 : :
677 : : !Initialize kpgx array related to tensors defined below
678 : : if ( nlang>=2 .or. choice==2 .or. choice==3 .or. choice==4 .or. choice==5&
679 [ - + - - : 91404 : & .or. choice==6 .or. choice==23) then
- - - - -
- - - -
- ]
680 [ - + ]: 91404 : if (nkpg>=3) then
681 [ # # ]: 0 : kpgx(1:nincpw,2)=kpg_k(ipw1+1:ipw1+nincpw,1)
682 [ # # ]: 0 : kpgx(1:nincpw,3)=kpg_k(ipw1+1:ipw1+nincpw,2)
683 [ # # ]: 0 : kpgx(1:nincpw,4)=kpg_k(ipw1+1:ipw1+nincpw,3)
684 : : else
685 : 91404 : ig=ipw1
686 [ + + ]: 15487040 : do ipw=1,nincpw
687 : 15395636 : kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig))
688 : 15395636 : kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig))
689 : 15395636 : kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig))
690 : 15487040 : ig=ig+1
691 : : end do
692 : : end if
693 : : end if
694 [ + + + + ]: 91404 : if (nlang>=3 .or. choice==3 .or. choice==6 .or. choice==23) then
695 : : ! Define (k+G) part of rank 2 symmetric tensor (6 components), l=2
696 : : ! Compressed storage is 11 22 33 32 31 21
697 [ + - ]: 49446 : if (nkpg>=9) then
698 [ # # ]: 0 : kpgx(1:nincpw,5) =kpg_k(ipw1+1:ipw1+nincpw,4)
699 [ # # ]: 0 : kpgx(1:nincpw,6) =kpg_k(ipw1+1:ipw1+nincpw,5)
700 [ # # ]: 0 : kpgx(1:nincpw,7) =kpg_k(ipw1+1:ipw1+nincpw,6)
701 [ # # ]: 0 : kpgx(1:nincpw,8) =kpg_k(ipw1+1:ipw1+nincpw,7)
702 [ # # ]: 0 : kpgx(1:nincpw,9) =kpg_k(ipw1+1:ipw1+nincpw,8)
703 [ # # ]: 0 : kpgx(1:nincpw,10)=kpg_k(ipw1+1:ipw1+nincpw,9)
704 : : else
705 [ + + ]: 8306224 : do ipw=1,nincpw
706 : 8256778 : kpgx(ipw, 5) = kpgx(ipw, 2)*kpgx(ipw, 2)
707 : 8256778 : kpgx(ipw, 6) = kpgx(ipw, 3)*kpgx(ipw, 3)
708 : 8256778 : kpgx(ipw, 7) = kpgx(ipw, 4)*kpgx(ipw, 4)
709 : 8256778 : kpgx(ipw, 8) = kpgx(ipw, 4)*kpgx(ipw, 3)
710 : 8256778 : kpgx(ipw, 9) = kpgx(ipw, 4)*kpgx(ipw, 2)
711 : 8306224 : kpgx(ipw,10) = kpgx(ipw, 3)*kpgx(ipw, 2)
712 : : end do
713 : : end if
714 : : end if
715 [ + + + + : 91404 : if (nlang>=4 .or. ((choice==3.or.choice==23) .and. nlang>=2) .or. choice==6) then
+ - + + ]
716 : : ! Define (k+G) part of rank 3 symmetric tensor (10 components), l=3
717 : : ! Compressed storage is 111 221 331 321 311 211 222 332 322 333
718 [ + + ]: 4004438 : do ipw=1,nincpw
719 : 3977660 : kpgx(ipw,11) = kpgx(ipw, 5)*kpgx(ipw, 2)
720 : 3977660 : kpgx(ipw,12) = kpgx(ipw, 6)*kpgx(ipw, 2)
721 : 3977660 : kpgx(ipw,13) = kpgx(ipw, 7)*kpgx(ipw, 2)
722 : 3977660 : kpgx(ipw,14) = kpgx(ipw, 8)*kpgx(ipw, 2)
723 : 3977660 : kpgx(ipw,15) = kpgx(ipw, 9)*kpgx(ipw, 2)
724 : 3977660 : kpgx(ipw,16) = kpgx(ipw,10)*kpgx(ipw, 2)
725 : 3977660 : kpgx(ipw,17) = kpgx(ipw, 6)*kpgx(ipw, 3)
726 : 3977660 : kpgx(ipw,18) = kpgx(ipw, 7)*kpgx(ipw, 3)
727 : 3977660 : kpgx(ipw,19) = kpgx(ipw, 8)*kpgx(ipw, 3)
728 : 4069064 : kpgx(ipw,20) = kpgx(ipw, 7)*kpgx(ipw, 4)
729 : : end do
730 : : end if
731 [ + + + + : 91404 : if (((choice==3.or.choice==23) .and. nlang>=3) .or. choice==6) then
+ + ]
732 : : ! Add additional tensors for strain gradients
733 : : ! Define (k+G) part of rank 4 symmetric tensor (15 components), l=2
734 : : ! Compressed storage is 1111 2211 3311 3211 3111 2111 2221 3321 3221
735 : : ! 3331 2222 3322 3222 3332 3333
736 [ + + ]: 569868 : do ipw=1,nincpw
737 : 566076 : kpgx(ipw,21) = kpgx(ipw, 5)*kpgx(ipw, 5)
738 : 566076 : kpgx(ipw,22) = kpgx(ipw, 6)*kpgx(ipw, 5)
739 : 566076 : kpgx(ipw,23) = kpgx(ipw, 7)*kpgx(ipw, 5)
740 : 566076 : kpgx(ipw,24) = kpgx(ipw, 8)*kpgx(ipw, 5)
741 : 566076 : kpgx(ipw,25) = kpgx(ipw, 9)*kpgx(ipw, 5)
742 : 566076 : kpgx(ipw,26) = kpgx(ipw,10)*kpgx(ipw, 5)
743 : 566076 : kpgx(ipw,27) = kpgx(ipw, 6)*kpgx(ipw,10)
744 : 566076 : kpgx(ipw,28) = kpgx(ipw, 7)*kpgx(ipw,10)
745 : 566076 : kpgx(ipw,29) = kpgx(ipw, 8)*kpgx(ipw,10)
746 : 566076 : kpgx(ipw,30) = kpgx(ipw, 7)*kpgx(ipw, 9)
747 : 566076 : kpgx(ipw,31) = kpgx(ipw, 6)*kpgx(ipw, 6)
748 : 566076 : kpgx(ipw,32) = kpgx(ipw, 7)*kpgx(ipw, 6)
749 : 566076 : kpgx(ipw,33) = kpgx(ipw, 8)*kpgx(ipw, 6)
750 : 566076 : kpgx(ipw,34) = kpgx(ipw, 7)*kpgx(ipw, 8)
751 : 657480 : kpgx(ipw,35) = kpgx(ipw, 7)*kpgx(ipw, 7)
752 : : end do
753 : : end if
754 [ + + + + : 91404 : if (((choice==3.or.choice==23) .and. nlang>=4) .or. (choice==6 .and. nlang>=2)) then
+ + - + ]
755 : : ! Define (k+G) part of rank 5 symmetric tensor (21 components), l=3
756 : : ! Compressed storage is 11111 22111 33111 32111 31111 21111
757 : : ! 22211 33211 32211 33311 22221 33221 32221 33321 33331
758 : : ! 22222 33222 32222 33322 33332 33333
759 [ + + ]: 551628 : do ipw=1,nincpw
760 : 547932 : kpgx(ipw,36) = kpgx(ipw,21)*kpgx(ipw, 2)
761 : 547932 : kpgx(ipw,37) = kpgx(ipw,22)*kpgx(ipw, 2)
762 : 547932 : kpgx(ipw,38) = kpgx(ipw,23)*kpgx(ipw, 2)
763 : 547932 : kpgx(ipw,39) = kpgx(ipw,24)*kpgx(ipw, 2)
764 : 547932 : kpgx(ipw,40) = kpgx(ipw,25)*kpgx(ipw, 2)
765 : 547932 : kpgx(ipw,41) = kpgx(ipw,26)*kpgx(ipw, 2)
766 : 547932 : kpgx(ipw,42) = kpgx(ipw,27)*kpgx(ipw, 2)
767 : 547932 : kpgx(ipw,43) = kpgx(ipw,28)*kpgx(ipw, 2)
768 : 547932 : kpgx(ipw,44) = kpgx(ipw,29)*kpgx(ipw, 2)
769 : 547932 : kpgx(ipw,45) = kpgx(ipw,30)*kpgx(ipw, 2)
770 : 547932 : kpgx(ipw,46) = kpgx(ipw,31)*kpgx(ipw, 2)
771 : 547932 : kpgx(ipw,47) = kpgx(ipw,32)*kpgx(ipw, 2)
772 : 547932 : kpgx(ipw,48) = kpgx(ipw,33)*kpgx(ipw, 2)
773 : 547932 : kpgx(ipw,49) = kpgx(ipw,34)*kpgx(ipw, 2)
774 : 547932 : kpgx(ipw,50) = kpgx(ipw,35)*kpgx(ipw, 2)
775 : 547932 : kpgx(ipw,51) = kpgx(ipw,31)*kpgx(ipw, 3)
776 : 547932 : kpgx(ipw,52) = kpgx(ipw,32)*kpgx(ipw, 3)
777 : 547932 : kpgx(ipw,53) = kpgx(ipw,33)*kpgx(ipw, 3)
778 : 547932 : kpgx(ipw,54) = kpgx(ipw,34)*kpgx(ipw, 3)
779 : 547932 : kpgx(ipw,55) = kpgx(ipw,35)*kpgx(ipw, 3)
780 : 639336 : kpgx(ipw,56) = kpgx(ipw,35)*kpgx(ipw, 4)
781 : : end do
782 : : end if
783 [ + + + + ]: 91404 : if (choice==6 .and. nlang>=3) then
784 : : ! Define (k+G) part of rank 6 symmetric tensor (28 components)
785 : : ! Compressed storage is
786 : : ! 111111 221111 331111 321111 311111 211111 222111 332111 322111
787 : : ! 333111 222211 332211 322211 333211 333311 222221 332221 322221
788 : : ! 333221 333321 333331 222222 332222 322222 333222 333322 333332
789 : : ! 333333
790 [ + + ]: 7164 : do ipw=1,nincpw
791 : 7116 : kpgx(ipw,57) = kpgx(ipw,36)*kpgx(ipw, 2)
792 : 7116 : kpgx(ipw,58) = kpgx(ipw,37)*kpgx(ipw, 2)
793 : 7116 : kpgx(ipw,59) = kpgx(ipw,38)*kpgx(ipw, 2)
794 : 7116 : kpgx(ipw,60) = kpgx(ipw,39)*kpgx(ipw, 2)
795 : 7116 : kpgx(ipw,61) = kpgx(ipw,40)*kpgx(ipw, 2)
796 : 7116 : kpgx(ipw,62) = kpgx(ipw,41)*kpgx(ipw, 2)
797 : 7116 : kpgx(ipw,63) = kpgx(ipw,42)*kpgx(ipw, 2)
798 : 7116 : kpgx(ipw,64) = kpgx(ipw,43)*kpgx(ipw, 2)
799 : 7116 : kpgx(ipw,65) = kpgx(ipw,44)*kpgx(ipw, 2)
800 : 7116 : kpgx(ipw,66) = kpgx(ipw,45)*kpgx(ipw, 2)
801 : 7116 : kpgx(ipw,67) = kpgx(ipw,46)*kpgx(ipw, 2)
802 : 7116 : kpgx(ipw,68) = kpgx(ipw,47)*kpgx(ipw, 2)
803 : 7116 : kpgx(ipw,69) = kpgx(ipw,48)*kpgx(ipw, 2)
804 : 7116 : kpgx(ipw,70) = kpgx(ipw,49)*kpgx(ipw, 2)
805 : 7116 : kpgx(ipw,71) = kpgx(ipw,50)*kpgx(ipw, 2)
806 : 7116 : kpgx(ipw,72) = kpgx(ipw,51)*kpgx(ipw, 2)
807 : 7116 : kpgx(ipw,73) = kpgx(ipw,52)*kpgx(ipw, 2)
808 : 7116 : kpgx(ipw,74) = kpgx(ipw,53)*kpgx(ipw, 2)
809 : 7116 : kpgx(ipw,75) = kpgx(ipw,54)*kpgx(ipw, 2)
810 : 7116 : kpgx(ipw,76) = kpgx(ipw,55)*kpgx(ipw, 2)
811 : 7116 : kpgx(ipw,77) = kpgx(ipw,56)*kpgx(ipw, 2)
812 : 7116 : kpgx(ipw,78) = kpgx(ipw,51)*kpgx(ipw, 3)
813 : 7116 : kpgx(ipw,79) = kpgx(ipw,52)*kpgx(ipw, 3)
814 : 7116 : kpgx(ipw,80) = kpgx(ipw,53)*kpgx(ipw, 3)
815 : 7116 : kpgx(ipw,81) = kpgx(ipw,54)*kpgx(ipw, 3)
816 : 7116 : kpgx(ipw,82) = kpgx(ipw,55)*kpgx(ipw, 3)
817 : 7116 : kpgx(ipw,83) = kpgx(ipw,56)*kpgx(ipw, 3)
818 : 7164 : kpgx(ipw,84) = kpgx(ipw,56)*kpgx(ipw, 4)
819 : : end do
820 : : end if
821 [ + + + + ]: 91404 : if (choice==6 .and. nlang==4) then
822 : : ! Define (k+G) part of rank 7 symmetric tensor (36 components)
823 : : ! Compressed storage is
824 : : ! 1111111 2211111 3311111 3211111 3111111 2111111 2221111 3321111 3221111
825 : : ! 3331111 2222111 3322111 3222111 3332111 3333111 2222211 3322211 3222211
826 : : ! 3332211 3333211 3333311 2222221 3322221 3222221 3332221 3333221 3333321
827 : : ! 3333331 2222222 3322222 3222222 3332222 3333222 3333322 3333332 3333333
828 [ + + ]: 7164 : do ipw=1,nincpw
829 : 7116 : kpgx(ipw,85) = kpgx(ipw,57)*kpgx(ipw, 2)
830 : 7116 : kpgx(ipw,86) = kpgx(ipw,58)*kpgx(ipw, 2)
831 : 7116 : kpgx(ipw,87) = kpgx(ipw,59)*kpgx(ipw, 2)
832 : 7116 : kpgx(ipw,88) = kpgx(ipw,60)*kpgx(ipw, 2)
833 : 7116 : kpgx(ipw,89) = kpgx(ipw,61)*kpgx(ipw, 2)
834 : 7116 : kpgx(ipw,90) = kpgx(ipw,62)*kpgx(ipw, 2)
835 : 7116 : kpgx(ipw,91) = kpgx(ipw,63)*kpgx(ipw, 2)
836 : 7116 : kpgx(ipw,92) = kpgx(ipw,64)*kpgx(ipw, 2)
837 : 7116 : kpgx(ipw,93) = kpgx(ipw,65)*kpgx(ipw, 2)
838 : 7116 : kpgx(ipw,94) = kpgx(ipw,66)*kpgx(ipw, 2)
839 : 7116 : kpgx(ipw,95) = kpgx(ipw,67)*kpgx(ipw, 2)
840 : 7116 : kpgx(ipw,96) = kpgx(ipw,68)*kpgx(ipw, 2)
841 : 7116 : kpgx(ipw,97) = kpgx(ipw,69)*kpgx(ipw, 2)
842 : 7116 : kpgx(ipw,98) = kpgx(ipw,70)*kpgx(ipw, 2)
843 : 7116 : kpgx(ipw,99) = kpgx(ipw,71)*kpgx(ipw, 2)
844 : 7116 : kpgx(ipw,100) = kpgx(ipw,72)*kpgx(ipw, 2)
845 : 7116 : kpgx(ipw,101) = kpgx(ipw,73)*kpgx(ipw, 2)
846 : 7116 : kpgx(ipw,102) = kpgx(ipw,74)*kpgx(ipw, 2)
847 : 7116 : kpgx(ipw,103) = kpgx(ipw,75)*kpgx(ipw, 2)
848 : 7116 : kpgx(ipw,104) = kpgx(ipw,76)*kpgx(ipw, 2)
849 : 7116 : kpgx(ipw,105) = kpgx(ipw,77)*kpgx(ipw, 2)
850 : 7116 : kpgx(ipw,106) = kpgx(ipw,78)*kpgx(ipw, 2)
851 : 7116 : kpgx(ipw,107) = kpgx(ipw,79)*kpgx(ipw, 2)
852 : 7116 : kpgx(ipw,108) = kpgx(ipw,80)*kpgx(ipw, 2)
853 : 7116 : kpgx(ipw,109) = kpgx(ipw,81)*kpgx(ipw, 2)
854 : 7116 : kpgx(ipw,110) = kpgx(ipw,82)*kpgx(ipw, 2)
855 : 7116 : kpgx(ipw,111) = kpgx(ipw,83)*kpgx(ipw, 2)
856 : 7116 : kpgx(ipw,112) = kpgx(ipw,84)*kpgx(ipw, 2)
857 : 7116 : kpgx(ipw,113) = kpgx(ipw,78)*kpgx(ipw, 3)
858 : 7116 : kpgx(ipw,114) = kpgx(ipw,79)*kpgx(ipw, 3)
859 : 7116 : kpgx(ipw,115) = kpgx(ipw,80)*kpgx(ipw, 3)
860 : 7116 : kpgx(ipw,116) = kpgx(ipw,81)*kpgx(ipw, 3)
861 : 7116 : kpgx(ipw,117) = kpgx(ipw,82)*kpgx(ipw, 3)
862 : 7116 : kpgx(ipw,118) = kpgx(ipw,83)*kpgx(ipw, 3)
863 : 7116 : kpgx(ipw,119) = kpgx(ipw,84)*kpgx(ipw, 3)
864 : 7164 : kpgx(ipw,120) = kpgx(ipw,84)*kpgx(ipw, 4)
865 : : end do
866 : : end if
867 : :
868 : : !*****************************************************************************
869 : : !
870 : : !Packing of composite projectors in ffkg
871 : :
872 : : iffkg=0
873 : :
874 : : !Treat composite projectors for the energy
875 : : iln0=0
876 [ + + ]: 457020 : do ilmn=1,lmnmax
877 : 365616 : iln=indlmn(5,ilmn,itypat)
878 [ + + ]: 457020 : if (iln>iln0) then
879 : 296802 : iln0=iln
880 : 296802 : ilang=1+indlmn(1,ilmn,itypat)
881 : 296802 : iproj=indlmn(3,ilmn,itypat)
882 [ + - ]: 296802 : if(iproj>0)then
883 : 296802 : ilang2=(ilang*(ilang+1))/2
884 : :
885 [ + + ]: 296802 : if(ilang==1)then
886 : : ! Treat s-component separately
887 : 136932 : ig=ipw1
888 : 136932 : iffkg=iffkg+1
889 [ + + ]: 24127092 : do ipw=1,nincpw
890 : 23990160 : ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)
891 : 24127092 : ig=ig+1
892 : : end do
893 : 136932 : parity(iffkg)=2
894 : : else
895 : : ! Treat other components (could be made faster by treating explicitely
896 : : ! each angular momentum)
897 [ + + ]: 937152 : do ii=1,ilang2
898 : : ! Get the starting address for the relevant tensor
899 : 777282 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
900 : 777282 : ig=ipw1
901 : 777282 : iffkg=iffkg+1
902 [ + + ]: 129846776 : do ipw=1,nincpw
903 : 129069494 : ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
904 : 129846776 : ig=ig+1
905 : : end do
906 [ + + ]: 777282 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
907 [ + + ]: 937152 : if(ilang==3)parity(iffkg)=2
908 : : end do
909 : : end if
910 : :
911 : : ! End condition if(iproj>0)
912 : : end if
913 : : ! End loop on ilang (ilmn)
914 : : end if
915 : : end do
916 : :
917 : : !This is the number of composite projectors for the energy
918 : 91404 : nffkge=iffkg
919 : :
920 : : !Second, treat forces : actually, this part could be rationalized,
921 : : !since the outcome is a multiplication by three of the number
922 : : !of composite projectors for the energy, while less should be needed
923 [ + + + + ]: 91404 : if((choice==2.or.choice==23) .and. ndgxdt/=1)then
924 [ + + ]: 3120 : do ii=1,nffkge
925 [ + + ]: 497004 : do ipw=1,nincpw
926 : 494172 : ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,2)
927 : 494172 : ffkg(ipw,iffkg+2)=ffkg(ipw,ii)*kpgx(ipw,3)
928 : 497004 : ffkg(ipw,iffkg+3)=ffkg(ipw,ii)*kpgx(ipw,4)
929 : : end do
930 : 2832 : parity(iffkg+1)=3-parity(ii)
931 : 2832 : parity(iffkg+2)=parity(iffkg+1)
932 : 2832 : parity(iffkg+3)=parity(iffkg+1)
933 : 3120 : iffkg=iffkg+3
934 : : end do
935 : : end if
936 : : !Note that the additional number of projectors for forces is 3*nffkge
937 : :
938 : : !Third, treat first-derivative of the non-local operator
939 : : !with respect to an atomic displacement in one direction :
940 [ + + + - ]: 91404 : if(choice==2 .and. ndgxdt==1)then
941 [ + + ]: 13248 : do ii=1,nffkge
942 [ + + ]: 1805328 : do ipw=1,nincpw
943 : 1805328 : ffkg(ipw,iffkg+1)=ffkg(ipw,ii)*kpgx(ipw,idir+1)
944 : : end do
945 : 12096 : parity(iffkg+1)=3-parity(ii)
946 : 13248 : iffkg=iffkg+1
947 : : end do
948 : : end if
949 : : !Note that the additional number of projectors for this case is nffkge
950 : :
951 : :
952 : : !Fourth, treat dynamical matrices : like forces, this part could be rationalized.
953 [ - + ]: 91404 : if(choice==4)then
954 [ # # ]: 0 : do ii=1,nffkge
955 [ # # ]: 0 : do ipw=1,nincpw
956 : 0 : kpg_x=kpgx(ipw,2) ; kpg_y=kpgx(ipw,3) ; kpg_z=kpgx(ipw,4)
957 : 0 : ffkg_now=ffkg(ipw,ii)
958 : 0 : ffkg(ipw,iffkg+1)=ffkg_now*kpg_x
959 : 0 : ffkg(ipw,iffkg+2)=ffkg_now*kpg_y
960 : 0 : ffkg(ipw,iffkg+3)=ffkg_now*kpg_z
961 : 0 : ffkg(ipw,iffkg+4)=ffkg_now*kpg_x*kpg_x
962 : 0 : ffkg(ipw,iffkg+5)=ffkg_now*kpg_y*kpg_y
963 : 0 : ffkg(ipw,iffkg+6)=ffkg_now*kpg_z*kpg_z
964 : 0 : ffkg(ipw,iffkg+7)=ffkg_now*kpg_z*kpg_y
965 : 0 : ffkg(ipw,iffkg+8)=ffkg_now*kpg_z*kpg_x
966 : 0 : ffkg(ipw,iffkg+9)=ffkg_now*kpg_y*kpg_x
967 : : end do
968 [ # # ]: 0 : parity(iffkg+1:iffkg+3)=3-parity(ii)
969 [ # # ]: 0 : parity(iffkg+4:iffkg+9)=parity(ii)
970 : 0 : iffkg=iffkg+9
971 : : end do
972 : : end if
973 : : !Note that the additional number of projectors for dynamical matrices is 9*nffkge
974 : :
975 : : !Treat composite projectors for the stress or 1st derivative contribution
976 : : !to frozen-wavefunction part of elastic tensor
977 : : !as well as, for ddk perturbation, the part that depend on ffnl(:,2,..)
978 [ + + ]: 91404 : if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then
979 : :
980 : : iln0=0
981 [ + + ]: 37440 : do ilmn=1,lmnmax
982 [ + + ]: 37440 : if (ispinor==indlmn(6,ilmn,itypat)) then
983 : 19008 : iln=indlmn(5,ilmn,itypat)
984 [ + - ]: 19008 : if (iln>iln0) then
985 : 19008 : iln0=iln
986 : 19008 : ilang=1+indlmn(1,ilmn,itypat)
987 : 19008 : iproj=indlmn(3,ilmn,itypat)
988 [ + - ]: 19008 : if(iproj>0)then
989 : : ! number of unique tensor components
990 [ + - ]: 19008 : if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2
991 [ - + ]: 19008 : if(choice==5)ilangx=(ilang*(ilang+1))/2
992 : :
993 [ + + ]: 237216 : do ii=1,ilangx
994 : : ! Get the starting address for the relevant tensor
995 : 218208 : if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
996 [ - + ]: 218208 : if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6
997 : 218208 : ig=ipw1
998 [ - + ]: 218208 : iffkg=iffkg+1
999 : : if(choice==3 .or. choice==6 .or. choice==23)then
1000 [ + + ]: 32837472 : do ipw=1,nincpw
1001 : 32619264 : ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
1002 : 32837472 : ig=ig+1
1003 : : end do
1004 : : else
1005 [ # # ]: 0 : do ipw=1,nincpw
1006 : : ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)*&
1007 : : & (kpgx(ipw,2)*gmet(1,idir)+ &
1008 : : & kpgx(ipw,3)*gmet(2,idir)+ &
1009 : 0 : & kpgx(ipw,4)*gmet(3,idir) )
1010 : 0 : ig=ig+1
1011 : : end do
1012 : : end if
1013 [ + + ]: 218208 : if(ilang==1 .or. ilang==3)parity(iffkg)=2
1014 [ + + ]: 218208 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
1015 [ - + ]: 237216 : if(choice==5)parity(iffkg)=3-parity(iffkg)
1016 : : end do
1017 : :
1018 : : ! End condition if(iproj>0)
1019 : : end if
1020 : :
1021 : : ! End loop on ilang (ilmn)
1022 : : end if
1023 : : end if
1024 : : end do
1025 : :
1026 : : ! End condition of stress
1027 : : end if
1028 : :
1029 : : !Treat composite projectors for the 2nd derivative wrt 2 strains
1030 : : !and wrt one strain and one atomic displacement (internal strain)
1031 : : !contributions to frozen-wavefunction part of (generalized) elastic tensor.
1032 : : !There are 3 sets on terms (in historical order):
1033 : : !first, terms with ffnl(:,3,...) and rank+4 tensors.
1034 : : !second, terms with ffnl(:,1,...) and rank+1 tensors.
1035 : : !third, terms with ffnl(:,2,...) and rank+3 tensors.
1036 : :
1037 [ + + ]: 91404 : if(choice==6)then
1038 : :
1039 : : iln0=0
1040 [ + + ]: 480 : do ilmn=1,lmnmax
1041 [ + + ]: 480 : if (ispinor==indlmn(6,ilmn,itypat)) then
1042 : 240 : iln=indlmn(5,ilmn,itypat)
1043 [ + - ]: 240 : if (iln>iln0) then
1044 : 240 : iln0=iln
1045 : 240 : ilang=1+indlmn(1,ilmn,itypat)
1046 : 240 : iproj=indlmn(3,ilmn,itypat)
1047 [ + - ]: 240 : if(iproj>0)then
1048 : : ! First set of terms
1049 : : ! number of unique tensor components
1050 : 240 : ilangx=((ilang+4)*(ilang+5))/2
1051 : :
1052 [ + + ]: 5760 : do ii=1,ilangx
1053 : : ! Get the starting address for the relevant tensor
1054 : 5520 : jj=ii+((ilang+3)*(ilang+4)*(ilang+5))/6
1055 : 5520 : ig=ipw1
1056 : 5520 : iffkg=iffkg+1
1057 [ + + ]: 823860 : do ipw=1,nincpw
1058 : 818340 : ffkg(ipw,iffkg)=ffnl(ig,3,ilmn,itypat)*kpgx(ipw,jj)
1059 : 823860 : ig=ig+1
1060 : : end do
1061 [ + + ]: 5520 : if(ilang==1 .or. ilang==3)parity(iffkg)=2
1062 [ + + ]: 5760 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
1063 : : end do
1064 : :
1065 : : ! Second set of terms
1066 : : ! number of unique tensor components
1067 : 240 : ilangx=((ilang+1)*(ilang+2))/2
1068 : :
1069 [ + + ]: 2016 : do ii=1,ilangx
1070 : : ! Get the starting address for the relevant tensor
1071 : 1776 : jj=ii+((ilang)*(ilang+1)*(ilang+2))/6
1072 : 1776 : ig=ipw1
1073 : 1776 : iffkg=iffkg+1
1074 [ + + ]: 265068 : do ipw=1,nincpw
1075 : 263292 : ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
1076 : 265068 : ig=ig+1
1077 : : end do
1078 [ + + ]: 1776 : if(ilang==1 .or. ilang==3)parity(iffkg)=1
1079 [ + + ]: 2016 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
1080 : : end do
1081 : :
1082 : : ! Third set of terms
1083 : : ! number of unique tensor components
1084 : 240 : ilangx=((ilang+3)*(ilang+4))/2
1085 : :
1086 [ + + ]: 4272 : do ii=1,ilangx
1087 : : ! Get the starting address for the relevant tensor
1088 : 4032 : jj=ii+((ilang+2)*(ilang+3)*(ilang+4))/6
1089 : 4032 : ig=ipw1
1090 : 4032 : iffkg=iffkg+1
1091 [ + + ]: 601776 : do ipw=1,nincpw
1092 : 597744 : ffkg(ipw,iffkg)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
1093 : 601776 : ig=ig+1
1094 : : end do
1095 [ + + ]: 4032 : if(ilang==1 .or. ilang==3)parity(iffkg)=1
1096 [ + + ]: 4272 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
1097 : : end do
1098 : :
1099 : : ! End condition if(iproj>0)
1100 : : end if
1101 : : ! End loop on ilang (ilmn)
1102 : : end if
1103 : : end if
1104 : : end do
1105 : :
1106 : : ! End condition of 2nd strain derivatives
1107 : : end if
1108 : :
1109 : : !For ddk perturbation, treat the part that depend on ffnl(:,1,..)
1110 : : !no contribution from s state
1111 [ + - - + ]: 91404 : if(nlang>=2 .and. choice==5)then
1112 : : iln0=0
1113 [ # # ]: 0 : do ilmn=1,lmnmax
1114 [ # # ]: 0 : if (ispinor==indlmn(6,ilmn,itypat)) then
1115 : 0 : iln=indlmn(5,ilmn,itypat)
1116 [ # # ]: 0 : if (iln>iln0) then
1117 : 0 : iln0=iln
1118 : 0 : ilang=1+indlmn(1,ilmn,itypat)
1119 [ # # ]: 0 : if (ilang>=2) then
1120 : 0 : iproj=indlmn(3,ilmn,itypat)
1121 [ # # ]: 0 : if(iproj>0)then
1122 : 0 : ilang2=(ilang*(ilang-1))/2
1123 : :
1124 [ # # ]: 0 : do ii=1,ilang2
1125 : : ! Get the starting address for the relevant tensor
1126 : 0 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
1127 : 0 : ig=ipw1
1128 : 0 : iffkg=iffkg+1
1129 [ # # ]: 0 : do ipw=1,nincpw
1130 : 0 : ffkg(ipw,iffkg)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
1131 : 0 : ig=ig+1
1132 : : end do
1133 [ # # ]: 0 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
1134 [ # # ]: 0 : if(ilang==3)parity(iffkg)=1
1135 : : end do
1136 : :
1137 : : ! End condition if(iproj>0)
1138 : : end if
1139 : : ! End loop on ilang>=2
1140 : : end if
1141 : : end if
1142 : : end if
1143 : : end do
1144 : : ! End condition of p,d or f state
1145 : : end if
1146 : :
1147 : 91404 : end subroutine mkffkg
1148 : : !!***
1149 : :
1150 : : end module m_mkffkg
1151 : : !!***
|
