Line data Source code
1 : !!****m* ABINIT/m_opernld_ylm
2 : !! NAME
3 : !! m_opernld_ylm
4 : !!
5 : !! FUNCTION
6 : !!
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MT)
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_opernld_ylm
23 :
24 : use defs_basis
25 : use m_errors
26 : use m_abicore
27 :
28 : implicit none
29 :
30 : private
31 : !!***
32 :
33 : public :: opernld_ylm
34 : !!***
35 :
36 : contains
37 : !!***
38 :
39 : !!****f* ABINIT/opernld_ylm
40 : !! NAME
41 : !! opernld_ylm
42 : !!
43 : !! FUNCTION
44 : !! * Operate with the non-local part of the hamiltonian,
45 : !! in order to get contributions to energy/forces/stress/dyn.matrix/elst tens.
46 : !! from projected scalars
47 : !! * Operate with the non-local projectors and the overlap matrix Sij
48 : !! in order to get contributions to <c|S|c>
49 : !! from projected scalars
50 : !!
51 : !! INPUTS
52 : !! choice=chooses possible output
53 : !! cplex=1 if <p_lmn|c> scalars are real (equivalent to istwfk>1)
54 : !! 2 if <p_lmn|c> scalars are complex
55 : !! cplex_fac=1 if gxfac scalars are real, 2 if gxfac scalars are complex
56 : !! d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)=2nd gradients of projected scalars
57 : !! dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)=gradients of projected scalars
58 : !! dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
59 : !! related to Vnl (NL operator)
60 : !! dgxdtfac_sij(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)=gradients of reduced projected scalars
61 : !! related to Sij (overlap)
62 : !! gx(cplex,nlmn,nincat,nspinor)= projected scalars
63 : !! gxfac(cplex_fac,nlmn,nincat,nspinor)= reduced projected scalars related to Vnl (NL operator)
64 : !! gxfac_sij(cplex,nlmn,nincat,nspinor)= reduced projected scalars related to Sij (overlap)
65 : !! ia3=gives the absolute number of the first atom in the subset presently treated
66 : !! natom=number of atoms in cell
67 : !! nd2gxdt=second dimension of d2gxdt
68 : !! ndgxdt=second dimension of dgxdt
69 : !! ndgxdtfac=second dimension of dgxdtfac
70 : !! nincat=number of atoms in the subset here treated
71 : !! nlmn=number of (l,m,n) numbers for current type of atom
72 : !! nnlout=dimension of enlout
73 : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
74 : !! paw_opt= define the nonlocal operator concerned with:
75 : !! paw_opt=0 : Norm-conserving Vnl (use of Kleinman-Bylander ener.)
76 : !! paw_opt=1 : PAW nonlocal part of H (use of Dij coeffs)
77 : !! paw_opt=2 : PAW: (Vnl-lambda.Sij) (Sij=overlap matrix)
78 : !! paw_opt=3 : PAW overlap matrix (Sij)
79 : !! paw_opt=4 : both PAW nonlocal part of H (Dij) and overlap matrix (Sij)
80 : !!
81 : !! OUTPUT
82 : !! (see side effects)
83 : !!
84 : !! SIDE EFFECTS
85 : !! --If (paw_opt==0, 1 or 2)
86 : !! enlout(nnlout)= contribution to the non-local part of the following properties:
87 : !! if choice=1 : enlout(1) -> the energy
88 : !! if choice=2 : enlout(3*natom) -> 1st deriv. of energy wrt atm. pos (forces)
89 : !! if choice=3 : enlout(6) -> 1st deriv. of energy wrt strain (stresses)
90 : !! if choice=4 : enlout(6*natom) -> 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
91 : !! if choice=23: enlout(6+3*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
92 : !! 1st deriv. of energy wrt strain (stresses)
93 : !! if choice=24: enlout(9*natom) -> 1st deriv. of energy wrt atm. pos (forces) and
94 : !! 2nd deriv. of energy wrt 2 atm. pos (dyn. mat.)
95 : !! if choice=5 : enlout(3) -> 1st deriv. of energy wrt k
96 : !! if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
97 : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
98 : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
99 : !! if choice=6 : enlout(36+18*natom) -> 2nd deriv. of energy wrt 2 strains (elast. tensor) and
100 : !! 2nd deriv. of energy wrt to atm. pos and strain (internal strain)
101 : !! if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
102 : !! if choice=81: enlout(18) -> 2nd deriv. of energy wrt k and right k
103 : !! --If (paw_opt==3)
104 : !! if choice=1 : enlout(1) -> contribution to <c|S|c> (note: not including <c|c>)
105 : !! if choice=2 : enlout(3*natom) -> contribution to <c|dS/d_atm.pos|c>
106 : !! if choice=53: enlout(3) -> 1st deriv. (twist) of energy wrt k
107 : !! if choice=54: enlout(18*natom) -> 2nd deriv. of energy wrt atm. pos and right k (Born eff. charge)
108 : !! if choice=55: enlout(36) -> 2nd deriv. of energy wrt strain and right k (piezoelastic tensor)
109 : !! if choice=8 : enlout(6) -> 2nd deriv. of energy wrt 2 k
110 : !! if choice=81: enlout(18) -> 2nd deriv. of energy wrt k and right k
111 : !! --If (paw_opt==4)
112 : !! not available
113 : !!
114 : !! NOTES
115 : !! Operate for one type of atom, and within this given type of atom,
116 : !! for a subset of at most nincat atoms.
117 : !!
118 : !! SOURCE
119 :
120 738502 : subroutine opernld_ylm(choice,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
121 738502 : & enlk,enlout,fnlk,gx,gxfac,gxfac_sij,ia3,natom,ndat_left,nd2gxdt,ndgxdt,&
122 : & ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt,strnlk,&
123 : & enlout_im)
124 :
125 : !Arguments ------------------------------------
126 : !scalars
127 : integer,intent(in) :: choice,cplex,cplex_fac,ia3,natom,nd2gxdt,ndgxdt
128 : integer,intent(in) :: ndgxdtfac,nincat,nlmn,nnlout,nspinor,paw_opt
129 : integer,intent(in) :: ndat_left
130 : real(dp),intent(inout) :: enlk
131 : !arrays
132 : real(dp),intent(in) :: d2gxdt(cplex,nd2gxdt,nlmn,nincat,nspinor)
133 : real(dp),intent(in) :: dgxdt(cplex,ndgxdt,nlmn,nincat,nspinor)
134 : real(dp),intent(in) :: dgxdtfac(cplex_fac,ndgxdtfac,nlmn,nincat,nspinor)
135 : real(dp),intent(in) :: dgxdtfac_sij(cplex,ndgxdtfac,nlmn,nincat,nspinor*(paw_opt/3))
136 : real(dp),intent(in) :: gx(cplex,nlmn,nincat,nspinor*ndat_left),gxfac(cplex_fac,nlmn,nincat,nspinor)
137 : real(dp),intent(in) :: gxfac_sij(cplex,nlmn,nincat,nspinor*(paw_opt/3))
138 : real(dp),intent(inout) :: ddkk(6),enlout(nnlout*ndat_left),fnlk(3*natom),strnlk(6)
139 : real(dp),intent(inout),optional :: enlout_im(nnlout*ndat_left)
140 :
141 : !Local variables-------------------------------
142 : !scalars
143 : integer :: ia,iashift,idat_left,ilmn,iplex,ishift,ispinor,mu,mua,mua1,mua2,mub,mushift,mut,muu
144 : integer :: nu,nushift
145 : real(dp) :: dummy
146 : !arrays
147 : integer,parameter :: alpha(6)=(/1,2,3,3,3,2/),beta(6)=(/1,2,3,2,1,1/)
148 : integer,parameter :: gamma(3,3)=reshape((/1,6,5,6,2,4,5,4,3/),(/3,3/))
149 : integer,parameter :: twist_dir(6)=(/2,3,3,1,1,2/)
150 1477004 : real(dp) :: d2gx(cplex),enlj(6*cplex),gxfacj(cplex)
151 738502 : real(dp),allocatable :: enljj(:)
152 738502 : complex(dp),allocatable :: cft(:,:), cfu(:,:)
153 : ! *************************************************************************
154 :
155 738502 : ABI_CHECK(cplex_fac>=cplex,'BUG: invalid cplex_fac<cplex!')
156 :
157 738502 : if (paw_opt==0.or.paw_opt==1.or.paw_opt==2) then
158 :
159 : ! ============== Accumulate the non-local energy ===============
160 473638 : if (choice==1) then
161 0 : if (present(enlout_im).and.cplex==2) then ! cplex=cplex_fac=2
162 0 : do idat_left=1,ndat_left
163 0 : do ispinor=1,nspinor
164 0 : do ia=1,nincat
165 0 : do ilmn=1,nlmn
166 0 : enlout (idat_left)=enlout (idat_left)+gxfac(1,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor)
167 0 : enlout (idat_left)=enlout (idat_left)+gxfac(2,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor)
168 0 : enlout_im(idat_left)=enlout_im(idat_left)+gxfac(2,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor)
169 0 : enlout_im(idat_left)=enlout_im(idat_left)-gxfac(1,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor)
170 : end do
171 : end do
172 : end do
173 : end do
174 0 : else if (present(enlout_im).and.cplex_fac==2) then ! cplex=1,cplex_fac=2
175 0 : do idat_left=1,ndat_left
176 0 : do ispinor=1,nspinor
177 0 : do ia=1,nincat
178 0 : do ilmn=1,nlmn
179 0 : enlout (idat_left)=enlout (idat_left)+gxfac(1,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor)
180 0 : enlout_im(idat_left)=enlout_im(idat_left)+gxfac(2,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor)
181 : end do
182 : end do
183 : end do
184 : end do
185 : else ! only the real part is needed or the imaginary part is zero
186 0 : do idat_left=1,ndat_left
187 0 : do ispinor=1,nspinor
188 0 : do ia=1,nincat
189 0 : do ilmn=1,nlmn
190 0 : do iplex=1,cplex
191 0 : enlout(idat_left)=enlout(idat_left)+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor+(idat_left-1)*nspinor)
192 : end do
193 : end do
194 : end do
195 : end do
196 : end do
197 : end if
198 : end if
199 :
200 : ! ============ Accumulate the forces contributions =============
201 473638 : if (choice==2.or.choice==23.or.choice==24) then
202 381305 : ishift=0;if (choice==23) ishift=6
203 828268 : do ispinor=1,nspinor
204 1950523 : do ia=1,nincat
205 4119688 : enlj(1:3)=zero
206 1029922 : iashift=3*(ia+ia3-2)+ishift
207 11176397 : do ilmn=1,nlmn
208 41615822 : do mu=1,3
209 : dummy = zero ! Dummy needed here to get the correct forces with intel -O3
210 80196027 : do iplex=1,cplex
211 : !enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
212 80196027 : dummy=dummy+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
213 : end do
214 40585900 : enlj(mu)=enlj(mu)+dummy
215 : end do
216 : end do
217 4566651 : enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3)
218 : end do
219 : end do
220 : end if
221 :
222 : ! ======== Accumulate the stress tensor contributions ==========
223 473638 : if (choice==3.or.choice==23) then
224 760606 : enlj(1:6)=zero
225 228414 : do ispinor=1,nspinor
226 411979 : do ia=1,nincat
227 2697174 : do ilmn=1,nlmn
228 7117731 : gxfacj(1:cplex)=gxfac(1:cplex,ilmn,ia,ispinor)
229 7117731 : do iplex=1,cplex
230 7117731 : enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor)
231 : end do
232 16940536 : do mu=1,6
233 45100239 : do iplex=1,cplex
234 42706386 : enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor)
235 : end do
236 : end do
237 : end do
238 : end do
239 : end do
240 760606 : enlout(1:6)=enlout(1:6)+two*enlj(1:6)
241 : end if
242 :
243 : ! ====== Accumulate the dynamical matrix contributions =========
244 473638 : if (choice==4.or.choice==24) then
245 42584 : ishift=0;if (choice==24) ishift=3*natom
246 85232 : do ispinor=1,nspinor
247 130920 : do ia=1,nincat
248 319816 : enlj(1:6)=zero
249 45688 : iashift=6*(ia+ia3-2)+ishift
250 514024 : do ilmn=1,nlmn
251 3324040 : do mu=1,6
252 2810016 : mua=alpha(mu);mub=beta(mu)
253 8898384 : do iplex=1,cplex
254 : enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
255 8430048 : & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,mua,ilmn,ia,ispinor)
256 : end do
257 : end do
258 : end do
259 362464 : enlout(iashift+1:iashift+6)=enlout(iashift+1:iashift+6)+two*enlj(1:6)
260 : end do
261 : end do
262 : end if
263 :
264 : ! ======== Accumulate the contributions of derivatives of E wrt to k ==========
265 473638 : if (choice==5) then
266 4608 : enlj(1:3)=zero
267 2304 : do ispinor=1,nspinor
268 3456 : do ia=1,nincat
269 2304 : if(cplex==2)then
270 10368 : do ilmn=1,nlmn
271 38016 : do mu=1,3
272 92160 : do iplex=1,cplex
273 82944 : enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor)
274 : end do
275 : end do
276 : end do
277 : ! If cplex=1, dgxdt is pure imaginary; thus there is no contribution
278 0 : else if (cplex_fac==2) then
279 0 : do ilmn=1,nlmn
280 0 : do mu=1,3
281 0 : enlj(mu)=enlj(mu)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
282 : end do
283 : end do
284 : end if
285 : end do
286 : end do
287 4608 : enlout(1:3)=enlout(1:3)+two*enlj(1:3)
288 : end if
289 :
290 : ! ======== Accumulate the contributions of partial derivatives of E wrt to k ==========
291 : ! Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k
292 473638 : if (choice==51.or.choice==52) then
293 0 : enlj(1:6)=zero
294 0 : do ispinor=1,nspinor
295 0 : do ia=1,nincat
296 0 : if(cplex==2)then
297 0 : do ilmn=1,nlmn
298 0 : do mu=1,3
299 : enlj(2*mu-1)=enlj(2*mu-1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) &
300 0 : & +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor)
301 : enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) &
302 0 : & -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
303 : end do
304 : end do
305 0 : else if (cplex_fac==2) then
306 0 : do ilmn=1,nlmn
307 0 : do mu=1,3
308 0 : enlj(2*mu-1)=enlj(2*mu-1)+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
309 0 : enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
310 : end do
311 : end do
312 0 : else if (cplex_fac==1) then
313 0 : do ilmn=1,nlmn
314 0 : do mu=1,3
315 0 : enlj(2*mu )=enlj(2*mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
316 : end do
317 : end do
318 : end if
319 : end do
320 : end do
321 0 : if (choice==52) then
322 0 : enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6)
323 : end if
324 0 : enlout(1:6)=enlout(1:6)+enlj(1:6)
325 : end if
326 :
327 : ! ======== Accumulate the contributions of twist derivatives of E wrt to k ==========
328 : ! accumulate <u|dp_i/dk_(idir+1)>D_ij<dp_j/dk(idir+2)|u>
329 473638 : if (choice==53) then
330 0 : enlj(:)=zero
331 0 : ABI_MALLOC(cft,(3,nlmn))
332 0 : ABI_MALLOC(cfu,(3,nlmn))
333 : ! If cplex=1, dgxdt is pure imaginary;
334 : ! If cplex_fac=1, dgxdtfac is pure imaginary;
335 0 : do ispinor=1,nspinor
336 0 : do ia=1,nincat
337 0 : if(cplex==2)then
338 0 : cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
339 : else
340 0 : cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
341 : end if
342 0 : if(cplex_fac==2)then
343 0 : cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor))
344 : else
345 0 : cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor))
346 : end if
347 0 : do ilmn=1,nlmn
348 0 : do mu=1,3
349 0 : mut = twist_dir(2*mu-1)
350 0 : muu = twist_dir(2*mu)
351 0 : if (cplex == 2) then
352 0 : enlj(2*mu-1) = enlj(2*mu-1) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
353 0 : enlj(2*mu) = enlj(2*mu) + aimag(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
354 : else
355 0 : enlj(mu) = enlj(mu) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
356 : end if
357 : end do ! end loop over mu=1,3
358 : end do ! end loop over ilmn states
359 : end do ! end loop over ia atoms
360 : end do ! end loop over ispinor
361 0 : do mu = 1, 3
362 0 : if (cplex == 2) then
363 0 : enlout(2*mu-1)=enlout(2*mu-1)+enlj(2*mu-1)
364 0 : enlout(2*mu) =enlout(2*mu) +enlj(2*mu)
365 : else
366 0 : enlout(mu)=enlout(mu)+enlj(mu)
367 : end if
368 : end do ! end loop over mu = 1, 3
369 0 : ABI_FREE(cft)
370 0 : ABI_FREE(cfu)
371 : end if
372 :
373 : ! ====== Accumulate the effective charges contributions =========
374 473638 : if (choice==54) then
375 0 : ABI_MALLOC(enljj,(18))
376 0 : do ispinor=1,nspinor
377 0 : do ia=1,nincat
378 0 : enljj(1:18)=zero
379 0 : iashift=18*(ia+ia3-2)
380 : ! If cplex=1, dgxdt is real for atm. pos, pure imaginary for k;
381 : ! If cplex_fac=1, dgxdtfac is pure imaginary for k;
382 0 : if(cplex==2.and.cplex_fac==2) then
383 0 : do ilmn=1,nlmn
384 : mu=1;nu=1
385 0 : do mua=1,3 ! atm. pos
386 0 : do mub=1,3 ! k
387 : enljj(nu)=enljj(nu) &
388 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
389 : & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) &
390 : & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) &
391 0 : & +gxfac(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor)
392 : enljj(nu+1)=enljj(nu+1) &
393 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) &
394 : & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
395 : & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) &
396 0 : & -gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
397 0 : mu=mu+1;nu=nu+2
398 : end do
399 : end do
400 : end do
401 0 : else if(cplex==1.and.cplex_fac==2)then
402 0 : do ilmn=1,nlmn
403 : mu=1;nu=1
404 0 : do mua=1,3 ! atm. pos
405 0 : do mub=1,3 ! k
406 : enljj(nu)=enljj(nu) &
407 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
408 0 : & +gxfac(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
409 : enljj(nu+1)=enljj(nu+1) &
410 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,3+mub,ilmn,ia,ispinor) &
411 0 : & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
412 0 : mu=mu+1;nu=nu+2
413 : end do
414 : end do
415 : end do
416 0 : else if(cplex==1.and.cplex_fac==1)then
417 0 : do ilmn=1,nlmn
418 : mu=1;nu=1
419 0 : do mua=1,3 ! atm. pos
420 0 : do mub=1,3 ! k
421 : enljj(nu+1)=enljj(nu+1) &
422 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,3+mub,ilmn,ia,ispinor) &
423 0 : & +gxfac(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
424 0 : mu=mu+1;nu=nu+2
425 : end do
426 : end do
427 : end do
428 : end if
429 0 : enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18)
430 : end do
431 : end do
432 0 : ABI_FREE(enljj)
433 : end if
434 :
435 : ! ====== Accumulate the piezoelectric tensor contributions =========
436 473638 : if (choice==55) then
437 0 : ABI_MALLOC(enljj,(36))
438 0 : do ispinor=1,nspinor
439 0 : do ia=1,nincat
440 0 : enljj(1:36)=zero;enlj(:)=zero
441 : ! If cplex=1, dgxdt is real for strain, pure imaginary for k;
442 : ! If cplex_fac=1, dgxdtfac is pure imaginary for k;
443 0 : if(cplex==2.and.cplex_fac==2) then
444 0 : do ilmn=1,nlmn
445 : ! First compute 2nd-derivative contribution
446 : mu=1
447 0 : do mua=1,6 ! strain (lambda,nu)
448 0 : mua1=alpha(mua) ! (nu)
449 0 : mua2=beta(mua) ! (lambda)
450 0 : do mub=1,3 ! k (mu)
451 0 : muu=3*(gamma(mua1,mub)-1)+mua2
452 0 : mut=3*(gamma(mua2,mub)-1)+mua1
453 : d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) &
454 0 : & +d2gxdt(1:cplex,mut,ilmn,ia,ispinor))
455 : enljj(mu)=enljj(mu) &
456 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
457 : & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) &
458 0 : & +gxfac(1,ilmn,ia,ispinor)*d2gx(1)+gxfac(2,ilmn,ia,ispinor)*d2gx(2)
459 : enljj(mu+1)=enljj(mu+1) &
460 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) &
461 : & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
462 0 : & +gxfac(1,ilmn,ia,ispinor)*d2gx(2)-gxfac(2,ilmn,ia,ispinor)*d2gx(1)
463 0 : mu=mu+2
464 : end do
465 : end do
466 : ! Then store 1st-derivative contribution
467 : mu=1
468 0 : do nu=1,3
469 : enlj(mu )=enlj(mu )+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) &
470 0 : & +gxfac(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor)
471 : enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) &
472 0 : & -gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
473 0 : mu=mu+2
474 : end do
475 : end do
476 0 : else if(cplex==1.and.cplex_fac==2)then
477 0 : do ilmn=1,nlmn
478 : ! First compute 2nd-derivative contribution
479 : mu=1
480 0 : do mua=1,6 ! strain (lambda,nu)
481 0 : mua1=alpha(mua) ! (nu)
482 0 : mua2=beta(mua) ! (lambda)
483 0 : do mub=1,3 ! k (mu)
484 0 : muu=3*(gamma(mua1,mub)-1)+mua2
485 0 : mut=3*(gamma(mua2,mub)-1)+mua1
486 0 : d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor))
487 : enljj(mu)=enljj(mu) &
488 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
489 0 : & +gxfac(2,ilmn,ia,ispinor)*d2gx(1)
490 : enljj(mu+1)=enljj(mu+1) &
491 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(2,6+mub,ilmn,ia,ispinor) &
492 0 : & +gxfac(1,ilmn,ia,ispinor)*d2gx(1)
493 0 : mu=mu+2
494 : end do
495 : end do
496 : ! Then store 1st-derivative contribution
497 : mu=1
498 0 : do nu=1,3
499 0 : enlj(mu )=enlj(mu )+gxfac(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
500 0 : enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
501 0 : mu=mu+2
502 : end do
503 : end do
504 0 : else if(cplex==1.and.cplex_fac==1)then
505 0 : do ilmn=1,nlmn
506 : mu=1
507 0 : do mua=1,6 ! strain (lambda,nu)
508 0 : mua1=alpha(mua) ! (nu)
509 0 : mua2=beta(mua) ! (lambda)
510 0 : do mub=1,3 ! k (mu)
511 0 : muu=3*(gamma(mua1,mub)-1)+mua2
512 0 : mut=3*(gamma(mua2,mub)-1)+mua1
513 0 : d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor))
514 : enljj(mu+1)=enljj(mu+1) &
515 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac(1,6+mub,ilmn,ia,ispinor) &
516 0 : & +gxfac(1,ilmn,ia,ispinor)*d2gx(1)
517 0 : mu=mu+2
518 : end do
519 : end do
520 : ! Then store 1st-derivative contribution
521 : mu=1
522 0 : do nu=1,3
523 0 : enlj(mu+1)=enlj(mu+1)+gxfac(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
524 0 : mu=mu+2
525 : end do
526 : end do
527 : end if
528 0 : enlout(1:36)=enlout(1:36)+enljj(1:36)
529 0 : ddkk(1:6)=ddkk(1:6)+enlj(1:6)
530 : end do
531 : end do
532 0 : ABI_FREE(enljj)
533 : end if
534 :
535 : ! ======= Accumulate the elastic tensor contributions ==========
536 473638 : if (choice==6) then
537 5760 : do ispinor=1,nspinor
538 11136 : do ia=1,nincat
539 5376 : iashift=3*(ia+ia3-2)
540 53312 : do ilmn=1,nlmn
541 135168 : do iplex=1,cplex
542 135168 : enlk=enlk+gxfac(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor)
543 : end do
544 180224 : enlj(1:3)=zero
545 180224 : do mu=1,3
546 450560 : do iplex=1,cplex
547 405504 : enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,6+mu,ilmn,ia,ispinor)
548 : end do
549 : end do
550 180224 : fnlk(iashift+1:iashift+3)=fnlk(iashift+1:iashift+3)+two*enlj(1:3)
551 315392 : enlj(1:6)=zero
552 315392 : do mu=1,6
553 856064 : do iplex=1,cplex
554 811008 : enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor)
555 : end do
556 : end do
557 315392 : strnlk(1:6)=strnlk(1:6)+two*enlj(1:6)
558 320768 : do mub=1,6
559 270336 : mushift=6*(mub-1);nushift=(3*natom+6)*(mub-1)
560 1892352 : do mua=1,6
561 1622016 : mu=mushift+mua;nu=nushift+mua
562 5136384 : do iplex=1,cplex
563 : enlout(nu)=enlout(nu)+two* &
564 : & (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
565 4866048 : & +dgxdtfac(iplex,mua,ilmn,ia,ispinor)*dgxdt(iplex,mub,ilmn,ia,ispinor))
566 : end do
567 : end do
568 270336 : mushift=36+3*(mub-1);nushift=6+iashift+(3*natom+6)*(mub-1)
569 1126400 : do mua=1,3
570 811008 : mu=mushift+mua;nu=nushift+mua
571 2703360 : do iplex=1,cplex
572 : enlout(nu)=enlout(nu)+two* &
573 : & (gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
574 2433024 : & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,6+mua,ilmn,ia,ispinor))
575 : end do
576 : end do
577 : end do
578 : end do
579 : end do
580 : end do
581 : end if
582 :
583 : ! ======== Accumulate the contributions of 2nd-derivatives of E wrt to k ==========
584 473638 : if (choice==8) then
585 2304 : ABI_MALLOC(cft,(3,nlmn))
586 1536 : ABI_MALLOC(cfu,(3,nlmn))
587 1536 : do ispinor=1,nspinor
588 2304 : do ia=1,nincat
589 5376 : enlj(1:6)=zero
590 6912 : do ilmn=1,nlmn
591 43776 : do mu=1,6
592 116736 : do iplex=1,cplex
593 110592 : enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)
594 : end do
595 : end do
596 : end do
597 : ! If cplex=1, dgxdt is pure imaginary;
598 : ! If cplex_fac=1, dgxdtfac is pure imaginary;
599 768 : if(cplex==2)then
600 25344 : cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
601 : else
602 0 : cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
603 : end if
604 768 : if(cplex_fac==2)then
605 25344 : cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor))
606 : else
607 0 : cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor))
608 : end if
609 6912 : do ilmn=1,nlmn
610 43776 : do mu=1,6
611 36864 : mua=alpha(mu);mub=beta(mu)
612 43008 : enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn))
613 : end do
614 : end do
615 6144 : enlout(1:6)=enlout(1:6)+two*enlj(1:6)
616 : end do
617 : end do
618 768 : ABI_FREE(cft)
619 768 : ABI_FREE(cfu)
620 : end if
621 :
622 : ! ======== Accumulate the contributions of partial 2nd-derivatives of E wrt to k ==========
623 : ! Full derivative wrt to k1, right derivative wrt to k2
624 473638 : if (choice==81) then
625 0 : ABI_MALLOC(cft,(3,nlmn))
626 0 : ABI_MALLOC(cfu,(6,nlmn))
627 0 : ABI_MALLOC(enljj,(18))
628 0 : do ispinor=1,nspinor
629 0 : do ia=1,nincat
630 0 : enljj(1:18)=zero
631 0 : if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real
632 0 : cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),gxfac(2,1:nlmn,ia,ispinor))
633 : else
634 0 : cft(1,1:nlmn)=cmplx(gxfac(1,1:nlmn,ia,ispinor),zero)
635 : end if
636 0 : if(cplex==2)then !If cplex=1, d2gxdt is pure real
637 0 : cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor))
638 : else
639 0 : cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero)
640 : end if
641 0 : do ilmn=1,nlmn
642 0 : do mu=1,3
643 0 : do nu=1,3
644 0 : muu=3*(mu-1)+nu ; mut=gamma(mu,nu)
645 0 : enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn))
646 0 : enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn))
647 : end do
648 : end do
649 : end do
650 0 : if(cplex==2)then !If cplex=1, dgxdt is pure imaginary
651 0 : cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
652 : else
653 0 : cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
654 : end if
655 0 : if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary
656 0 : cfu(1:3,1:nlmn)=cmplx(dgxdtfac(1,1:3,1:nlmn,ia,ispinor),dgxdtfac(2,1:3,1:nlmn,ia,ispinor))
657 : else
658 0 : cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac(1,1:3,1:nlmn,ia,ispinor))
659 : end if
660 0 : do ilmn=1,nlmn
661 0 : do mu=1,3
662 0 : do nu=1,3
663 0 : muu=3*(mu-1)+nu
664 0 : enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
665 0 : enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
666 : end do
667 : end do
668 : end do
669 0 : enlout(1:18)=enlout(1:18)+enljj(1:18)
670 : end do
671 : end do
672 0 : ABI_FREE(cft)
673 0 : ABI_FREE(cfu)
674 0 : ABI_FREE(enljj)
675 : end if
676 :
677 : end if
678 :
679 738502 : if (paw_opt==3) then
680 :
681 : ! ============== Accumulate contribution to <c|S|c> ===============
682 264864 : if (choice==1) then
683 220960 : if (present(enlout_im).and.cplex==2) then ! cplex=2
684 2217120 : do idat_left=1,ndat_left
685 4213280 : do ispinor=1,nspinor
686 8549120 : do ia=1,nincat
687 54296960 : do ilmn=1,nlmn
688 47744000 : enlout (idat_left)=enlout (idat_left)+gxfac_sij(1,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor)
689 47744000 : enlout (idat_left)=enlout (idat_left)+gxfac_sij(2,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor)
690 47744000 : enlout_im(idat_left)=enlout_im(idat_left)+gxfac_sij(2,ilmn,ia,ispinor)*gx(1,ilmn,ia,ispinor+(idat_left-1)*nspinor)
691 52300800 : enlout_im(idat_left)=enlout_im(idat_left)-gxfac_sij(1,ilmn,ia,ispinor)*gx(2,ilmn,ia,ispinor+(idat_left-1)*nspinor)
692 : end do
693 : end do
694 : end do
695 : end do
696 : else ! only the real part is needed or the imaginary part is zero
697 0 : do idat_left=1,ndat_left
698 0 : do ispinor=1,nspinor
699 0 : do ia=1,nincat
700 0 : do ilmn=1,nlmn
701 0 : do iplex=1,cplex
702 : enlout(idat_left)=enlout(idat_left)+&
703 0 : & gxfac_sij(iplex,ilmn,ia,ispinor)*gx(iplex,ilmn,ia,ispinor+(idat_left-1)*nspinor)
704 : end do
705 : end do
706 : end do
707 : end do
708 : end do
709 : end if
710 : end if
711 :
712 : ! ============== Accumulate contribution to <c|dS/d_atm_pos|c> ===============
713 264864 : if (choice==2.or.choice==23) then
714 13824 : ishift=0;if (choice==23) ishift=6
715 27648 : do ispinor=1,nspinor
716 41472 : do ia=1,nincat
717 55296 : enlj(1:3)=zero
718 13824 : iashift=3*(ia+ia3-2)
719 124416 : do ilmn=1,nlmn
720 456192 : do mu=1,3
721 1105920 : do iplex=1,cplex
722 995328 : enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
723 : end do
724 : end do
725 : end do
726 69120 : enlout(iashift+1:iashift+3)=enlout(iashift+1:iashift+3)+two*enlj(1:3)
727 : end do
728 : end do
729 : end if
730 :
731 : ! ============== Accumulate contribution to <c|dS/d_strain|c> ===============
732 264864 : if (choice==3.or.choice==23) then
733 0 : enlj(1:6)=zero
734 0 : do ispinor=1,nspinor
735 0 : do ia=1,nincat
736 0 : do ilmn=1,nlmn
737 0 : gxfacj(1:cplex)=gxfac_sij(1:cplex,ilmn,ia,ispinor)
738 0 : do iplex=1,cplex
739 0 : enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor)
740 : end do
741 0 : do mu=1,6
742 0 : do iplex=1,cplex
743 0 : enlj(mu)=enlj(mu)+gxfacj(iplex)*dgxdt(iplex,mu,ilmn,ia,ispinor)
744 : end do
745 : end do
746 : end do
747 : end do
748 : end do
749 0 : enlout(1:6)=enlout(1:6)+two*enlj(1:6)
750 : end if
751 :
752 : ! ======== Accumulate the contributions of derivatives of <c|S|c> wrt to k ==========
753 264864 : if (choice==5) then
754 4608 : enlj(1:3)=zero
755 : ! If cplex=1, gxfac is real and dgxdt is pure imaginary; thus there is no contribution
756 1152 : if(cplex==2)then
757 2304 : do ispinor=1,nspinor
758 3456 : do ia=1,nincat
759 11520 : do ilmn=1,nlmn
760 38016 : do mu=1,3
761 92160 : do iplex=1,cplex
762 82944 : enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu,ilmn,ia,ispinor)
763 : end do
764 : end do
765 : end do
766 : end do
767 : end do
768 : end if
769 4608 : enlout(1:3)=enlout(1:3)+two*enlj(1:3)
770 : end if
771 :
772 : ! ====== Accumulate the contributions of left or right derivatives of <c|S|c> wrt to k ==========
773 : ! Choice 51: right derivative wrt to k ; Choice 52: left derivative wrt to k
774 264864 : if (choice==51.or.choice==52) then
775 0 : enlj(1:6)=zero
776 0 : do ispinor=1,nspinor
777 0 : do ia=1,nincat
778 0 : if(cplex==2)then
779 0 : do ilmn=1,nlmn
780 0 : do mu=1,3
781 : enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor) &
782 0 : & +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor)
783 : enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,mu,ilmn,ia,ispinor) &
784 0 : & -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
785 : end do
786 : end do
787 0 : else if (cplex_fac==2) then
788 0 : do ilmn=1,nlmn
789 0 : do mu=1,3
790 0 : enlj(2*mu-1)=enlj(2*mu-1)+gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
791 0 : enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
792 : end do
793 : end do
794 0 : else if (cplex_fac==1) then
795 0 : do ilmn=1,nlmn
796 0 : do mu=1,3
797 0 : enlj(2*mu )=enlj(2*mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,mu,ilmn,ia,ispinor)
798 : end do
799 : end do
800 : end if
801 : end do
802 : end do
803 0 : if (choice==52) then
804 0 : enlj(2)=-enlj(2);enlj(4)=-enlj(4);enlj(6)=-enlj(6)
805 : end if
806 0 : enlout(1:6)=enlout(1:6)+enlj(1:6)
807 : end if
808 :
809 : ! ====== Accumulate the contributions of twist derivatives of <c|S|c> wrt to k ==========
810 : ! Choice 53: <u|dp_i/dk_(idir+1)>Sij<dp_j/dk_(idir+2)|u>
811 264864 : if (choice==53) then
812 0 : enlj(:)=zero
813 0 : ABI_MALLOC(cft,(3,nlmn))
814 0 : ABI_MALLOC(cfu,(3,nlmn))
815 : ! If cplex=1, dgxdt is pure imaginary;
816 : ! If cplex_fac=1, dgxdtfac is pure imaginary;
817 0 : do ispinor=1,nspinor
818 0 : do ia=1,nincat
819 0 : if(cplex==2)then
820 0 : cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
821 : else
822 0 : cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
823 : end if
824 0 : if(cplex_fac==2)then
825 0 : cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor))
826 : else
827 0 : cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor))
828 : end if
829 0 : do ilmn=1,nlmn
830 0 : do mu=1,3
831 0 : mut = twist_dir(2*mu-1)
832 0 : muu = twist_dir(2*mu)
833 0 : if (cplex == 2) then
834 0 : enlj(2*mu-1) = enlj(2*mu-1) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
835 0 : enlj(2*mu) = enlj(2*mu) + aimag(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
836 : else
837 0 : enlj(mu) = enlj(mu) + real(conjg(cft(mut,ilmn))*cfu(muu,ilmn))
838 : end if
839 : end do ! end loop over mu=1,3
840 : end do ! end loop over ilmn states
841 : end do ! end loop over ia atoms
842 : end do ! end loop over ispinor
843 0 : do mu = 1, 3
844 0 : if (cplex == 2) then
845 0 : enlout(2*mu-1)=enlout(2*mu-1)+enlj(2*mu-1)
846 0 : enlout(2*mu) =enlout(2*mu) +enlj(2*mu)
847 : else
848 0 : enlout(mu)=enlout(mu)+enlj(mu)
849 : end if
850 : end do ! end loop over mu = 1, 3
851 0 : ABI_FREE(cft)
852 0 : ABI_FREE(cfu)
853 : end if
854 :
855 : ! ====== Accumulate contribution to <c|d2S/d_atm_pos d_left_k|c> =========
856 264864 : if (choice==54) then
857 22624 : ABI_MALLOC(enljj,(18))
858 45248 : do ispinor=1,nspinor
859 68672 : do ia=1,nincat
860 445056 : enljj(1:18)=zero
861 23424 : iashift=18*(ia+ia3-2)
862 23424 : if(cplex==2) then
863 220736 : do ilmn=1,nlmn
864 : mu=1;nu=1
865 812672 : do mua=1,3 ! atm. pos
866 2565056 : do mub=1,3 ! k
867 : enljj(nu)=enljj(nu) &
868 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) &
869 : & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) &
870 : & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor) &
871 1775808 : & +gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor)
872 :
873 : enljj(nu+1)=enljj(nu+1) &
874 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,3+mub,ilmn,ia,ispinor) &
875 : & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) &
876 : & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(2,mu,ilmn,ia,ispinor) &
877 1775808 : & -gxfac_sij(2,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
878 :
879 2367744 : mu=mu+1;nu=nu+2
880 : end do
881 : end do
882 : end do
883 : ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k
884 : else
885 0 : do ilmn=1,nlmn
886 : mu=1;nu=1
887 0 : do mua=1,3 ! atm. pos
888 0 : do mub=1,3 ! k
889 : enljj(nu+1)=enljj(nu+1) &
890 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,3+mub,ilmn,ia,ispinor) &
891 0 : & +gxfac_sij(1,ilmn,ia,ispinor)*d2gxdt(1,mu,ilmn,ia,ispinor)
892 0 : mu=mu+1;nu=nu+2
893 : end do
894 : end do
895 : end do
896 : end if
897 467680 : enlout(iashift+1:iashift+18)=enlout(iashift+1:iashift+18)+enljj(1:18)
898 : end do
899 : end do
900 22624 : ABI_FREE(enljj)
901 : end if
902 :
903 : ! ====== Accumulate contribution to <c|d2S/d_dstrain d_right_k|c> =========
904 264864 : if (choice==55) then
905 5152 : ABI_MALLOC(enljj,(36))
906 10304 : do ispinor=1,nspinor
907 16096 : do ia=1,nincat
908 289600 : enljj(1:36)=zero;enlj(:)=zero
909 : ! If cplex=1, dgxdt is real for strain, pure imaginary for k;
910 : ! If cplex_fac=1, dgxdtfac is pure imaginary for k;
911 5792 : if(cplex==2.and.cplex_fac==2) then
912 56288 : do ilmn=1,nlmn
913 : ! First compute 2nd-derivative contribution
914 : mu=1
915 353472 : do mua=1,6 ! strain (lambda,nu)
916 302976 : mua1=alpha(mua) ! (nu)
917 302976 : mua2=beta(mua) ! (lambda)
918 1262400 : do mub=1,3 ! k (mu)
919 908928 : muu=3*(gamma(mua1,mub)-1)+mua2
920 908928 : mut=3*(gamma(mua2,mub)-1)+mua1
921 : d2gx(1:cplex)=half*(d2gxdt(1:cplex,muu,ilmn,ia,ispinor) &
922 2726784 : & +d2gxdt(1:cplex,mut,ilmn,ia,ispinor))
923 : enljj(mu)=enljj(mu) &
924 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) &
925 : & +dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) &
926 908928 : & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1)+gxfac_sij(2,ilmn,ia,ispinor)*d2gx(2)
927 : enljj(mu+1)=enljj(mu+1) &
928 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(2,6+mub,ilmn,ia,ispinor) &
929 : & -dgxdt(2,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) &
930 908928 : & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(2)-gxfac_sij(2,ilmn,ia,ispinor)*d2gx(1)
931 1211904 : mu=mu+2
932 : end do
933 : end do
934 : ! Then store 1st-derivative contribution
935 : mu=1
936 207776 : do nu=1,3
937 : enlj(mu )=enlj(mu )+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor) &
938 151488 : & +gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor)
939 : enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(2,6+nu,ilmn,ia,ispinor) &
940 151488 : & -gxfac_sij(2,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
941 201984 : mu=mu+2
942 : end do
943 : end do
944 : ! If cplex=1, dgxdt, d2gxdt and dgxdtfac_sij are real for atm. pos, pure imaginary for k
945 : else
946 0 : do ilmn=1,nlmn
947 : mu=1
948 0 : do mua=1,6 ! strain (lambda,nu)
949 0 : mua1=alpha(mua) ! (nu)
950 0 : mua2=beta(mua) ! (lambda)
951 0 : do mub=1,3 ! k (mu)
952 0 : muu=3*(gamma(mua1,mub)-1)+mua2
953 0 : mut=3*(gamma(mua2,mub)-1)+mua1
954 0 : d2gx(1)=half*(d2gxdt(1,muu,ilmn,ia,ispinor)+d2gxdt(1,mut,ilmn,ia,ispinor))
955 : enljj(mu+1)=enljj(mu+1) &
956 : & +dgxdt(1,mua,ilmn,ia,ispinor)*dgxdtfac_sij(1,6+mub,ilmn,ia,ispinor) &
957 0 : & +gxfac_sij(1,ilmn,ia,ispinor)*d2gx(1)
958 0 : mu=mu+2
959 : end do
960 : end do
961 : ! Then store 1st-derivative contribution
962 : mu=1
963 0 : do nu=1,3
964 0 : enlj(mu+1)=enlj(mu+1)+gxfac_sij(1,ilmn,ia,ispinor)*dgxdt(1,6+nu,ilmn,ia,ispinor)
965 0 : mu=mu+2
966 : end do
967 : end do
968 : end if
969 214304 : enlout(1:36)=enlout(1:36)+enljj(1:36)
970 45696 : ddkk(1:6)=ddkk(1:6)+enlj(1:6)
971 : end do
972 : end do
973 5152 : ABI_FREE(enljj)
974 : end if
975 :
976 : ! ====== Accumulate contribution to <c|d2S/d_k d_k|c> =========
977 264864 : if (choice==8) then
978 0 : ABI_MALLOC(cft,(3,nlmn))
979 0 : ABI_MALLOC(cfu,(3,nlmn))
980 0 : do ispinor=1,nspinor
981 0 : do ia=1,nincat
982 0 : enlj(1:6)=zero
983 0 : do ilmn=1,nlmn
984 0 : do mu=1,6
985 0 : do iplex=1,cplex
986 0 : enlj(mu)=enlj(mu)+gxfac_sij(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)
987 : end do
988 : end do
989 : end do
990 : ! If cplex=1, dgxdt is pure imaginary, dgxdtfac_sij is pure imaginary;
991 0 : if(cplex==2)then
992 0 : cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
993 0 : cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor))
994 : else
995 0 : cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
996 0 : cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor))
997 : end if
998 0 : do ilmn=1,nlmn
999 0 : do mu=1,6
1000 0 : mua=alpha(mu);mub=beta(mu)
1001 0 : enlj(mu)=enlj(mu)+real(conjg(cfu(mub,ilmn))*cft(mua,ilmn))
1002 : end do
1003 : end do
1004 0 : enlout(1:6)=enlout(1:6)+two*enlj(1:6)
1005 : end do
1006 : end do
1007 0 : ABI_FREE(cft)
1008 0 : ABI_FREE(cfu)
1009 : end if
1010 :
1011 : ! ====== Accumulate contribution to <c|d/d_k[d(right)S/d_k]|c> =========
1012 : ! Full derivative wrt to k1, right derivative wrt to k2
1013 264864 : if (choice==81) then
1014 3456 : ABI_MALLOC(cft,(3,nlmn))
1015 3456 : ABI_MALLOC(cfu,(6,nlmn))
1016 1152 : ABI_MALLOC(enljj,(18))
1017 2304 : do ispinor=1,nspinor
1018 3456 : do ia=1,nincat
1019 21888 : enljj(1:18)=zero
1020 1152 : if(cplex_fac==2)then !If cplex_fac=1, gxfac is pure real
1021 10368 : cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),gxfac_sij(2,1:nlmn,ia,ispinor))
1022 : else
1023 0 : cft(1,1:nlmn)=cmplx(gxfac_sij(1,1:nlmn,ia,ispinor),zero)
1024 : end if
1025 1152 : if(cplex==2)then !If cplex=1, d2gxdt is pure real
1026 65664 : cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),d2gxdt(2,1:6,1:nlmn,ia,ispinor))
1027 : else
1028 0 : cfu(1:6,1:nlmn)=cmplx(d2gxdt(1,1:6,1:nlmn,ia,ispinor),zero)
1029 : end if
1030 10368 : do ilmn=1,nlmn
1031 38016 : do mu=1,3
1032 119808 : do nu=1,3
1033 82944 : muu=3*(mu-1)+nu ; mut=gamma(mu,nu)
1034 82944 : enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(1,ilmn))*cfu(mut,ilmn))
1035 110592 : enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(1,ilmn))*cfu(mut,ilmn))
1036 : end do
1037 : end do
1038 : end do
1039 1152 : if(cplex==2)then !If cplex=1, dgxdt is pure imaginary
1040 38016 : cft(1:3,1:nlmn)=cmplx(dgxdt(1,1:3,1:nlmn,ia,ispinor),dgxdt(2,1:3,1:nlmn,ia,ispinor))
1041 : else
1042 0 : cft(1:3,1:nlmn)=cmplx(zero,dgxdt(1,1:3,1:nlmn,ia,ispinor))
1043 : end if
1044 1152 : if(cplex_fac==2)then !If cplex_fac=1, dgxdtfac is pure imaginary
1045 38016 : cfu(1:3,1:nlmn)=cmplx(dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor),dgxdtfac_sij(2,1:3,1:nlmn,ia,ispinor))
1046 : else
1047 0 : cfu(1:3,1:nlmn)=cmplx(zero,dgxdtfac_sij(1,1:3,1:nlmn,ia,ispinor))
1048 : end if
1049 10368 : do ilmn=1,nlmn
1050 38016 : do mu=1,3
1051 119808 : do nu=1,3
1052 82944 : muu=3*(mu-1)+nu
1053 82944 : enljj(2*muu-1)=enljj(2*muu-1)+ real(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
1054 110592 : enljj(2*muu )=enljj(2*muu )+aimag(conjg(cft(mu,ilmn))*cfu(nu,ilmn))
1055 : end do
1056 : end do
1057 : end do
1058 23040 : enlout(1:18)=enlout(1:18)+enljj(1:18)
1059 : end do
1060 : end do
1061 1152 : ABI_FREE(cft)
1062 1152 : ABI_FREE(cfu)
1063 1152 : ABI_FREE(enljj)
1064 : end if
1065 :
1066 : end if
1067 :
1068 738502 : end subroutine opernld_ylm
1069 : !!***
1070 :
1071 : end module m_opernld_ylm
1072 : !!***
|
