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 183818278 : subroutine dfpt_mkffkg(choice,ffkg,ffnl,gmet,idir,indlmn,ipw1,ispinor,itypat,&
86 183818278 : & kg_k,kpg_k,kpgx,kpt,lmnmax,mblkpw,ndgxdt,nffkg,nffnl,nincpw,nkpg,nlang,&
87 183818278 : & 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 183818278 : jj=0;ilangx=0
109 :
110 : !This will be useless after all the modifications have been done
111 27086708953 : do ipw=1,nincpw
112 27086708953 : 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 183818278 : & .or. choice==6 .or. choice==23) then
118 163609786 : 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 163609786 : ig=ipw1
124 24416698703 : do ipw=1,nincpw
125 24253088917 : kpgx(ipw,2)=kpt(1)+dble(kg_k(1,ig))
126 24253088917 : kpgx(ipw,3)=kpt(2)+dble(kg_k(2,ig))
127 24253088917 : kpgx(ipw,4)=kpt(3)+dble(kg_k(3,ig))
128 24436907195 : ig=ig+1
129 : end do
130 : end if
131 : end if
132 183818278 : 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 95660788 : 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 13948772665 : do ipw=1,nincpw
144 13853111877 : kpgx(ipw, 5) = kpgx(ipw, 2)*kpgx(ipw, 2)
145 13853111877 : kpgx(ipw, 6) = kpgx(ipw, 3)*kpgx(ipw, 3)
146 13853111877 : kpgx(ipw, 7) = kpgx(ipw, 4)*kpgx(ipw, 4)
147 13853111877 : kpgx(ipw, 8) = kpgx(ipw, 4)*kpgx(ipw, 3)
148 13853111877 : kpgx(ipw, 9) = kpgx(ipw, 4)*kpgx(ipw, 2)
149 13948772665 : kpgx(ipw,10) = kpgx(ipw, 3)*kpgx(ipw, 2)
150 : end do
151 : end if
152 : end if
153 183818278 : 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 487277992 : do ipw=1,nincpw
157 483366666 : kpgx(ipw,11) = kpgx(ipw, 5)*kpgx(ipw, 2)
158 483366666 : kpgx(ipw,12) = kpgx(ipw, 6)*kpgx(ipw, 2)
159 483366666 : kpgx(ipw,13) = kpgx(ipw, 7)*kpgx(ipw, 2)
160 483366666 : kpgx(ipw,14) = kpgx(ipw, 8)*kpgx(ipw, 2)
161 483366666 : kpgx(ipw,15) = kpgx(ipw, 9)*kpgx(ipw, 2)
162 483366666 : kpgx(ipw,16) = kpgx(ipw,10)*kpgx(ipw, 2)
163 483366666 : kpgx(ipw,17) = kpgx(ipw, 6)*kpgx(ipw, 3)
164 483366666 : kpgx(ipw,18) = kpgx(ipw, 7)*kpgx(ipw, 3)
165 483366666 : kpgx(ipw,19) = kpgx(ipw, 8)*kpgx(ipw, 3)
166 667184944 : kpgx(ipw,20) = kpgx(ipw, 7)*kpgx(ipw, 4)
167 : end do
168 : end if
169 183818278 : 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 171536330 : do ipw=1,nincpw
175 169672386 : kpgx(ipw,21) = kpgx(ipw, 5)*kpgx(ipw, 5)
176 169672386 : kpgx(ipw,22) = kpgx(ipw, 6)*kpgx(ipw, 5)
177 169672386 : kpgx(ipw,23) = kpgx(ipw, 7)*kpgx(ipw, 5)
178 169672386 : kpgx(ipw,24) = kpgx(ipw, 8)*kpgx(ipw, 5)
179 169672386 : kpgx(ipw,25) = kpgx(ipw, 9)*kpgx(ipw, 5)
180 169672386 : kpgx(ipw,26) = kpgx(ipw,10)*kpgx(ipw, 5)
181 169672386 : kpgx(ipw,27) = kpgx(ipw, 6)*kpgx(ipw,10)
182 169672386 : kpgx(ipw,28) = kpgx(ipw, 7)*kpgx(ipw,10)
183 169672386 : kpgx(ipw,29) = kpgx(ipw, 8)*kpgx(ipw,10)
184 169672386 : kpgx(ipw,30) = kpgx(ipw, 7)*kpgx(ipw, 9)
185 169672386 : kpgx(ipw,31) = kpgx(ipw, 6)*kpgx(ipw, 6)
186 169672386 : kpgx(ipw,32) = kpgx(ipw, 7)*kpgx(ipw, 6)
187 169672386 : kpgx(ipw,33) = kpgx(ipw, 8)*kpgx(ipw, 6)
188 169672386 : kpgx(ipw,34) = kpgx(ipw, 7)*kpgx(ipw, 8)
189 353490664 : kpgx(ipw,35) = kpgx(ipw, 7)*kpgx(ipw, 7)
190 : end do
191 : end if
192 183818278 : 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 194221754 : kpgx(ipw,56) = kpgx(ipw,35)*kpgx(ipw, 4)
219 : end do
220 : end if
221 183818278 : 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 183818278 : 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 1218379489 : do ilmn=1,lmnmax
315 1034561211 : iln=indlmn(5,ilmn,itypat)
316 1218379489 : if (iln>iln0) then
317 782906663 : iln0=iln
318 782906663 : ilang=1+indlmn(1,ilmn,itypat)
319 782906663 : iproj=indlmn(3,ilmn,itypat)
320 782906663 : if(iproj>0)then
321 782906663 : ilang2=(ilang*(ilang+1))/2
322 :
323 782906663 : if(ilang==1)then
324 : ! Treat s-component separately
325 292741458 : ig=ipw1
326 292741458 : iffkg=iffkg+1
327 44453821405 : do ipw=1,nincpw
328 44161079947 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)
329 44453821405 : ig=ig+1
330 : end do
331 292741458 : parity(iffkg)=2
332 : else
333 : ! Treat other components (could be made faster by treating explicitely
334 : ! each angular momentum)
335 2505981302 : do ii=1,ilang2
336 : ! Get the starting address for the relevant tensor
337 2015816097 : jj=ii+((ilang-1)*ilang*(ilang+1))/6
338 2015816097 : ig=ipw1
339 2015816097 : iffkg=iffkg+1
340 >29817*10^7 : do ipw=1,nincpw
341 >29615*10^7 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
342 >29817*10^7 : ig=ig+1
343 : end do
344 2015816097 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
345 2505981302 : 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 183818278 : 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 183818278 : if((choice==2.or.choice==23) .and. ndgxdt/=1)then
363 12157563 : do ii=1,nffkge
364 1751357957 : do ipw=1,nincpw
365 1739884833 : ffkg(iffkg+1,ipw)=ffkg(ii,ipw)*kpgx(ipw,2)
366 1739884833 : ffkg(iffkg+2,ipw)=ffkg(ii,ipw)*kpgx(ipw,3)
367 1751357957 : ffkg(iffkg+3,ipw)=ffkg(ii,ipw)*kpgx(ipw,4)
368 : end do
369 11473124 : parity(iffkg+1)=3-parity(ii)
370 11473124 : parity(iffkg+2)=parity(iffkg+1)
371 11473124 : parity(iffkg+3)=parity(iffkg+1)
372 12157563 : 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 183818278 : 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 183818278 : 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 183818278 : if(choice==3 .or. choice==5 .or. choice==6 .or. choice==23)then
418 :
419 : iln0=0
420 34368187 : do ilmn=1,lmnmax
421 34368187 : if (ispinor==indlmn(6,ilmn,itypat)) then
422 24638987 : iln=indlmn(5,ilmn,itypat)
423 24638987 : if (iln>iln0) then
424 24638987 : iln0=iln
425 24638987 : ilang=1+indlmn(1,ilmn,itypat)
426 24638987 : iproj=indlmn(3,ilmn,itypat)
427 24638987 : if(iproj>0)then
428 : ! number of unique tensor components
429 24638987 : if(choice==3 .or. choice==6 .or. choice==23)ilangx=((ilang+2)*(ilang+3))/2
430 24638987 : if(choice==5)ilangx=(ilang*(ilang+1))/2
431 :
432 185726656 : do ii=1,ilangx
433 : ! Get the starting address for the relevant tensor
434 122731083 : if(choice==3 .or. choice==6 .or. choice==23)jj=ii+((ilang+1)*(ilang+2)*(ilang+3))/6
435 161087669 : if(choice==5)jj=ii+((ilang-1)*ilang*(ilang+1))/6
436 161087669 : ig=ipw1
437 161087669 : iffkg=iffkg+1
438 : if(choice==3 .or. choice==6 .or. choice==23)then
439 11648646404 : do ipw=1,nincpw
440 11525915321 : ffkg(iffkg,ipw)=ffnl(ig,2,ilmn,itypat)*kpgx(ipw,jj)
441 11648646404 : ig=ig+1
442 : end do
443 : else
444 3782945300 : 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 3744588714 : & kpgx(ipw,4)*gmet(3,idir) )
449 3782945300 : ig=ig+1
450 : end do
451 : end if
452 161087669 : if(ilang==1 .or. ilang==3)parity(iffkg)=2
453 161087669 : if(ilang==2 .or. ilang==4)parity(iffkg)=1
454 185726656 : 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 183818278 : 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 183818278 : if(nlang>=2 .and. choice==5)then
555 : iln0=0
556 17087660 : do ilmn=1,lmnmax
557 17087660 : if (ispinor==indlmn(6,ilmn,itypat)) then
558 11908224 : iln=indlmn(5,ilmn,itypat)
559 11908224 : if (iln>iln0) then
560 11908224 : iln0=iln
561 11908224 : ilang=1+indlmn(1,ilmn,itypat)
562 11908224 : if (ilang>=2) then
563 7830966 : iproj=indlmn(3,ilmn,itypat)
564 7830966 : if(iproj>0)then
565 7830966 : ilang2=(ilang*(ilang-1))/2
566 :
567 22234852 : do ii=1,ilang2
568 : ! Get the starting address for the relevant tensor
569 14403886 : jj=ii+((ilang-2)*(ilang-1)*ilang)/6
570 14403886 : ig=ipw1
571 14403886 : iffkg=iffkg+1
572 1386787770 : do ipw=1,nincpw
573 1372383884 : ffkg(iffkg,ipw)=ffnl(ig,1,ilmn,itypat)*kpgx(ipw,jj)
574 1386787770 : ig=ig+1
575 : end do
576 14403886 : if(ilang==2 .or. ilang==4)parity(iffkg)=2
577 22234852 : 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 183818278 : 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 : !!***
|
