Branch data 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 : 768799 : subroutine opernld_ylm(choice,cplex,cplex_fac,ddkk,dgxdt,dgxdtfac,dgxdtfac_sij,d2gxdt,&
121 : 768799 : & 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 : 1537598 : real(dp) :: d2gx(cplex),enlj(6*cplex),gxfacj(cplex)
151 : 768799 : real(dp),allocatable :: enljj(:)
152 : 768799 : complex(dp),allocatable :: cft(:,:), cfu(:,:)
153 : : ! *************************************************************************
154 : :
155 [ - + ]: 768799 : ABI_CHECK(cplex_fac>=cplex,'BUG: invalid cplex_fac<cplex!')
156 : :
157 [ + + ]: 768799 : if (paw_opt==0.or.paw_opt==1.or.paw_opt==2) then
158 : :
159 : : ! ============== Accumulate the non-local energy ===============
160 [ - + ]: 503935 : 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 [ + + - + ]: 503935 : if (choice==2.or.choice==23.or.choice==24) then
202 [ + + ]: 405298 : ishift=0;if (choice==23) ishift=6
203 [ + + ]: 870083 : do ispinor=1,nspinor
204 [ + + ]: 1996826 : do ia=1,nincat
205 [ + + ]: 4112424 : enlj(1:3)=zero
206 : 1028106 : iashift=3*(ia+ia3-2)+ishift
207 [ + + ]: 10864183 : do ilmn=1,nlmn
208 [ + + ]: 40372414 : do mu=1,3
209 : : dummy = zero ! Dummy needed here to get the correct forces with intel -O3
210 [ + + ]: 77331987 : do iplex=1,cplex
211 : : !enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
212 : 77331987 : dummy=dummy+gxfac(iplex,ilmn,ia,ispinor)*dgxdt(iplex,mu+ishift,ilmn,ia,ispinor)
213 : : end do
214 : 39344308 : enlj(mu)=enlj(mu)+dummy
215 : : end do
216 : : end do
217 [ + + ]: 4577209 : 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 [ + + ]: 503935 : if (choice==3.or.choice==23) then
224 [ + + ]: 802445 : enlj(1:6)=zero
225 [ + + ]: 243573 : do ispinor=1,nspinor
226 [ + + ]: 442602 : do ia=1,nincat
227 [ + + ]: 3138142 : do ilmn=1,nlmn
228 [ + + ]: 8356907 : gxfacj(1:cplex)=gxfac(1:cplex,ilmn,ia,ispinor)
229 [ + + ]: 8356907 : do iplex=1,cplex
230 : 8356907 : enlk=enlk+gxfacj(iplex)*gx(iplex,ilmn,ia,ispinor)
231 : : end do
232 [ + + ]: 19870254 : do mu=1,6
233 [ + + ]: 52951617 : do iplex=1,cplex
234 : 50141442 : 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 [ + + ]: 802445 : enlout(1:6)=enlout(1:6)+two*enlj(1:6)
241 : : end if
242 : :
243 : : ! ====== Accumulate the dynamical matrix contributions =========
244 [ + + ]: 503935 : if (choice==4.or.choice==24) then
245 [ - + ]: 43544 : ishift=0;if (choice==24) ishift=3*natom
246 [ + + ]: 87152 : do ispinor=1,nspinor
247 [ + + ]: 134760 : do ia=1,nincat
248 [ + + ]: 333256 : enlj(1:6)=zero
249 : 47608 : iashift=6*(ia+ia3-2)+ishift
250 [ + + ]: 536104 : do ilmn=1,nlmn
251 [ + + ]: 3467080 : do mu=1,6
252 : 2930976 : mua=alpha(mu);mub=beta(mu)
253 [ + + ]: 9281424 : do iplex=1,cplex
254 : : enlj(mu)=enlj(mu)+gxfac(iplex,ilmn,ia,ispinor)*d2gxdt(iplex,mu,ilmn,ia,ispinor)&
255 : 8792928 : & +dgxdtfac(iplex,mub,ilmn,ia,ispinor)*dgxdt(iplex,mua,ilmn,ia,ispinor)
256 : : end do
257 : : end do
258 : : end do
259 [ + + ]: 376864 : 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 [ + + ]: 503935 : 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 [ - + ]: 503935 : 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 [ - + ]: 503935 : 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 [ - + ]: 503935 : 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 [ - + ]: 503935 : 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 [ + + ]: 503935 : 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 [ + + ]: 503935 : 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 [ - + ]: 503935 : 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 [ + + ]: 768799 : 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 : 768799 : end subroutine opernld_ylm
1069 : : !!***
1070 : :
1071 : : end module m_opernld_ylm
1072 : : !!***
|
