TABLE OF CONTENTS


ABINIT/elpolariz [ Functions ]

[ Top ] [ Functions ]

NAME

 elpolariz

FUNCTION

 Calculate corrections to total energy from polarising
 electric field with or without Berry phases (berryopt keyword)

INPUTS

 atindx1(natom)=index table for atoms, inverse of atindx (see gstate.f)
 cg(2,mcg)=planewave coefficients of wavefunctions
 cprj(natom,mcprj*usecprj)=<p_lmn|Cnk> coefficients for each WF |Cnk>
                           and each |p_lmn> non-local projector
 dtfil <type(datafiles_type)>=variables related to files
 dtset <type(dataset_type)>=all input variables in this dataset
 gprimd(3,3)=reciprocal space dimensional primitive translations
 hdr <type(hdr_type)>=the header of wf, den and pot files
 kg(3,mpw*mkmem)=reduced planewave coordinates
 mband=maximum number of bands
 mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
 mkmem=number of k points treated by this node.
 mpi_enreg=information about MPI parallelization
 mpw=maximum dimensioned size of npw
 my_natom=number of atoms treated by current processor
 natom=number of atoms in cell
 nattyp(ntypat)= # atoms of each type.
 nkpt=number of k points
 npwarr(nkpt)=number of planewaves in basis at this k point
 nsppol=1 for unpolarized, 2 for spin-polarized
 ntypat=number of types of atoms in unit cell
 nkpt=number of k-points
 option = 1: compute Berryphase polarization
          2: compute finite difference expression of the ddk
          3: compute polarization & ddk
 pawrhoij(my_natom*usepaw) <type(pawrhoij_type)> atomic occupancies
 pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
 pel_cg(3) = reduced coordinates of the electronic polarization (a. u.)
             computed in the SCF loop
 pelev(3)= expectation value polarization term (PAW only) in cartesian coordinates
 psps <type(pseudopotential_type)>=variables related to pseudopotentials
 pwind(pwind_alloc,2,3) = array used to compute
           the overlap matrix smat between k-points (see initberry.f)
 pwind_alloc = first dimension of pwind
 pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations
 rprimd(3,3)=dimensional real space primitive translations (bohr)
 ucvol=unit cell volume in bohr**3.
 usecprj=1 if cprj datastructure has been allocated
 xred(3,natom)=reduced atomic coordinates

OUTPUT

  (see side effects)

SIDE EFFECTS

 dtefield <type(efield_type)> = variables related to Berry phase
       and electric field calculations (see initberry.f).
       In case berryopt = 4/6/7/14/16/17, the overlap matrices computed
       in this routine are stored in dtefield%smat in order
       to be used in the electric field calculation.
 enefield=field energy
 etotal=total energy, might be correct by improved polarization computation
 pel(3) = reduced coordinates of the electronic polarization (a. u.)
 pion(3)= reduced coordinates of the ionic polarization (a. u.)

PARENTS

      afterscfloop

CHILDREN

      berryphase,berryphase_new,metric,uderiv,wrtout

SOURCE

