1#include "f_defs.h"
2
3!-----------------------------------------------------------------------
4subroutine genwf_kernel(crys,gvec,kg,kgq,syms, &
5  wfnc,wfnv,nspin,ik,ic,iv,indexq,xct,intwfnv,intwfnc)
6!
7!     input: crys, gvec, kg,  syms types
8!            nspin  (number of spins)
9!            ik     (index of k-point in full BZ)
10!            ib     (index of ck block)
11!     output: wfnc   (conduction states at point ik, block ib)
12!             wfnv   (valence states at point ik)
13!
14!-----------------------------------------------------------------------
15
16  use global_m
17  use blas_m
18  use gmap_m
19  use input_utils_m
20  use sort_m
21  use misc_m
22  implicit none
23
24  type (crystal), intent(in) :: crys
25  type (gspace), intent(in) :: gvec
26  type (grid), intent(in) :: kg
27  type (grid), intent(in) :: kgq
28  type (symmetry), intent(in) :: syms
29  type (wavefunction), intent(out) :: wfnc,wfnv
30  integer, intent(in) :: nspin,ik,ic,iv
31  integer, intent(in) :: indexq(kg%nf)
32  type (int_wavefunction), intent(in) :: intwfnc,intwfnv
33  type (xctinfo), intent(in) :: xct
34
35  character :: errmsg*100
36  integer :: cnb,cng,cns,cnsp,vng,vns,vnsp
37  integer :: ig,ii,jj,kk,info
38  integer :: ivown,icown,ikown,invband,incband,ijk
39  real(DP) :: xnorm,qk(3)
40
41  integer, save :: ikold=0
42  integer, save :: ikold2=0
43  integer, save :: vngold=0
44  integer, save :: cngold=0
45  integer, save :: ifirst=0
46
47  integer, allocatable :: isorti(:)
48  integer, allocatable :: ind(:),isort(:)
49  integer, allocatable :: indq(:),isortq(:)
50  integer, allocatable, save :: ind_old(:),isort_old(:)
51  integer, allocatable, save :: indq_old(:),isortq_old(:)
52  integer, allocatable, save :: ind_old2(:),isort_old2(:)
53  integer, allocatable, save :: indq_old2(:),isortq_old2(:)
54  real(DP), allocatable :: ekin(:)
55  SCALAR, allocatable :: ccg(:,:,:),vcg(:,:,:),ph(:),phq(:)
56  SCALAR, allocatable, save :: ph_old(:),phq_old(:)
57  SCALAR, allocatable, save :: ph_old2(:),phq_old2(:)
58
59  PUSH_SUB(genwf_kernel)
60
61!-----------------------------------------------------------------------
62! Deal with the valence wavefunctions
63
64!      write(6,*) peinf%inode, 'in genwf',iv,ic,ik
65
66  if (xct%qflag .eq. 0) then
67  else
68    ivown = peinf%ipev(peinf%inode+1,iv,kg%indr(ik))
69    ikown = peinf%ipek(peinf%inode+1,kg%indr(ik))
70  endif
71
72  vng = intwfnv%ng(ikown)
73  vns = intwfnv%nspin
74  vnsp = intwfnv%nspinor
75  if (xct%ivpar .eq. 1) then
76    SAFE_ALLOCATE(vcg, (vng,1,vns*vnsp))
77  else
78    SAFE_ALLOCATE(vcg, (vng,xct%nvb_co,vns*vnsp))
79  endif
80
81  SAFE_ALLOCATE(indq, (vng))
82  SAFE_ALLOCATE(phq, (vng))
83  SAFE_ALLOCATE(isortq, (gvec%ng))
84
85  if (ik .ne. ikold .and. ik .ne. ikold2 .and. ifirst .ne. 0) then
86    SAFE_DEALLOCATE(indq_old2)
87    SAFE_DEALLOCATE(phq_old2)
88    SAFE_DEALLOCATE(isortq_old2)
89    SAFE_ALLOCATE(indq_old2, (vngold))
90    SAFE_ALLOCATE(phq_old2, (vngold))
91    SAFE_ALLOCATE(isortq_old2, (gvec%ng))
92
93    indq_old2 = indq_old
94    phq_old2 = phq_old
95    isortq_old2 = isortq_old
96
97    SAFE_DEALLOCATE(indq_old)
98    SAFE_DEALLOCATE(phq_old)
99    SAFE_DEALLOCATE(isortq_old)
100    SAFE_ALLOCATE(indq_old, (vng))
101    SAFE_ALLOCATE(phq_old, (vng))
102    SAFE_ALLOCATE(isortq_old, (gvec%ng))
103  endif
104  if (ifirst .eq. 0) then
105    SAFE_ALLOCATE(indq_old, (vng))
106    SAFE_ALLOCATE(phq_old, (vng))
107    SAFE_ALLOCATE(isortq_old, (gvec%ng))
108    SAFE_ALLOCATE(indq_old2, (vng))
109    SAFE_ALLOCATE(phq_old2, (vng))
110    SAFE_ALLOCATE(isortq_old2, (gvec%ng))
111  endif
112
113! Initalize parameters in variable wfnv
114
115  wfnv%ng=vng
116  wfnv%nband=1
117  wfnv%nspin=vns
118  wfnv%nspinor=vnsp
119  if (intwfnv%nspin.ne.nspin) then
120    write(errmsg,'(2a,2i2)') 'spin number mismatch: ', nspin, vns
121    call die(errmsg, only_root_writes = .true.)
122  endif
123
124  if (xct%ivpar .eq. 1) then
125    SAFE_ALLOCATE(wfnv%cg, (wfnv%ng,1,wfnv%nspin*wfnv%nspinor))
126    vcg(1:vng,1,:) = intwfnv%cg(1:vng,ivown,:)
127  else
128    SAFE_ALLOCATE(wfnv%cg, (wfnv%ng,xct%nvb_co,wfnv%nspin*wfnv%nspinor))
129    vcg(1:vng,:,:) = intwfnv%cg(1:vng,ivown:ivown+xct%nvb_co-1,:)
130  endif
131  SAFE_ALLOCATE(wfnv%isort, (gvec%ng))
132
133  isortq(:) = intwfnv%isort(:,ikown)
134
135! Compute inverse index array of Fourier components around rk-kpoint
136
137  if (ik .ne. ikold .and. ik .ne. ikold2) then
138    SAFE_ALLOCATE(isorti, (gvec%ng))
139    isorti(:)=0
140    do ii=1,wfnv%ng
141      isorti(isortq(ii))=ii
142    enddo
143  endif
144
145! Compute index array of Fourier components around fk-kpoint
146
147  if (ik .ne. ikold .and. ik .ne. ikold2) then
148    if (peinf%verb_debug) then
149      write(6,*) 'ikneikold',peinf%inode,ik,ikold
150    endif
151    SAFE_ALLOCATE(ekin, (gvec%ng))
152    if (xct%qflag .eq. 0) then
153      qk(1:3)=kgq%f(1:3,indexq(ik))
154    elseif (xct%qflag .eq. 2) then
155      qk(1:3)=kg%f(1:3,indexq(ik))
156    else
157      qk(1:3)=kg%f(1:3,ik)
158    endif
159    call kinetic_energies(gvec, crys%bdot, ekin, qvec = qk)
160    call sortrx(gvec%ng,ekin,isortq)
161    SAFE_DEALLOCATE(ekin)
162    isortq_old(:)=isortq(:)
163  elseif (ik .eq. ikold) then
164    isortq(:)=isortq_old(:)
165  else
166    isortq(:)=isortq_old2(:)
167  endif
168
169! Find ind and ph relating wavefunctions in fk to rk-kpoint
170
171  if (ik .ne. ikold .and. ik .ne. ikold2) then
172    indq=0
173    phq=ZERO
174    if (xct%qflag .eq. 0) then
175    else
176      call gmap(gvec,syms,wfnv%ng,kg%itran(ik), &
177        kg%kg0(:,ik),isortq,isorti,indq,phq,.true.)
178    endif
179    SAFE_DEALLOCATE(isorti)
180    indq_old(:)=indq(:)
181    phq_old(:)=phq(:)
182  else if (ik .eq. ikold) then
183    indq(:)=indq_old(:)
184    phq(:)=phq_old(:)
185  else
186    indq(:)=indq_old2(:)
187    phq(:)=phq_old2(:)
188  endif
189
190! Compute and renormalize valence wavefunctions in fk-kpoint
191
192  if (xct%ivpar .eq. 1) invband = 1
193  if (xct%ivpar .eq. 0) invband = xct%nvb_co
194
195  do ijk = 1, invband
196    do kk=1,wfnv%nspin
197      do jj=1,wfnv%nspinor
198        do ii=1,wfnv%ng
199          if (indq(ii) .gt. 0) then
200            wfnv%cg(ii,ijk,kk*jj)=phq(ii)*vcg(indq(ii),ijk,kk*jj)
201          else
202            wfnv%cg(ii,ijk,kk*jj) = ZERO
203          endif
204        enddo !ii
205      enddo ! jj
206      call checknorm('wfnv%cg',ijk,ik,wfnv%ng,kk,wfnv%nspinor,wfnv%cg(1:wfnv%ng,ijk,:))
207    enddo ! kk
208  enddo ! ijk
209  vcg=wfnv%cg
210
211
212  wfnv%cg=vcg
213  wfnv%isort=isortq
214
215
216  SAFE_DEALLOCATE(vcg)
217  SAFE_DEALLOCATE(indq)
218  SAFE_DEALLOCATE(phq)
219  SAFE_DEALLOCATE(isortq)
220
221!-----------------------------------------------------------------------
222! Deal with the conduction wavefunctions
223
224  icown = peinf%ipec(peinf%inode+1,ic,kg%indr(ik))
225  ikown = peinf%ipek(peinf%inode+1,kg%indr(ik))
226
227  cng = intwfnc%ng(ikown)
228  cnb = 1
229  cns = intwfnc%nspin
230  cnsp = intwfnc%nspinor
231
232  if (xct%icpar .eq. 1) then
233    SAFE_ALLOCATE(ccg, (cng,1,cns*cnsp))
234  else
235    SAFE_ALLOCATE(ccg, (cng,xct%ncb_co,cns*cnsp))
236  endif
237
238  SAFE_ALLOCATE(ind, (cng))
239  SAFE_ALLOCATE(ph, (cng))
240  SAFE_ALLOCATE(isort, (gvec%ng))
241
242  if (ik .ne. ikold .and. ik .ne. ikold2 .and. ifirst .ne. 0) then
243    SAFE_DEALLOCATE(ind_old2)
244    SAFE_DEALLOCATE(ph_old2)
245    SAFE_DEALLOCATE(isort_old2)
246    SAFE_ALLOCATE(ind_old2, (cngold))
247    SAFE_ALLOCATE(ph_old2, (cngold))
248    SAFE_ALLOCATE(isort_old2, (gvec%ng))
249
250    ind_old2 = ind_old
251    ph_old2 = ph_old
252    isort_old2 = isort_old
253
254    SAFE_DEALLOCATE(ind_old)
255    SAFE_DEALLOCATE(ph_old)
256    SAFE_DEALLOCATE(isort_old)
257    SAFE_ALLOCATE(ind_old, (cng))
258    SAFE_ALLOCATE(ph_old, (cng))
259    SAFE_ALLOCATE(isort_old, (gvec%ng))
260  endif
261  if (ifirst .eq. 0) then
262    SAFE_ALLOCATE(ind_old, (cng))
263    SAFE_ALLOCATE(ph_old, (cng))
264    SAFE_ALLOCATE(isort_old, (gvec%ng))
265    SAFE_ALLOCATE(ind_old2, (cng))
266    SAFE_ALLOCATE(ph_old2, (cng))
267    SAFE_ALLOCATE(isort_old2, (gvec%ng))
268  endif
269
270  wfnc%ng=cng
271  wfnc%nband=1
272  wfnc%nspin=cns
273  wfnc%nspinor=cnsp
274
275  if (cns.ne.nspin) then
276    write(errmsg,'(2a,2i2)') 'spin number mismatch: ', nspin, cns
277    call die(errmsg, only_root_writes = .true.)
278  endif
279
280  SAFE_ALLOCATE(wfnc%isort, (gvec%ng))
281
282  isort(:)=intwfnc%isort(:,ikown)
283  if (xct%icpar .eq. 1) then
284    SAFE_ALLOCATE(wfnc%cg, (wfnc%ng,1,wfnc%nspin*wfnc%nspinor))
285    ccg(1:cng,1,:) = intwfnc%cg(1:cng,icown,:)
286  else
287    SAFE_ALLOCATE(wfnc%cg, (wfnc%ng,xct%ncb_co,wfnc%nspin*wfnc%nspinor))
288    ccg(1:cng,:,:) = intwfnc%cg(1:cng,icown:icown+xct%ncb_co-1,:)
289  endif
290
291! JRD: Below is now necessary because kg might be different from kgq
292! Compute inverse index array of Fourier components around rk-kpoint
293
294  if (ik .ne. ikold .and. ik .ne. ikold2) then
295    SAFE_ALLOCATE(isorti, (gvec%ng))
296    isorti(:)=0
297    do ii=1,wfnc%ng
298      isorti(isort(ii))=ii
299    enddo
300  endif
301
302! Compute index array of Fourier components around fk-kpoint
303
304  if (ik .ne. ikold .and. ik .ne. ikold2) then
305    SAFE_ALLOCATE(ekin, (gvec%ng))
306    call kinetic_energies(gvec, crys%bdot, ekin, qvec = kg%f(:, ik))
307    call sortrx(gvec%ng,ekin,isort)
308    SAFE_DEALLOCATE(ekin)
309    isort_old(:) = isort(:)
310  else if (ik .eq. ikold) then
311    isort(:)=isort_old(:)
312  else
313    isort(:)=isort_old2(:)
314  endif
315
316! Find ind and ph relating wavefunctions in fk to rk-kpoint
317
318  if (ik .ne. ikold .and. ik .ne. ikold2) then
319    ind=0
320    ph=ZERO
321    call gmap(gvec,syms,wfnc%ng,kg%itran(ik), &
322      kg%kg0(:,ik),isort,isorti,ind,ph,.true.)
323    SAFE_DEALLOCATE(isorti)
324    ind_old(:)=ind(:)
325    ph_old(:)=ph(:)
326  else if (ik .eq. ikold) then
327    ind(:)=ind_old(:)
328    ph(:)=ph_old(:)
329  else
330    ind(:)=ind_old2(:)
331    ph(:)=ph_old2(:)
332  endif
333
334! Compute and renormalize conduction wavefunctions
335
336  if (xct%icpar .eq. 1) incband =1
337  if (xct%icpar .eq. 0) incband = xct%ncb_co
338
339  do ijk =1, incband
340    do kk=1,wfnc%nspin
341      do jj=1,wfnc%nspinor
342        do ii=1,wfnc%ng
343          if (ind(ii) .gt. 0) then
344            wfnc%cg(ii,ijk,kk*jj)=ph(ii)*ccg(ind(ii),ijk,kk*jj)
345          else
346            wfnc%cg(ii,ijk,kk*jj)=ZERO
347          endif
348        enddo
349      enddo ! jj
350      call checknorm('wfnc%cg',ijk,ik,wfnc%ng,kk,wfnc%nspinor,wfnc%cg(1:wfnc%ng,ijk,:))
351    enddo ! kk
352  enddo ! ijk
353  ccg=wfnc%cg
354
355
356  wfnc%cg=ccg
357  wfnc%isort=isort
358
359  SAFE_DEALLOCATE(ccg)
360  SAFE_DEALLOCATE(ind)
361  SAFE_DEALLOCATE(ph)
362  SAFE_DEALLOCATE(isort)
363
364  if (ik .ne. ikold .and. ik .ne. ikold2) then
365    ikold2=ikold
366    ikold=ik
367    cngold=cng
368    vngold=vng
369  endif
370
371  ifirst=-1
372
373  POP_SUB(genwf_kernel)
374
375  return
376end subroutine genwf_kernel
377