Line data Source code
1 : !!****m* ABINIT/m_dfpt_fef
2 : !! NAME
3 : !! m_dfpt_fef
4 : !!
5 : !! FUNCTION
6 : !! Response calculations in finite electric field.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2004-2026 ABINIT group (XW).
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_dfpt_fef
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_efield
28 : use m_dtset
29 :
30 : use defs_abitypes, only : MPI_type
31 : use m_kg, only : kpgio
32 : use m_cgtools, only : overlap_g
33 : use m_xmpi
34 : use m_mpinfo
35 :
36 : implicit none
37 :
38 : private
39 : !!***
40 :
41 : public :: dfptff_initberry
42 : public :: dfptff_gradberry
43 : public :: dfptff_gbefd
44 : public :: dfptff_edie
45 : public :: dfptff_ebp
46 : public :: dfptff_die
47 : public :: dfptff_bec
48 : public :: qmatrix
49 : !!***
50 :
51 : contains
52 : !!***
53 :
54 : !!****f* ABINIT/dfptff_initberry
55 : !! NAME
56 : !! dfptff_initberry
57 : !!
58 : !! FUNCTION
59 : !! Initialization of response calculations in finite electric field.
60 : !!
61 : !! INPUTS
62 : !! dtset <type(dataset_type)> = all input variables in this dataset
63 : !! gmet(3,3) = reciprocal space metric tensor in bohr**-2
64 : !! kg(3,mpw*mkmem_rbz) = reduced (integer) coordinates of G vecs in basis sphere
65 : !! kg1(3,mpw1*mkmem_rbz) = reduced (integer) coordinates of G vecs for response wfs
66 : !! mband = maximum number of bands
67 : !! mkmem_rbz = maximum number of k-points in core memory
68 : !! mpw = maximum number of plane waves
69 : !! mpw1 = maximum number of plane waves for response wavefunctions
70 : !! nkpt = number of k points
71 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
72 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
73 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
74 : !! occ(mband*nkpt*nsppol) = occup number for each band at each k point
75 : !! rprimd(3,3) = dimensional primitive vectors
76 : !!
77 : !! OUTPUT
78 : !! dtefield=variables related to response Berry-phase calculation
79 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
80 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
81 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
82 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
83 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
84 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
85 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
86 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
87 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
88 : !!
89 : !! NOTES
90 : !! this duplicates in part initberry for the init of the dtefield - should be made
91 : !! into a common constructor in m_dtefield or somethin
92 : !!
93 : !! SOURCE
94 :
95 9 : subroutine dfptff_initberry(dtefield,dtset,gmet,kg,kg1,mband,mkmem_rbz,mpi_enreg,&
96 9 : & mpw,mpw1,nkpt,npwarr,npwar1,nsppol,occ,pwindall,rprimd)
97 :
98 : !Arguments ----------------------------------------
99 : !scalars
100 : integer,intent(in) :: mband,mkmem_rbz,mpw,mpw1,nkpt,nsppol
101 : type(MPI_type),intent(inout) :: mpi_enreg
102 : type(dataset_type),intent(in) :: dtset
103 : type(efield_type),intent(inout) :: dtefield !vz_i needs efield2
104 : !arrays
105 : integer,intent(in) :: kg(3,mpw*mkmem_rbz),kg1(3,mpw1*mkmem_rbz),npwar1(nkpt)
106 : integer,intent(in) :: npwarr(nkpt)
107 : integer,intent(out) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
108 : real(dp),intent(in) :: gmet(3,3),occ(mband*nkpt*nsppol),rprimd(3,3)
109 :
110 : !Local variables ----------------------------------
111 : !scalars
112 : integer :: flag,iband,icg,idir,ifor,ikg,ikg1,ikpt,ikpt1,ikpt2,ikstr
113 : integer :: index,ipw,isppol,istr,iunmark,jpw,nband_k,mband_occ_k,nkstr,npw_k
114 : integer :: npw_k1,orig
115 : real(dp) :: ecut_eff,occ_val
116 : character(len=500) :: message
117 : !arrays
118 : integer :: dg(3)
119 9 : integer,allocatable :: kg_tmp(:,:),kpt_mark(:),npwarr_tmp(:),npwtot(:)
120 : real(dp) :: diffk(3),dk(3)
121 9 : real(dp),allocatable :: kpt1(:,:)
122 :
123 : ! *************************************************************************
124 :
125 : !-----------------------------------------------------------
126 : !---------------- initilize dtefield //---------------------
127 : !-----------------------------------------------------------
128 :
129 : !Initialization of efield_type variables
130 36 : dtefield%efield_dot(:) = zero
131 117 : dtefield%dkvecs(:,:) = zero
132 9 : dtefield%maxnstr = 0 ; dtefield%maxnkstr = 0
133 63 : dtefield%nstr(:) = 0 ; dtefield%nkstr(:) = 0
134 36 : ABI_MALLOC(dtefield%ikpt_dk,(nkpt,9,3))
135 36 : ABI_MALLOC(dtefield%cgindex,(nkpt,nsppol*2))
136 27 : ABI_MALLOC(dtefield%kgindex,(nkpt))
137 2223 : dtefield%ikpt_dk(:,:,:) = 0
138 171 : dtefield%cgindex(:,:) = 0
139 9 : dtefield%mband_occ = 0
140 27 : ABI_MALLOC(dtefield%nband_occ,(nsppol))
141 18 : dtefield%nband_occ = 0
142 74556 : pwindall(:,:,:) = 0
143 :
144 : !Compute the number of occupied bands and check that--------
145 : !it is the same for each k-point----------------------------
146 :
147 9 : occ_val = two/(dtset%nsppol*one)
148 :
149 9 : index = 0
150 18 : do isppol = 1, nsppol
151 90 : do ikpt = 1, nkpt
152 :
153 72 : mband_occ_k = 0
154 72 : nband_k = dtset%nband(ikpt)
155 :
156 360 : do iband = 1, nband_k
157 288 : index = index + 1
158 360 : if (abs(occ(index) - occ_val) < tol8) mband_occ_k = mband_occ_k + 1
159 : end do
160 :
161 81 : if (ikpt > 1) then
162 63 : if (dtefield%nband_occ(isppol) /= mband_occ_k) then
163 0 : message = ' The number of valence bands is not the same for every k-point for present spin'
164 0 : ABI_ERROR(message)
165 : end if
166 : else
167 9 : dtefield%mband_occ = max(dtefield%mband_occ,mband_occ_k)
168 9 : dtefield%nband_occ(isppol) = mband_occ_k
169 : end if
170 :
171 : end do
172 : end do
173 :
174 : !DEBUG
175 : !write(std_out,*)'dfptff_initberry:nkpt',nkpt
176 : !ENDDEBUG
177 :
178 : !Compute the location of zero-order wavefunction --------------
179 : icg = 0
180 18 : do isppol = 1, nsppol
181 90 : do ikpt = 1, nkpt
182 :
183 72 : dtefield%cgindex(ikpt,isppol) = icg
184 72 : nband_k = dtset%nband(ikpt)
185 72 : npw_k = npwarr(ikpt)
186 : icg = icg + dtset%nspinor*npw_k*&
187 81 : & proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,1,mpi_enreg%me_band)
188 :
189 : end do
190 : end do
191 :
192 : !Compute the location of kg --------------
193 9 : ikg = 0
194 81 : do ikpt = 1, nkpt
195 :
196 72 : dtefield%kgindex(ikpt) = ikg
197 72 : npw_k = npwarr(ikpt)
198 81 : ikg = ikg + npw_k
199 :
200 : end do
201 :
202 : !Compute the location of first-order wavefunction -------------------
203 : icg = 0
204 18 : do isppol = 1, nsppol
205 90 : do ikpt = 1, nkpt
206 :
207 72 : dtefield%cgindex(ikpt,isppol+nsppol) = icg
208 72 : nband_k = dtset%nband(ikpt)
209 72 : npw_k = npwar1(ikpt)
210 : icg = icg + dtset%nspinor*npw_k* &
211 81 : & proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,1,mpi_enreg%me_band)
212 :
213 : end do
214 : end do
215 :
216 : !Compute the reciprocal lattice coordinates of the electric field-----
217 :
218 36 : dtefield%efield_dot(1) = dot_product(dtset%efield(:),rprimd(:,1))
219 36 : dtefield%efield_dot(2) = dot_product(dtset%efield(:),rprimd(:,2))
220 36 : dtefield%efield_dot(3) = dot_product(dtset%efield(:),rprimd(:,3))
221 :
222 9 : write(message,'(a,a,a,a,3(2x,f16.9),a)')ch10,&
223 9 : & ' initberry: Reciprocal lattice coordinates of the electric field',ch10,&
224 18 : & ' efield_dot(1:3) = ',dtefield%efield_dot(1:3),ch10
225 9 : call wrtout(std_out,message,'COLL')
226 :
227 : !find the related k points to every k point in full BZ-----------------
228 : !TODO: import hash tables for k-points, to make lookup fast and check for missing
229 : ! matches
230 :
231 : !loop over three reciprocal directions
232 36 : do idir = 1, 3
233 :
234 27 : if (dtset%rfdir(idir) == 1) then
235 :
236 : ! Compute dk(:), the vector between a k-point and its nearest
237 : ! neighbour along the direction idir
238 :
239 27 : dk(:) = zero
240 27 : dk(idir) = 1000_dp ! initialize with a big number
241 216 : do ikpt = 2, nkpt
242 756 : diffk(:) = abs(dtset%kptns(:,ikpt) - dtset%kptns(:,1))
243 189 : if ((diffk(1) < dk(1)+tol8).and.(diffk(2) < dk(2)+tol8).and.&
244 54 : & (diffk(3) < dk(3)+tol8)) dk(:) = diffk(:)
245 : end do
246 108 : dtefield%dkvecs(:,idir) = dk(:)
247 :
248 : ! For each k point, find k_prim such that k_prim= k + dk mod(G)
249 : ! where G is a vector of the reciprocal lattice
250 :
251 243 : do ikpt = 1, nkpt
252 :
253 : ! First: k + dk
254 972 : do ikpt1 = 1, nkpt
255 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:))
256 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
257 216 : dtefield%ikpt_dk(ikpt,1,idir) = ikpt1
258 216 : exit
259 : end if
260 : end do
261 :
262 : ! Second: k - dk
263 972 : do ikpt1 = 1, nkpt
264 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:))
265 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
266 216 : dtefield%ikpt_dk(ikpt,2,idir) = ikpt1
267 216 : exit
268 : end if
269 : end do
270 :
271 : ! new
272 : ! 3rd: k + (n+1)dk
273 :
274 972 : do ikpt1 = 1, nkpt
275 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) - dtset%qptn(:))
276 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
277 216 : dtefield%ikpt_dk(ikpt,3,idir) = ikpt1
278 216 : exit
279 : end if
280 : end do
281 :
282 : ! 6th: k - (n+1)dk
283 :
284 972 : do ikpt1 = 1, nkpt
285 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) + dtset%qptn(:))
286 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
287 216 : dtefield%ikpt_dk(ikpt,6,idir) = ikpt1
288 216 : exit
289 : end if
290 : end do
291 :
292 : ! 4th: k + (n-1)dk
293 :
294 972 : do ikpt1 = 1, nkpt
295 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dk(:) - dtset%qptn(:))
296 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
297 216 : dtefield%ikpt_dk(ikpt,4,idir) = ikpt1
298 216 : exit
299 : end if
300 : end do
301 :
302 : ! 5th: k - (n-1)dk
303 :
304 972 : do ikpt1 = 1, nkpt
305 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dk(:) + dtset%qptn(:))
306 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
307 216 : dtefield%ikpt_dk(ikpt,5,idir) = ikpt1
308 216 : exit
309 : end if
310 : end do
311 :
312 : ! 7th: k+n dk
313 :
314 972 : do ikpt1 = 1, nkpt
315 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) - dtset%qptn(:))
316 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
317 216 : dtefield%ikpt_dk(ikpt,7,idir) = ikpt1
318 216 : exit
319 : end if
320 : end do
321 :
322 : ! 8th: k-n dk
323 :
324 972 : do ikpt1 = 1, nkpt
325 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) - dtset%kptns(:,ikpt) + dtset%qptn(:))
326 3888 : if(sum(abs(diffk(:) - nint(diffk(:)))) < 3*tol8) then
327 216 : dtefield%ikpt_dk(ikpt,8,idir) = ikpt1
328 216 : exit
329 : end if
330 : end do
331 :
332 : ! 9th: -k
333 :
334 999 : do ikpt1 = 1, nkpt
335 3888 : diffk(:) = abs(dtset%kptns(:,ikpt1) + dtset%kptns(:,ikpt))
336 3888 : if(sum(diffk(:)) < 3*tol8) then
337 216 : dtefield%ikpt_dk(ikpt,9,idir) = ikpt1
338 216 : exit
339 : end if
340 : end do
341 :
342 : end do ! ikpt
343 :
344 :
345 :
346 : ! Find the string length, starting from k point 1
347 : ! (all strings must have the same number of points)
348 :
349 27 : nkstr = 1
350 27 : ikpt1 = 1
351 54 : do ikpt = 1, nkpt
352 54 : ikpt1 = dtefield%ikpt_dk(ikpt1,1,idir)
353 54 : if (ikpt1 == 1) exit
354 54 : nkstr = nkstr + 1
355 : end do
356 :
357 : ! Check that the string length is a divisor of nkpt
358 27 : if(mod(nkpt,nkstr) /= 0) then
359 0 : write(message,'(a,i0,a,i0)')' The string length = ',nkstr,', is not a divisor of nkpt =',nkpt
360 0 : ABI_BUG(message)
361 : end if
362 27 : dtefield%nkstr(idir) = nkstr
363 27 : dtefield%nstr(idir) = nkpt/nkstr
364 :
365 : end if ! dtset%rfdir(idir) == 1
366 :
367 : write(message,'(a,i1,a,i3,a,i3)')&
368 27 : & ' dfptff_initberry: for direction ',idir,', nkstr = ',dtefield%nkstr(idir),&
369 54 : & ', nstr = ',dtefield%nstr(idir)
370 36 : call wrtout(std_out,message,'COLL')
371 :
372 : end do ! close loop over idir
373 :
374 36 : dtefield%maxnstr = maxval(dtefield%nstr(:))
375 36 : dtefield%maxnkstr = maxval(dtefield%nkstr(:))
376 45 : ABI_MALLOC(dtefield%idxkstr,(dtefield%maxnkstr,dtefield%maxnstr,3))
377 360 : dtefield%idxkstr(:,:,:) = 0
378 :
379 :
380 : !Build the different strings------------------------------------------
381 :
382 18 : ABI_MALLOC(kpt_mark,(nkpt))
383 36 : do idir = 1, 3
384 :
385 36 : if (dtset%rfdir(idir) == 1) then
386 :
387 243 : iunmark = 1
388 243 : kpt_mark(:) = 0
389 135 : do istr = 1, dtefield%nstr(idir)
390 :
391 234 : do while(kpt_mark(iunmark) /= 0)
392 126 : iunmark = iunmark + 1
393 : end do
394 108 : dtefield%idxkstr(1,istr,idir) = iunmark
395 108 : kpt_mark(iunmark) = 1
396 243 : do ikstr = 2, dtefield%nkstr(idir)
397 108 : ikpt1 = dtefield%idxkstr(ikstr-1,istr,idir)
398 108 : ikpt2 = dtefield%ikpt_dk(ikpt1,1,idir)
399 108 : dtefield%idxkstr(ikstr,istr,idir) = ikpt2
400 216 : kpt_mark(ikpt2) = 1
401 : end do
402 :
403 : end do ! istr
404 :
405 : end if ! rfdir(idir) == 1
406 :
407 : end do ! close loop over idir
408 :
409 9 : ABI_FREE(kpt_mark)
410 :
411 :
412 : !Build the array pwindall that is needed to compute the different overlap matrices
413 : !at k +- dk
414 :
415 27 : ABI_MALLOC(kg_tmp,(3,max(mpw,mpw1)*mkmem_rbz))
416 27 : ABI_MALLOC(kpt1,(3,nkpt))
417 18 : ABI_MALLOC(npwarr_tmp,(nkpt))
418 18 : ABI_MALLOC(npwtot,(nkpt))
419 9 : ecut_eff = dtset%ecut*(dtset%dilatmx)**2
420 :
421 36 : do idir = 1, 3
422 :
423 36 : if (dtset%rfdir(idir) == 1) then
424 :
425 108 : dk(:) = dtefield%dkvecs(:,idir)
426 :
427 81 : do ifor = 1, 2
428 :
429 135 : if (ifor == 2) dk(:) = -1_dp*dk(:)
430 :
431 : ! Build the array kpt1 = kptns + dk
432 : ! all k-poins of kpt1 must be in the same BZ as those of kptns
433 1782 : kpt1(:,:) = zero
434 486 : do ikpt = 1, nkpt
435 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
436 1782 : kpt1(:,ikpt) = dtset%kptns(:,ikpt1)
437 : end do ! close loop over ikpt
438 :
439 : ! Set up the basis sphere of plane waves at kpt1
440 74358 : kg_tmp(:,:) = 0
441 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
442 : & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,&
443 54 : & npwarr_tmp,npwtot,dtset%nsppol)
444 :
445 54 : ikg = 0 ; ikg1 = 0
446 513 : do ikpt = 1, nkpt
447 :
448 432 : nband_k = dtset%nband(ikpt)
449 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
450 :
451 :
452 : dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
453 1728 : & dtset%kptns(:,ikpt1))
454 :
455 432 : flag = 0; orig = 1
456 432 : if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
457 :
458 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
459 432 : npw_k = npwarr(ikpt)
460 432 : npw_k1 = npwarr_tmp(ikpt)
461 :
462 18360 : do ipw = 1, npw_k
463 275328 : do jpw = orig, npw_k1
464 : if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
465 271440 : (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
466 3456 : (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
467 14472 : pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor,idir) = jpw
468 25164 : if (flag == 0) orig = jpw
469 : exit
470 : end if
471 : end do
472 : end do
473 :
474 432 : ikg = ikg + npw_k
475 486 : ikg1 = ikg1 + npw_k1
476 :
477 : end do ! close loop over ikpt
478 :
479 : end do ! close loop over ifor
480 :
481 : end if ! rfdir(idir) == 1
482 :
483 : end do ! close loop over idir
484 :
485 : !------------------------------------------------------------------------------
486 : !<u_q|u_q>
487 : !at k +- dk
488 :
489 36 : do idir = 1, 3
490 :
491 36 : if (dtset%rfdir(idir) == 1) then
492 :
493 108 : dk(:) = dtefield%dkvecs(:,idir)
494 :
495 81 : do ifor = 1, 2
496 :
497 135 : if (ifor == 2) dk(:) = -1_dp*dk(:)
498 :
499 : ! Build the array kpt1 = kptns + qptn + dk
500 : ! all k-poins of kpt1 must be in the same BZ as those of kptns
501 1782 : kpt1(:,:) = zero
502 486 : do ikpt = 1, nkpt
503 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
504 1782 : kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+dtset%qptn(:)
505 : end do ! close loop over ikpt
506 :
507 : ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1
508 74358 : kg_tmp(:,:) = 0
509 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
510 : & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,&
511 54 : & npwarr_tmp,npwtot,dtset%nsppol)
512 :
513 :
514 54 : ikg = 0 ; ikg1 = 0
515 513 : do ikpt = 1, nkpt
516 :
517 432 : nband_k = dtset%nband(ikpt)
518 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
519 :
520 432 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle
521 :
522 : dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
523 1728 : & dtset%kptns(:,ikpt1))
524 :
525 432 : flag = 0; orig = 1
526 432 : if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
527 :
528 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor,idir)
529 432 : npw_k = npwar1(ikpt)
530 432 : npw_k1 = npwarr_tmp(ikpt)
531 :
532 18360 : do ipw = 1, npw_k
533 275328 : do jpw = orig, npw_k1
534 : if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
535 271440 : (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
536 3456 : (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
537 14472 : pwindall((ikpt-1)*max(mpw,mpw1) +ipw,ifor+2,idir) = jpw
538 25164 : if (flag == 0) orig = jpw
539 : exit
540 : end if
541 : end do
542 : end do
543 :
544 432 : ikg = ikg + npw_k
545 486 : ikg1 = ikg1 + npw_k1
546 :
547 : end do ! close loop over ikpt
548 : end do ! close loop over ifor
549 : end if ! rfdir(idir) == 1
550 : end do ! close loop over idir
551 :
552 :
553 : !---------------------------------------------------------------------------
554 :
555 36 : do idir = 1, 3
556 :
557 36 : if (dtset%rfdir(idir) == 1) then
558 :
559 108 : dk(:) = dtset%qptn(:) + dtefield%dkvecs(:,idir)
560 :
561 81 : do ifor = 1, 2
562 :
563 135 : if (ifor == 2) dk(:) = dtset%qptn(:) - dtefield%dkvecs(:,idir)
564 :
565 : ! Build the array kpt1 = kptns + qptn + dk
566 : ! all k-poins of kpt1 must be in the same BZ as those of kptns
567 1782 : kpt1(:,:) = zero
568 486 : do ikpt = 1, nkpt
569 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir)
570 1782 : kpt1(:,ikpt) = dtset%kptns(:,ikpt1)
571 : end do ! close loop over ikpt
572 :
573 : ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1
574 74358 : kg_tmp(:,:) = 0
575 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
576 : & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw,&
577 54 : & npwarr_tmp,npwtot,dtset%nsppol)
578 :
579 54 : ikg = 0 ; ikg1 = 0
580 513 : do ikpt = 1, nkpt
581 :
582 432 : nband_k = dtset%nband(ikpt)
583 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir)
584 :
585 432 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle
586 :
587 : dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
588 1728 : & dtset%kptns(:,ikpt1))
589 :
590 432 : flag = 0; orig = 1
591 432 : if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
592 :
593 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor+2,idir)
594 432 : npw_k = npwar1(ikpt)
595 432 : npw_k1 = npwarr_tmp(ikpt)
596 :
597 18360 : do ipw = 1, npw_k
598 275328 : do jpw = orig, npw_k1
599 : if ((kg1(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
600 271440 : (kg1(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
601 3456 : (kg1(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
602 14472 : pwindall((ikpt-1)*max(mpw,mpw1)+ipw,ifor+4,idir) = jpw
603 25164 : if (flag == 0) orig = jpw
604 : exit
605 : end if
606 : end do
607 : end do
608 432 : ikg = ikg + npw_k
609 486 : ikg1 = ikg1 + npw_k1
610 :
611 : end do ! close loop over ikpt
612 : end do ! close loop over ifor
613 : end if ! rfdir(idir) == 1
614 : end do ! close loop over idir===============================================================
615 :
616 :
617 : !Build the array pwind3 that is needed to compute the overlap matrices===============================
618 :
619 36 : do idir = 1, 3
620 :
621 36 : if (dtset%rfdir(idir) == 1) then
622 :
623 108 : dk(:) = - dtset%qptn(:) + dtefield%dkvecs(:,idir)
624 :
625 81 : do ifor = 1, 2
626 :
627 135 : if (ifor == 2) dk(:) = - dtset%qptn(:) - dtefield%dkvecs(:,idir)
628 :
629 : ! Build the array kpt1 = kptns + qptn + dk
630 : ! all k-poins of kpt1 must be in the same BZ as those of kptns
631 1782 : kpt1(:,:) = zero
632 486 : do ikpt = 1, nkpt
633 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir)
634 1782 : kpt1(:,ikpt) = dtset%kptns(:,ikpt1)+ dtset%qptn(:)
635 : end do ! close loop over ikpt
636 :
637 : ! Set UP THE BASIS SPHERE OF PLANE waves at kpt1
638 74358 : kg_tmp(:,:) = 0
639 : call kpgio(ecut_eff,dtset%exchn2n3d,gmet,dtset%istwfk,kg_tmp,&
640 : & kpt1,mkmem_rbz,dtset%nband,nkpt,'PERS',mpi_enreg,mpw1,&
641 54 : & npwarr_tmp,npwtot,dtset%nsppol)
642 :
643 54 : ikg = 0 ; ikg1 = 0
644 513 : do ikpt = 1, nkpt
645 :
646 432 : nband_k = dtset%nband(ikpt)
647 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir)
648 :
649 432 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,1,mpi_enreg%me)) cycle
650 :
651 : dg(:) = nint(dtset%kptns(:,ikpt) + dk(:) + tol10 - &
652 1728 : & dtset%kptns(:,ikpt1))
653 :
654 432 : flag = 0; orig = 1
655 432 : if (dg(1)*dg(1) + dg(2)*dg(2) + dg(3)*dg(3) > 0) flag = 1
656 :
657 432 : ikpt1 = dtefield%ikpt_dk(ikpt,ifor+4,idir)
658 432 : npw_k = npwarr(ikpt)
659 432 : npw_k1 = npwarr_tmp(ikpt)
660 :
661 18360 : do ipw = 1, npw_k
662 275328 : do jpw = orig, npw_k1
663 : if ((kg(1,ikg + ipw) == kg_tmp(1,ikg1 + jpw) - dg(1)).and. &
664 271440 : (kg(2,ikg + ipw) == kg_tmp(2,ikg1 + jpw) - dg(2)).and. &
665 3456 : (kg(3,ikg + ipw) == kg_tmp(3,ikg1 + jpw) - dg(3))) then
666 14472 : pwindall((ikpt-1)*max(mpw,mpw1) + ipw,ifor+6,idir) = jpw
667 25164 : if (flag == 0) orig = jpw
668 : exit
669 : end if
670 : end do
671 : end do
672 432 : ikg = ikg + npw_k
673 486 : ikg1 = ikg1 + npw_k1
674 :
675 : end do ! close loop over ikpt
676 : end do ! close loop over ifor
677 : end if ! rfdir(idir) == 1
678 : end do ! close loop over idir====================================================================
679 :
680 9 : ABI_FREE(kg_tmp)
681 9 : ABI_FREE(kpt1)
682 9 : ABI_FREE(npwarr_tmp)
683 9 : ABI_FREE(npwtot)
684 :
685 9 : end subroutine dfptff_initberry
686 : !!***
687 :
688 : !!****f* ABINIT/dfptff_gradberry
689 : !! NAME
690 : !! dfptff_gradberry
691 : !!
692 : !! FUNCTION
693 : !! Calculation of the gradient of Berry-phase term in finite electric field.
694 : !!
695 : !! COPYRIGHT
696 : !! Copyright (C) 2004-2026 ABINIT group (XW).
697 : !! This file is distributed under the terms of the
698 : !! GNU General Public License, see ~abinit/COPYING
699 : !! or http://www.gnu.org/copyleft/gpl.txt .
700 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
701 : !!
702 : !! INPUTS
703 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
704 : !! cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
705 : !! RF wavefunctions at k,q.
706 : !! dtefield = variables related to finite electric field calculation
707 : !! ikpt = the index of the current k point
708 : !! isppol=1 for unpolarized, 2 for spin-polarized
709 : !! mband = maximum number of bands
710 : !! mkmem_rbz = maximum number of k-points in core memory
711 : !! mpw = maximum number of plane waves
712 : !! mpw1 = maximum number of plane waves for response wavefunctions
713 : !! nkpt = number of k points
714 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
715 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
716 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
717 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
718 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
719 : !! inverse of the overlap matrix
720 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
721 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
722 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
723 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
724 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
725 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
726 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
727 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
728 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
729 : !!
730 : !! OUTPUT
731 : !! grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
732 : !!
733 : !! SOURCE
734 :
735 488 : subroutine dfptff_gradberry(cg,cg1,dtefield,grad_berry,ikpt,isppol,&
736 : & mband,mband_mem,mpw,mpw1,mkmem_rbz,mk1mem,&
737 : & mpi_enreg,nkpt,&
738 488 : & npwarr,npwar1,nspinor,nsppol,qmat,pwindall)
739 :
740 : !Arguments ----------------------------------------
741 : !scalars
742 : integer,intent(in) :: ikpt,isppol,mband,mk1mem,mkmem_rbz,mpw,mpw1,nkpt,nspinor
743 : integer,intent(in) :: mband_mem
744 : integer,intent(in) :: nsppol
745 : type(efield_type),intent(in) :: dtefield
746 : type(MPI_type),intent(in) :: mpi_enreg
747 : !arrays
748 : integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
749 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
750 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol)
751 : real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)
752 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
753 : !TODO MJV: grad_berry should also be dimensioned with mband_mem
754 : real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ)
755 :
756 : !Local variables -------------------------
757 : !scalars
758 : integer :: iband,icg,icg1,idir,ikpt1
759 : integer :: iband_me, jband_me, me_band, ierr
760 : integer :: ikpt1m,ikptn,ikptnm,ikptnp1,ipw,jband,jpw,kband
761 : integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
762 : real(dp) :: doti,dotr,fac,wfi,wfr
763 : !arrays
764 976 : integer :: band_procs(mband)
765 488 : integer,allocatable :: pwind_tmp(:)
766 : real(dp) :: z1(2),z2(2)
767 488 : real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:)
768 : real(dp),allocatable :: vect2(:,:)
769 :
770 : ! *************************************************************************
771 :
772 488 : mpw_tmp=max(mpw,mpw1)
773 1464 : ABI_MALLOC(vect1,(2,0:mpw_tmp))
774 976 : ABI_MALLOC(vect2,(2,0:mpw_tmp))
775 1952 : ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
776 1464 : ABI_MALLOC(pwind_tmp,(mpw_tmp))
777 1464 : ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
778 1464 : ABI_MALLOC(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ))
779 2928 : vect1(:,0) = zero ; vect2(:,0) = zero
780 25864 : s1mat(:,:,:)=zero
781 254248 : grad_berry(:,:,:) = zero
782 :
783 488 : me_band = mpi_enreg%me_band
784 488 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,mband,me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
785 :
786 1952 : do idir=1,3
787 : fac = dtefield%efield_dot(idir)*dble(nkpt)/&
788 1464 : & (dble(dtefield%nstr(idir))*four_pi)
789 :
790 : ! prepare
791 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
792 1464 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
793 1464 : npw_k1 = npwar1(ikpt)
794 1464 : npw_k2 = npwar1(ikpt1)
795 62220 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
796 :
797 : !TODO: this algorithm is very inefficient: the looping is in the wrong order.
798 : ! should be possible to store a temp vector and make it into a BLAS call...
799 : ! basically boils down to an internal sum over iband
800 : ! grad_berry(ipw,jband) = fac * wf(ipw,:) * qmat(:,jband)
801 62220 : do ipw = 1, npw_k1
802 60756 : jpw = pwind_tmp(ipw)
803 :
804 62220 : if (jpw > 0) then
805 : iband_me = 0
806 245220 : do iband = 1, dtefield%mband_occ
807 196176 : if (band_procs(iband) /= me_band) cycle
808 196176 : iband_me = iband_me + 1
809 196176 : wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
810 196176 : wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
811 :
812 1029924 : do jband = 1, dtefield%mband_occ
813 :
814 : grad_berry(1,ipw,jband) = &
815 : & grad_berry(1,ipw,jband) + &
816 784704 : & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi
817 :
818 : grad_berry(2,ipw,jband) = &
819 : & grad_berry(2,ipw,jband) + &
820 980880 : & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr
821 :
822 : end do
823 : end do
824 : end if
825 :
826 : end do
827 :
828 : ! compute <u^(0)_{k_j+n}|u^(1)_{k_j+1,q}> matrix----------------------------------------------------
829 :
830 : ! prepare to calculate overlap matrix
831 1464 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
832 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
833 1464 : icg = dtefield%cgindex(ikptn,isppol)
834 1464 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
835 1464 : npw_k1 = npwarr(ikptn)
836 : npw_k2 = npwar1(ikpt1)
837 62220 : pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir)
838 :
839 :
840 7320 : vect1(:,0) = zero ; vect2(:,0) = zero
841 7320 : do jband = 1, dtefield%mband_occ
842 : vect2(:,1:npw_k2) = &
843 734928 : & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
844 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
845 :
846 30744 : do iband = 1, dtefield%mband_occ
847 :
848 23424 : pwmin = (iband-1)*npw_k1*nspinor
849 23424 : pwmax = pwmin + npw_k1*nspinor
850 : vect1(:,1:npw_k1) = &
851 2939712 : & cg(:,icg + 1 + pwmin:icg + pwmax)
852 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
853 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
854 23424 : & vect1,vect2)
855 23424 : s1mat(1,iband,jband) = dotr
856 29280 : s1mat(2,iband,jband) = doti
857 :
858 : end do ! iband
859 : end do !jband
860 :
861 : ! compute <u^(0)_{-k_j+1}|u^(1)_{-k_j+n},q> matrix--------------------
862 :
863 : ! prepare to calculate overlap matrix
864 1464 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
865 1464 : ikptnm= dtefield%ikpt_dk(ikptn,9,idir)
866 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
867 1464 : ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir)
868 :
869 1464 : icg = dtefield%cgindex(ikpt1m,isppol)
870 1464 : icg1 = dtefield%cgindex(ikptnm,isppol+nsppol)
871 1464 : npw_k1 = npwarr(ikpt1m)
872 1464 : npw_k2 = npwar1(ikptnm)
873 62220 : pwind_tmp(1:npw_k1) = pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,7,idir)
874 :
875 8784 : vect1(:,0) = zero ; vect2(:,0) = zero
876 7320 : do jband = 1, dtefield%mband_occ
877 : vect2(:,1:npw_k2) = &
878 734928 : & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
879 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
880 :
881 30744 : do iband = 1, dtefield%mband_occ
882 :
883 23424 : pwmin = (iband-1)*npw_k1*nspinor
884 23424 : pwmax = pwmin + npw_k1*nspinor
885 : vect1(:,1:npw_k1) = &
886 2939712 : & cg(:,icg + 1 + pwmin:icg + pwmax)
887 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
888 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
889 23424 : & vect1,vect2)
890 :
891 23424 : s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
892 29280 : s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
893 :
894 : end do ! iband
895 : end do !jband
896 :
897 77592 : Amat(:,:,:)=zero
898 :
899 : ! calculate Amat
900 7320 : do iband=1, dtefield%mband_occ
901 30744 : do jband=1, dtefield%mband_occ
902 122976 : do kband=1, dtefield%mband_occ
903 : Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*&
904 : & qmat(1,kband,jband,ikpt,1,idir)&
905 93696 : & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)
906 : Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*&
907 : & qmat(2,kband,jband,ikpt,1,idir)&
908 117120 : & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)
909 : end do
910 : end do
911 : end do
912 :
913 77592 : Bmat(:,:,:)=zero
914 :
915 : ! calculate Bmat
916 1464 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
917 7320 : do iband=1, dtefield%mband_occ
918 30744 : do jband=1, dtefield%mband_occ
919 122976 : do kband=1, dtefield%mband_occ
920 : Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*&
921 : & qmat(1,jband,iband,ikptn,1,idir)&
922 93696 : & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)
923 : Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*&
924 : & qmat(2,jband,iband,ikptn,1,idir)&
925 117120 : & + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)
926 : end do
927 : end do
928 : end do
929 :
930 : ! calc. the second term of gradient------------------------------
931 :
932 : ! preparation
933 :
934 1464 : ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir)
935 1464 : icg = dtefield%cgindex(ikptnp1,isppol)
936 1464 : npw_k1 = npwar1(ikpt)
937 1464 : npw_k2 = npwarr(ikptnp1)
938 62220 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
939 :
940 : z1(:) = zero
941 : z2(:) = zero
942 :
943 62220 : do ipw = 1, npw_k1
944 :
945 60756 : jpw = pwind_tmp(ipw)
946 :
947 62220 : if (jpw > 0) then
948 :
949 : iband_me = 0
950 245220 : do iband = 1, dtefield%mband_occ
951 196176 : if (band_procs(iband) /= me_band) cycle
952 196176 : iband_me = iband_me + 1
953 196176 : wfr = cg(1,icg + (iband_me - 1)*npw_k2*nspinor + jpw)
954 196176 : wfi = cg(2,icg + (iband_me - 1)*npw_k2*nspinor + jpw)
955 :
956 1029924 : do jband=1, dtefield%mband_occ
957 :
958 : grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) &
959 784704 : & - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
960 : grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) &
961 980880 : & - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
962 :
963 : end do
964 : end do
965 : end if
966 : end do
967 :
968 : ! Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
969 :
970 8784 : vect1(:,0) = zero ; vect2(:,0) = zero
971 :
972 : ! prepare
973 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
974 1464 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
975 : npw_k1 = npwar1(ikpt)
976 1464 : npw_k2 = npwar1(ikpt1)
977 62220 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir)
978 :
979 : ! write(std_out,*)'dfpt_cgwf:pwind_tmp',pwind_tmp
980 : ! stop
981 :
982 62220 : do ipw = 1, npw_k1
983 :
984 60756 : jpw = pwind_tmp(ipw)
985 :
986 62220 : if (jpw > 0) then
987 : iband_me = 0
988 245220 : do iband = 1, dtefield%mband_occ
989 196176 : if (band_procs(iband) /= me_band) cycle
990 196176 : iband_me = iband_me + 1
991 196176 : wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
992 196176 : wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
993 :
994 1029924 : do jband = 1, dtefield%mband_occ
995 :
996 : grad_berry(1,ipw,jband) = &
997 : & grad_berry(1,ipw,jband) - &
998 784704 : & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi
999 :
1000 : grad_berry(2,ipw,jband) = &
1001 : & grad_berry(2,ipw,jband) - &
1002 980880 : & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr
1003 :
1004 : end do
1005 : end do
1006 : end if
1007 : end do
1008 :
1009 :
1010 : ! compute <u^(0)_{k_j+n}|u^(1)_{k_j-1,q}> matrix----------------------------------------------------
1011 :
1012 : ! prepare to calculate overlap matrix
1013 1464 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1014 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1015 1464 : icg = dtefield%cgindex(ikptn,isppol)
1016 1464 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1017 1464 : npw_k1 = npwarr(ikptn)
1018 : npw_k2 = npwar1(ikpt1)
1019 62220 : pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir)
1020 :
1021 8784 : vect1(:,0) = zero ; vect2(:,0) = zero
1022 : jband_me = 0
1023 7320 : do jband = 1, dtefield%mband_occ
1024 5856 : if (band_procs(jband) == me_band) then
1025 5856 : jband_me = jband_me + 1
1026 : vect2(:,1:npw_k2) = &
1027 734928 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1028 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1029 : end if
1030 5856 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1031 :
1032 5856 : iband_me = 0
1033 30744 : do iband = 1, dtefield%mband_occ
1034 23424 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1035 23424 : iband_me = iband_me + 1
1036 23424 : pwmin = (iband-1)*npw_k1*nspinor
1037 23424 : pwmax = pwmin + npw_k1*nspinor
1038 : vect1(:,1:npw_k1) = &
1039 2939712 : & cg(:,icg + 1 + pwmin:icg + pwmax)
1040 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1041 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1042 23424 : & vect1,vect2)
1043 23424 : s1mat(1,iband,jband) = dotr
1044 29280 : s1mat(2,iband,jband) = doti
1045 :
1046 : end do ! iband
1047 : end do !jband
1048 :
1049 : ! compute <u^(0)_{-k_j-1}|u^(1)_{-k_j+n,q}> matrix-----------------------------------------------------
1050 :
1051 : ! prepare to calculate overlap matrix
1052 1464 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1053 1464 : ikptnm= dtefield%ikpt_dk(ikptn,9,idir)
1054 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1055 1464 : ikpt1m= dtefield%ikpt_dk(ikpt1,9,idir)
1056 1464 : icg = dtefield%cgindex(ikpt1m,isppol)
1057 1464 : icg1 = dtefield%cgindex(ikptnm,isppol+nsppol)
1058 1464 : npw_k1 = npwarr(ikpt1m)
1059 1464 : npw_k2 = npwar1(ikptnm)
1060 62220 : pwind_tmp(1:npw_k1) =pwindall((ikpt1m-1)*mpw_tmp+1:(ikpt1m-1)*mpw_tmp+npw_k1,8,idir)
1061 :
1062 :
1063 :
1064 7320 : vect1(:,0) = zero ; vect2(:,0) = zero
1065 : jband_me = 0
1066 7320 : do jband = 1, dtefield%mband_occ
1067 5856 : if (band_procs(jband) == me_band) then
1068 5856 : jband_me = jband_me + 1
1069 : vect2(:,1:npw_k2) = &
1070 734928 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1071 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1072 : end if
1073 5856 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1074 :
1075 5856 : iband_me = 0
1076 30744 : do iband = 1, dtefield%mband_occ
1077 23424 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1078 23424 : iband_me = iband_me + 1
1079 23424 : pwmin = (iband_me-1)*npw_k1*nspinor
1080 23424 : pwmax = pwmin + npw_k1*nspinor
1081 : vect1(:,1:npw_k1) = &
1082 2939712 : & cg(:,icg + 1 + pwmin:icg + pwmax)
1083 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1084 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1085 23424 : & vect1,vect2)
1086 23424 : s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
1087 29280 : s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
1088 :
1089 : end do ! iband
1090 : end do !jband
1091 : ! accumulate all iband values on each proc
1092 1464 : call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr)
1093 :
1094 77592 : Amat(:,:,:)=zero
1095 :
1096 : ! calculate Amat: s1mat is complete so sum over all bands here
1097 : ! TODO: could reduce one loop over bands
1098 : ! TODO: replace this with a BLAS call: A = s1*q in complex numbers
1099 7320 : do iband=1, dtefield%mband_occ
1100 30744 : do jband=1, dtefield%mband_occ
1101 122976 : do kband=1, dtefield%mband_occ
1102 : Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*&
1103 : & qmat(1,kband,jband,ikpt,2,idir)&
1104 93696 : & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)
1105 : Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*&
1106 : & qmat(2,kband,jband,ikpt,2,idir)&
1107 117120 : & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)
1108 : end do
1109 : end do
1110 : end do
1111 :
1112 77592 : Bmat(:,:,:)=zero
1113 :
1114 : ! calculate Bmat: as above, all matrices A B q are band-complete at this stage
1115 : ! TODO: could reduce one loop over bands
1116 : ! TODO: replace this with a BLAS call: B = q*A in complex numbers
1117 1464 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1118 7320 : do iband=1, dtefield%mband_occ
1119 30744 : do jband=1, dtefield%mband_occ
1120 122976 : do kband=1, dtefield%mband_occ
1121 : Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*&
1122 : & qmat(1,jband,iband,ikptn,2,idir)&
1123 93696 : & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)
1124 : Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*&
1125 : & qmat(2,jband,iband,ikptn,2,idir)&
1126 117120 : + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)
1127 : end do
1128 : end do
1129 : end do
1130 :
1131 : ! calc. the second term of gradient------------------------------
1132 :
1133 : ! preparation
1134 :
1135 1464 : ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir)
1136 1464 : icg = dtefield%cgindex(ikptnp1,isppol)
1137 1464 : npw_k1 = npwar1(ikpt)
1138 1464 : npw_k2 = npwarr(ikptnp1)
1139 62220 : pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1140 :
1141 : z1(:) = zero
1142 : z2(:) = zero
1143 64172 : do ipw = 1, npw_k1
1144 60756 : jpw = pwind_tmp(ipw)
1145 62220 : if (jpw > 0) then
1146 : iband_me = 0
1147 245220 : do iband = 1, dtefield%mband_occ
1148 196176 : if (band_procs(iband) /= me_band) cycle
1149 196176 : iband_me = iband_me + 1
1150 196176 : wfr = cg(1,icg + (iband_me - 1)*npw_k2*nspinor + jpw)
1151 196176 : wfi = cg(2,icg + (iband_me - 1)*npw_k2*nspinor + jpw)
1152 :
1153 1029924 : do jband=1, dtefield%mband_occ
1154 784704 : grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
1155 980880 : grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
1156 : end do
1157 : end do
1158 : end if
1159 : end do
1160 :
1161 : end do !idir
1162 :
1163 488 : call xmpi_sum(grad_berry,mpi_enreg%comm_band,ierr) ! sum over iband for all previous loops
1164 :
1165 488 : ABI_FREE(vect1)
1166 488 : ABI_FREE(vect2)
1167 488 : ABI_FREE(s1mat)
1168 488 : ABI_FREE(Amat)
1169 488 : ABI_FREE(Bmat)
1170 488 : ABI_FREE(pwind_tmp)
1171 :
1172 488 : end subroutine dfptff_gradberry
1173 : !!***
1174 :
1175 : !!****f* ABINIT/dfptff_gbefd
1176 : !! NAME
1177 : !! dfptff_gbefd
1178 : !!
1179 : !! FUNCTION
1180 : !! calculate the gradient of the second order \Omega E \cdot P
1181 : !! term, Eq.(23) in PRB 75, 115116(2007) [[cite:Wang2007]].
1182 : !!
1183 : !! INPUTS
1184 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
1185 : !! cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
1186 : !! RF wavefunctions at k,q.
1187 : !! dtefield = variables related to response Berry-phase calculation
1188 : !! ikpt = the index of the current k point
1189 : !! isppol = the index of the spin component
1190 : !! mband = maximum number of bands
1191 : !! mband_mem = maximum number of bands on this cpu
1192 : !! mkmem_rbz = maximum number of k-points in core memory
1193 : !! mpi_enreg = parallel distribution datastructure
1194 : !! mpw = maximum number of plane waves
1195 : !! mpw1 = maximum number of plane waves for response wavefunctions
1196 : !! nkpt = number of k points
1197 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
1198 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
1199 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
1200 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
1201 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
1202 : !! inverse of the overlap matrix
1203 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
1204 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
1205 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
1206 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
1207 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
1208 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
1209 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
1210 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
1211 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
1212 : !!
1213 : !! OUTPUT
1214 : !! grad_berry = the gradient of the Berry phase term
1215 : !!
1216 : !! SOURCE
1217 :
1218 288 : subroutine dfptff_gbefd(cg,cg1,dtefield,grad_berry,idir_efield,ikpt,isppol,&
1219 : & mband,mband_mem,mpw,mpw1,mkmem_rbz,mk1mem,&
1220 288 : & mpi_enreg,nkpt,npwarr,npwar1,nspinor,nsppol,qmat,pwindall,rprimd)
1221 :
1222 : !Arguments ----------------------------------------
1223 : !scalars
1224 : integer,intent(in) :: idir_efield,ikpt,isppol,mband,mk1mem,mkmem_rbz,mpw,mpw1,nkpt
1225 : integer,intent(in) :: mband_mem
1226 : integer,intent(in) :: nspinor,nsppol
1227 : type(efield_type),intent(in) :: dtefield
1228 : type(MPI_type),intent(in) :: mpi_enreg
1229 : !arrays
1230 : integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
1231 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
1232 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol)
1233 : real(dp),intent(in) :: cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol)
1234 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
1235 : real(dp),intent(in) :: rprimd(3,3)
1236 : !TODO MJV: this array is still npw*nband, and should be parallelized as well at some point...
1237 : real(dp),intent(out) :: grad_berry(2,mpw1,dtefield%mband_occ)
1238 :
1239 : !Local variables -------------------------
1240 : !scalars
1241 : integer :: iband,icg,icg1,idir,ikpt1
1242 : integer :: iband_me, jband_me, ierr
1243 : integer :: ikptn,ikptnp1,ipw,jband,jpw,kband
1244 : integer :: mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
1245 : real(dp) :: doti,dotr,fac,wfi,wfr
1246 : !arrays
1247 288 : integer,allocatable :: pwind_tmp(:)
1248 576 : integer :: band_procs(mband)
1249 : real(dp) :: z1(2),z2(2)
1250 288 : real(dp),allocatable :: Amat(:,:,:),Bmat(:,:,:),s1mat(:,:,:),vect1(:,:)
1251 : real(dp),allocatable :: vect2(:,:)
1252 :
1253 : ! *************************************************************************
1254 :
1255 288 : mpw_tmp=max(mpw,mpw1)
1256 864 : ABI_MALLOC(vect1,(2,0:mpw_tmp))
1257 576 : ABI_MALLOC(vect2,(2,0:mpw_tmp))
1258 1152 : ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
1259 864 : ABI_MALLOC(pwind_tmp,(mpw_tmp))
1260 864 : ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
1261 864 : ABI_MALLOC(Bmat,(2,dtefield%mband_occ,dtefield%mband_occ))
1262 1440 : vect1(:,0) = zero ; vect2(:,0) = zero
1263 15264 : s1mat(:,:,:)=zero
1264 150048 : grad_berry(:,:,:) = zero
1265 :
1266 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,mband,&
1267 288 : & mpi_enreg%me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
1268 :
1269 1152 : do idir=1,3
1270 : fac = dtefield%efield_dot(idir)*dble(nkpt)/&
1271 864 : & (dble(dtefield%nstr(idir))*four_pi)
1272 :
1273 : ! prepare
1274 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1275 864 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1276 864 : npw_k1 = npwar1(ikpt)
1277 864 : npw_k2 = npwar1(ikpt1)
1278 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
1279 :
1280 36720 : do ipw = 1, npw_k1
1281 35856 : jpw = pwind_tmp(ipw)
1282 36720 : if (jpw > 0) then
1283 28944 : iband_me = 0
1284 144720 : do iband = 1, dtefield%mband_occ
1285 115776 : if (proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,iband,iband,isppol,mpi_enreg%me_band)) cycle
1286 115776 : iband_me = iband_me + 1
1287 115776 : wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
1288 115776 : wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
1289 :
1290 607824 : do jband = 1, dtefield%mband_occ
1291 : grad_berry(1,ipw,jband) = &
1292 : & grad_berry(1,ipw,jband) + &
1293 463104 : & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi
1294 :
1295 : grad_berry(2,ipw,jband) = &
1296 : & grad_berry(2,ipw,jband) + &
1297 578880 : & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr
1298 : end do
1299 : end do
1300 : end if
1301 : end do
1302 : ! accumulate full sum over iband below at the end
1303 :
1304 : ! compute <u^(0)_{k_j}|u^(1)_{k_j+1}> matrix----------------------------------------------------
1305 :
1306 : ! prepare to calculate overlap matrix
1307 864 : ikptn = ikpt
1308 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1309 864 : icg = dtefield%cgindex(ikptn,isppol)
1310 864 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1311 864 : npw_k1 = npwarr(ikptn)
1312 : npw_k2 = npwar1(ikpt1)
1313 36720 : pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,7,idir)
1314 :
1315 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1316 : jband_me = 0
1317 4320 : do jband = 1, dtefield%mband_occ
1318 3456 : if (band_procs(jband) == mpi_enreg%me_band) then
1319 3456 : jband_me = jband_me + 1
1320 : vect2(:,1:npw_k2) = &
1321 433728 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1322 : end if
1323 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1324 :
1325 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1326 3456 : iband_me = 0
1327 21600 : do iband = 1, dtefield%mband_occ
1328 13824 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1329 13824 : iband_me = iband_me + 1
1330 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1331 13824 : pwmax = pwmin + npw_k1*nspinor
1332 1734912 : vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax)
1333 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1334 13824 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,vect1,vect2)
1335 13824 : s1mat(1,iband,jband) = dotr
1336 17280 : s1mat(2,iband,jband) = doti
1337 : end do ! iband
1338 : end do !jband
1339 :
1340 : ! compute <u^(1)_{k_j}|u^(0)_{k_j+1}> matrix-----------------------------------------------------
1341 :
1342 : ! prepare to calculate overlap matrix
1343 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1344 864 : icg = dtefield%cgindex(ikpt,isppol+nsppol)
1345 864 : icg1 = dtefield%cgindex(ikpt1,isppol)
1346 864 : npw_k1 = npwar1(ikpt)
1347 864 : npw_k2 = npwarr(ikpt1)
1348 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
1349 :
1350 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1351 : jband_me = 0
1352 4320 : do jband = 1, dtefield%mband_occ
1353 3456 : if (band_procs(jband) == mpi_enreg%me_band) then
1354 3456 : jband_me = jband_me + 1
1355 : vect2(:,1:npw_k2) = &
1356 433728 : & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1357 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1358 : end if
1359 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1360 :
1361 3456 : iband_me = 0
1362 18144 : do iband = 1, dtefield%mband_occ
1363 13824 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1364 13824 : iband_me = iband_me + 1
1365 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1366 13824 : pwmax = pwmin + npw_k1*nspinor
1367 : vect1(:,1:npw_k1) = &
1368 1734912 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
1369 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1370 13824 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,vect1,vect2)
1371 13824 : s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
1372 17280 : s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
1373 : end do ! iband
1374 : end do !jband
1375 : ! accumulate all iband values on each proc
1376 864 : call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr)
1377 :
1378 45792 : Amat(:,:,:)=zero
1379 :
1380 : ! calculate Amat
1381 4320 : do iband=1, dtefield%mband_occ
1382 18144 : do jband=1, dtefield%mband_occ
1383 72576 : do kband=1, dtefield%mband_occ
1384 : Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)&
1385 55296 : & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)
1386 : Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*qmat(2,kband,jband,ikpt,1,idir)&
1387 69120 : & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,1,idir)
1388 : end do
1389 : end do
1390 : end do
1391 :
1392 45792 : Bmat(:,:,:)=zero
1393 :
1394 : ! calculate Bmat
1395 864 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1396 4320 : do iband=1, dtefield%mband_occ
1397 18144 : do jband=1, dtefield%mband_occ
1398 72576 : do kband=1, dtefield%mband_occ
1399 : Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)&
1400 55296 : & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)
1401 : Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*qmat(2,jband,iband,ikptn,1,idir)&
1402 69120 : & + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,1,idir)
1403 : end do
1404 : end do
1405 : end do
1406 :
1407 : ! calc. the second term of gradient------------------------------
1408 :
1409 : ! preparation
1410 :
1411 864 : ikptnp1 = dtefield%ikpt_dk(ikpt,3,idir)
1412 864 : icg = dtefield%cgindex(ikptnp1,isppol)
1413 : npw_k1 = npwar1(ikpt)
1414 864 : npw_k2 = npwarr(ikptnp1)
1415 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
1416 :
1417 : z1(:) = zero
1418 : z2(:) = zero
1419 :
1420 36720 : do ipw = 1, npw_k1
1421 35856 : jpw = pwind_tmp(ipw)
1422 36720 : if (jpw > 0) then
1423 28944 : iband_me = 0
1424 144720 : do iband = 1, dtefield%mband_occ
1425 115776 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1426 115776 : iband_me = iband_me + 1
1427 115776 : wfr = cg(1,icg + (iband_me - 1)*npw_k2*nspinor + jpw)
1428 115776 : wfi = cg(2,icg + (iband_me - 1)*npw_k2*nspinor + jpw)
1429 :
1430 607824 : do jband=1, dtefield%mband_occ
1431 463104 : grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) - fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
1432 578880 : grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) - fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
1433 : end do
1434 : end do
1435 : end if
1436 : end do
1437 :
1438 : ! Second part of gradient of Berry phase++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1439 :
1440 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1441 :
1442 : ! prepare
1443 864 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1444 864 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1445 : npw_k1 = npwar1(ikpt)
1446 864 : npw_k2 = npwar1(ikpt1)
1447 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,4,idir)
1448 :
1449 36720 : do ipw = 1, npw_k1
1450 35856 : jpw = pwind_tmp(ipw)
1451 36720 : if (jpw > 0) then
1452 28944 : iband_me = 0
1453 144720 : do iband = 1, dtefield%mband_occ
1454 115776 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1455 115776 : iband_me = iband_me + 1
1456 115776 : wfr = cg1(1,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
1457 115776 : wfi = cg1(2,icg1 + (iband_me - 1)*npw_k2*nspinor + jpw)
1458 :
1459 607824 : do jband = 1, dtefield%mband_occ
1460 : grad_berry(1,ipw,jband) = &
1461 : & grad_berry(1,ipw,jband) - &
1462 463104 : & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi
1463 :
1464 : grad_berry(2,ipw,jband) = &
1465 : & grad_berry(2,ipw,jband) - &
1466 578880 : & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr
1467 : end do
1468 : end do
1469 : end if
1470 : end do
1471 :
1472 : ! compute <u^(0)_{k_j}|u^(1)_{k_j-1}> matrix----------------------------------------------------
1473 :
1474 : ! prepare to calculate overlap matrix
1475 864 : ikptn = ikpt
1476 864 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1477 864 : icg = dtefield%cgindex(ikptn,isppol)
1478 864 : icg1 = dtefield%cgindex(ikpt1,isppol+nsppol)
1479 864 : npw_k1 = npwarr(ikptn)
1480 : npw_k2 = npwar1(ikpt1)
1481 36720 : pwind_tmp(1:npw_k1) =pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,8,idir)
1482 :
1483 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1484 : jband_me = 0
1485 4320 : do jband = 1, dtefield%mband_occ
1486 3456 : if (band_procs(jband) == mpi_enreg%me_band) then
1487 3456 : jband_me = jband_me + 1
1488 : vect2(:,1:npw_k2) = &
1489 433728 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1490 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1491 : end if
1492 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1493 :
1494 :
1495 3456 : iband_me = 0
1496 18144 : do iband = 1, dtefield%mband_occ
1497 13824 : if (band_procs(iband) /= mpi_enreg%me_band) cycle
1498 13824 : iband_me = iband_me + 1
1499 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1500 13824 : pwmax = pwmin + npw_k1*nspinor
1501 : vect1(:,1:npw_k1) = &
1502 1734912 : & cg(:,icg + 1 + pwmin:icg + pwmax)
1503 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1504 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1505 13824 : & vect1,vect2)
1506 13824 : s1mat(1,iband,jband) = dotr
1507 17280 : s1mat(2,iband,jband) = doti
1508 : end do ! iband
1509 : end do !jband
1510 :
1511 : ! compute <u^(1)_{k_j}|u^(0)_{k_j-1}> matrix-----------------------------------------------------
1512 :
1513 : ! prepare to calculate overlap matrix
1514 864 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1515 864 : icg = dtefield%cgindex(ikpt,isppol+nsppol)
1516 864 : icg1 = dtefield%cgindex(ikpt1,isppol)
1517 : npw_k1 = npwarr(ikpt)
1518 864 : npw_k2 = npwar1(ikpt1)
1519 36720 : pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1520 :
1521 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1522 : jband_me = 0
1523 4320 : do jband = 1, dtefield%mband_occ
1524 3456 : if (band_procs(jband) == mpi_enreg%me_band) then
1525 3456 : jband_me = jband_me + 1
1526 : vect2(:,1:npw_k2) = &
1527 433728 : & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1528 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1529 : end if
1530 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1531 :
1532 18144 : do iband = 1, dtefield%mband_occ
1533 13824 : pwmin = (iband-1)*npw_k1*nspinor
1534 13824 : pwmax = pwmin + npw_k1*nspinor
1535 : vect1(:,1:npw_k1) = &
1536 1734912 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
1537 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1538 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1539 13824 : & vect1,vect2)
1540 13824 : s1mat(1,jband,iband) = s1mat(1,jband,iband) + dotr
1541 17280 : s1mat(2,jband,iband) = s1mat(2,jband,iband) + doti
1542 : end do ! iband
1543 : end do !jband
1544 :
1545 45792 : Amat(:,:,:)=zero
1546 :
1547 : ! calculate Amat
1548 4320 : do iband=1, dtefield%mband_occ
1549 18144 : do jband=1, dtefield%mband_occ
1550 72576 : do kband=1, dtefield%mband_occ
1551 : Amat(1,iband,jband) = Amat(1,iband,jband) + s1mat(1,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)&
1552 55296 : & - s1mat(2,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)
1553 : Amat(2,iband,jband) = Amat(2,iband,jband) + s1mat(1,iband,kband)*qmat(2,kband,jband,ikpt,2,idir)&
1554 69120 : & + s1mat(2,iband,kband)*qmat(1,kband,jband,ikpt,2,idir)
1555 : end do
1556 : end do
1557 : end do
1558 :
1559 45792 : Bmat(:,:,:)=zero
1560 :
1561 : ! calculate Bmat
1562 864 : ikptn = dtefield%ikpt_dk(ikpt,7,idir)
1563 4320 : do iband=1, dtefield%mband_occ
1564 18144 : do jband=1, dtefield%mband_occ
1565 72576 : do kband=1, dtefield%mband_occ
1566 : Bmat(1,jband,kband) = Bmat(1,jband,kband) + Amat(1,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)&
1567 55296 : & - Amat(2,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)
1568 : Bmat(2,jband,kband) = Bmat(2,jband,kband) + Amat(1,iband,kband)*qmat(2,jband,iband,ikptn,2,idir)&
1569 69120 : + Amat(2,iband,kband)*qmat(1,jband,iband,ikptn,2,idir)
1570 : end do
1571 : end do
1572 : end do
1573 :
1574 : ! calc. the second term of gradient------------------------------
1575 :
1576 : ! preparation
1577 :
1578 864 : ikptnp1 = dtefield%ikpt_dk(ikpt,4,idir)
1579 864 : icg = dtefield%cgindex(ikptnp1,isppol)
1580 864 : npw_k1 = npwar1(ikpt)
1581 864 : npw_k2 = npwarr(ikptnp1)
1582 36720 : pwind_tmp(1:npw_k1) =pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1583 : z1(:) = zero
1584 : z2(:) = zero
1585 37872 : do ipw = 1, npw_k1
1586 35856 : jpw = pwind_tmp(ipw)
1587 36720 : if (jpw > 0) then
1588 144720 : do iband = 1, dtefield%mband_occ
1589 115776 : wfr = cg(1,icg + (iband - 1)*npw_k2*nspinor + jpw)
1590 115776 : wfi = cg(2,icg + (iband - 1)*npw_k2*nspinor + jpw)
1591 607824 : do jband=1, dtefield%mband_occ
1592 463104 : grad_berry(1,ipw,jband) = grad_berry(1,ipw,jband) + fac*(Bmat(1,iband,jband)*wfr - Bmat(2,iband,jband)*wfi)
1593 578880 : grad_berry(2,ipw,jband) = grad_berry(2,ipw,jband) + fac*(Bmat(1,iband,jband)*wfi + Bmat(2,iband,jband)*wfr)
1594 : end do
1595 : end do
1596 : end if
1597 : end do
1598 :
1599 : end do !idir
1600 :
1601 : !!----------------------------------------third part of gradient------------------------------------------------------
1602 1152 : do idir=1,3
1603 : fac = rprimd(idir_efield,idir)*dble(nkpt)/&
1604 864 : & (dble(dtefield%nstr(idir))*four_pi)
1605 :
1606 : ! prepare
1607 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1608 864 : icg1 = dtefield%cgindex(ikpt1,isppol)
1609 864 : npw_k1 = npwar1(ikpt)
1610 864 : npw_k2 = npwarr(ikpt1)
1611 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
1612 36720 : do ipw = 1, npw_k1
1613 35856 : jpw = pwind_tmp(ipw)
1614 36720 : if (jpw > 0) then
1615 144720 : do iband = 1, dtefield%mband_occ
1616 115776 : wfr = cg(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1617 115776 : wfi = cg(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1618 607824 : do jband = 1, dtefield%mband_occ
1619 : grad_berry(1,ipw,jband) = &
1620 : & grad_berry(1,ipw,jband) + &
1621 463104 : & fac*qmat(1,iband,jband,ikpt,1,idir)*wfr - fac*qmat(2,iband,jband,ikpt,1,idir)*wfi
1622 : grad_berry(2,ipw,jband) = &
1623 : & grad_berry(2,ipw,jband) + &
1624 578880 : & fac*qmat(1,iband,jband,ikpt,1,idir)*wfi + fac*qmat(2,iband,jband,ikpt,1,idir)*wfr
1625 : end do
1626 : end do
1627 : end if
1628 : end do
1629 :
1630 : ! prepare
1631 864 : ikpt1 = dtefield%ikpt_dk(ikpt,2,idir)
1632 864 : icg1 = dtefield%cgindex(ikpt1,isppol)
1633 : npw_k1 = npwar1(ikpt)
1634 864 : npw_k2 = npwarr(ikpt1)
1635 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,6,idir)
1636 :
1637 37008 : do ipw = 1, npw_k1
1638 35856 : jpw = pwind_tmp(ipw)
1639 36720 : if (jpw > 0) then
1640 144720 : do iband = 1, dtefield%mband_occ
1641 115776 : wfr = cg(1,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1642 115776 : wfi = cg(2,icg1 + (iband - 1)*npw_k2*nspinor + jpw)
1643 607824 : do jband = 1, dtefield%mband_occ
1644 : grad_berry(1,ipw,jband) = &
1645 : & grad_berry(1,ipw,jband) - &
1646 463104 : & fac*qmat(1,iband,jband,ikpt,2,idir)*wfr + fac*qmat(2,iband,jband,ikpt,2,idir)*wfi
1647 :
1648 : grad_berry(2,ipw,jband) = &
1649 : & grad_berry(2,ipw,jband) - &
1650 578880 : & fac*qmat(1,iband,jband,ikpt,2,idir)*wfi - fac*qmat(2,iband,jband,ikpt,2,idir)*wfr
1651 :
1652 : end do
1653 : end do
1654 : end if
1655 : end do
1656 :
1657 : end do !idir
1658 :
1659 : ! accumulate full sum over iband in each jband entry
1660 288 : call xmpi_sum(grad_berry,mpi_enreg%comm_band,ierr)
1661 :
1662 288 : ABI_FREE(vect1)
1663 288 : ABI_FREE(vect2)
1664 288 : ABI_FREE(s1mat)
1665 288 : ABI_FREE(Amat)
1666 288 : ABI_FREE(Bmat)
1667 288 : ABI_FREE(pwind_tmp)
1668 :
1669 288 : end subroutine dfptff_gbefd
1670 : !!***
1671 :
1672 : !!****f* ABINIT/dfptff_edie
1673 : !! NAME
1674 : !! dfptff_edie
1675 : !!
1676 : !! FUNCTION
1677 : !! calculate the second order energy from the contribution of \Omega E \cdot P
1678 : !! term, Eq.(6) in PRB 75, 115116(2007) [[cite:Wang2007]].
1679 : !!
1680 : !! INPUTS
1681 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
1682 : !! cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
1683 : !! RF wavefunctions at k,q.
1684 : !! dtefield = variables related to response Berry-phase calculation
1685 : !! ikpt = the index of the current k point
1686 : !! isppol = the index of the spin component
1687 : !! mband = maximum number of bands
1688 : !! mkmem_rbz = maximum number of k-points in core memory
1689 : !! mpw = maximum number of plane waves
1690 : !! mpw1 = maximum number of plane waves for response wavefunctions
1691 : !! nkpt = number of k points
1692 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
1693 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
1694 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
1695 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
1696 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
1697 : !! inverse of the overlap matrix
1698 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
1699 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
1700 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
1701 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
1702 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
1703 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
1704 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
1705 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
1706 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
1707 : !!
1708 : !! OUTPUT
1709 : !! eberry = the energy of the Berry phase term
1710 : !!
1711 : !! SOURCE
1712 :
1713 36 : subroutine dfptff_edie(cg,cg1,dtefield,eberry,idir_efield,mband,mband_mem,mkmem_rbz,&
1714 36 : & mpi_enreg,mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
1715 :
1716 : !Arguments ----------------------------------------
1717 : !scalars
1718 : integer,intent(in) :: idir_efield,mband,mkmem_rbz,mpw,mpw1,nkpt,nspinor,nsppol
1719 : integer,intent(in) :: mband_mem
1720 : real(dp),intent(out) :: eberry
1721 : type(efield_type),intent(in) :: dtefield
1722 : type(MPI_type),intent(in) :: mpi_enreg
1723 : !arrays
1724 : integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
1725 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
1726 : real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol)
1727 : real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol)
1728 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
1729 : real(dp),intent(in) :: rprimd(3,3)
1730 :
1731 : !Local variables ----------------------------------
1732 : !scalars
1733 : integer :: iband,icg,icg1,idir
1734 : integer :: jband_me, iband_me, me_band, ierr
1735 : integer :: ikpt,ikpt1,ikptn,ikptnm
1736 : integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
1737 : real(dp) :: doti,dotr,e0,fac
1738 : !arrays
1739 36 : integer,allocatable :: pwind_tmp(:)
1740 72 : integer :: band_procs(mband)
1741 : real(dp) :: z1(2)
1742 36 : real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:)
1743 :
1744 : ! *************************************************************************
1745 :
1746 : !calculate 4 matrices -----------------------------
1747 36 : mpw_tmp=max(mpw,mpw1)
1748 180 : ABI_MALLOC(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4))
1749 108 : ABI_MALLOC(vect1,(2,0:mpw_tmp))
1750 72 : ABI_MALLOC(vect2,(2,0:mpw_tmp))
1751 108 : ABI_MALLOC(pwind_tmp,(mpw_tmp))
1752 144 : ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
1753 216 : vect1(:,0) = zero ; vect2(:,0) = zero
1754 36 : eberry=zero
1755 :
1756 36 : me_band = mpi_enreg%me_band
1757 :
1758 324 : do ikpt=1,nkpt
1759 :
1760 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,&
1761 288 : & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
1762 :
1763 1188 : do idir=1,3
1764 : fac = dtefield%efield_dot(idir)/&
1765 864 : & (dble(dtefield%nstr(idir))*four_pi)
1766 :
1767 : ! compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix----------------------------------------------------
1768 :
1769 : ! prepare to calculate overlap matrix
1770 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1771 864 : icg = dtefield%cgindex(ikpt,1+nsppol)
1772 864 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1773 864 : npw_k1 = npwar1(ikpt)
1774 864 : npw_k2 = npwar1(ikpt1)
1775 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
1776 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1777 : jband_me = 0
1778 4320 : do jband = 1, dtefield%mband_occ
1779 3456 : if (band_procs(jband) == me_band) then
1780 3456 : jband_me = jband_me + 1
1781 : vect2(:,1:npw_k2) = &
1782 433728 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
1783 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1784 : end if
1785 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1786 :
1787 3456 : iband_me = 0
1788 18144 : do iband = 1, dtefield%mband_occ
1789 13824 : if (band_procs(iband) /= me_band) cycle
1790 13824 : iband_me = iband_me + 1
1791 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1792 13824 : pwmax = pwmin + npw_k1*nspinor
1793 : vect1(:,1:npw_k1) = &
1794 1734912 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
1795 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1796 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1797 13824 : & vect1,vect2)
1798 13824 : umat(1,iband,jband,1) = dotr
1799 17280 : umat(2,iband,jband,1) = doti
1800 : end do ! iband
1801 : end do !jband
1802 :
1803 : ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix----------------------------------------------------
1804 :
1805 : ! prepare to calculate overlap matrix
1806 864 : ikpt1 = dtefield%ikpt_dk(ikpt,5,idir)
1807 864 : icg = dtefield%cgindex(ikpt,1)
1808 864 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1809 864 : npw_k1 = npwarr(ikpt)
1810 864 : npw_k2 = npwar1(ikpt1)
1811 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
1812 :
1813 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1814 : jband_me = 0
1815 4320 : do jband = 1, dtefield%mband_occ
1816 3456 : if (band_procs(jband) == me_band) then
1817 3456 : jband_me = jband_me + 1
1818 : vect2(:,1:npw_k2) = &
1819 433728 : & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1820 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1821 : end if
1822 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1823 :
1824 3456 : iband_me = 0
1825 18144 : do iband = 1, dtefield%mband_occ
1826 13824 : if (band_procs(iband) /= me_band) cycle
1827 13824 : iband_me = iband_me + 1
1828 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1829 13824 : pwmax = pwmin + npw_k1*nspinor
1830 : vect1(:,1:npw_k1) = &
1831 1734912 : & cg(:,icg + 1 + pwmin:icg + pwmax)
1832 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1833 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1834 13824 : & vect1,vect2)
1835 13824 : umat(1,iband,jband,2) = dotr
1836 17280 : umat(2,iband,jband,2) = doti
1837 : end do ! iband
1838 : end do !jband
1839 :
1840 : ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix----------------------------------------------------
1841 :
1842 : ! prepare to calculate overlap matrix
1843 864 : ikptn = dtefield%ikpt_dk(ikpt,8,idir)
1844 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1845 864 : icg = dtefield%cgindex(ikptn,1+nsppol)
1846 864 : icg1 = dtefield%cgindex(ikpt1,1)
1847 864 : npw_k1 = npwar1(ikptn)
1848 864 : npw_k2 = npwarr(ikpt1)
1849 36720 : pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir)
1850 :
1851 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1852 : jband_me = 0
1853 4320 : do jband = 1, dtefield%mband_occ
1854 3456 : if (band_procs(jband) == me_band) then
1855 3456 : jband_me = jband_me + 1
1856 : vect2(:,1:npw_k2) = &
1857 433728 : & cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1858 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1859 : end if
1860 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1861 :
1862 3456 : iband_me = 0
1863 18144 : do iband = 1, dtefield%mband_occ
1864 13824 : if (band_procs(iband) /= me_band) cycle
1865 13824 : iband_me = iband_me + 1
1866 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1867 13824 : pwmax = pwmin + npw_k1*nspinor
1868 : vect1(:,1:npw_k1) = &
1869 1734912 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
1870 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1871 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1872 13824 : & vect1,vect2)
1873 13824 : umat(1,iband,jband,3) = dotr
1874 17280 : umat(2,iband,jband,3) = doti
1875 : end do ! iband
1876 : end do !jband
1877 :
1878 : ! compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix----------------------------------------------------
1879 :
1880 : ! prepare to calculate overlap matrix
1881 864 : ikptn = dtefield%ikpt_dk(ikpt,5,idir)
1882 864 : ikptnm = dtefield%ikpt_dk(ikptn,9,idir)
1883 864 : ikpt1 = dtefield%ikpt_dk(ikpt,9,idir)
1884 864 : icg = dtefield%cgindex(ikptnm,1)
1885 864 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
1886 864 : npw_k1 = npwarr(ikptnm)
1887 864 : npw_k2 = npwar1(ikpt1)
1888 36720 : pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir)
1889 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1890 : jband_me = 0
1891 4320 : do jband = 1, dtefield%mband_occ
1892 3456 : if (band_procs(jband) == me_band) then
1893 3456 : jband_me = jband_me + 1
1894 : vect2(:,1:npw_k2) = &
1895 433728 : & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1896 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1897 : end if
1898 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1899 :
1900 3456 : iband_me = 0
1901 18144 : do iband = 1, dtefield%mband_occ
1902 13824 : if (band_procs(iband) /= me_band) cycle
1903 13824 : iband_me = iband_me + 1
1904 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1905 13824 : pwmax = pwmin + npw_k1*nspinor
1906 : vect1(:,1:npw_k1) = &
1907 1734912 : & cg(:,icg + 1 + pwmin:icg + pwmax)
1908 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1909 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1910 13824 : & vect1,vect2)
1911 13824 : umat(1,iband,jband,4) = dotr
1912 17280 : umat(2,iband,jband,4) = doti
1913 : end do ! iband
1914 : end do !jband
1915 :
1916 : ! sum up over all iband and all procs
1917 864 : call xmpi_sum(umat,mpi_enreg%comm_band,ierr)
1918 :
1919 : ! sum over the whole------------------------------------------------------------
1920 :
1921 864 : e0=zero
1922 4320 : do iband=1,dtefield%mband_occ
1923 18144 : do jband=1,dtefield%mband_occ
1924 : e0 = e0 + 2_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)&
1925 17280 : & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir))
1926 : end do
1927 : end do
1928 864 : eberry = eberry - e0*fac
1929 864 : e0=zero
1930 864 : ikptn=dtefield%ikpt_dk(ikpt,8,idir)
1931 45792 : Amat(:,:,:)=zero
1932 : ! calculate Amat
1933 4320 : do iband=1, dtefield%mband_occ
1934 18144 : do jband=1, dtefield%mband_occ
1935 72576 : do kband=1, dtefield%mband_occ
1936 : Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)&
1937 55296 : & - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)
1938 : Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)&
1939 69120 : & + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)
1940 : end do
1941 : end do
1942 : end do
1943 :
1944 4320 : do iband=1, dtefield%mband_occ
1945 18144 : do jband=1, dtefield%mband_occ
1946 72576 : do kband=1, dtefield%mband_occ
1947 : z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir)&
1948 55296 : & - (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir)
1949 : z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*qmat(2,jband,kband,ikptn,1,idir)&
1950 55296 : & + (umat(2,jband,iband,4)+umat(2,iband,jband,2))*qmat(1,jband,kband,ikptn,1,idir)
1951 :
1952 69120 : e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband))
1953 : end do
1954 : end do
1955 : end do
1956 :
1957 864 : eberry = eberry - e0*fac
1958 :
1959 : ! !---------------------------------last part---------------------------------------------
1960 :
1961 : fac = rprimd(idir_efield,idir)/&
1962 864 : & (dble(dtefield%nstr(idir))*two_pi)
1963 :
1964 : ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix----------------------------------------------------
1965 :
1966 : ! prepare to calculate overlap matrix
1967 864 : ikptn = dtefield%ikpt_dk(ikpt,8,idir)
1968 864 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
1969 864 : icg = dtefield%cgindex(ikptn,1+nsppol)
1970 864 : icg1 = dtefield%cgindex(ikpt1,1)
1971 864 : npw_k1 = npwar1(ikptn)
1972 864 : npw_k2 = npwarr(ikpt1)
1973 36720 : pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir)
1974 :
1975 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
1976 : jband_me = 0
1977 4320 : do jband = 1, dtefield%mband_occ
1978 3456 : if (band_procs(jband) == me_band) then
1979 3456 : jband_me = jband_me + 1
1980 : vect2(:,1:npw_k2) = &
1981 433728 : & cg(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
1982 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
1983 : end if
1984 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
1985 :
1986 3456 : iband_me = 0
1987 18144 : do iband = 1, dtefield%mband_occ
1988 13824 : if (band_procs(iband) /= me_band) cycle
1989 13824 : iband_me = iband_me + 1
1990 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
1991 13824 : pwmax = pwmin + npw_k1*nspinor
1992 : vect1(:,1:npw_k1) = &
1993 1734912 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
1994 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
1995 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
1996 13824 : & vect1,vect2)
1997 13824 : umat(1,iband,jband,1) = dotr
1998 17280 : umat(2,iband,jband,1) = doti
1999 : end do ! iband
2000 : end do !jband
2001 :
2002 : ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix----------------------------------------------------
2003 :
2004 : ! prepare to calculate overlap matrix
2005 864 : ikpt1 = dtefield%ikpt_dk(ikpt,5,idir)
2006 864 : icg = dtefield%cgindex(ikpt,1)
2007 864 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2008 864 : npw_k1 = npwarr(ikpt)
2009 864 : npw_k2 = npwar1(ikpt1)
2010 36720 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2011 4320 : vect1(:,0) = zero ; vect2(:,0) = zero
2012 : jband_me = 0
2013 4320 : do jband = 1, dtefield%mband_occ
2014 3456 : if (band_procs(jband) == me_band) then
2015 3456 : jband_me = jband_me + 1
2016 : vect2(:,1:npw_k2) = &
2017 433728 : & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2018 19008 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2019 : end if
2020 3456 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2021 :
2022 3456 : iband_me = 0
2023 18144 : do iband = 1, dtefield%mband_occ
2024 13824 : if (band_procs(iband) /= me_band) cycle
2025 13824 : iband_me = iband_me + 1
2026 13824 : pwmin = (iband_me-1)*npw_k1*nspinor
2027 13824 : pwmax = pwmin + npw_k1*nspinor
2028 : vect1(:,1:npw_k1) = &
2029 1734912 : & cg(:,icg + 1 + pwmin:icg + pwmax)
2030 76032 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2031 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2032 13824 : & vect1,vect2)
2033 13824 : umat(1,iband,jband,1) = umat(1,iband,jband,1) + dotr
2034 17280 : umat(2,iband,jband,1) = umat(2,iband,jband,1) + doti
2035 : end do ! iband
2036 : end do !jband
2037 :
2038 : ! sum up over all iband and all procs
2039 864 : call xmpi_sum(umat,mpi_enreg%comm_band,ierr)
2040 :
2041 864 : e0=zero
2042 :
2043 4320 : do iband=1,dtefield%mband_occ
2044 18144 : do jband=1,dtefield%mband_occ
2045 : e0 = e0 + (umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)&
2046 17280 : & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir))
2047 : end do
2048 : end do
2049 :
2050 2880 : eberry = eberry - e0*fac
2051 :
2052 : end do !end idir
2053 : end do !end ikpt
2054 :
2055 36 : ABI_FREE(umat)
2056 36 : ABI_FREE(vect1)
2057 36 : ABI_FREE(vect2)
2058 36 : ABI_FREE(pwind_tmp)
2059 36 : ABI_FREE(Amat)
2060 :
2061 36 : end subroutine dfptff_edie
2062 : !!***
2063 :
2064 : !!****f* ABINIT/dfptff_ebp
2065 : !! NAME
2066 : !! dfptff_ebp
2067 : !!
2068 : !! FUNCTION
2069 : !! calculation of the energy from the term \Omega E \cdot P
2070 : !!
2071 : !! INPUTS
2072 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
2073 : !! cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
2074 : !! RF wavefunctions at k,q.
2075 : !! dtefield = variables related to response Berry-phase calculation
2076 : !! ikpt = the index of the current k point
2077 : !! isppol = the index of the spin component
2078 : !! mband = maximum number of bands
2079 : !! mkmem_rbz = maximum number of k-points in core memory
2080 : !! mpw = maximum number of plane waves
2081 : !! mpw1 = maximum number of plane waves for response wavefunctions
2082 : !! nkpt = number of k points
2083 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
2084 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
2085 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
2086 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
2087 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
2088 : !! inverse of the overlap matrix
2089 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
2090 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
2091 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
2092 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
2093 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
2094 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
2095 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
2096 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
2097 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
2098 : !!
2099 : !! OUTPUT
2100 : !! grad_berry(2,mpw1,dtefield%mband_occ) = the gradient of the Berry phase term
2101 : !!
2102 : !! SOURCE
2103 :
2104 61 : subroutine dfptff_ebp(cg,cg1,dtefield,eberry,mband,mband_mem,mkmem_rbz,&
2105 61 : & mpi_enreg,mpw,mpw1,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat)
2106 :
2107 : !Arguments ----------------------------------------
2108 : !scalars
2109 : integer,intent(in) :: mband,mkmem_rbz,mpw,mpw1,nkpt,nspinor,nsppol
2110 : integer,intent(in) :: mband_mem
2111 : real(dp),intent(out) :: eberry
2112 : type(efield_type),intent(in) :: dtefield
2113 : type(MPI_type),intent(in) :: mpi_enreg
2114 : !arrays
2115 : integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
2116 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
2117 : real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol)
2118 : real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol)
2119 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2120 :
2121 : !Local variables ----------------------------------
2122 : !scalars
2123 : integer :: iband,icg,icg1,idir
2124 : integer :: jband_me, iband_me, me_band, ierr
2125 : integer :: ikpt,ikpt1,ikptn,ikptnm
2126 : integer :: jband,kband,mpw_tmp,npw_k1,npw_k2,pwmax,pwmin
2127 : real(dp) :: doti,dotr,e0,fac
2128 : !arrays
2129 122 : integer :: band_procs(mband)
2130 61 : integer,allocatable :: pwind_tmp(:)
2131 : real(dp) :: z1(2)
2132 61 : real(dp),allocatable :: Amat(:,:,:),umat(:,:,:,:),vect1(:,:),vect2(:,:)
2133 :
2134 : ! *************************************************************************
2135 :
2136 : !calculate 4 matrices -----------------------------
2137 61 : mpw_tmp=max(mpw,mpw1)
2138 305 : ABI_MALLOC(umat,(2,dtefield%mband_occ,dtefield%mband_occ,4))
2139 183 : ABI_MALLOC(vect1,(2,0:mpw_tmp))
2140 122 : ABI_MALLOC(vect2,(2,0:mpw_tmp))
2141 183 : ABI_MALLOC(pwind_tmp,(mpw_tmp))
2142 244 : ABI_MALLOC(Amat,(2,dtefield%mband_occ,dtefield%mband_occ))
2143 366 : vect1(:,0) = zero ; vect2(:,0) = zero
2144 61 : eberry=zero
2145 :
2146 61 : me_band = mpi_enreg%me_band
2147 :
2148 : !TODO no info on sppol here - I default to isppol 1
2149 549 : do ikpt=1,nkpt
2150 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,&
2151 488 : & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
2152 :
2153 :
2154 2013 : do idir=1,3
2155 :
2156 : fac = dtefield%efield_dot(idir)/&
2157 1464 : & (dble(dtefield%nstr(idir))*four_pi)
2158 :
2159 : ! compute <u^(1)_{k_j,q}|u^(1)_{k_j+1,q}> matrix---------
2160 :
2161 : ! prepare to calculate overlap matrix
2162 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2163 1464 : icg = dtefield%cgindex(ikpt,1+nsppol)
2164 1464 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2165 1464 : npw_k1 = npwar1(ikpt)
2166 1464 : npw_k2 = npwar1(ikpt1)
2167 62220 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,3,idir)
2168 7320 : vect1(:,0) = zero ; vect2(:,0) = zero
2169 : jband_me = 0
2170 7320 : do jband = 1, dtefield%mband_occ
2171 5856 : if (band_procs(jband) == me_band) then
2172 5856 : jband_me = jband_me + 1
2173 734928 : vect2(:,1:npw_k2) = cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2174 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2175 : end if
2176 5856 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2177 :
2178 :
2179 5856 : iband_me = 0
2180 30744 : do iband = 1, dtefield%mband_occ
2181 23424 : if (band_procs(iband) /= me_band) cycle
2182 23424 : iband_me = iband_me + 1
2183 :
2184 23424 : pwmin = (iband_me-1)*npw_k1*nspinor
2185 23424 : pwmax = pwmin + npw_k1*nspinor
2186 : vect1(:,1:npw_k1) = &
2187 2939712 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
2188 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2189 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2190 23424 : & vect1,vect2)
2191 23424 : umat(1,iband,jband,1) = dotr
2192 29280 : umat(2,iband,jband,1) = doti
2193 :
2194 : end do ! iband
2195 : end do !jband
2196 :
2197 : ! compute <u^(0)_{k_j}|u^(1)_{k_j-n+1,q}> matrix----------------------------------------------------
2198 :
2199 : ! prepare to calculate overlap matrix
2200 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,5,idir)
2201 1464 : icg = dtefield%cgindex(ikpt,1)
2202 1464 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2203 1464 : npw_k1 = npwarr(ikpt)
2204 1464 : npw_k2 = npwar1(ikpt1)
2205 62220 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2206 :
2207 7320 : vect1(:,0) = zero ; vect2(:,0) = zero
2208 : jband_me = 0
2209 7320 : do jband = 1, dtefield%mband_occ
2210 5856 : if (band_procs(jband) == me_band) then
2211 5856 : jband_me = jband_me + 1
2212 : vect2(:,1:npw_k2) = &
2213 734928 : & cg1(:,icg1 + 1 + (jband-1)*npw_k2*nspinor:icg1 + jband*npw_k2*nspinor)
2214 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2215 : end if
2216 5856 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2217 :
2218 5856 : iband_me = 0
2219 30744 : do iband = 1, dtefield%mband_occ
2220 23424 : if (band_procs(iband) /= me_band) cycle
2221 23424 : iband_me = iband_me + 1
2222 :
2223 23424 : pwmin = (iband_me-1)*npw_k1*nspinor
2224 23424 : pwmax = pwmin + npw_k1*nspinor
2225 : vect1(:,1:npw_k1) = &
2226 2939712 : & cg(:,icg + 1 + pwmin:icg + pwmax)
2227 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2228 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2229 23424 : & vect1,vect2)
2230 23424 : umat(1,iband,jband,2) = dotr
2231 29280 : umat(2,iband,jband,2) = doti
2232 :
2233 : end do ! iband
2234 : end do !jband
2235 :
2236 : ! compute <u^(1)_{k_j-n,q}|u^(0)_{k_j+1}> matrix----------------------------------------------------
2237 :
2238 : ! prepare to calculate overlap matrix
2239 1464 : ikptn = dtefield%ikpt_dk(ikpt,8,idir)
2240 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2241 1464 : icg = dtefield%cgindex(ikptn,1+nsppol)
2242 1464 : icg1 = dtefield%cgindex(ikpt1,1)
2243 1464 : npw_k1 = npwar1(ikptn)
2244 1464 : npw_k2 = npwarr(ikpt1)
2245 62220 : pwind_tmp(1:npw_k1) = pwindall((ikptn-1)*mpw_tmp+1:(ikptn-1)*mpw_tmp+npw_k1,5,idir)
2246 :
2247 7320 : vect1(:,0) = zero ; vect2(:,0) = zero
2248 : jband_me = 0
2249 7320 : do jband = 1, dtefield%mband_occ
2250 5856 : if (band_procs(jband) == me_band) then
2251 5856 : jband_me = jband_me + 1
2252 : vect2(:,1:npw_k2) = &
2253 734928 : & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2254 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2255 : end if
2256 5856 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2257 :
2258 5856 : iband_me = 0
2259 30744 : do iband = 1, dtefield%mband_occ
2260 23424 : if (band_procs(iband) /= me_band) cycle
2261 23424 : iband_me = iband_me + 1
2262 :
2263 23424 : pwmin = (iband_me-1)*npw_k1*nspinor
2264 23424 : pwmax = pwmin + npw_k1*nspinor
2265 : vect1(:,1:npw_k1) = &
2266 2939712 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
2267 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2268 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2269 23424 : & vect1,vect2)
2270 23424 : umat(1,iband,jband,3) = dotr
2271 29280 : umat(2,iband,jband,3) = doti
2272 :
2273 : end do ! iband
2274 : end do !jband
2275 :
2276 : ! compute <u^(0)_{-k_j-n+1}|u^(1)_{-k_j,q}> matrix----------------------------------------------------
2277 :
2278 : ! prepare to calculate overlap matrix
2279 1464 : ikptn = dtefield%ikpt_dk(ikpt,5,idir)
2280 1464 : ikptnm = dtefield%ikpt_dk(ikptn,9,idir)
2281 1464 : ikpt1 = dtefield%ikpt_dk(ikpt,9,idir)
2282 1464 : icg = dtefield%cgindex(ikptnm,1)
2283 1464 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2284 1464 : npw_k1 = npwarr(ikptnm)
2285 1464 : npw_k2 = npwar1(ikpt1)
2286 62220 : pwind_tmp(1:npw_k1) = pwindall((ikptnm-1)*mpw_tmp+1:(ikptnm-1)*mpw_tmp+npw_k1,7,idir)
2287 :
2288 7320 : vect1(:,0) = zero ; vect2(:,0) = zero
2289 : jband_me = 0
2290 7320 : do jband = 1, dtefield%mband_occ
2291 5856 : if (band_procs(jband) == me_band) then
2292 5856 : jband_me = jband_me + 1
2293 : vect2(:,1:npw_k2) = &
2294 734928 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2295 32208 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2296 : end if
2297 5856 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2298 :
2299 5856 : iband_me = 0
2300 30744 : do iband = 1, dtefield%mband_occ
2301 23424 : if (band_procs(iband) /= me_band) cycle
2302 23424 : iband_me = iband_me + 1
2303 :
2304 23424 : pwmin = (iband_me-1)*npw_k1*nspinor
2305 23424 : pwmax = pwmin + npw_k1*nspinor
2306 : vect1(:,1:npw_k1) = &
2307 2939712 : & cg(:,icg + 1 + pwmin:icg + pwmax)
2308 128832 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2309 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2310 23424 : & vect1,vect2)
2311 23424 : umat(1,iband,jband,4) = dotr
2312 29280 : umat(2,iband,jband,4) = doti
2313 :
2314 : end do ! iband
2315 : end do !jband
2316 :
2317 : ! recompose full umat on each proc, summing over iband
2318 1464 : call xmpi_sum(umat,mpi_enreg%comm_band,ierr)
2319 :
2320 : ! sum over the whole------------------------------------------------------------
2321 :
2322 1464 : e0=zero
2323 7320 : do iband=1,dtefield%mband_occ
2324 30744 : do jband=1,dtefield%mband_occ
2325 : e0 = e0 + 4_dp*(umat(1,iband,jband,1)*qmat(2,jband,iband,ikpt,1,idir)&
2326 29280 : & + umat(2,iband,jband,1)*qmat(1,jband,iband,ikpt,1,idir))
2327 :
2328 : end do
2329 : end do
2330 :
2331 1464 : eberry = eberry - e0*fac
2332 :
2333 1464 : e0=zero
2334 :
2335 1464 : ikptn=dtefield%ikpt_dk(ikpt,8,idir)
2336 :
2337 77592 : Amat(:,:,:)=zero
2338 :
2339 : ! calculate Amat
2340 7320 : do iband=1, dtefield%mband_occ
2341 30744 : do jband=1, dtefield%mband_occ
2342 122976 : do kband=1, dtefield%mband_occ
2343 : Amat(1,iband,jband) = Amat(1,iband,jband) + (umat(1,iband,kband,3))*&
2344 : & qmat(1,kband,jband,ikpt,1,idir)&
2345 93696 : & - (umat(2,iband,kband,3))*qmat(2,kband,jband,ikpt,1,idir)
2346 : Amat(2,iband,jband) = Amat(2,iband,jband) + (umat(1,iband,kband,3))*&
2347 : & qmat(2,kband,jband,ikpt,1,idir)&
2348 117120 : & + (umat(2,iband,kband,3))*qmat(1,kband,jband,ikpt,1,idir)
2349 : end do
2350 : end do
2351 : end do
2352 :
2353 7320 : do iband=1, dtefield%mband_occ
2354 30744 : do jband=1, dtefield%mband_occ
2355 122976 : do kband=1, dtefield%mband_occ
2356 :
2357 : z1(1) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*&
2358 : & qmat(1,jband,kband,ikptn,1,idir)&
2359 : & - (umat(2,jband,iband,4)+umat(2,iband,jband,2))*&
2360 93696 : & qmat(2,jband,kband,ikptn,1,idir)
2361 : z1(2) = (umat(1,jband,iband,4)+umat(1,iband,jband,2))*&
2362 : & qmat(2,jband,kband,ikptn,1,idir)&
2363 : & + (umat(2,jband,iband,4)+umat(2,iband,jband,2))*&
2364 93696 : & qmat(1,jband,kband,ikptn,1,idir)
2365 :
2366 117120 : e0 = e0 - 4_dp*(z1(1)*Amat(2,kband,iband)+z1(2)*Amat(1,kband,iband))
2367 :
2368 : end do
2369 : end do
2370 : end do
2371 :
2372 3416 : eberry = eberry - e0*fac
2373 :
2374 : end do !end idir
2375 : end do !end ikpt
2376 :
2377 61 : ABI_FREE(umat)
2378 61 : ABI_FREE(vect1)
2379 61 : ABI_FREE(vect2)
2380 61 : ABI_FREE(pwind_tmp)
2381 61 : ABI_FREE(Amat)
2382 :
2383 61 : end subroutine dfptff_ebp
2384 : !!***
2385 :
2386 : !!****f* ABINIT/dfptff_die
2387 : !! NAME
2388 : !! dfptff_die
2389 : !!
2390 : !! FUNCTION
2391 : !! calculate electric susceptibility tensor in Eq.(28) in PRB 75, 115116(2007) [[cite:Wang2007]].
2392 : !!
2393 : !! INPUTS
2394 : !! cg(2,mpw*nspinor*mband*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
2395 : !! cg1(2,mpw1*nspinor*mband*mk1mem*nsppol) = pw coefficients of
2396 : !! RF wavefunctions at k,q.
2397 : !! dtefield = variables related to response Berry-phase calculation
2398 : !! idirpert = the current coloumn of the dielectric permittivity tensor
2399 : !! mband = maximum number of bands
2400 : !! mkmem_rbz = maximum number of k-points in core memory
2401 : !! mpw = maximum number of plane waves
2402 : !! mpw1 = maximum number of plane waves for response wavefunctions
2403 : !! nkpt = number of k points
2404 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
2405 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
2406 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
2407 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
2408 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
2409 : !! inverse of the overlap matrix
2410 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
2411 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
2412 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
2413 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
2414 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
2415 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
2416 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
2417 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
2418 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
2419 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2420 : !!
2421 : !! OUTPUT
2422 : !! diet = electric susceptibility tensor
2423 : !!
2424 : !! SOURCE
2425 :
2426 3 : subroutine dfptff_die(cg,cg1,dtefield,d2lo,idirpert,ipert,mband,mband_mem,mkmem_rbz,&
2427 3 : & mpi_enreg,mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
2428 :
2429 : !Arguments ----------------------------------------
2430 : !scalars
2431 : integer,intent(in) :: idirpert,ipert,mband,mkmem_rbz,mpert,mpw,mpw1,nkpt,nspinor
2432 : integer,intent(in) :: mband_mem
2433 : integer,intent(in) :: nsppol
2434 : type(efield_type),intent(in) :: dtefield
2435 : type(MPI_type),intent(in) :: mpi_enreg
2436 : !arrays
2437 : integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
2438 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
2439 : real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol)
2440 : real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol)
2441 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2442 : real(dp),intent(in) :: rprimd(3,3)
2443 : real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i
2444 :
2445 : !Local variables ----------------------------------
2446 : !scalars
2447 : integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1
2448 : integer :: jband_me, iband_me, me_band, ierr
2449 : integer :: npw_k2,pwmax,pwmin
2450 : real(dp) :: doti,dotr,e0,fac
2451 : !arrays
2452 6 : integer :: band_procs(mband)
2453 3 : integer,allocatable :: pwind_tmp(:)
2454 : real(dp) :: edir(3)
2455 3 : real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:)
2456 :
2457 : ! *************************************************************************
2458 :
2459 : !calculate s1 matrices -----------------------------
2460 3 : mpw_tmp=max(mpw,mpw1)
2461 12 : ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
2462 9 : ABI_MALLOC(vect1,(2,0:mpw_tmp))
2463 6 : ABI_MALLOC(vect2,(2,0:mpw_tmp))
2464 9 : ABI_MALLOC(pwind_tmp,(mpw_tmp))
2465 18 : vect1(:,0) = zero ; vect2(:,0) = zero
2466 :
2467 3 : edir(:)=zero
2468 :
2469 3 : me_band = mpi_enreg%me_band
2470 :
2471 : !TODO no info on sppol here - I default to isppol 1
2472 27 : do ikpt=1,nkpt
2473 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,&
2474 24 : & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
2475 99 : do idir=1,3
2476 : ! compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ----------------------------------------
2477 :
2478 : ! prepare to calculate overlap matrix
2479 72 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2480 72 : icg = dtefield%cgindex(ikpt,1)
2481 72 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2482 72 : npw_k1 = npwarr(ikpt)
2483 72 : npw_k2 = npwar1(ikpt1)
2484 3060 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2485 :
2486 360 : vect1(:,0) = zero ; vect2(:,0) = zero
2487 : jband_me = 0
2488 360 : do jband = 1, dtefield%mband_occ
2489 288 : if (band_procs(jband) == me_band) then
2490 288 : jband_me = jband_me + 1
2491 : vect2(:,1:npw_k2) = &
2492 36144 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2493 1584 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2494 : end if
2495 288 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2496 :
2497 288 : iband_me = 0
2498 1512 : do iband = 1, dtefield%mband_occ
2499 1152 : if (band_procs(iband) /= me_band) cycle
2500 1152 : iband_me = iband_me + 1
2501 1152 : pwmin = (iband_me-1)*npw_k1*nspinor
2502 1152 : pwmax = pwmin + npw_k1*nspinor
2503 : vect1(:,1:npw_k1) = &
2504 144576 : & cg(:,icg + 1 + pwmin:icg + pwmax)
2505 6336 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2506 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2507 1152 : & vect1,vect2)
2508 1152 : s1mat(1,iband,jband) = dotr
2509 1440 : s1mat(2,iband,jband) = doti
2510 : end do ! iband
2511 : end do !jband
2512 :
2513 : ! compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 -------------------------------------
2514 :
2515 : ! prepare to calculate overlap matrix
2516 72 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2517 72 : icg = dtefield%cgindex(ikpt,1+nsppol)
2518 72 : icg1 = dtefield%cgindex(ikpt1,1)
2519 72 : npw_k1 = npwar1(ikpt)
2520 72 : npw_k2 = npwarr(ikpt1)
2521 3060 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
2522 360 : vect1(:,0) = zero ; vect2(:,0) = zero
2523 : jband_me = 0
2524 360 : do jband = 1, dtefield%mband_occ
2525 288 : if (band_procs(jband) == me_band) then
2526 288 : jband_me = jband_me + 1
2527 : vect2(:,1:npw_k2) = &
2528 36144 : & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2529 1584 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2530 : end if
2531 288 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2532 :
2533 288 : iband_me = 0
2534 1512 : do iband = 1, dtefield%mband_occ
2535 1152 : if (band_procs(iband) /= me_band) cycle
2536 1152 : iband_me = iband_me + 1
2537 1152 : pwmin = (iband_me-1)*npw_k1*nspinor
2538 1152 : pwmax = pwmin + npw_k1*nspinor
2539 : vect1(:,1:npw_k1) = &
2540 144576 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
2541 6336 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2542 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2543 1152 : & vect1,vect2)
2544 1152 : s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr
2545 1440 : s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti
2546 : end do ! iband
2547 : end do !jband
2548 :
2549 : ! recompose full s1mat on each proc, summing over iband
2550 72 : call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr)
2551 :
2552 : ! sum over the whole------------------------------------------------------------
2553 :
2554 72 : e0=zero
2555 :
2556 360 : do iband=1,dtefield%mband_occ
2557 1512 : do jband=1,dtefield%mband_occ
2558 : e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)&
2559 1440 : & + s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir))
2560 :
2561 : end do
2562 : end do
2563 :
2564 384 : do ialpha=1,3
2565 : fac = rprimd(ialpha,idir)/&
2566 216 : & (dble(dtefield%nstr(idir))*pi)
2567 288 : edir(ialpha)=edir(ialpha)+ e0*fac
2568 : end do
2569 :
2570 : end do !idir
2571 : end do !ikpt
2572 :
2573 12 : d2lo(1,1:3,ipert,idirpert,ipert)=edir(:)
2574 :
2575 3 : ABI_FREE(s1mat)
2576 3 : ABI_FREE(vect1)
2577 3 : ABI_FREE(vect2)
2578 3 : ABI_FREE(pwind_tmp)
2579 :
2580 3 : end subroutine dfptff_die
2581 : !!***
2582 :
2583 : !!****f* ABINIT/dfptff_bec
2584 : !! NAME
2585 : !! dfptff_bec
2586 : !!
2587 : !! FUNCTION
2588 : !! calculate Born effective charge tensor in Eq.(33) in PRB 75, 115116(2007) [[cite:Wang2007]].
2589 : !!
2590 : !! INPUTS
2591 : !! cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
2592 : !! cg1(2,mpw1*nspinor*mband_mem*mk1mem*nsppol) = pw coefficients of
2593 : !! RF wavefunctions at k,q.
2594 : !! dtefield = variables related to response Berry-phase calculation
2595 : !! idirpert = the current coloumn of the dielectric permittivity tensor
2596 : !! mband = maximum number of bands
2597 : !! mband_mem = maximum number of bands on this cpu
2598 : !! mkmem_rbz = maximum number of k-points in core memory
2599 : !! mpw = maximum number of plane waves
2600 : !! mpw1 = maximum number of plane waves for response wavefunctions
2601 : !! nkpt = number of k points
2602 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
2603 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
2604 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
2605 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
2606 : !! qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3) =
2607 : !! inverse of the overlap matrix
2608 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
2609 : !! pwindall(:,1,:) <- <u^(0)_i|u^(0)_i+1>
2610 : !! pwindall(:,2,:) <- <u^(0)_i|u^(0)_i-1>
2611 : !! pwindall(:,3,:) <- <u^(1)_i|u^(1)_i+1>
2612 : !! pwindall(:,4,:) <- <u^(1)_i|u^(1)_i-1>
2613 : !! pwindall(:,5,:) <- <u^(1)_i|u^(0)_i+n+1>
2614 : !! pwindall(:,6,:) <- <u^(1)_i|u^(0)_i+n-1>
2615 : !! pwindall(:,7,:) <- <u^(0)_i|u^(1)_i-n+1>
2616 : !! pwindall(:,8,:) <- <u^(0)_i|u^(1)_i-n-1>
2617 : !! rprimd(3,3)=dimensional primitive translations in real space (bohr)
2618 : !!
2619 : !! OUTPUT
2620 : !! d2lo(1,1:3,natom+5,1:3,1:natom) = Born effective charge tensor
2621 : !!
2622 : !! SOURCE
2623 :
2624 6 : subroutine dfptff_bec(cg,cg1,dtefield,natom,d2lo,idirpert,ipert,mband,mband_mem,mkmem_rbz,&
2625 6 : & mpi_enreg,mpw,mpw1,mpert,nkpt,npwarr,npwar1,nsppol,nspinor,pwindall,qmat,rprimd)
2626 :
2627 : !Arguments ----------------------------------------
2628 : !scalars
2629 : integer,intent(in) :: idirpert,ipert,mband,mkmem_rbz,mpert,mpw,mpw1,natom,nkpt
2630 : integer,intent(in) :: mband_mem
2631 : integer,intent(in) :: nspinor,nsppol
2632 : type(efield_type),intent(in) :: dtefield
2633 : type(MPI_type),intent(in) :: mpi_enreg
2634 : !arrays
2635 : integer,intent(in) :: npwar1(nkpt),npwarr(nkpt)
2636 : integer,intent(in) :: pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
2637 : real(dp),intent(in) :: cg(2,mpw*mband_mem*mkmem_rbz*nspinor*nsppol)
2638 : real(dp),intent(in) :: cg1(2,mpw1*mband_mem*mkmem_rbz*nspinor*nsppol)
2639 : real(dp),intent(in) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2640 : real(dp),intent(in) :: rprimd(3,3)
2641 : real(dp),intent(inout) :: d2lo(2,3,mpert,3,mpert) !vz_i
2642 :
2643 : !Local variables ----------------------------------
2644 : !scalars
2645 : integer :: ialpha,iband,icg,icg1,idir,ikpt,ikpt1,jband,mpw_tmp,npw_k1
2646 : integer :: npw_k2,pwmax,pwmin
2647 : integer :: me_band,ierr, iband_me, jband_me
2648 : real(dp) :: doti,dotr,e0,fac
2649 : !arrays
2650 6 : integer,allocatable :: pwind_tmp(:)
2651 12 : integer :: band_procs(mband)
2652 : real(dp) :: edir(3)
2653 6 : real(dp),allocatable :: s1mat(:,:,:),vect1(:,:),vect2(:,:)
2654 :
2655 : ! *************************************************************************
2656 :
2657 : !calculate s1 matrices -----------------------------
2658 6 : mpw_tmp=max(mpw,mpw1)
2659 24 : ABI_MALLOC(s1mat,(2,dtefield%mband_occ,dtefield%mband_occ))
2660 18 : ABI_MALLOC(vect1,(2,0:mpw_tmp))
2661 12 : ABI_MALLOC(vect2,(2,0:mpw_tmp))
2662 18 : ABI_MALLOC(pwind_tmp,(mpw_tmp))
2663 36 : vect1(:,0) = zero ; vect2(:,0) = zero
2664 :
2665 6 : me_band = mpi_enreg%me_band
2666 :
2667 6 : edir(:)=zero
2668 :
2669 : !TODO MJV: where is the loop over sppol? proc_distrb_band_procs chooses isppol 1 below
2670 54 : do ikpt=1,nkpt
2671 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,1,mband,&
2672 48 : & me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
2673 :
2674 198 : do idir=1,3
2675 :
2676 : ! compute <u^(0)_{k_j}|u^(1)_{k_j+1,q}> matrix--- q=0 ----------------------------------------
2677 :
2678 : ! prepare to calculate overlap matrix
2679 144 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2680 144 : icg = dtefield%cgindex(ikpt,1)
2681 144 : icg1 = dtefield%cgindex(ikpt1,1+nsppol)
2682 144 : npw_k1 = npwarr(ikpt)
2683 144 : npw_k2 = npwar1(ikpt1)
2684 6120 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,7,idir)
2685 720 : vect1(:,0) = zero ; vect2(:,0) = zero
2686 144 : jband_me = 0
2687 7632 : s1mat = zero
2688 720 : do jband = 1, dtefield%mband_occ
2689 576 : if (band_procs(jband) == me_band) then
2690 576 : jband_me = jband_me + 1
2691 : vect2(:,1:npw_k2) = &
2692 72288 : & cg1(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2693 : end if
2694 576 : call xmpi_bcast(vect2,band_procs(jband), mpi_enreg%comm_band,ierr)
2695 :
2696 : ! now everyone has vect2 for present jband
2697 3168 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2698 :
2699 : iband_me = 0
2700 3600 : do iband = 1, dtefield%mband_occ
2701 2304 : if (band_procs(iband) /= me_band) cycle
2702 2304 : iband_me = iband_me + 1
2703 :
2704 2304 : pwmin = (iband_me-1)*npw_k1*nspinor
2705 2304 : pwmax = pwmin + npw_k1*nspinor
2706 : vect1(:,1:npw_k1) = &
2707 289152 : & cg(:,icg + 1 + pwmin:icg + pwmax)
2708 12672 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2709 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2710 2304 : & vect1,vect2)
2711 2304 : s1mat(1,iband,jband) = dotr
2712 2880 : s1mat(2,iband,jband) = doti
2713 : end do ! iband
2714 : end do !jband
2715 : ! for the moment s1mat is only filled for iband on current cpu
2716 :
2717 :
2718 : ! compute <u^(1)_{k_j,q}|u^(0)_{k_j+1}> matrix-- q=0 -------------------------------------
2719 :
2720 : ! prepare to calculate overlap matrix
2721 144 : ikpt1 = dtefield%ikpt_dk(ikpt,1,idir)
2722 144 : icg = dtefield%cgindex(ikpt,1+nsppol)
2723 144 : icg1 = dtefield%cgindex(ikpt1,1)
2724 144 : npw_k1 = npwar1(ikpt)
2725 144 : npw_k2 = npwarr(ikpt1)
2726 6120 : pwind_tmp(1:npw_k1) = pwindall((ikpt-1)*mpw_tmp+1:(ikpt-1)*mpw_tmp+npw_k1,5,idir)
2727 720 : vect1(:,0) = zero ; vect2(:,0) = zero
2728 : jband_me = 0
2729 720 : do jband = 1, dtefield%mband_occ
2730 576 : if (band_procs(jband) == me_band) then
2731 576 : jband_me = jband_me + 1
2732 : vect2(:,1:npw_k2) = &
2733 72288 : & cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2734 : end if
2735 576 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2736 :
2737 3168 : if (npw_k2 < mpw_tmp) vect2(:,npw_k2+1:mpw_tmp) = zero
2738 : iband_me = 0
2739 3600 : do iband = 1, dtefield%mband_occ
2740 2304 : if (band_procs(iband) /= me_band) cycle
2741 2304 : iband_me = iband_me + 1
2742 :
2743 2304 : pwmin = (iband_me-1)*npw_k1*nspinor
2744 2304 : pwmax = pwmin + npw_k1*nspinor
2745 : vect1(:,1:npw_k1) = &
2746 289152 : & cg1(:,icg + 1 + pwmin:icg + pwmax)
2747 12672 : if (npw_k1 < mpw_tmp) vect1(:,npw_k1+1:mpw_tmp) = zero
2748 : call overlap_g(doti,dotr,mpw_tmp,npw_k1,npw_k2,nspinor,pwind_tmp,&
2749 2304 : & vect1,vect2)
2750 :
2751 2304 : s1mat(1,iband,jband) = s1mat(1,iband,jband) + dotr
2752 2880 : s1mat(2,iband,jband) = s1mat(2,iband,jband) + doti
2753 : end do ! iband
2754 : end do !jband
2755 :
2756 : ! recompose full s1mat on each proc
2757 144 : call xmpi_sum(s1mat,mpi_enreg%comm_band,ierr)
2758 :
2759 : ! sum over the whole------------------------------------------------------------
2760 :
2761 144 : e0=zero
2762 :
2763 720 : do iband=1,dtefield%mband_occ
2764 3024 : do jband=1,dtefield%mband_occ
2765 : e0 = e0 + (s1mat(1,iband,jband)*qmat(2,jband,iband,ikpt,1,idir)&
2766 2880 : & + s1mat(2,iband,jband)*qmat(1,jband,iband,ikpt,1,idir))
2767 : end do
2768 : end do
2769 :
2770 768 : do ialpha=1,3
2771 : fac = rprimd(ialpha,idir)/&
2772 432 : & (dble(dtefield%nstr(idir))*pi)
2773 :
2774 576 : edir(ialpha)=edir(ialpha)+ e0*fac
2775 : end do
2776 :
2777 : end do ! idir
2778 : end do ! ikpt
2779 :
2780 24 : d2lo(1,1:3,natom+2,idirpert,ipert)=edir(:)
2781 :
2782 6 : ABI_FREE(s1mat)
2783 6 : ABI_FREE(vect1)
2784 6 : ABI_FREE(vect2)
2785 6 : ABI_FREE(pwind_tmp)
2786 :
2787 6 : end subroutine dfptff_bec
2788 : !!***
2789 :
2790 : !!****f* ABINIT/qmatrix
2791 : !! NAME
2792 : !! qmatrix
2793 : !!
2794 : !! FUNCTION
2795 : !! calculation of the inverse of the overlap matrix
2796 : !!
2797 : !! INPUTS
2798 : !! cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol) = planewave coefficients of wavefunctions
2799 : !! RF wavefunctions at k,q.
2800 : !! dtefield = variables related to response Berry-phase calculation
2801 : !! ikpt = the index of the current k point
2802 : !! mband = maximum number of bands
2803 : !! mband_mem = maximum number of bands on this cpu
2804 : !! mkmem_rbz = maximum number of k-points in core memory
2805 : !! mpw = maximum number of plane waves
2806 : !! mpw1 = maximum number of plane waves for response wavefunctions
2807 : !! nkpt = number of k points
2808 : !! npwarr(nkpt) = number of planewaves in basis and boundary at this k point
2809 : !! npwar1(nkpt) = number of planewaves in basis and boundary for response wfs
2810 : !! nspinor = 1 for scalar wfs, 2 for spinor wfs
2811 : !! nsppol = 1 for unpolarized, 2 for spin-polarized
2812 : !! pwindall(max(mpw,mpw1)*mkmem_rbz,8,3) = array used to compute the overlap matrices
2813 : !!
2814 : !! OUTPUT
2815 : !! qmat(2,dtefield%nband_occ,dtefield%nband_occ,nkpt,2,3) = inverse of the overlap matrix
2816 : !!
2817 : !! SOURCE
2818 :
2819 9 : subroutine qmatrix(cg,dtefield,qmat,mpi_enreg,mpw,mpw1,mkmem_rbz,mband,mband_mem,npwarr,nkpt,nspinor,nsppol,pwindall)
2820 :
2821 : use m_hide_lapack, only : dzgedi, dzgefa
2822 :
2823 : !Arguments ----------------------------------------
2824 : !scalars
2825 : integer,intent(in) :: mband,mkmem_rbz,mpw,mpw1,nkpt,nspinor,nsppol
2826 : integer,intent(in) :: mband_mem
2827 : type(efield_type),intent(in) :: dtefield
2828 : type(MPI_type),intent(in) :: mpi_enreg
2829 : !arrays
2830 : integer,intent(in) :: npwarr(nkpt),pwindall(max(mpw,mpw1)*mkmem_rbz,8,3)
2831 : real(dp),intent(in) :: cg(2,mpw*nspinor*mband_mem*mkmem_rbz*nsppol)
2832 : real(dp),intent(out) :: qmat(2,dtefield%mband_occ,dtefield%mband_occ,nkpt,2,3)
2833 :
2834 : !Local variables -------------------------
2835 : !scalars
2836 : integer :: iband,icg,icg1,idir,ifor,ikpt,ikpt2,info,jband,job
2837 : integer :: iband_me, jband_me, me_band, ierr
2838 : integer :: npw_k1,npw_k2,pwmax,pwmin
2839 : integer :: isppol
2840 : real(dp) :: doti,dotr
2841 : !arrays
2842 9 : integer,allocatable :: ipvt(:),pwind_k(:)
2843 : real(dp) :: det(2,2)
2844 18 : integer :: band_procs(mband)
2845 9 : real(dp),allocatable :: sinv(:,:,:),smat_k(:,:,:),vect1(:,:),vect2(:,:)
2846 9 : real(dp),allocatable :: zgwork(:,:)
2847 :
2848 : ! *************************************************************************
2849 :
2850 27 : ABI_MALLOC(ipvt,(dtefield%mband_occ))
2851 36 : ABI_MALLOC(sinv,(2,dtefield%mband_occ,dtefield%mband_occ))
2852 27 : ABI_MALLOC(zgwork,(2,dtefield%mband_occ))
2853 27 : ABI_MALLOC(vect1,(2,0:mpw))
2854 18 : ABI_MALLOC(vect2,(2,0:mpw))
2855 27 : ABI_MALLOC(smat_k,(2,dtefield%mband_occ,dtefield%mband_occ))
2856 27 : ABI_MALLOC(pwind_k,(max(mpw,mpw1)))
2857 54 : vect1(:,0) = zero ; vect2(:,0) = zero
2858 :
2859 9 : job = 11
2860 :
2861 9 : me_band = mpi_enreg%me_band
2862 :
2863 : !loop over k points
2864 18 : do isppol = 1, nsppol
2865 90 : do ikpt = 1, nkpt
2866 72 : call proc_distrb_band(band_procs,mpi_enreg%proc_distrb,ikpt,isppol,mband,me_band,mpi_enreg%me_kpt,mpi_enreg%comm_band)
2867 :
2868 72 : npw_k1 = npwarr(ikpt)
2869 72 : icg = dtefield%cgindex(ikpt,1)
2870 297 : do idir = 1, 3
2871 720 : do ifor = 1, 2
2872 :
2873 432 : ikpt2 = dtefield%ikpt_dk(ikpt,ifor,idir)
2874 432 : npw_k2 = npwarr(ikpt2)
2875 432 : icg1 = dtefield%cgindex(ikpt2,1)
2876 18360 : pwind_k(1:npw_k1) = pwindall((ikpt-1)*max(mpw,mpw1)+1:(ikpt-1)*max(mpw,mpw1)+npw_k1,ifor,idir)
2877 :
2878 22896 : smat_k = zero
2879 432 : jband_me = 0
2880 2160 : do jband = 1, dtefield%nband_occ(isppol)
2881 1728 : if (band_procs(jband) == me_band) then
2882 1728 : jband_me = jband_me + 1
2883 216864 : vect2(:,1:npw_k2) = cg(:,icg1 + 1 + (jband_me-1)*npw_k2*nspinor:icg1 + jband_me*npw_k2*nspinor)
2884 : end if
2885 1728 : call xmpi_bcast(vect2,band_procs(jband),mpi_enreg%comm_band,ierr)
2886 :
2887 9504 : if (npw_k2 < mpw) vect2(:,npw_k2+1:mpw) = zero
2888 :
2889 1728 : iband_me = 0
2890 10800 : do iband = 1, dtefield%nband_occ(isppol)
2891 6912 : if (band_procs(iband) /= me_band) cycle
2892 6912 : iband_me = iband_me + 1
2893 :
2894 6912 : pwmin = (iband_me-1)*npw_k1*nspinor
2895 6912 : pwmax = pwmin + npw_k1*nspinor
2896 867456 : vect1(:,1:npw_k1) = cg(:,icg + 1 + pwmin:icg + pwmax)
2897 38016 : if (npw_k1 < mpw) vect1(:,npw_k1+1:mpw) = zero
2898 6912 : call overlap_g(doti,dotr,mpw,npw_k1,npw_k2,nspinor,pwind_k,vect1,vect2)
2899 6912 : smat_k(1,iband,jband) = dotr
2900 8640 : smat_k(2,iband,jband) = doti
2901 : end do ! iband
2902 : end do !jband
2903 :
2904 : ! recompose full s1mat on each proc
2905 432 : call xmpi_sum(smat_k,mpi_enreg%comm_band,ierr)
2906 :
2907 22896 : sinv(:,:,:) = smat_k(:,:,:)
2908 :
2909 432 : call dzgefa(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,info)
2910 432 : call dzgedi(sinv,dtefield%mband_occ,dtefield%nband_occ(isppol),ipvt,det,zgwork,job)
2911 :
2912 23544 : qmat(:,:,:,ikpt,ifor,idir) = sinv(:,:,:)
2913 : end do
2914 : end do
2915 : end do !end loop over k
2916 : end do
2917 :
2918 9 : ABI_FREE(ipvt)
2919 9 : ABI_FREE(sinv)
2920 9 : ABI_FREE(zgwork)
2921 9 : ABI_FREE(vect1)
2922 9 : ABI_FREE(vect2)
2923 9 : ABI_FREE(smat_k)
2924 9 : ABI_FREE(pwind_k)
2925 :
2926 9 : end subroutine qmatrix
2927 : !!***
2928 :
2929 : end module m_dfpt_fef
2930 : !!***
|