131 subroutine elpolariz(atindx1,cg,cprj,dtefield,dtfil,dtset,etotal,enefield,gprimd,hdr,&
132 & kg,mband,mcg,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,nattyp,nkpt,&
133 & npwarr,nsppol,ntypat,pawrhoij,pawtab,&
134 & pel,pel_cg,pelev,pion,psps,pwind,pwind_alloc,&
135 & pwnsfac,rprimd,ucvol,usecprj,xred)
136 
137 
138 !This section has been created automatically by the script Abilint (TD).
139 !Do not modify the following lines by hand.
140 #undef ABI_FUNC
141 #define ABI_FUNC 'elpolariz'
142 !End of the abilint section
143 
144  implicit none
145 
146 !Arguments ------------------------------------
147 !scalars
148  integer,intent(in) :: mband,mcg,mcprj,mkmem,mpw,my_natom,natom,nkpt,nsppol,ntypat
149  integer,intent(in) :: pwind_alloc,usecprj
150  real(dp),intent(in) :: ucvol
151  real(dp),intent(inout) :: enefield,etotal
152  type(MPI_type),intent(in) :: mpi_enreg
153  type(datafiles_type),intent(in) :: dtfil
154  type(dataset_type),intent(inout) :: dtset
155  type(efield_type),intent(inout) :: dtefield
156  type(hdr_type),intent(inout) :: hdr
157  type(pseudopotential_type),intent(in) :: psps
158 !arrays
159  integer,intent(in) :: atindx1(natom),kg(3,mpw*mkmem),nattyp(ntypat)
160  integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3)
161  real(dp),intent(in) :: cg(2,mcg),gprimd(3,3)
162  real(dp),intent(in) :: pel_cg(3),pwnsfac(2,pwind_alloc),rprimd(3,3)
163  real(dp),intent(inout) :: pel(3),pelev(3),pion(3),xred(3,natom)
164  type(pawcprj_type),intent(in) :: cprj(natom,mcprj*usecprj)
165  type(pawrhoij_type), intent(in) :: pawrhoij(my_natom*psps%usepaw)
166  type(pawtab_type),intent(in) :: pawtab(dtset%ntypat*psps%usepaw)
167 
168 !Local variables-------------------------------
169 !scalars
170  integer :: my_nspinor,option,unit_out,iir,jjr,kkr
171  real(dp) :: pdif_mod,eenth,ucvol_local
172  character(len=500) :: message
173 !arrays
174  real(dp) :: gmet(3,3),gprimdlc(3,3),pdif(3),ptot(3),red_ptot(3),rmet(3,3)
175 !! ptot(3) = total polarization (not reduced) REC
176 !! red_ptot(3) = internal reduced total polarization REC
177  real(dp) ::  A(3,3),A1(3,3),A_new(3,3),efield_new(3)
178 
179 ! *************************************************************************
180 
181  DBG_ENTER("COLL")
182 
183  if (usecprj==0.and.psps%usepaw==1) then
184    write (message,'(3a)')&
185 &   'cprj datastructure must be allocated !',ch10,&
186 &   'Action: change pawusecp input keyword.'
187    MSG_ERROR(message)
188  end if
189 
190  my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor)
191 
192  if(dtset%berryopt>0 .and. dtset%berryopt/=4 .and. dtset%berryopt/=6 .and. dtset%berryopt/=7 .and. &
193 & dtset%berryopt/=14 .and. dtset%berryopt/=16 .and. dtset%berryopt/=17)then  !!HONG
194 
195    if (dtset%berryopt==1 .or. dtset%berryopt==3) then
196      call berryphase(atindx1,dtset%bdberry,cg,gprimd,dtset%istwfk,&
197 &     dtset%kberry,kg,dtset%kptns,dtset%kptopt,dtset%kptrlatt,&
198 &     mband,mcg,mkmem,mpw,natom,nattyp,dtset%nband,dtset%nberry,npwarr,&
199 &     my_nspinor,nsppol,psps%ntypat,nkpt,rprimd,ucvol,&
200 &     xred,psps%ziontypat)
201    end if
202 
203    if (dtset%berryopt==2 .or. dtset%berryopt==3) then
204      call uderiv(dtset%bdberry,cg,gprimd,hdr,dtset%istwfk,&
205 &     dtset%kberry,kg,dtset%kptns,dtset%kptopt,&
206 &     dtset%kptrlatt,mband,mcg,mkmem,mpi_enreg,mpw,&
207 &     natom,dtset%nband,dtset%nberry,npwarr,my_nspinor,nsppol,&
208 &     nkpt,dtfil%unddk,dtfil%fnameabo_1wf)
209    end if
210 
211  else if(dtset%berryopt<0 .or. dtset%berryopt==4 .or. dtset%berryopt==6 .or. dtset%berryopt==7 .or.  &
212 &   dtset%berryopt==14 .or. dtset%berryopt==16 .or. dtset%berryopt==17)then   !!HONG
213 
214    select case (dtset%berryopt)
215    case (-5)
216      option = 2
217    case (-3)
218      option = 3
219    case (-2)
220      option = 2
221    case (-1)
222      option = 1
223    case (4)
224      option = 1
225      pel(:) = zero
226      pelev(:) = zero
227    case (6)                !!HONG
228      option = 1
229      pel(:) = zero
230      pelev(:) = zero
231    case (7)                !!HONG
232      option = 1
233      pel(:) = zero
234      pelev(:) = zero
235    case (14)                !!HONG
236      option = 1
237      pel(:) = zero
238      pelev(:) = zero
239    case (16)                !!HONG
240      option = 1
241      pel(:) = zero
242      pelev(:) = zero
243    case (17)                !!HONG
244      option = 1
245      pel(:) = zero
246      pelev(:) = zero
247    end select
248 
249    unit_out = ab_out
250    call berryphase_new(atindx1,cg,cprj,dtefield,dtfil,dtset,psps,&
251 &   gprimd,hdr,psps%indlmn,kg,&
252 &   psps%lmnmax,mband,mcg,mcprj,mkmem,mpi_enreg,mpw,my_natom,natom,npwarr,&
253 &   nsppol,psps%ntypat,nkpt,option,pawrhoij,&
254 &   pawtab,pel,pelev,pion,ptot,red_ptot,pwind,&                            !!REC
255 &  pwind_alloc,pwnsfac,rprimd,dtset%typat,ucvol,&
256 &   unit_out,usecprj,psps%usepaw,xred,psps%ziontypat)
257 
258    dtefield%red_ptot1(:)=red_ptot(:)
259 
260    if (dtset%berryopt == 4 .or. dtset%berryopt == 6 .or. dtset%berryopt == 7 .or.  &
261 &   dtset%berryopt == 14 .or. dtset%berryopt == 16 .or. dtset%berryopt == 17 ) then   !!HONG
262 
263 !    Check if pel has the same value as pel_cg
264 !    if (psps%usepaw == 1) pel(:) = pel(:) + pelev(:) ! add on-site term for PAW
265 !    if (psps%usepaw == 1) red_ptot(:) = red_ptot(:) + pelev(:) ! add on-site term for PAW  !! REC
266 !    above line suppressed because in the PAW case, pel already includes all on-site
267 !    terms and pelev should not be added in additionally. We are computing pelev separately for
268 !    reporting purposes only.
269 !    13 June 2012 J Zwanziger
270 
271      pdif(:) = pel_cg(:) - pel(:)
272      pdif_mod = pdif(1)**2 + pdif(2)**2 + pdif(3)**2
273 
274      if (pdif_mod > tol8) then
275        write(message,'(11(a),e16.9)')ch10,&
276 &       ' scfcv (electric field calculation) : WARNING -',ch10,&
277 &       '   The difference between pel (electronic Berry phase updated ',ch10,&
278 &       '   at each SCF cycle)',ch10,&
279 &       '   and pel_cg (electronic Berryphase computed using the ',&
280 &       'berryphase routine) is',ch10,&
281 &       '   pdif_mod = ',pdif_mod
282        call wrtout(std_out,message,'COLL')
283        write(message,'(a,6(a,e16.9,a))') ch10,&
284 &       'pel_cg(1) = ',pel_cg(1),ch10,&
285 &       'pel_cg(2) = ',pel_cg(2),ch10,&
286 &       'pel_cg(3) = ',pel_cg(3),ch10,&
287 &       'pel(1) = ',pel(1),ch10,&
288 &       'pel(2) = ',pel(2),ch10,&
289 &       'pel(3) = ',pel(3),ch10
290        MSG_ERROR(message)
291      end if
292 
293 !    Use this (more accurate) value of P to recompute enefield
294      if (dtset%berryopt == 4 .or. dtset%berryopt == 14 ) then             !!HONG
295        etotal = etotal - enefield
296 
297        enefield = -dot_product(dtset%red_efieldbar,red_ptot)
298        call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
299        eenth = zero
300        do iir=1,3
301          do jjr=1,3
302            eenth= eenth+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)         !! HONG g^{-1})_ij ebar_i ebar_j
303          end do
304        end do
305        eenth=-1_dp*(ucvol_local/(8.0d0*pi))*eenth
306        enefield=enefield+eenth
307 
308        etotal = etotal + enefield
309 
310        write(message,'(a,a)')ch10,&
311 &       ' Stress tensor under a constant electric field:'
312        call wrtout(std_out,message,'COLL')
313        call wrtout(ab_out,message,'COLL')
314 
315      end if
316 
317 !    ! In finite D-field case, turn it into internal energy    !!HONG
318      if (dtset%berryopt == 6 .or. dtset%berryopt == 16 )  then
319        etotal = etotal - enefield
320 
321        enefield=zero
322        call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
323        do iir=1,3
324          do jjr=1,3
325            enefield= enefield+gmet(iir,jjr)*dtset%red_efieldbar(iir)*dtset%red_efieldbar(jjr)         !! HONG g^{-1})_ij ebar_i ebar_j
326          end do
327        end do
328        enefield= ucvol_local/(8.0d0*pi)*enefield
329 
330        etotal = etotal + enefield
331 
332        write(message,'(a,a)')ch10,&
333 &       ' Stress tensor under a constant electric displacement field:'
334        call wrtout(std_out,message,'COLL')
335        call wrtout(ab_out,message,'COLL')
336 
337      end if
338 
339 !    HONG  calculate internal energy and electric enthalpy for mixed BC case.
340      if ( dtset%berryopt == 17 ) then
341        etotal = etotal - enefield
342        enefield = zero
343 
344        call metric(gmet,gprimdlc,-1,rmet,rprimd,ucvol_local)
345        A(:,:)=(4*pi/ucvol_local)*rmet(:,:)
346        A1(:,:)=A(:,:)
347        A_new(:,:)=A(:,:)
348        efield_new(:)=dtset%red_efield(:)
349        eenth = zero
350 
351        do kkr=1,3
352          if (dtset%jfielddir(kkr)==1) then    ! fixed ebar direction
353 !          step 1 add -ebar*p
354            eenth=eenth - dtset%red_efieldbar(kkr)*red_ptot(kkr)
355 
356 !          step 2  chang to e_new (change e to ebar)
357            efield_new(kkr)=dtset%red_efieldbar(kkr)
358 
359 !          step 3  chang matrix A to A1
360 
361            do iir=1,3
362              do jjr=1,3
363                if (iir==kkr .and. jjr==kkr) A1(iir,jjr)=-1.0/A(kkr,kkr)
364                if ((iir==kkr .and. jjr/=kkr) .or.  (iir/=kkr .and.  jjr==kkr)) &
365 &               A1(iir,jjr)=-1.0*A(iir,jjr)/A(kkr,kkr)
366                if (iir/=kkr .and. jjr/=kkr) A1(iir,jjr)=A(iir,jjr)-A(iir,kkr)*A(kkr,jjr)/A(kkr,kkr)
367              end do
368            end do
369 
370            A(:,:)=A1(:,:)
371            A_new(:,:)=A1(:,:)
372          end if
373 
374        end do  ! end fo kkr
375 
376 
377        do iir=1,3
378          do jjr=1,3
379            eenth= eenth+(1/2.0)*A_new(iir,jjr)*efield_new(iir)*efield_new(jjr)
380          end do
381        end do
382 
383        enefield=eenth
384        etotal = etotal + enefield
385 
386        write(message,'(a,a)')ch10,&
387 &       ' Stress tensor under a constant (mixed) electric and electric displacement field:'
388        call wrtout(std_out,message,'COLL')
389        call wrtout(ab_out,message,'COLL')
390 
391      end if   ! berryopt==17
392 
393 
394 !    MVeithen: to clarify
395 !    Which stress tensor should be used in structural optimizations?
396 !    The one at constant electric field or at constant potential drop.
397 !    write(message,'(a,a)')ch10,&
398 !    &     ' Stress tensor imposing a constant electric field:'
399 !    call wrtout(std_out,message,'COLL')
400 !    call wrtout(ab_out,message,'COLL')
401 
402    end if ! dtset%berryopt == 4/6/7/14/16/17
403 
404  end if ! dtset%berryopt>0 or dtset%berryopt/=4/6/7/14/16/17
405 
406  DBG_EXIT("COLL")
407 
408 end subroutine elpolariz

