Branch data 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 : : !!***
|
