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