ABINIT/m_elpolariz [ Modules ]

[ Top ] [ Modules ]

NAME

 m_elpolariz

FUNCTION

COPYRIGHT

  Copyright (C) 2005-2018 ABINIT group (XG, NSAI, MKV)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_elpolariz
28 
29  use defs_basis
30  use defs_datatypes
31  use defs_abitypes
32  use m_abicore
33  use m_errors
34  use m_efield
35  use m_xmpi
36  use m_wffile
37  use m_hdr
38 
39  use m_geometry, only : metric
40  use m_symtk,    only : matr3inv
41  use m_hide_lapack,  only : dzgedi, dzgefa
42  use m_rwwf,     only : rwwf
43  use m_pawtab,   only : pawtab_type
44  use m_pawrhoij, only : pawrhoij_type
45  use m_pawcprj,  only : pawcprj_type
46  use m_berryphase, only : berryphase
47  use m_berryphase_new, only : berryphase_new
48 
49  implicit none
50 
51  private

ABINIT/uderiv [ Functions ]

[ Top ] [ Functions ]

NAME

 uderiv

FUNCTION

 This routine is called computes the derivative of
 ground-state wavefunctions with respect to k (du/dk) by finite differencing
 on neighbouring k points
 Work for nsppol=1 or 2, but only accept nspinor=1,

INPUTS

  bdberry(4)=band limits for Berry phase contributions (or du/dk)
   spin up and spin down (bdberry(3:4) is irrelevant when nsppol=1)
  cg(2,mcg)=planewave coefficients of wavefunctions
  gprimd(3,3)=reciprocal space dimensional primitive translations
  hdr <type(hdr_type)>=the header of wf, den and pot files
  istwfk(nkpt_)=input option parameter that describes the storage of wfs
  kberry(3,20)= different delta k for Berry phases(or du/dk),
   in unit of kptrlatt only kberry(1:3,1:nberry) is relevant
  kg(3,mpw*mkmem)=reduced planewave coordinates
  kpt_(3,nkpt_)=reduced coordinates of k points generated by ABINIT,
               kpt_ samples half the BZ if time-reversal symetrie is used
  kptopt=2 when time-reversal symmetry is used
  kptrlatt(3,3)=k-point lattice specification
  mband=maximum number of bands
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mkmem=number of k points treated by this node.
  mpi_enreg=information about MPI parallelization
  mpw=maximum dimensioned size of npw
  natom=number of atoms in cell
  nband(nkpt*nsppol)=number of bands at each k point, for each polarization
  nberry=number of Berry phases(or du/dk) to be computed
  nkpt=number of k points
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspinor=number of spinorial components of the wavefunctions (on current proc)
  nsppol=1 for unpolarized, 2 for spin-polarized
  unddk=unit number for ddk file

