Branch data Line data Source code
1 : : !!****m* ABINIT/m_elpolariz
2 : : !! NAME
3 : : !! m_elpolariz
4 : : !!
5 : : !! FUNCTION
6 : : !!
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2005-2026 ABINIT group (XG, NSAI, MKV)
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_elpolariz
23 : :
24 : : use defs_basis
25 : : use m_abicore
26 : : use m_errors
27 : : use m_efield
28 : : use m_xmpi
29 : : use m_wffile
30 : : use m_hdr
31 : : use m_dtset
32 : : use m_dtfil
33 : :
34 : : use defs_datatypes, only : pseudopotential_type
35 : : use defs_abitypes, only : MPI_type
36 : : use m_geometry, only : metric
37 : : use m_matrix, only : matr3inv
38 : : use m_hide_lapack, only : dzgedi, dzgefa
39 : : use m_rwwf, only : rwwf
40 : : use m_pawtab, only : pawtab_type
41 : : use m_pawrhoij, only : pawrhoij_type
42 : : use m_pawcprj, only : pawcprj_type
43 : : use m_berryphase, only : berryphase
44 : : use m_berryphase_new, only : berryphase_new
45 : :
46 : : implicit none
47 : :
48 : : private
49 : : !!***
50 : :
51 : : public :: elpolariz
52 : : !!***
53 : :
54 : : contains
55 : : !!***
56 : :
57 : : !!****f* ABINIT/elpolariz
58 : : !! NAME
59 : : !! elpolariz
60 : : !!
61 : : !! FUNCTION
62 : : !! Calculate corrections to total energy from polarising
63 : : !! electric field with or without Berry phases (berryopt keyword)
64 : : !!
65 : : !! INPUTS
66 : : !! atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
67 : : !! cg(2,mcg)=planewave coefficients of wavefunctions
68 : : !! cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk>
69 : : !! and each |p_lmn> non-local projector
70 : : !! dtfil <type(datafiles_type)>=variables related to files
71 : : !! dtset <type(dataset_type)>=all input variables in this dataset
72 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
73 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
74 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates
75 : : !! mband=maximum number of bands
76 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
77 : : !! mkmem=number of k points treated by this node.
78 : : !! mpi_enreg=information about MPI parallelization
79 : : !! mpw=maximum dimensioned size of npw
80 : : !! my_natom=number of atoms treated by current processor
81 : : !! natom=number of atoms in cell
82 : : !! nattyp(ntypat)= # atoms of each type.
83 : : !! nkpt=number of k points
84 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
85 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
86 : : !! ntypat=number of types of atoms in unit cell
87 : : !! nkpt=number of k-points
88 : : !! option = 1: compute Berryphase polarization
89 : : !! 2: compute finite difference expression of the ddk
90 : : !! 3: compute polarization & ddk
91 : : !! pawrhoij(my_natom*usepaw) <type(pawrhoij_type)> atomic occupancies
92 : : !! pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
93 : : !! pel_cg(3) = reduced coordinates of the electronic polarization (a. u.)
94 : : !! computed in the SCF loop
95 : : !! pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates
96 : : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
97 : : !! pwind(pwind_alloc,2,3) = array used to compute
98 : : !! the overlap matrix smat between k-points (see initberry.f)
99 : : !! pwind_alloc = first dimension of pwind
100 : : !! pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
101 : : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
102 : : !! ucvol=unit cell volume in bohr**3.
103 : : !! usecprj=1 if cprj datastructure has been allocated
104 : : !! xred(3,natom)=reduced atomic coordinates
105 : : !!
106 : : !! OUTPUT
107 : : !! (see side effects)
108 : : !!
109 : : !! SIDE EFFECTS
110 : : !! dtefield <type(efield_type)> = variables related to Berry phase
111 : : !! and electric field calculations (see initberry.f).
112 : : !! In case berryopt = 4/6/7/14/16/17, the overlap matrices computed
113 : : !! in this routine are stored in dtefield%smat in order
114 : : !! to be used in the electric field calculation.
115 : : !! enefield=field energy
116 : : !! etotal=total energy, might be correct by improved polarization computation
117 : : !! pel(3) = reduced coordinates of the electronic polarization (a. u.)
118 : : !! pion(3)= reduced coordinates of the ionic polarization (a. u.)
119 : : !!
120 : : !! SOURCE
121 : :
122 : 85 : subroutine elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,enefield,gprimd,hdr,&
123 : 85 : & kg,mband,mcg,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,nkpt,&
124 : 85 : & npwarr,nsppol,ntypat,pawrhoij,pawtab,&
125 : 85 : & pel,pel_cg,pelev,pion,psps,pwind,pwind_alloc,&
126 : 85 : & pwnsfac,rprimd,ucvol,usecprj,xred)
127 : :
128 : : !Arguments ------------------------------------
129 : : !scalars
130 : : integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,my_natom,natom,nkpt,nsppol,ntypat
131 : : integer,intent(in) :: pwind_alloc,usecprj
132 : : real(dp),intent(in) :: ucvol
133 : : real(dp),intent(inout) :: enefield,etotal
134 : : type(MPI_type),intent(in) :: mpi_enreg
135 : : type(datafiles_type),intent(in) :: dtfil
136 : : type(dataset_type),intent(inout) :: dtset
137 : : type(efield_type),intent(inout) :: dtefield
138 : : type(hdr_type),intent(inout) :: hdr
139 : : type(pseudopotential_type),intent(in) :: psps
140 : : !arrays
141 : : integer,intent(in) :: atindx1(natom),kg(3,mpw*mkmem),nattyp(ntypat)
142 : : integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
143 : : real(dp),intent(in) :: cg(2,mcg),gprimd(3,3)
144 : : real(dp),intent(in) :: pel_cg(3),pwnsfac(2,pwind_alloc),rprimd(3,3)
145 : : real(dp),intent(inout) :: pel(3),pelev(3),pion(3),xred(3,natom)
146 : : type(pawcprj_type),intent(in) :: cprj(natom,mcprj*usecprj)
147 : : type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*psps%usepaw)
148 : : type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
149 : :
150 : : !Local variables-------------------------------
151 : : !scalars
152 : : integer :: mcg1_3,my_nspinor,option,unit_out,iir,jjr,kkr
153 : : real(dp) :: pdif_mod,eenth,ucvol_local
154 : : logical :: save_cg1_3
155 : : character(len=500) :: message
156 : : !arrays
157 : : real(dp) :: gmet(3,3),gprimdlc(3,3),pdif(3),ptot(3),red_ptot(3),rmet(3,3)
158 : : !! ptot(3) = total polarization (not reduced) REC
159 : : !! red_ptot(3) = internal reduced total polarization REC
160 : : real(dp) :: A(3,3),A1(3,3),A_new(3,3),efield_new(3)
161 : 85 : real(dp),allocatable :: cg1_3(:,:,:)
162 : :
163 : : ! *************************************************************************
164 : :
165 : : DBG_ENTER("COLL")
166 : :
167 [ + + - + ]: 85 : if (usecprj==0.and.psps%usepaw==1) then
168 : : write (message,'(3a)')&
169 : 0 : & 'cprj datastructure must be allocated !',ch10,&
170 : 0 : & 'Action: change pawusecp input keyword.'
171 : 0 : ABI_ERROR(message)
172 : : end if
173 : :
174 : 85 : my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
175 : :
176 : : if(dtset%berryopt>0 .and. dtset%berryopt/=4 .and. dtset%berryopt/=6 .and. dtset%berryopt/=7 .and. &
177 [ - + ]: 85 : & dtset%berryopt/=14 .and. dtset%berryopt/=16 .and. dtset%berryopt/=17)then !!HONG
178 : :
179 : : if (dtset%berryopt==1 .or. dtset%berryopt==3) then
180 : : call berryphase(atindx1,dtset%bdberry,cg,gprimd,dtset%istwfk,&
181 : : & dtset%kberry,kg,dtset%kptns,dtset%kptopt,dtset%kptrlatt,&
182 : : & mband,mcg,mkmem,mpw,natom,nattyp,dtset%nband,dtset%nberry,npwarr,&
183 : : & my_nspinor,nsppol,psps%ntypat,nkpt,rprimd,ucvol,&
184 : 0 : & xred,psps%ziontypat)
185 : : end if
186 : :
187 [ # # ]: 0 : if (dtset%berryopt==2 .or. dtset%berryopt==3) then
188 : : call uderiv(dtset%bdberry,cg,gprimd,hdr,dtset%istwfk,&
189 : : & dtset%kberry,kg,dtset%kptns,dtset%kptopt,&
190 : : & dtset%kptrlatt,mband,mcg,mkmem,mpi_enreg,mpw,&
191 : : & natom,dtset%nband,dtset%nberry,npwarr,my_nspinor,nsppol,&
192 : 0 : & nkpt,dtfil%unddk,dtfil%fnameabo_1wf)
193 : : end if
194 : :
195 : : else if(dtset%berryopt<0 .or. dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or. &
196 [ - + ]: 85 : & dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)then !!HONG
197 : :
198 : 0 : select case (dtset%berryopt)
199 : : case (-5)
200 : 0 : option = 2
201 : : case (-3)
202 : 1 : option = 3
203 : : case (-2)
204 : 3 : option = 2
205 : : case (-1)
206 : 56 : option = 1
207 : : case (4)
208 : 25 : option = 1
209 : 25 : pel(:) = zero
210 : 25 : pelev(:) = zero
211 : : case (6) !!HONG
212 : 0 : option = 1
213 : 0 : pel(:) = zero
214 : 0 : pelev(:) = zero
215 : : case (7) !!HONG
216 : 0 : option = 1
217 : 0 : pel(:) = zero
218 : 0 : pelev(:) = zero
219 : : case (14) !!HONG
220 : 0 : option = 1
221 : 0 : pel(:) = zero
222 : 0 : pelev(:) = zero
223 : : case (16) !!HONG
224 : 0 : option = 1
225 : 0 : pel(:) = zero
226 : 0 : pelev(:) = zero
227 : : case (17) !!HONG
228 : 0 : option = 1
229 : 0 : pel(:) = zero
230 [ - + + + : 85 : pelev(:) = zero
+ - - - -
- - ]
231 : : end select
232 : :
233 : 85 : unit_out = ab_out
234 : :
235 : : ! by default berryphase_new writes ddk wavefunctions to disk
236 : : ! with save_cg1_3 set to true it can return them in memory
237 : : ! this feature is used at present when orbmag is also called,
238 : : ! as in afterscfloop
239 : 85 : save_cg1_3 = .FALSE.
240 : 85 : mcg1_3 = 0
241 [ - + ]: 85 : ABI_MALLOC(cg1_3,(2,mcg1_3,3))
242 : : call berryphase_new(atindx1,cg,cg1_3,cprj,dtefield,dtfil,dtset,psps,&
243 : : & gprimd,hdr,psps%indlmn,kg,&
244 : : & psps%lmnmax,mband,mcg,mcg1_3,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,&
245 : : & nsppol,psps%ntypat,nkpt,option,pawrhoij,&
246 : : & pawtab,pel,pelev,pion,ptot,red_ptot,pwind,& !!REC
247 : : & pwind_alloc,pwnsfac,rprimd,save_cg1_3,dtset%typat,ucvol,&
248 : 85 : & unit_out,usecprj,psps%usepaw,xred,psps%ziontypat)
249 [ - + ]: 85 : ABI_FREE(cg1_3)
250 : :
251 [ + + ]: 340 : dtefield%red_ptot1(:)=red_ptot(:)
252 : :
253 : : if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or. &
254 [ + + ]: 85 : & dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt == 17 ) then !!HONG
255 : :
256 : : ! Check if pel has the same value as pel_cg
257 : : ! if (psps%usepaw == 1) pel(:) = pel(:) + pelev(:) ! add on-site term for PAW
258 : : ! if (psps%usepaw == 1) red_ptot(:) = red_ptot(:) + pelev(:) ! add on-site term for PAW !! REC
259 : : ! above line suppressed because in the PAW case, pel already includes all on-site
260 : : ! terms and pelev should not be added in additionally. We are computing pelev separately for
261 : : ! reporting purposes only.
262 : : ! 13 June 2012 J Zwanziger
263 : :
264 [ + + ]: 100 : pdif(:) = pel_cg(:) - pel(:)
265 : 25 : pdif_mod = pdif(1)**2 + pdif(2)**2 + pdif(3)**2
266 : :
267 [ - + ]: 25 : if (pdif_mod > tol8) then
268 : 0 : write(message,'(11(a),e16.9)')ch10,&
269 : 0 : & ' scfcv (electric field calculation) : WARNING -',ch10,&
270 : 0 : & ' The difference between pel (electronic Berry phase updated ',ch10,&
271 : 0 : & ' at each SCF cycle)',ch10,&
272 : 0 : & ' and pel_cg (electronic Berryphase computed using the ',&
273 : 0 : & 'berryphase routine) is',ch10,&
274 : 0 : & ' pdif_mod = ',pdif_mod
275 : 0 : call wrtout(std_out,message,'COLL')
276 : 0 : write(message,'(a,6(a,e16.9,a))') ch10,&
277 : 0 : & 'pel_cg(1) = ',pel_cg(1),ch10,&
278 : 0 : & 'pel_cg(2) = ',pel_cg(2),ch10,&
279 : 0 : & 'pel_cg(3) = ',pel_cg(3),ch10,&
280 : 0 : & 'pel(1) = ',pel(1),ch10,&
281 : 0 : & 'pel(2) = ',pel(2),ch10,&
282 : 0 : & 'pel(3) = ',pel(3),ch10
283 : 0 : ABI_ERROR(message)
284 : : end if
285 : :
286 : : ! Use this (more accurate) value of P to recompute enefield
287 [ + - ]: 25 : if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then !!HONG
288 : 25 : etotal = etotal - enefield
289 : :
290 [ + + ]: 100 : enefield = -dot_product(dtset%red_efieldbar,red_ptot)
291 : 25 : call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
292 : 25 : eenth = zero
293 [ + + ]: 100 : do iir=1,3
294 [ + + ]: 325 : do jjr=1,3
295 : 300 : eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j
296 : : end do
297 : : end do
298 : 25 : eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth
299 : 25 : enefield=enefield+eenth
300 : :
301 : 25 : etotal = etotal + enefield
302 : :
303 : 25 : write(message,'(a,a)')ch10,&
304 : 50 : & ' Stress tensor under a constant electric field:'
305 : 25 : call wrtout(std_out,message,'COLL')
306 : 25 : call wrtout(ab_out,message,'COLL')
307 : :
308 : : end if
309 : :
310 : : ! ! In finite D-field case, turn it into internal energy !!HONG
311 [ - + ]: 25 : if (dtset%berryopt == 6 .or. dtset%berryopt == 16 ) then
312 : 0 : etotal = etotal - enefield
313 : :
314 : 0 : enefield=zero
315 : 0 : call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
316 [ # # ]: 0 : do iir=1,3
317 [ # # ]: 0 : do jjr=1,3
318 : 0 : enefield= enefield+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr) !! HONG g^{-1})_ij ebar_i ebar_j
319 : : end do
320 : : end do
321 : 0 : enefield= ucvol_local/(8.0d0*pi)*enefield
322 : :
323 : 0 : etotal = etotal + enefield
324 : :
325 : 0 : write(message,'(a,a)')ch10,&
326 : 0 : & ' Stress tensor under a constant electric displacement field:'
327 : 0 : call wrtout(std_out,message,'COLL')
328 : 0 : call wrtout(ab_out,message,'COLL')
329 : :
330 : : end if
331 : :
332 : : ! HONG calculate internal energy and electric enthalpy for mixed BC case.
333 [ - + ]: 25 : if ( dtset%berryopt == 17 ) then
334 : 0 : etotal = etotal - enefield
335 : 0 : enefield = zero
336 : :
337 : 0 : call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
338 [ # # # # ]: 0 : A(:,:)=(4*pi/ucvol_local)*rmet(:,:)
339 : 0 : A1(:,:)=A(:,:)
340 : 0 : A_new(:,:)=A(:,:)
341 [ # # ]: 0 : efield_new(:)=dtset%red_efield(:)
342 : : eenth = zero
343 : :
344 [ # # ]: 0 : do kkr=1,3
345 [ # # ]: 0 : if (dtset%jfielddir(kkr)==1) then ! fixed ebar direction
346 : : ! step 1 add -ebar*p
347 : 0 : eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr)
348 : :
349 : : ! step 2 chang to e_new (change e to ebar)
350 : 0 : efield_new(kkr)=dtset%red_efieldbar(kkr)
351 : :
352 : : ! step 3 chang matrix A to A1
353 : :
354 [ # # ]: 0 : do iir=1,3
355 [ # # ]: 0 : do jjr=1,3
356 [ # # ]: 0 : if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
357 [ # # ]: 0 : if ((iir==kkr .and. jjr/=kkr) .or. (iir/=kkr .and. jjr==kkr)) &
358 : 0 : & A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
359 [ # # ]: 0 : if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
360 : : end do
361 : : end do
362 : :
363 : 0 : A(:,:)=A1(:,:)
364 : 0 : A_new(:,:)=A1(:,:)
365 : : end if
366 : :
367 : : end do ! end fo kkr
368 : :
369 : :
370 [ # # ]: 0 : do iir=1,3
371 [ # # ]: 0 : do jjr=1,3
372 : 0 : eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
373 : : end do
374 : : end do
375 : :
376 : 0 : enefield=eenth
377 : 0 : etotal = etotal + enefield
378 : :
379 : 0 : write(message,'(a,a)')ch10,&
380 : 0 : & ' Stress tensor under a constant (mixed) electric and electric displacement field:'
381 : 0 : call wrtout(std_out,message,'COLL')
382 : 0 : call wrtout(ab_out,message,'COLL')
383 : :
384 : : end if ! berryopt==17
385 : :
386 : :
387 : : ! MVeithen: to clarify
388 : : ! Which stress tensor should be used in structural optimizations?
389 : : ! The one at constant electric field or at constant potential drop.
390 : : ! write(message,'(a,a)')ch10,&
391 : : ! & ' Stress tensor imposing a constant electric field:'
392 : : ! call wrtout(std_out,message,'COLL')
393 : : ! call wrtout(ab_out,message,'COLL')
394 : :
395 : : end if ! dtset%berryopt == 4/6/7/14/16/17
396 : :
397 : : end if ! dtset%berryopt>0 or dtset%berryopt/=4/6/7/14/16/17
398 : :
399 : : DBG_EXIT("COLL")
400 : :
401 : 85 : end subroutine elpolariz
402 : : !!***
403 : :
404 : : !!****f* ABINIT/uderiv
405 : : !! NAME
406 : : !! uderiv
407 : : !!
408 : : !! FUNCTION
409 : : !! This routine is called computes the derivative of
410 : : !! ground-state wavefunctions with respect to k (du/dk) by finite differencing
411 : : !! on neighbouring k points
412 : : !! Work for nsppol=1 or 2, but only accept nspinor=1,
413 : : !!
414 : : !! INPUTS
415 : : !! bdberry(4)=band limits for Berry phase contributions (or du/dk)
416 : : !! spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1)
417 : : !! cg(2,mcg)=planewave coefficients of wavefunctions
418 : : !! gprimd(3,3)=reciprocal space dimensional primitive translations
419 : : !! hdr <type(hdr_type)>=the header of wf, den and pot files
420 : : !! istwfk(nkpt_)=input option parameter that describes the storage of wfs
421 : : !! kberry(3,20)= different delta k for Berry phases(or du/dk),
422 : : !! in unit of kptrlatt only kberry(1:3,1:nberry) is relevant
423 : : !! kg(3,mpw*mkmem)=reduced planewave coordinates
424 : : !! kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT,
425 : : !! kpt_ samples half the BZ if time-reversal symetrie is used
426 : : !! kptopt=2 when time-reversal symmetry is used
427 : : !! kptrlatt(3,3)=k-point lattice specification
428 : : !! mband=maximum number of bands
429 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
430 : : !! mkmem=number of k points treated by this node.
431 : : !! mpi_enreg=information about MPI parallelization
432 : : !! mpw=maximum dimensioned size of npw
433 : : !! natom=number of atoms in cell
434 : : !! nband(nkpt*nsppol)=number of bands at each k point, for each polarization
435 : : !! nberry=number of Berry phases(or du/dk) to be computed
436 : : !! nkpt=number of k points
437 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
438 : : !! nspinor=number of spinorial components of the wavefunctions (on current proc)
439 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
440 : : !! unddk=unit number for ddk file
441 : : !!
442 : : !! OUTPUT
443 : : !! (the ddk wavefunctions are written on disk)
444 : : !!
445 : : !! SIDE EFFECTS
446 : : !!
447 : : !! TODO
448 : : !! Cleaning, checking for rules
449 : : !! Should allow for time-reversal symmetry (istwfk)
450 : : !! WARNING : the use of nspinor is completely erroneous
451 : : !!
452 : : !! NOTES
453 : : !! Local Variables:
454 : : !! cmatrix(:,:,:)= overlap matrix of size maxband*maxband
455 : : !! cg_index(:,:,:)= unpacked cg index array for specific band,
456 : : !! k point and polarization.
457 : : !! det(2,2)= intermediate output of Lapack routine zgedi.f
458 : : !! dk(3)= step taken to the next k mesh point along the kberry direction
459 : : !! gpard(3)= dimensionalreciprocal lattice vector G along which the
460 : : !! polarization is computed
461 : : !! kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
462 : : !! planewave and k point
463 : : !! kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
464 : : !! kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT
465 : : !! code
466 : : !! kpt_flag(ikpt) gives the indices of the k-point
467 : : !! related to ikpt by time revers
468 : : !! maxband/minband= control the minimum and maximum band calculated in the
469 : : !! overlap matrix
470 : : !! npw_k= npwarr(ikpt), number of planewaves in basis at this k point
471 : : !! shift_g_2(nkpt,nkpt)= .true. if the k point should be shifted by a G vector;
472 : : !! .false. if not
473 : : !! tr(2)=variable that changes k to -k
474 : : !! G to -G
475 : : !! $c_g$ to $c_g^*$ when time-reversal symetrie is used
476 : : !!
477 : : !! SOURCE
478 : :
479 : 0 : subroutine uderiv(bdberry,cg,gprimd,hdr,istwfk,kberry,kg,kpt_,kptopt,kptrlatt,&
480 : 0 : & mband,mcg,mkmem,mpi_enreg,mpw,natom,nband,nberry,npwarr,nspinor,nsppol,nkpt_,&
481 : : & unddk,fnameabo_1wf)
482 : :
483 : : !Arguments ------------------------------------
484 : : !scalars
485 : : integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_,nspinor
486 : : integer,intent(in) :: nsppol,unddk
487 : : type(MPI_type),intent(in) :: mpi_enreg
488 : : type(hdr_type),intent(inout) :: hdr
489 : : !arrays
490 : : integer,intent(in) :: bdberry(4),istwfk(nkpt_),kberry(3,20),kg(3,mpw*mkmem)
491 : : integer,intent(in) :: kptrlatt(3,3),nband(nkpt_*nsppol),npwarr(nkpt_)
492 : : real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3)
493 : : real(dp),intent(in) :: kpt_(3,nkpt_)
494 : : character(len=fnlen),intent(in) :: fnameabo_1wf
495 : :
496 : : !Local variables -------------------------
497 : : !scalars
498 : : integer,parameter :: master=0
499 : : integer :: iomode,band_in,cg_index_iband,fform,flag1
500 : : integer :: formeig,iband,iberry,icg,idir,ierr,ifor,ii,ikpt,ikpt2,ikpt_
501 : : integer :: index,index1,info,ipert,ipw,isgn,isppol,jband,jj,jkpt,jkpt_
502 : : integer :: maxband,mcg_disk,me,minband,nband_diff,nband_k
503 : : integer :: nkpt,npw_k,pertcase,rdwr,read_k,spaceComm
504 : : integer :: tim_rwwf
505 : : real(dp) :: gmod,twodk
506 : : character(len=500) :: message
507 : : character(len=fnlen) :: fiwf1o
508 : : type(wffile_type) :: wffddk
509 : : !arrays
510 : 0 : integer :: kg_jl(0,0,0),kpt_flag(2*nkpt_)
511 : 0 : integer,allocatable :: cg_index(:,:,:),ikpt_dk(:,:),ipvt(:)
512 : 0 : integer,allocatable :: kg_kpt(:,:,:)
513 : : real(dp) :: det(2,2),diffk(3),diffk2(3),dk(3),gpard(3),klattice(3,3)
514 : : real(dp) :: kptrlattr(3,3),tr(2)
515 : : real(dp) :: cg_disk(0,0,0)
516 : 0 : real(dp),allocatable :: cmatrix(:,:,:),dudk(:,:)
517 : 0 : real(dp),allocatable :: eig_dum_2(:),kpt(:,:)
518 : 0 : real(dp),allocatable :: occ_dum_2(:),phi(:,:,:),u_tilde(:,:,:,:),zgwork(:,:)
519 : 0 : logical,allocatable :: shift_g_2(:,:)
520 : :
521 : : ! *************************************************************************
522 : :
523 [ # # ]: 0 : if(min(2,(1+mpi_enreg%paral_spinor)*nspinor)==2)then
524 : 0 : ABI_ERROR('uderiv: does not yet work for nspinor=2')
525 : : end if
526 : :
527 [ # # # # ]: 0 : if(maxval(istwfk(:))/=1)then
528 : : write(message,'(3a)')&
529 : 0 : & 'Sorry, this routine does not work yet with istwfk/=1.',ch10,&
530 : 0 : & 'This should have been tested previously ...'
531 : 0 : ABI_BUG(message)
532 : : end if
533 : :
534 [ # # ]: 0 : if (kptopt==3) then
535 : 0 : nkpt = nkpt_
536 [ # # # # : 0 : ABI_MALLOC(kpt,(3,nkpt))
# # # # #
# # # ]
537 [ # # # # ]: 0 : kpt(:,:)=kpt_(:,:)
538 [ # # ]: 0 : else if (kptopt==2) then
539 : 0 : nkpt = nkpt_*2
540 [ # # # # : 0 : ABI_MALLOC(kpt,(3,nkpt))
# # # # #
# # # ]
541 [ # # ]: 0 : do ikpt = 1,nkpt/2
542 : 0 : kpt_flag(ikpt) = 0
543 [ # # ]: 0 : kpt(:,ikpt)=kpt_(:,ikpt)
544 : : end do
545 : 0 : index = 0
546 [ # # ]: 0 : do ikpt = (nkpt/2+1),nkpt
547 : 0 : flag1 = 0
548 [ # # ]: 0 : do jkpt = 1, nkpt/2
549 : : if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.&
550 : : & (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))&
551 : : & .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.&
552 : : & (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))&
553 [ # # # # : 0 : & .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.&
# # # # #
# ]
554 [ # # ]: 0 : & (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then
555 : 0 : flag1 = 1
556 : 0 : index = index + 1
557 : : exit
558 : : end if
559 : : end do
560 : 0 : if (flag1==0) then
561 : 0 : kpt_flag(ikpt-index)=ikpt-nkpt/2
562 [ # # ]: 0 : kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2)
563 : : end if
564 : : end do
565 : 0 : nkpt = nkpt - index
566 : : end if
567 : :
568 : : !DEBUG
569 : : !write(101,*) 'beginning write kpt'
570 : : !do ikpt=1,nkpt
571 : : !write(101,*) kpt(:,ikpt)
572 : : !end do
573 : : !ENDDEBUG
574 : :
575 [ # # # # : 0 : ABI_MALLOC(shift_g_2,(nkpt,nkpt))
# # # # #
# # # # #
# # ]
576 : :
577 : : !Compute primitive vectors of the k point lattice
578 : : !Copy to real(dp)
579 [ # # # # ]: 0 : kptrlattr(:,:)=kptrlatt(:,:)
580 : : !Go to reciprocal space (in reduced coordinates)
581 : 0 : call matr3inv(kptrlattr,klattice)
582 : :
583 [ # # ]: 0 : do iberry=1,nberry
584 : :
585 : : ! **************************************************************************
586 : : ! Determine the appended index for ddk 1WF files
587 : :
588 [ # # ]: 0 : do idir=1,3
589 [ # # ]: 0 : if (kberry(idir,iberry) ==1) then
590 : 0 : ipert=natom+1
591 : 0 : pertcase=idir+(ipert-1)*3
592 : : end if
593 : : end do
594 : :
595 : : ! open ddk 1WF file
596 : 0 : formeig=1
597 : :
598 : 0 : call appdig(pertcase,fnameabo_1wf,fiwf1o)
599 : : !call wfk_open_read(wfk, fiwf1o, formeig, iomode, unddk, spaceComm)
600 : :
601 : 0 : spaceComm=xmpi_comm_self; me=0 ; iomode=IO_MODE_FORTRAN
602 : 0 : call WffOpen(iomode,spaceComm,fiwf1o,ierr,wffddk,master,me,unddk)
603 : :
604 : 0 : rdwr=2 ; fform=2
605 : 0 : call hdr_io(fform,hdr,rdwr,wffddk)
606 : :
607 : : ! Define offsets, in case of MPI I/O
608 : 0 : call xdefineOff(formeig,wffddk,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_)
609 : :
610 : : ! *****************************************************************************
611 : : ! Calculate dimensional recip lattice vector along which P is calculated
612 : : ! dk = step to the nearest k point along that direction
613 : : ! in reduced coordinates
614 : :
615 : : dk(:)=dble(kberry(1,iberry))*klattice(:,1)+&
616 : : & dble(kberry(2,iberry))*klattice(:,2)+&
617 [ # # ]: 0 : & dble(kberry(3,iberry))*klattice(:,3)
618 : :
619 [ # # ]: 0 : do idir=1,3
620 [ # # ]: 0 : if (dk(idir)/=0) then
621 : 0 : twodk=2*dk(idir)
622 : : end if
623 : : end do
624 : :
625 [ # # ]: 0 : gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3)
626 : : gmod=sqrt(dot_product(gpard,gpard))
627 : :
628 : : ! ******************************************************************************
629 : : ! Select the k grid points along the direction to compute dudk
630 : : ! dk = step to the nearest k point along that direction
631 : :
632 : : ! For each k point, find k_prim such that k_prim= k + dk mod(G)
633 : : ! where G is a vector of the reciprocal lattice
634 [ # # # # : 0 : ABI_MALLOC(ikpt_dk,(2,nkpt))
# # # # #
# # # ]
635 [ # # # # ]: 0 : ikpt_dk(1:2,1:nkpt)=0
636 [ # # # # ]: 0 : shift_g_2(:,:)= .false.
637 : :
638 [ # # ]: 0 : do ikpt=1,nkpt
639 [ # # ]: 0 : do ikpt2=1,nkpt
640 [ # # ]: 0 : diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:))
641 [ # # ]: 0 : diffk2(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)+dk(:))
642 [ # # # # ]: 0 : if (sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then
643 : 0 : ikpt_dk(1,ikpt)=ikpt2
644 [ # # # # ]: 0 : if(sum(diffk(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
645 : : end if
646 [ # # # # ]: 0 : if (sum(abs(diffk2(:)-nint(diffk2(:))))<3*tol8)then
647 : 0 : ikpt_dk(2,ikpt)=ikpt2
648 [ # # # # ]: 0 : if(sum(diffk2(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
649 : : end if
650 : : end do
651 : : end do
652 : :
653 : 0 : write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
654 : 0 : & ' Computing the derivative for reciprocal vector:',ch10,&
655 : 0 : & dk(:),' (in reduced coordinates)',ch10,&
656 : 0 : & gpard(1:3),' (in cartesian coordinates - atomic units)'
657 : 0 : call wrtout(ab_out,message,'COLL')
658 : 0 : call wrtout(std_out,message,'COLL')
659 : :
660 [ # # ]: 0 : if(nsppol==1)then
661 : : write(message, '(a,i5,a,i5)')&
662 : 0 : & ' From band number',bdberry(1),' to band number',bdberry(2)
663 : : else
664 : : write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')&
665 : 0 : & ' From band number',bdberry(1),' to band number',bdberry(2),' for spin up,',&
666 : 0 : & ch10,&
667 : 0 : & ' from band number',bdberry(3),' to band number',bdberry(4),' for spin down.'
668 : : end if
669 : 0 : call wrtout(ab_out,message,'COLL')
670 : 0 : call wrtout(std_out,message,'COLL')
671 : :
672 : : ! *****************************************************************************
673 [ # # # # : 0 : ABI_MALLOC(dudk,(2,mpw*nspinor*mband*nsppol))
# # # # #
# # # ]
674 [ # # # # : 0 : ABI_MALLOC(eig_dum_2,((2*mband)**formeig*mband))
# # # # #
# # # ]
675 [ # # # # : 0 : ABI_MALLOC(occ_dum_2,((2*mband)**formeig*mband))
# # # # ]
676 [ # # # # ]: 0 : dudk(1:2,:)=0.0_dp
677 [ # # ]: 0 : eig_dum_2=0.0_dp
678 [ # # ]: 0 : occ_dum_2=0.0_dp
679 : :
680 [ # # ]: 0 : if (mkmem/=0) then
681 : :
682 : : ! Find the location of each wavefunction
683 : :
684 [ # # # # : 0 : ABI_MALLOC(cg_index,(mband,nkpt_,nsppol))
# # # # #
# # # # #
# # # # #
# ]
685 : 0 : icg = 0
686 [ # # ]: 0 : do isppol=1,nsppol
687 [ # # ]: 0 : do ikpt=1,nkpt_
688 : 0 : nband_k=nband(ikpt+(isppol-1)*nkpt_)
689 : 0 : npw_k=npwarr(ikpt)
690 [ # # ]: 0 : do iband=1,nband_k
691 : 0 : cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg
692 : : end do
693 : 0 : icg=icg+npw_k*nspinor*nband(ikpt)
694 : : end do
695 : : end do
696 : :
697 : : ! Find the planewave vectors for each k point
698 : : ! SHOULD BE REMOVED WHEN ANOTHER INDEXING TECHNIQUE WILL BE USED FOR kg
699 [ # # # # : 0 : ABI_MALLOC(kg_kpt,(3,mpw*nspinor,nkpt_))
# # # # #
# # # # #
# # ]
700 [ # # # # : 0 : kg_kpt(:,:,:) = 0
# # ]
701 : 0 : index1 = 0
702 [ # # ]: 0 : do ikpt=1,nkpt_
703 : 0 : npw_k=npwarr(ikpt)
704 [ # # ]: 0 : do ipw=1,npw_k*nspinor
705 [ # # ]: 0 : kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1)
706 : : end do
707 : 0 : index1=index1+npw_k*nspinor
708 : : end do
709 : : end if
710 : :
711 : : ! *************************************************************************
712 : : ! Loop over spins
713 [ # # ]: 0 : do isppol=1,nsppol
714 : :
715 : 0 : minband=bdberry(2*isppol-1)
716 : 0 : maxband=bdberry(2*isppol)
717 : :
718 [ # # ]: 0 : if(minband<1)then
719 : 0 : write(message,'(a,i0,a)')' The band limit minband= ',minband,', is lower than 0.'
720 : 0 : ABI_BUG(message)
721 : : end if
722 : :
723 [ # # ]: 0 : if(maxband<1)then
724 : 0 : write(message,'(a,i0,a)')' The band limit maxband= ',maxband,', is lower than 0.'
725 : 0 : ABI_BUG(message)
726 : : end if
727 : :
728 [ # # ]: 0 : if(maxband<minband)then
729 : 0 : write(message,'(a,i0,a,i0)')' maxband= ',maxband,', is lower than minband= ',minband
730 : 0 : ABI_BUG(message)
731 : : end if
732 : :
733 : : ! Loop over k points
734 [ # # ]: 0 : do ikpt_=1,nkpt_
735 : :
736 : 0 : read_k = 0
737 : :
738 : 0 : ikpt=ikpt_
739 : 0 : tr(1) = 1.0_dp
740 : :
741 [ # # ]: 0 : if (kptopt==2) then
742 : : if (read_k == 0) then
743 [ # # ]: 0 : if (kpt_flag(ikpt_)/=0) then
744 : 0 : tr(1) = -1.0_dp
745 : 0 : ikpt= kpt_flag(ikpt_)
746 : : end if
747 : : else !read_k
748 : : if (kpt_flag(ikpt_)/=0) then
749 : : tr(-1*read_k+3) = -1.0_dp
750 : : ikpt= kpt_flag(ikpt_)
751 : : end if
752 : : end if !read_k
753 : : end if !kptopt
754 : :
755 : 0 : nband_k=nband(ikpt+(isppol-1)*nkpt_)
756 : :
757 [ # # ]: 0 : if(nband_k<maxband)then
758 : 0 : write(message,'(a,i0,a,i0)')' maxband=',maxband,', is larger than nband(i,isppol)=',nband_k
759 : 0 : ABI_BUG(message)
760 : : end if
761 : :
762 : 0 : npw_k=npwarr(ikpt)
763 : :
764 [ # # # # : 0 : ABI_MALLOC(u_tilde,(2,npw_k,maxband,2))
# # # # #
# # # # #
# # # # ]
765 [ # # # # : 0 : u_tilde(1:2,1:npw_k,1:maxband,1:2)=0.0_dp
# # # # ]
766 : :
767 : : ! ifor = 1,2 represents forward and backward neighbouring k points of ikpt
768 : : ! respectively along dk direction
769 : :
770 [ # # ]: 0 : do ifor=1,2
771 : :
772 [ # # # # : 0 : ABI_MALLOC(phi,(2,mpw,mband))
# # # # #
# # # # #
# # ]
773 [ # # # # : 0 : ABI_MALLOC(cmatrix,(2,maxband,maxband))
# # # # #
# # # # #
# # ]
774 [ # # # # : 0 : phi(1:2,1:mpw,1:mband)=0.0_dp; cmatrix(1:2,1:maxband,1:maxband)=0.0_dp
# # # # #
# # # ]
775 : :
776 : 0 : isgn=(-1)**ifor
777 : 0 : jkpt_= ikpt_dk(ifor,ikpt_)
778 : :
779 : 0 : tr(2) = 1.0_dp
780 : :
781 : 0 : jkpt=jkpt_
782 : :
783 [ # # ]: 0 : if (kptopt==2) then
784 [ # # ]: 0 : if (read_k == 0) then
785 [ # # ]: 0 : if (kpt_flag(jkpt_)/=0) then
786 : 0 : tr(2) = -1.0_dp
787 : 0 : jkpt= kpt_flag(jkpt_)
788 : : end if
789 : : else !read_k
790 [ # # ]: 0 : if (kpt_flag(jkpt_)/=0) then
791 : 0 : tr(read_k) = -1.0_dp
792 : 0 : jkpt= kpt_flag(jkpt_)
793 : : end if
794 : : end if !read_k
795 : : end if !kptopt
796 : :
797 [ # # ]: 0 : if (ifor==1) read_k = 2
798 : :
799 : 0 : jj = read_k
800 : 0 : ii = -1*read_k+3
801 : :
802 : : call waveformat(cg,cg_disk,cg_index,phi,dk,ii,ikpt,&
803 : : & ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,&
804 : 0 : & minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr)
805 : :
806 : : ! Compute the overlap matrix <u_k|u_k+b>
807 : :
808 [ # # ]: 0 : do iband=minband,maxband
809 : 0 : cg_index_iband=cg_index(iband,ikpt,isppol)
810 [ # # ]: 0 : do jband=minband,maxband
811 [ # # ]: 0 : do ipw=1,npwarr(ikpt)
812 : : cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
813 : : & cg(1,ipw+cg_index_iband)*phi(1,ipw,jband)+&
814 : 0 : & tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)
815 : :
816 : : cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
817 : : & cg(1,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)-&
818 : 0 : & tr(ii)*cg(2,ipw+cg_index_iband)*phi(1,ipw,jband)
819 : : end do
820 : : end do
821 : : end do
822 : :
823 : : ! Compute the inverse of cmatrix(1:2,minband:maxband, minband:maxband)
824 : :
825 : 0 : band_in = maxband - minband + 1
826 [ # # # # : 0 : ABI_MALLOC(ipvt,(maxband))
# # # # #
# # # ]
827 [ # # # # : 0 : ABI_MALLOC(zgwork,(2,1:maxband))
# # # # #
# # # ]
828 : :
829 : : ! Last argument of zgedi means calculate inverse only
830 : 0 : call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
831 : 0 : call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01)
832 : :
833 [ # # ]: 0 : ABI_FREE(zgwork)
834 [ # # ]: 0 : ABI_FREE(ipvt)
835 : :
836 : : ! Compute the product of Inverse overlap matrix with the wavefunction
837 : :
838 [ # # ]: 0 : do iband=minband,maxband
839 [ # # ]: 0 : do ipw=1,npwarr(ikpt)
840 : : u_tilde(1,ipw,iband,ifor)= &
841 : : & dot_product(cmatrix(1,minband:maxband,iband),&
842 : : & phi(1,ipw,minband:maxband))-&
843 : : & dot_product(cmatrix(2,minband:maxband,iband),&
844 [ # # # # ]: 0 : & tr(jj)*phi(2,ipw,minband:maxband))
845 : : u_tilde(2,ipw,iband,ifor)= &
846 : : & dot_product(cmatrix(1,minband:maxband,iband),&
847 : : & tr(jj)*phi(2,ipw,minband:maxband))+&
848 : : & dot_product(cmatrix(2,minband:maxband,iband),&
849 [ # # # # ]: 0 : & phi(1,ipw,minband:maxband))
850 : : end do
851 : : end do
852 [ # # ]: 0 : ABI_FREE(cmatrix)
853 [ # # ]: 0 : ABI_FREE(phi)
854 : :
855 : : end do !ifor
856 : :
857 : : ! Compute dudk for ikpt
858 : :
859 : : npw_k=npwarr(ikpt)
860 : :
861 [ # # ]: 0 : do iband=minband,maxband
862 : :
863 : 0 : icg=(iband-minband)*npw_k
864 : :
865 : : dudk(1,1+icg:npw_k+icg)=(u_tilde(1,1:npw_k,iband,1)-&
866 [ # # ]: 0 : & u_tilde(1,1:npw_k,iband,2))/twodk
867 : :
868 : : dudk(2,1+icg:npw_k+icg)=(u_tilde(2,1:npw_k,iband,1)-&
869 [ # # ]: 0 : & u_tilde(2,1:npw_k,iband,2))/twodk
870 : :
871 : : end do
872 : :
873 : 0 : tim_rwwf=0
874 : 0 : mcg_disk=mpw*nspinor*mband
875 : 0 : nband_diff=maxband-minband+1
876 : : call rwwf(dudk,eig_dum_2,formeig,0,0,ikpt,isppol,kg_kpt(:,:,ikpt),&
877 : : & mband,mcg_disk,mpi_enreg,nband_diff,nband_diff,&
878 : 0 : & npw_k,nspinor,occ_dum_2,2,1,tim_rwwf,wffddk)
879 : :
880 : : !call wfk_read_band_block(wfk, band_block, ikpt, isppol, sc_mode,
881 : : ! kg_k=kg_kpt(:,:,ikpt), cg_k=dudk, eig_k=eig_dum, occ_k=occ_dum)
882 : :
883 [ # # ]: 0 : ABI_FREE(u_tilde)
884 : :
885 : : end do !ikpt
886 : : end do !isppol
887 : :
888 [ # # ]: 0 : ABI_FREE(eig_dum_2)
889 [ # # ]: 0 : ABI_FREE(occ_dum_2)
890 [ # # ]: 0 : ABI_FREE(dudk)
891 : :
892 : 0 : call WffClose(wffddk,ierr)
893 : : !call wfk_close(wfk)
894 : :
895 [ # # ]: 0 : ABI_FREE(kg_kpt)
896 [ # # ]: 0 : ABI_FREE(cg_index)
897 [ # # ]: 0 : ABI_FREE(ikpt_dk)
898 : :
899 : : end do ! iberry
900 : :
901 [ # # ]: 0 : ABI_FREE(shift_g_2)
902 [ # # ]: 0 : ABI_FREE(kpt)
903 : :
904 : 0 : write(std_out,*) 'uderiv: exit '
905 : :
906 : 0 : end subroutine uderiv
907 : : !!***
908 : :
909 : :
910 : : !!****f* ABINIT/waveformat
911 : : !! NAME
912 : : !! waveformat
913 : : !!
914 : : !! FUNCTION
915 : : !! This routine is to find the matched pairs of plane waves between
916 : : !! two neighbouring k points and load a new pw coefficients array cg_new
917 : : !! Was written first by Na Sai (thanks), but unfortunately without
918 : : !! any comment ...
919 : : !!
920 : : !! INPUTS
921 : : !! cg(2,mpw*nspinor*mband*mkmem*nsppol)= input planewave coefficients, in case mkmem/=0
922 : : !! cg_disk(2,mpw*nspinor*mband,2)= input planewave coefficients, in case mkmem==0
923 : : !! cg_index(mband,nkpt_,nsppol)=index of wavefunction iband,ikpt,isppol in the array cg.
924 : : !! dk(3)= step taken to the next k mesh point along the kberry direction (see also isgn)
925 : : !! ii=(to be documented)
926 : : !! ikpt=index of the first k-point in the reduced Brillouin zone
927 : : !! ikpt_=index of the first k-point in the full Brillouin zone
928 : : !! isgn=1 if dk(3) is connecting the k-points (ikpt_ and jkpt)
929 : : !! =-1 if -dk(3) is connecting the k-points
930 : : !! isppol=1 if spin-up, =2 if spin-down
931 : : !! jj=(to be documented)
932 : : !! jkpt=index of the second k-point in the reduced Brillouin zone
933 : : !! jkpt_=index of the second k-point in the full Brillouin zone
934 : : !! kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
935 : : !! planewave and k point
936 : : !! kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
937 : : !! kg_jl(3,mpw,2)=(to be documented)
938 : : !! maxband/minband= control the minimum and maximum band calculated in the
939 : : !! overlap matrix
940 : : !! mband=maximum number of bands (dimension of several cg* arrays)
941 : : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
942 : : !! mcg_disk=size of wave-functions array (cg_disk) =mpw*nspinor*mband
943 : : !! mkmem= if 0, the wavefunctions are input in cg_disk, otherwise in cg
944 : : !! mpw=maximum number of planewaves (dimension of several cg* arrays)
945 : : !! nkpt=number of k points (full Brillouin zone !?!)
946 : : !! nkpt_=number of k points (reduced Brillouin zone !?!)
947 : : !! npwarr(nkpt)=number of planewaves in basis at this k point
948 : : !! nspinor=number of spinorial components of the wavefunctions
949 : : !! nsppol=1 for unpolarized, 2 for spin-polarized
950 : : !! shift_g_2(nkpt,nkpt)=non-zero if a G vector along kberry is needed to connect k points
951 : : !! tr(2)=variable that changes k to -k
952 : : !! G to -G
953 : : !! $c_g$ to $c_g^*$ when time-reversal symetrie is used
954 : : !!
955 : : !! OUTPUT
956 : : !! cg_new(2,mpw,maxband)=planewave coefficients transferred onto the
957 : : !! set of planewaves at k
958 : : !!
959 : : !! SIDE EFFECTS
960 : : !!
961 : : !! SOURCE
962 : :
963 : 0 : subroutine waveformat(cg,cg_disk,cg_index,cg_new,dk,ii,ikpt,&
964 : 0 : & ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,&
965 : 0 : & minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr)
966 : :
967 : : !Arguments ------------------------------------
968 : : !scalars
969 : : integer,intent(in) :: ii,ikpt,ikpt_,isgn,isppol,jj,jkpt,jkpt_,maxband,mband,mcg,mcg_disk
970 : : integer,intent(in) :: minband,mkmem,mpw,nkpt,nkpt_,nspinor,nsppol
971 : : !arrays
972 : : integer,intent(in) :: cg_index(mband,nkpt_,nsppol),kg_jl(3,mpw,2)
973 : : integer,intent(in) :: kg_kpt(3,mpw*nspinor,nkpt_),npwarr(nkpt_)
974 : : real(dp),intent(in) :: cg(2,mcg)
975 : : real(dp),intent(in) :: cg_disk(2,mcg_disk,2),dk(3),kpt(3,nkpt),tr(2)
976 : : real(dp),intent(out) :: cg_new(2,mpw,maxband)
977 : : logical,intent(in) :: shift_g_2(nkpt,nkpt)
978 : :
979 : : !Local variables -------------------------
980 : : !scalars
981 : : integer :: cg_index_iband,iband,ipw,jpw,nomatch,npw_k
982 : : logical :: found_match
983 : : !arrays
984 : : integer :: dg(3)
985 : :
986 : : ! ***********************************************************************
987 : :
988 : 0 : npw_k=npwarr(ikpt)
989 : :
990 : :
991 : 0 : nomatch=0
992 : :
993 : : !If there is no shift of G-vector between ikpt_ and jkpt_
994 [ # # ]: 0 : if(shift_g_2(ikpt_,jkpt_) .eqv. .false.) then
995 : :
996 : : ! DEBUG
997 : : ! write(111,*)'pair', ikpt_,jkpt_,'noshift'
998 : : ! ENDDEBUG
999 : :
1000 : : ! If the original wavefunction is contained in cg_disk
1001 [ # # ]: 0 : if(mkmem==0) then
1002 : :
1003 [ # # ]: 0 : do ipw=1,npw_k
1004 : :
1005 : 0 : found_match = .false.
1006 : :
1007 [ # # ]: 0 : do jpw=1,npwarr(jkpt)
1008 [ # # # # ]: 0 : if (sum(abs(tr(ii)*kg_jl(:,ipw,ii)-tr(jj)*kg_jl(:,jpw,jj)))<3*tol8)then
1009 [ # # ]: 0 : do iband=minband, maxband
1010 : 0 : cg_index_iband=(iband-1)*npwarr(jkpt)
1011 [ # # ]: 0 : cg_new(1:2,ipw,iband)=cg_disk(1:2,jpw+cg_index_iband,jj)
1012 : : end do
1013 : : found_match = .true.
1014 : : exit
1015 : : end if
1016 : : end do
1017 : :
1018 [ # # ]: 0 : if (found_match .eqv. .false.) then
1019 [ # # ]: 0 : do iband=minband,maxband
1020 [ # # ]: 0 : cg_new(1:2,ipw,iband)=zero
1021 : : end do
1022 : 0 : nomatch = nomatch + 1
1023 : : end if
1024 : :
1025 : : end do
1026 : :
1027 : : ! Here, the wavefunctions are contained in cg
1028 : : else
1029 : :
1030 [ # # ]: 0 : do ipw=1,npw_k
1031 : :
1032 : 0 : found_match = .false.
1033 : :
1034 [ # # ]: 0 : do jpw=1,npwarr(jkpt)
1035 [ # # # # ]: 0 : if (sum(abs(tr(ii)*kg_kpt(:,ipw,ikpt)-tr(jj)*kg_kpt(:,jpw,jkpt)))<3*tol8)then
1036 [ # # ]: 0 : do iband=minband, maxband
1037 : 0 : cg_index_iband=cg_index(iband,jkpt,isppol)
1038 [ # # ]: 0 : cg_new(1:2,ipw,iband)=cg(1:2,jpw+cg_index_iband)
1039 : : end do
1040 : : found_match = .true.
1041 : : exit
1042 : : end if
1043 : : end do
1044 : :
1045 [ # # ]: 0 : if (found_match .eqv. .false.) then
1046 [ # # ]: 0 : do iband=minband,maxband
1047 [ # # ]: 0 : cg_new(1:2,ipw,iband)=(0.0_dp,0.0_dp)
1048 : : end do
1049 : 0 : nomatch = nomatch + 1
1050 : : end if
1051 : : end do
1052 : :
1053 : : end if
1054 : :
1055 : : ! DEBUG
1056 : : ! write(111,*) 'normal pair nomatch=',nomatch
1057 : : ! ENDDEBUG
1058 : :
1059 : : ! If there is a G-vector shift between ikpt_ and jkpt_
1060 : : else
1061 : :
1062 : : ! DEBUG
1063 : : ! write(111,*) 'pair',ikpt_,jkpt_,' need shift'
1064 : : ! ENDDEBUG
1065 : :
1066 [ # # ]: 0 : dg(:) = -1*nint(tr(jj)*kpt(:,jkpt)-tr(ii)*kpt(:,ikpt)+isgn*dk(:))
1067 : :
1068 : : ! If the original wavefunction is contained in cg_disk
1069 [ # # ]: 0 : if(mkmem==0) then
1070 : :
1071 [ # # ]: 0 : do ipw=1,npw_k
1072 : :
1073 : 0 : found_match = .false.
1074 : :
1075 [ # # ]: 0 : do jpw=1,npwarr(jkpt)
1076 [ # # # # ]: 0 : if (sum(abs(tr(ii)*kg_jl(:,ipw,ii)-(tr(jj)*kg_jl(:,jpw,jj)-&
1077 : 0 : & dg(:))))<3*tol8)then
1078 : :
1079 [ # # ]: 0 : do iband=minband, maxband
1080 : 0 : cg_index_iband=(iband-1)*npwarr(jkpt)
1081 [ # # ]: 0 : cg_new(1:2,ipw,iband)=cg_disk(1:2,jpw+cg_index_iband,jj)
1082 : : end do
1083 : : found_match = .true.
1084 : : exit
1085 : : end if
1086 : : end do
1087 : :
1088 [ # # ]: 0 : if (found_match .eqv. .false.) then
1089 [ # # ]: 0 : do iband=minband,maxband
1090 [ # # ]: 0 : cg_new(1:2,ipw,iband)=(0.0_dp,0.0_dp)
1091 : : end do
1092 : 0 : nomatch = nomatch + 1
1093 : : end if
1094 : : end do
1095 : :
1096 : : ! Here, the wavefunctions are contained in cg
1097 : : else
1098 : :
1099 [ # # ]: 0 : do ipw=1,npw_k
1100 : :
1101 : 0 : found_match = .false.
1102 : :
1103 [ # # ]: 0 : do jpw=1,npwarr(jkpt)
1104 [ # # # # ]: 0 : if (sum(abs(tr(ii)*kg_kpt(:,ipw,ikpt)-(tr(jj)*kg_kpt(:,jpw,jkpt)-&
1105 : 0 : & dg(:))))<3*tol8)then
1106 [ # # ]: 0 : do iband=minband, maxband
1107 : 0 : cg_index_iband=cg_index(iband,jkpt,isppol)
1108 [ # # ]: 0 : cg_new(1:2,ipw,iband)=cg(1:2,jpw+cg_index_iband)
1109 : : end do
1110 : : found_match = .true.
1111 : : exit
1112 : : end if
1113 : : end do
1114 : :
1115 [ # # ]: 0 : if (found_match .eqv. .false.) then
1116 [ # # ]: 0 : do iband=minband,maxband
1117 [ # # ]: 0 : cg_new(1:2,ipw,iband)=zero
1118 : : end do
1119 : 0 : nomatch = nomatch + 1
1120 : : end if
1121 : : end do
1122 : :
1123 : : end if
1124 : :
1125 : : ! DEBUG
1126 : : ! write(111,*) 'special pair nomatch=',nomatch
1127 : : ! ENDDEBUG
1128 : :
1129 : : end if
1130 : :
1131 : 0 : end subroutine waveformat
1132 : : !!***
1133 : :
1134 : : end module m_elpolariz
1135 : : !!***
|