OUTPUT

  (the ddk wavefunctions are written on disk)

SIDE EFFECTS

TODO

  Cleaning, checking for rules
  Should allow for time-reversal symmetry (istwfk)
  WARNING : the use of nspinor is completely erroneous

NOTES

 Local Variables:
  cmatrix(:,:,:)= overlap matrix of size maxband*maxband
  cg_index(:,:,:)= unpacked cg index array for specific band,
   k point and polarization.
  det(2,2)= intermediate output of Lapack routine zgedi.f
  dk(3)= step taken to the next k mesh point along the kberry direction
  gpard(3)= dimensionalreciprocal lattice vector G along which the
          polarization is computed
  kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
          planewave and k point
  kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
  kpt_flag(nkpt)=kpt_flag(ikpt)=0 when the wf was generated by the ABINIT
                 code
                 kpt_flag(ikpt) gives the indices of the k-point
                 related to ikpt by time revers
  maxband/minband= control the minimum and maximum band calculated in the
           overlap matrix
  npw_k= npwarr(ikpt), number of planewaves in basis at this k point
  shift_g_2(nkpt,nkpt)= .true. if the k point should be shifted by a G vector;
  .false. if not
  tr(2)=variable that changes k to -k
                              G to -G
                     $c_g$ to $c_g^*$ when time-reversal symetrie is used

PARENTS

      elpolariz

CHILDREN

      appdig,dzgedi,dzgefa,hdr_io,matr3inv,rwwf,waveformat,wffclose,wffopen
      wrtout,xdefineoff

SOURCE

493 subroutine uderiv(bdberry,cg,gprimd,hdr,istwfk,kberry,kg,kpt_,kptopt,kptrlatt,&
494 & mband,mcg,mkmem,mpi_enreg,mpw,natom,nband,nberry,npwarr,nspinor,nsppol,nkpt_,&
495 & unddk,fnameabo_1wf)
496 
497 
498 !This section has been created automatically by the script Abilint (TD).
499 !Do not modify the following lines by hand.
500 #undef ABI_FUNC
501 #define ABI_FUNC 'uderiv'
502 !End of the abilint section
503 
504  implicit none
505 
506 !Arguments ------------------------------------
507 !scalars
508  integer,intent(in) :: kptopt,mband,mcg,mkmem,mpw,natom,nberry,nkpt_,nspinor
509  integer,intent(in) :: nsppol,unddk
510  type(MPI_type),intent(in) :: mpi_enreg
511  type(hdr_type),intent(inout) :: hdr
512 !arrays
513  integer,intent(in) :: bdberry(4),istwfk(nkpt_),kberry(3,20),kg(3,mpw*mkmem)
514  integer,intent(in) :: kptrlatt(3,3),nband(nkpt_*nsppol),npwarr(nkpt_)
515  real(dp),intent(in) :: cg(2,mcg),gprimd(1:3,1:3)
516  real(dp),intent(in) :: kpt_(3,nkpt_)
517  character(len=fnlen),intent(in) :: fnameabo_1wf
518 
519 !Local variables -------------------------
520 !scalars
521  integer,parameter :: master=0
522  integer :: iomode,band_in,cg_index_iband,fform,flag1
523  integer :: formeig,iband,iberry,icg,idir,ierr,ifor,ii,ikpt,ikpt2,ikpt_
524  integer :: index,index1,info,ipert,ipw,isgn,isppol,jband,jj,jkpt,jkpt_
525  integer :: maxband,mcg_disk,me,minband,nband_diff,nband_k
526  integer :: nkpt,npw_k,pertcase,rdwr,read_k,spaceComm
527  integer :: tim_rwwf
528  real(dp) :: gmod,twodk
529  character(len=500) :: message
530  character(len=fnlen) :: fiwf1o
531  type(wffile_type) :: wffddk
532 !arrays
533  integer :: kg_jl(0,0,0),kpt_flag(2*nkpt_)
534  integer,allocatable :: cg_index(:,:,:),ikpt_dk(:,:),ipvt(:)
535  integer,allocatable :: kg_kpt(:,:,:)
536  real(dp) :: det(2,2),diffk(3),diffk2(3),dk(3),gpard(3),klattice(3,3)
537  real(dp) :: kptrlattr(3,3),tr(2)
538  real(dp) :: cg_disk(0,0,0)
539  real(dp),allocatable :: cmatrix(:,:,:),dudk(:,:)
540  real(dp),allocatable :: eig_dum_2(:),kpt(:,:)
541  real(dp),allocatable :: occ_dum_2(:),phi(:,:,:),u_tilde(:,:,:,:),zgwork(:,:)
542  logical,allocatable :: shift_g_2(:,:)
543 
544 ! *************************************************************************
545 
546  if(min(2,(1+mpi_enreg%paral_spinor)*nspinor)==2)then
547    MSG_ERROR('uderiv: does not yet work for nspinor=2')
548  end if
549 
550  if(maxval(istwfk(:))/=1)then
551    write(message,'(3a)')&
552 &   'Sorry, this routine does not work yet with istwfk/=1.',ch10,&
553 &   'This should have been tested previously ...'
554    MSG_BUG(message)
555  end if
556 
557  if (kptopt==3) then
558    nkpt = nkpt_
559    ABI_ALLOCATE(kpt,(3,nkpt))
560    kpt(:,:)=kpt_(:,:)
561  else if (kptopt==2) then
562    nkpt = nkpt_*2
563    ABI_ALLOCATE(kpt,(3,nkpt))
564    do ikpt = 1,nkpt/2
565      kpt_flag(ikpt) = 0
566      kpt(:,ikpt)=kpt_(:,ikpt)
567    end do
568    index = 0
569    do ikpt = (nkpt/2+1),nkpt
570      flag1 = 0
571      do jkpt = 1, nkpt/2
572        if (((abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt))<1.0d-8).or.&
573 &       (abs(1-abs(kpt_(1,ikpt-nkpt/2)+kpt_(1,jkpt)))<1.0d-8))&
574 &       .and.((abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt))<1.0d-8).or.&
575 &       (abs(1-abs(kpt_(2,ikpt-nkpt/2)+kpt_(2,jkpt)))<1.0d-8))&
576 &       .and.((abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt))<1.0d-8).or.&
577 &       (abs(1-abs(kpt_(3,ikpt-nkpt/2)+kpt_(3,jkpt)))<1.0d-8))) then
578          flag1 = 1
579          index = index + 1
580          exit
581        end if
582      end do
583      if (flag1==0) then
584        kpt_flag(ikpt-index)=ikpt-nkpt/2
585        kpt(:,ikpt-index)=-kpt_(:,ikpt-nkpt/2)
586      end if
587    end do
588    nkpt = nkpt - index
589  end if
590 
591 !DEBUG
592 !write(101,*) 'beginning write kpt'
593 !do ikpt=1,nkpt
594 !write(101,*) kpt(:,ikpt)
595 !end do
596 !ENDDEBUG
597 
598  ABI_ALLOCATE(shift_g_2,(nkpt,nkpt))
599 
600 !Compute primitive vectors of the k point lattice
601 !Copy to real(dp)
602  kptrlattr(:,:)=kptrlatt(:,:)
603 !Go to reciprocal space (in reduced coordinates)
604  call matr3inv(kptrlattr,klattice)
605 
606  do iberry=1,nberry
607 
608 !  **************************************************************************
609 !  Determine the appended index for ddk 1WF files
610 
611    do idir=1,3
612      if (kberry(idir,iberry) ==1) then
613        ipert=natom+1
614        pertcase=idir+(ipert-1)*3
615      end if
616    end do
617 
618 !  open ddk 1WF file
619    formeig=1
620 
621    call appdig(pertcase,fnameabo_1wf,fiwf1o)
622    !call wfk_open_read(wfk, fiwf1o, formeig, iomode, unddk, spaceComm)
623 
624    spaceComm=xmpi_comm_self; me=0 ; iomode=IO_MODE_FORTRAN
625    call WffOpen(iomode,spaceComm,fiwf1o,ierr,wffddk,master,me,unddk)
626 
627    rdwr=2 ; fform=2
628    call hdr_io(fform,hdr,rdwr,wffddk)
629 
630 !  Define offsets, in case of MPI I/O
631    call xdefineOff(formeig,wffddk,mpi_enreg,nband,npwarr,nspinor,nsppol,nkpt_)
632 
633 !  *****************************************************************************
634 !  Calculate dimensional recip lattice vector along which P is calculated
635 !  dk =  step to the nearest k point along that direction
636 !  in reduced coordinates
637 
638    dk(:)=dble(kberry(1,iberry))*klattice(:,1)+&
639 &   dble(kberry(2,iberry))*klattice(:,2)+&
640 &   dble(kberry(3,iberry))*klattice(:,3)
641 
642    do idir=1,3
643      if (dk(idir)/=0) then
644        twodk=2*dk(idir)
645      end if
646    end do
647 
648    gpard(:)=dk(1)*gprimd(:,1)+dk(2)*gprimd(:,2)+dk(3)*gprimd(:,3)
649    gmod=sqrt(dot_product(gpard,gpard))
650 
651 !  ******************************************************************************
652 !  Select the k grid  points along the direction to compute dudk
653 !  dk =  step to the nearest k point along that direction
654 
655 !  For each k point, find k_prim such that k_prim= k + dk mod(G)
656 !  where G is a vector of the reciprocal lattice
657    ABI_ALLOCATE(ikpt_dk,(2,nkpt))
658    ikpt_dk(1:2,1:nkpt)=0
659    shift_g_2(:,:)= .false.
660 
661    do ikpt=1,nkpt
662      do ikpt2=1,nkpt
663        diffk(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)-dk(:))
664        diffk2(:)=abs(kpt(:,ikpt2)-kpt(:,ikpt)+dk(:))
665        if (sum(abs(diffk(:)-nint(diffk(:))))<3*tol8)then
666          ikpt_dk(1,ikpt)=ikpt2
667          if(sum(diffk(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
668        end if
669        if (sum(abs(diffk2(:)-nint(diffk2(:))))<3*tol8)then
670          ikpt_dk(2,ikpt)=ikpt2
671          if(sum(diffk2(:))>=3*tol8) shift_g_2(ikpt,ikpt2) = .true.
672        end if
673      end do
674    end do
675 
676    write(message,'(a,a,a,3f9.5,a,a,3f9.5,a)')ch10,&
677 &   ' Computing the derivative for reciprocal vector:',ch10,&
678 &   dk(:),' (in reduced coordinates)',ch10,&
679 &   gpard(1:3),' (in cartesian coordinates - atomic units)'
680    call wrtout(ab_out,message,'COLL')
681    call wrtout(std_out,message,'COLL')
682 
683    if(nsppol==1)then
684      write(message, '(a,i5,a,i5)')&
685 &     ' From band number',bdberry(1),'  to band number',bdberry(2)
686    else
687      write(message, '(a,i5,a,i5,a,a,a,i5,a,i5,a)')&
688 &     ' From band number',bdberry(1),'  to band number',bdberry(2),' for spin up,',&
689 &     ch10,&
690 &     ' from band number',bdberry(3),'  to band number',bdberry(4),' for spin down.'
691    end if
692    call wrtout(ab_out,message,'COLL')
693    call wrtout(std_out,message,'COLL')
694 
695 !  *****************************************************************************
696    ABI_ALLOCATE(dudk,(2,mpw*nspinor*mband*nsppol))
697    ABI_ALLOCATE(eig_dum_2,((2*mband)**formeig*mband))
698    ABI_ALLOCATE(occ_dum_2,((2*mband)**formeig*mband))
699    dudk(1:2,:)=0.0_dp
700    eig_dum_2=0.0_dp
701    occ_dum_2=0.0_dp
702 
703    if (mkmem/=0) then
704 
705 !    Find the location of each wavefunction
706 
707      ABI_ALLOCATE(cg_index,(mband,nkpt_,nsppol))
708      icg = 0
709      do isppol=1,nsppol
710        do ikpt=1,nkpt_
711          nband_k=nband(ikpt+(isppol-1)*nkpt_)
712          npw_k=npwarr(ikpt)
713          do iband=1,nband_k
714            cg_index(iband,ikpt,isppol)=(iband-1)*npw_k*nspinor+icg
715          end do
716          icg=icg+npw_k*nspinor*nband(ikpt)
717        end do
718      end do
719 
720 !    Find the planewave vectors for each k point
721 !    SHOULD BE REMOVED WHEN ANOTHER INDEXING TECHNIQUE WILL BE USED FOR kg
722      ABI_ALLOCATE(kg_kpt,(3,mpw*nspinor,nkpt_))
723      kg_kpt(:,:,:) = 0
724      index1 = 0
725      do ikpt=1,nkpt_
726        npw_k=npwarr(ikpt)
727        do ipw=1,npw_k*nspinor
728          kg_kpt(1:3,ipw,ikpt)=kg(1:3,ipw+index1)
729        end do
730        index1=index1+npw_k*nspinor
731      end do
732    end if
733 
734 !  *************************************************************************
735 !  Loop over spins
736    do isppol=1,nsppol
737 
738      minband=bdberry(2*isppol-1)
739      maxband=bdberry(2*isppol)
740 
741      if(minband<1)then
742        write(message,'(a,i0,a)')'  The band limit minband= ',minband,', is lower than 0.'
743        MSG_BUG(message)
744      end if
745 
746      if(maxband<1)then
747        write(message,'(a,i0,a)')' The band limit maxband= ',maxband,', is lower than 0.'
748        MSG_BUG(message)
749      end if
750 
751      if(maxband<minband)then
752        write(message,'(a,i0,a,i0)')' maxband= ',maxband,', is lower than minband= ',minband
753        MSG_BUG(message)
754      end if
755 
756 !    Loop over k points
757      do ikpt_=1,nkpt_
758 
759        read_k = 0
760 
761        ikpt=ikpt_
762        tr(1) = 1.0_dp
763 
764        if (kptopt==2) then
765          if (read_k == 0) then
766            if (kpt_flag(ikpt_)/=0) then
767              tr(1) = -1.0_dp
768              ikpt= kpt_flag(ikpt_)
769            end if
770          else           !read_k
771            if (kpt_flag(ikpt_)/=0) then
772              tr(-1*read_k+3) = -1.0_dp
773              ikpt= kpt_flag(ikpt_)
774            end if
775          end if       !read_k
776        end if           !kptopt
777 
778        nband_k=nband(ikpt+(isppol-1)*nkpt_)
779 
780        if(nband_k<maxband)then
781          write(message,'(a,i0,a,i0)')'  maxband=',maxband,', is larger than nband(i,isppol)=',nband_k
782          MSG_BUG(message)
783        end if
784 
785        npw_k=npwarr(ikpt)
786 
787        ABI_ALLOCATE(u_tilde,(2,npw_k,maxband,2))
788        u_tilde(1:2,1:npw_k,1:maxband,1:2)=0.0_dp
789 
790 !      ifor = 1,2 represents forward and backward neighbouring k points of ikpt
791 !      respectively along dk direction
792 
793        do ifor=1,2
794 
795          ABI_ALLOCATE(phi,(2,mpw,mband))
796          ABI_ALLOCATE(cmatrix,(2,maxband,maxband))
797          phi(1:2,1:mpw,1:mband)=0.0_dp; cmatrix(1:2,1:maxband,1:maxband)=0.0_dp
798 
799          isgn=(-1)**ifor
800          jkpt_= ikpt_dk(ifor,ikpt_)
801 
802          tr(2) = 1.0_dp
803 
804          jkpt=jkpt_
805 
806          if (kptopt==2) then
807            if (read_k == 0) then
808              if (kpt_flag(jkpt_)/=0) then
809                tr(2) = -1.0_dp
810                jkpt= kpt_flag(jkpt_)
811              end if
812            else           !read_k
813              if (kpt_flag(jkpt_)/=0) then
814                tr(read_k) = -1.0_dp
815                jkpt= kpt_flag(jkpt_)
816              end if
817            end if       !read_k
818          end if           !kptopt
819 
820          if (ifor==1) read_k = 2
821 
822          jj = read_k
823          ii = -1*read_k+3
824 
825          call waveformat(cg,cg_disk,cg_index,phi,dk,ii,ikpt,&
826 &         ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,&
827 &         minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr)
828 
829 !        Compute the overlap matrix <u_k|u_k+b>
830 
831          do iband=minband,maxband
832            cg_index_iband=cg_index(iband,ikpt,isppol)
833            do jband=minband,maxband
834              do ipw=1,npwarr(ikpt)
835                cmatrix(1,iband,jband)=cmatrix(1,iband,jband)+&
836 &               cg(1,ipw+cg_index_iband)*phi(1,ipw,jband)+&
837 &               tr(ii)*cg(2,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)
838 
839                cmatrix(2,iband,jband)=cmatrix(2,iband,jband)+&
840 &               cg(1,ipw+cg_index_iband)*tr(jj)*phi(2,ipw,jband)-&
841 &               tr(ii)*cg(2,ipw+cg_index_iband)*phi(1,ipw,jband)
842              end do
843            end do
844          end do
845 
846 !        Compute the inverse of cmatrix(1:2,minband:maxband, minband:maxband)
847 
848          band_in = maxband - minband + 1
849          ABI_ALLOCATE(ipvt,(maxband))
850          ABI_ALLOCATE(zgwork,(2,1:maxband))
851 
852 !        Last argument of zgedi means calculate inverse only
853          call dzgefa(cmatrix(1,minband,minband),maxband, band_in,ipvt,info)
854          call dzgedi(cmatrix(1,minband,minband),maxband, band_in,ipvt,det,zgwork,01)
855 
856          ABI_DEALLOCATE(zgwork)
857          ABI_DEALLOCATE(ipvt)
858 
859 !        Compute the product of Inverse overlap matrix with the wavefunction
860 
861          do iband=minband,maxband
862            do ipw=1,npwarr(ikpt)
863              u_tilde(1,ipw,iband,ifor)= &
864 &             dot_product(cmatrix(1,minband:maxband,iband),&
865 &             phi(1,ipw,minband:maxband))-&
866 &             dot_product(cmatrix(2,minband:maxband,iband),&
867 &             tr(jj)*phi(2,ipw,minband:maxband))
868              u_tilde(2,ipw,iband,ifor)= &
869 &             dot_product(cmatrix(1,minband:maxband,iband),&
870 &             tr(jj)*phi(2,ipw,minband:maxband))+&
871 &             dot_product(cmatrix(2,minband:maxband,iband),&
872 &             phi(1,ipw,minband:maxband))
873            end do
874          end do
875          ABI_DEALLOCATE(cmatrix)
876          ABI_DEALLOCATE(phi)
877 
878        end do !ifor
879 
880 !      Compute dudk for ikpt
881 
882        npw_k=npwarr(ikpt)
883 
884        do iband=minband,maxband
885 
886          icg=(iband-minband)*npw_k
887 
888          dudk(1,1+icg:npw_k+icg)=(u_tilde(1,1:npw_k,iband,1)-&
889 &         u_tilde(1,1:npw_k,iband,2))/twodk
890 
891          dudk(2,1+icg:npw_k+icg)=(u_tilde(2,1:npw_k,iband,1)-&
892 &         u_tilde(2,1:npw_k,iband,2))/twodk
893 
894        end do
895 
896        tim_rwwf=0
897        mcg_disk=mpw*nspinor*mband
898        nband_diff=maxband-minband+1
899        call rwwf(dudk,eig_dum_2,formeig,0,0,ikpt,isppol,kg_kpt(:,:,ikpt),&
900 &       mband,mcg_disk,mpi_enreg,nband_diff,nband_diff,&
901 &       npw_k,nspinor,occ_dum_2,2,1,tim_rwwf,wffddk)
902 
903        !call wfk_read_band_block(wfk, band_block, ikpt, isppol, sc_mode,
904        !  kg_k=kg_kpt(:,:,ikpt), cg_k=dudk, eig_k=eig_dum, occ_k=occ_dum)
905 
906        ABI_DEALLOCATE(u_tilde)
907 
908      end do !ikpt
909    end do  !isppol
910 
911    ABI_DEALLOCATE(eig_dum_2)
912    ABI_DEALLOCATE(occ_dum_2)
913    ABI_DEALLOCATE(dudk)
914 
915    call WffClose(wffddk,ierr)
916    !call wfk_close(wfk)
917 
918    ABI_DEALLOCATE(kg_kpt)
919    ABI_DEALLOCATE(cg_index)
920    ABI_DEALLOCATE(ikpt_dk)
921 
922  end do ! iberry
923 
924  ABI_DEALLOCATE(shift_g_2)
925  ABI_DEALLOCATE(kpt)
926 
927  write(std_out,*) 'uderiv:  exit '
928 
929 end subroutine uderiv

ABINIT/waveformat [ Functions ]

[ Top ] [ Functions ]

NAME

 waveformat

FUNCTION

 This routine is to find the matched pairs of plane waves between
 two neighbouring k points and load a new pw coefficients array cg_new
 Was written first by Na Sai (thanks), but unfortunately without
 any comment ...

INPUTS

  cg(2,mpw*nspinor*mband*mkmem*nsppol)= input planewave coefficients, in case mkmem/=0
  cg_disk(2,mpw*nspinor*mband,2)= input planewave coefficients, in case mkmem==0
  cg_index(mband,nkpt_,nsppol)=index of wavefunction iband,ikpt,isppol in the array cg.
  dk(3)= step taken to the next k mesh point along the kberry direction (see also isgn)
  ii=(to be documented)
  ikpt=index of the first k-point in the reduced Brillouin zone
  ikpt_=index of the first k-point in the full Brillouin zone
  isgn=1 if dk(3) is connecting the k-points (ikpt_ and jkpt)
      =-1 if -dk(3) is connecting the k-points
  isppol=1 if spin-up, =2 if spin-down
  jj=(to be documented)
  jkpt=index of the second k-point in the reduced Brillouin zone
  jkpt_=index of the second k-point in the full Brillouin zone
  kg_kpt(:,:,:)= unpacked reduced planewave coordinates with subscript of
          planewave and k point
  kpt(3,nkpt)=reduced coordinates of k-point grid that samples the whole BZ
  kg_jl(3,mpw,2)=(to be documented)
  maxband/minband= control the minimum and maximum band calculated in the
           overlap matrix
  mband=maximum number of bands (dimension of several cg* arrays)
  mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
  mcg_disk=size of wave-functions array (cg_disk) =mpw*nspinor*mband
  mkmem= if 0, the wavefunctions are input in cg_disk, otherwise in cg
  mpw=maximum number of planewaves (dimension of several cg* arrays)
  nkpt=number of k points (full Brillouin zone !?!)
  nkpt_=number of k points (reduced Brillouin zone !?!)
  npwarr(nkpt)=number of planewaves in basis at this k point
  nspinor=number of spinorial components of the wavefunctions
  nsppol=1 for unpolarized, 2 for spin-polarized
  shift_g_2(nkpt,nkpt)=non-zero if a G vector along kberry is needed to connect k points
  tr(2)=variable that changes k to -k
                              G to -G
                     $c_g$ to $c_g^*$ when time-reversal symetrie is used

OUTPUT

  cg_new(2,mpw,maxband)=planewave coefficients transferred onto the
   set of planewaves at k

SIDE EFFECTS

PARENTS

      uderiv

CHILDREN

SOURCE

 991 subroutine waveformat(cg,cg_disk,cg_index,cg_new,dk,ii,ikpt,&
 992 & ikpt_,isgn,isppol,jj,jkpt,jkpt_,kg_kpt,kpt,kg_jl,maxband,mband,mcg,mcg_disk,&
 993 & minband,mkmem,mpw,nkpt,nkpt_,npwarr,nsppol,nspinor,shift_g_2,tr)
 994 
 995 
 996 !This section has been created automatically by the script Abilint (TD).
 997 !Do not modify the following lines by hand.
 998 #undef ABI_FUNC
 999 #define ABI_FUNC 'waveformat'
1000 !End of the abilint section
1001 
1002  implicit none
1003 
1004 !Arguments ------------------------------------
1005 !scalars
1006  integer,intent(in) :: ii,ikpt,ikpt_,isgn,isppol,jj,jkpt,jkpt_,maxband,mband,mcg,mcg_disk
1007  integer,intent(in) :: minband,mkmem,mpw,nkpt,nkpt_,nspinor,nsppol
1008 !arrays
1009  integer,intent(in) :: cg_index(mband,nkpt_,nsppol),kg_jl(3,mpw,2)
1010  integer,intent(in) :: kg_kpt(3,mpw*nspinor,nkpt_),npwarr(nkpt_)
1011  real(dp),intent(in) :: cg(2,mcg)
1012  real(dp),intent(in) :: cg_disk(2,mcg_disk,2),dk(3),kpt(3,nkpt),tr(2)
1013  real(dp),intent(out) :: cg_new(2,mpw,maxband)
1014  logical,intent(in) :: shift_g_2(nkpt,nkpt)
1015 
1016 !Local variables -------------------------
1017 !scalars
1018  integer :: cg_index_iband,iband,ipw,jpw,nomatch,npw_k
1019  logical :: found_match
1020 !arrays
1021  integer :: dg(3)
1022 
1023 ! ***********************************************************************
1024 
1025  npw_k=npwarr(ikpt)
1026 
1027 
1028  nomatch=0
1029 
1030 !If there is no shift of G-vector between ikpt_ and jkpt_
1031  if(shift_g_2(ikpt_,jkpt_) .eqv. .false.) then
1032 
1033 !  DEBUG
1034 !  write(111,*)'pair', ikpt_,jkpt_,'noshift'
1035 !  ENDDEBUG
1036 
1037 !  If the original wavefunction is contained in cg_disk
1038    if(mkmem==0) then
1039 
1040      do ipw=1,npw_k
1041 
1042        found_match = .false.
1043 
1044        do jpw=1,npwarr(jkpt)
1045          if (sum(abs(tr(ii)*kg_jl(:,ipw,ii)-tr(jj)*kg_jl(:,jpw,jj)))<3*tol8)then
1046            do iband=minband, maxband
1047              cg_index_iband=(iband-1)*npwarr(jkpt)
1048              cg_new(1:2,ipw,iband)=cg_disk(1:2,jpw+cg_index_iband,jj)
1049            end do
1050            found_match = .true.
1051            exit
1052          end if
1053        end do
1054 
1055        if (found_match .eqv. .false.) then
1056          do iband=minband,maxband
1057            cg_new(1:2,ipw,iband)=zero
1058          end do
1059          nomatch = nomatch + 1
1060        end if
1061 
1062      end do
1063 
1064 !    Here, the wavefunctions are contained in cg
1065    else
1066 
1067      do ipw=1,npw_k
1068 
1069        found_match = .false.
1070 
1071        do jpw=1,npwarr(jkpt)
1072          if (sum(abs(tr(ii)*kg_kpt(:,ipw,ikpt)-tr(jj)*kg_kpt(:,jpw,jkpt)))<3*tol8)then
1073            do iband=minband, maxband
1074              cg_index_iband=cg_index(iband,jkpt,isppol)
1075              cg_new(1:2,ipw,iband)=cg(1:2,jpw+cg_index_iband)
1076            end do
1077            found_match = .true.
1078            exit
1079          end if
1080        end do
1081 
1082        if (found_match .eqv. .false.) then
1083          do iband=minband,maxband
1084            cg_new(1:2,ipw,iband)=(0.0_dp,0.0_dp)
1085          end do
1086          nomatch = nomatch + 1
1087        end if
1088      end do
1089 
1090    end if
1091 
1092 !  DEBUG
1093 !  write(111,*) 'normal pair nomatch=',nomatch
1094 !  ENDDEBUG
1095 
1096 !  If there is a G-vector shift between ikpt_ and jkpt_
1097  else
1098 
1099 !  DEBUG
1100 !  write(111,*) 'pair',ikpt_,jkpt_,' need shift'
1101 !  ENDDEBUG
1102 
1103    dg(:) = -1*nint(tr(jj)*kpt(:,jkpt)-tr(ii)*kpt(:,ikpt)+isgn*dk(:))
1104 
1105 !  If the original wavefunction is contained in cg_disk
1106    if(mkmem==0) then
1107 
1108      do ipw=1,npw_k
1109 
1110        found_match = .false.
1111 
1112        do jpw=1,npwarr(jkpt)
1113          if (sum(abs(tr(ii)*kg_jl(:,ipw,ii)-(tr(jj)*kg_jl(:,jpw,jj)-&
1114 &         dg(:))))<3*tol8)then
1115 
1116            do iband=minband, maxband
1117              cg_index_iband=(iband-1)*npwarr(jkpt)
1118              cg_new(1:2,ipw,iband)=cg_disk(1:2,jpw+cg_index_iband,jj)
1119            end do
1120            found_match = .true.
1121            exit
1122          end if
1123        end do
1124 
1125        if (found_match .eqv. .false.) then
1126          do iband=minband,maxband
1127            cg_new(1:2,ipw,iband)=(0.0_dp,0.0_dp)
1128          end do
1129          nomatch = nomatch + 1
1130        end if
1131      end do
1132 
1133 !    Here, the wavefunctions are contained in cg
1134    else
1135 
1136      do ipw=1,npw_k
1137 
1138        found_match = .false.
1139 
1140        do jpw=1,npwarr(jkpt)
1141          if (sum(abs(tr(ii)*kg_kpt(:,ipw,ikpt)-(tr(jj)*kg_kpt(:,jpw,jkpt)-&
1142 &         dg(:))))<3*tol8)then
1143            do iband=minband, maxband
1144              cg_index_iband=cg_index(iband,jkpt,isppol)
1145              cg_new(1:2,ipw,iband)=cg(1:2,jpw+cg_index_iband)
1146            end do
1147            found_match = .true.
1148            exit
1149          end if
1150        end do
1151 
1152        if (found_match .eqv. .false.) then
1153          do iband=minband,maxband
1154            cg_new(1:2,ipw,iband)=zero
1155          end do
1156          nomatch = nomatch + 1
1157        end if
1158      end do
1159 
1160    end if
1161 
1162 !  DEBUG
1163 !  write(111,*) 'special pair nomatch=',nomatch
1164 !  ENDDEBUG
1165 
1166  end if
1167 
1168 end subroutine waveformat