TABLE OF CONTENTS


ABINIT/interpolate_gkk [ Functions ]

[ Top ] [ Functions ]

NAME

 interpolate_gkk

FUNCTION

 This routine interpolates the gkk matrices for all q vectors
  between points on the full kpt_phon grid.

COPYRIGHT

 Copyright (C) 2004-2018 ABINIT group (MVer)
 This file is distributed under the terms of the
 GNU General Public Licence, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

INPUTS

   elph_ds = elphon datastructure with data and dimensions
   kpt_phon = coordinates of all kpoints close to the FS

OUTPUT

   elph_ds = modified gkq

NOTES

  inspired to some extent by epcouple.f from the DecAFT package by J. Kay Dewhurst
  most inputs taken from mkifc.f
  in anaddb set ifcflag 1 such that the IFC are calculated in atmfrc prior to calling elphon

PARENTS

      get_all_gkk2

CHILDREN

      ftgkk,ifc_fourq,wrap2_pmhalf,zhpev

SOURCE

 38 #if defined HAVE_CONFIG_H
 39 #include "config.h"
 40 #endif
 41 
 42 #include "abi_common.h"
 43 
 44 
 45 subroutine interpolate_gkk(crystal,ifc,elph_ds,kpt_phon)
 46 
 47  use defs_basis
 48  use defs_elphon
 49  use m_profiling_abi
 50  use m_errors
 51 
 52  use m_numeric_tools,   only : wrap2_pmhalf
 53  use m_crystal,         only : crystal_t
 54  use m_ifc,             only : ifc_type, ifc_fourq
 55 
 56 !This section has been created automatically by the script Abilint (TD).
 57 !Do not modify the following lines by hand.
 58 #undef ABI_FUNC
 59 #define ABI_FUNC 'interpolate_gkk'
 60  use interfaces_77_ddb, except_this_one => interpolate_gkk
 61 !End of the abilint section
 62 
 63  implicit none
 64 
 65 !Arguments ------------------------------------
 66 !scalars
 67  type(crystal_t),intent(in) :: crystal
 68  type(ifc_type),intent(in) :: ifc
 69  type(elph_type),intent(inout) :: elph_ds
 70 !arrays
 71  real(dp),intent(in) :: kpt_phon(3,elph_ds%k_phon%nkpt)
 72 
 73 !Local variables-------------------------------
 74   ! output variables for dfpt_phfrq
 75 ! variables for zhpev
 76 ! variables for phonon interpolation
 77 !scalars
 78  integer :: i1,i2,ikpt_phon2,iFSqpt,ib1,ib2,ier,ii
 79  integer :: iost,isppol,qtor,natom
 80  integer :: sz1,sz2,sz3,sz4,unit_gkkp
 81  real(dp) :: qphnrm,res
 82  !character(len=500) :: msg
 83 !arrays
 84  real(dp) :: gprim(3,3)
 85  real(dp) :: displ(2,elph_ds%nbranch,elph_ds%nbranch),eigval(3*crystal%natom)
 86  real(dp) :: eigvec(3*3*crystal%natom*3*crystal%natom)
 87  real(dp) :: pheigvec(2*elph_ds%nbranch*elph_ds%nbranch)
 88  real(dp) :: phfrq_tmp(elph_ds%nbranch),qphon(3),redkpt(3)
 89  real(dp),allocatable :: gkk2_diag_tmp(:,:,:,:),gkk2_tmp(:,:,:,:,:,:,:)
 90  real(dp),allocatable :: matrx(:,:),zhpev1(:,:)
 91  real(dp),allocatable :: zhpev2(:)
 92 
 93 ! *************************************************************************
 94 
 95 !
 96 !NOTE: mjv 18/5/2008 reverted to old style of ftgkk with all kpt done together.
 97 !may want to modify this later to use the new cleaner format with 1 FT at a
 98 !time.
 99 !
100  write(std_out,*) 'interpolate_gkk : enter'
101 
102  natom = crystal%natom
103  gprim = ifc%gprim
104 
105  if (elph_ds%nsppol /= 1) then
106    MSG_ERROR("interpolate_gkk not coded with nsppol>1 yet")
107  end if
108  isppol = 1
109 
110 
111 !------------------------------------------------------
112 !complete dynamical matrices for all qpts between points
113 !on full kpt grid (interpolation from IFC)
114 !------------------------------------------------------
115 
116  sz1=elph_ds%ngkkband;sz2=elph_ds%nbranch
117  sz3=elph_ds%k_phon%nkpt;sz4=elph_ds%nFSband
118 !allocate (gkk_tmp(2,sz1,sz1,sz2,sz2,1,1))
119 !DEBUG
120 !allocate (gkk_tmp_full(2,sz1,sz1,sz2,elph_ds%nFSband,sz3))
121 !allocate (gkk_tmp_full(2,s2,sz4,sz4,sz3))
122 !ENDDEBUG
123  ABI_ALLOCATE(gkk2_tmp,(2,sz1,sz1,sz2,sz2,sz3,1))
124  ABI_ALLOCATE(gkk2_diag_tmp,(sz1,sz1,sz2,sz3))
125  ABI_ALLOCATE(zhpev1,(2,2*3*natom-1))
126  ABI_ALLOCATE(zhpev2,(3*3*natom-2))
127  ABI_ALLOCATE(matrx,(2,(3*natom*(3*natom+1))/2))
128 
129  qphnrm = one
130 !in this part use the inverse Fourier transform to get 1 (arbitrary) qpt at a
131 !time
132  ii = 0
133  qtor = 0
134  unit_gkkp = 150
135  open (unit=unit_gkkp,file='gkkp_file_ascii',form='formatted',status='unknown',iostat=iost)
136  if (iost /= 0) then
137    MSG_ERROR("error opening gkkpfile as new")
138  end if
139 
140 !loop over all FS pairs.
141 !do ikpt1=1,elph_ds%k_phon%nkptirr
142 !do iFSqpt=1,elph_ds%k_phon%nkpt
143 
144 !
145 !this should run through the sparse mesh of 2x2x2 kpoints
146 !
147  do iFSqpt=1,elph_ds%k_phon%nkpt
148    res = 2.0_dp*(kpt_phon(1,iFSqpt)+one)
149    if (abs(res-int(res)) > tol10) cycle
150    res = 2.0_dp*(kpt_phon(2,iFSqpt)+one)
151    if (abs(res-int(res)) > tol10) cycle
152    res = 2.0_dp*(kpt_phon(3,iFSqpt)+one)
153    if (abs(res-int(res)) > tol10) cycle
154 
155 !  do ikpt1=1,1
156 !  
157 !  NOTE: should be very easy to parallelize!
158 !  
159 !  write(std_out,*) ' interpolate_gkk : ikpt1 = ',ikpt1, ' / ', elph_ds%k_phon%nkptirr
160    write(std_out,*) ' interpolate_gkk : ikpt1 = ',iFSqpt, ' / ', elph_ds%k_phon%nkpt
161 
162 !  DEBUG
163 !  write(std_out,*) ' interpolate_gkk : Warning debug version'
164 !  cycle
165 !  ENDDEBUG
166 
167    gkk2_tmp(:,:,:,:,:,:,:) = zero
168 
169 !  qphon = 1 - 2    ie.  1 = 2+qphon
170    qphon(:) = kpt_phon(:,iFSqpt)
171 
172 !  shouldnt be necessary here, but oh well
173    call wrap2_pmhalf(qphon(1),redkpt(1),res)
174    call wrap2_pmhalf(qphon(2),redkpt(2),res)
175    call wrap2_pmhalf(qphon(3),redkpt(3),res)
176 
177    qphon(:) = redkpt(:)
178    redkpt(1) = qphon(1)*gprim(1,1)+qphon(2)*gprim(1,2)+qphon(3)*gprim(1,3)
179    redkpt(2) = qphon(1)*gprim(2,1)+qphon(2)*gprim(2,2)+qphon(3)*gprim(2,3)
180    redkpt(3) = qphon(1)*gprim(3,1)+qphon(2)*gprim(3,2)+qphon(3)*gprim(3,3)
181    write (unit_gkkp,*) 'qp= ', redkpt
182 
183    call ifc_fourq(ifc,crystal,qphon,phfrq_tmp,displ,out_eigvec=pheigvec)
184    write (unit_gkkp,*) phfrq_tmp(:)*Ha_cmm1
185 
186    ii = ii+1
187 !  if(ii > 0 .and. ii < 1000) write(std_out,'(a,i5,3E16.6,2x)') &
188 !  &   ' wrote phfrq_tmp for time ', ii, phfrq_tmp
189 !  end if
190 
191 !  phonon eigenvectors are in eigvec
192 !  real and imaginary parts
193 !  phonon displacements = eigvec/sqrt(M_i) are in displ
194 !  real and imaginary parts
195 
196 !  DEBUG
197 !  test: uniform phonon frequency
198 !  phfrq_tmp(:) = 0.0001_dp
199 !  ENDDEBUG
200 
201 !  FT gamma matrices for all kpt_phon points, and
202 !  for qpoint = qphon(:) = kpt_phon(ikpt_phon)
203 
204    call ftgkk(ifc%wghatm,gkk2_tmp,elph_ds%gkk_rpt,elph_ds%gkqwrite,&
205 &   elph_ds%gkk_rptwrite,gprim,1,&
206 &   natom,elph_ds%k_phon%nkpt,elph_ds%ngkkband,elph_ds%k_phon%nkpt,1,ifc%nrpt,elph_ds%nsppol,&
207 &   qtor,ifc%rpt,qphon,elph_ds%unit_gkk_rpt,elph_ds%unitgkq)
208 
209 !  NOTE: Normally the eigenvectors of the gkk2_tmp should be the same as eigvec
210 
211 !  Diagonalize gamma matrices at qpoint (complex matrix) for all kpt_phon.
212 !  Copied from dfpt_phfrq
213    do ikpt_phon2=1,elph_ds%k_phon%nkpt
214      res = 8.0_dp*(kpt_phon(1,ikpt_phon2)+one)
215      if (abs(res-int(res)) > tol10) cycle
216      res = 8.0_dp*(kpt_phon(2,ikpt_phon2)+one)
217      if (abs(res-int(res)) > tol10) cycle
218      res = 8.0_dp*(kpt_phon(3,ikpt_phon2)+one)
219      if (abs(res-int(res)) > tol10) cycle
220 
221      write (unit_gkkp,*) 'kp= ', kpt_phon(:,ikpt_phon2)
222 
223      do ib1=1,elph_ds%ngkkband
224        do ib2=1,elph_ds%ngkkband
225          ier=0
226          ii=1
227          do i2=1,3*natom
228            do i1=1,i2
229              matrx(1,ii)=gkk2_tmp(1,ib1,ib2,i1,i2,ikpt_phon2,1)
230              matrx(2,ii)=gkk2_tmp(2,ib1,ib2,i1,i2,ikpt_phon2,1)
231              ii=ii+1
232            end do
233          end do
234          call ZHPEV ('N','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,&
235 &         zhpev2,ier)
236 
237          gkk2_diag_tmp(ib2,ib1,:,ikpt_phon2) = eigval(:)
238          do i1=1,3*natom
239            write (unit_gkkp,*) elph_ds%minFSband-1+ib1,elph_ds%minFSband-1+ib2,i1,&
240 &           eigval(i1)
241          end do
242        end do
243      end do
244    end do
245 
246    if (elph_ds%gkk2write == 1) then
247      write(std_out,*) 'WARNING COMMENTED WRITE TO BINARY FILE!!!'
248 !    write (elph_ds%unit_gkk2,REC=iFSqpt) gkk2_diag_tmp(:,:,:,:)
249      write(std_out,'(a,i4,4(2E16.6,2x))') ' gkk2 loop ', &
250 &     iFSqpt,gkk2_diag_tmp(1,1,:,1:2),gkk2_diag_tmp(1,1,:,elph_ds%k_phon%nkpt-1:elph_ds%k_phon%nkpt)
251 !    &    ikpt1,gkk2_tmp(:,1,1,1,1,1:2),gkk2_tmp(:,1,1,elph_ds%k_phon%nkpt-1:elph_ds%k_phon%nkpt)
252    else if (elph_ds%gkk2write == 0) then
253      elph_ds%gkk2(:,:,:,:,iFSqpt,isppol) = gkk2_diag_tmp(:,:,:,:)
254 !    elph_ds%gkk2(:,:,:,:,ikpt1) = gkk2_tmp
255      write(std_out,*) ' interpolate_gkk : gkk2(b=1,b=1,:,kpt=1,iFSqpt) = '
256      write(std_out,*) gkk2_diag_tmp(1,1,:,1)
257    end if
258 
259  end do
260 !end do on iFSqpt
261 
262  ABI_DEALLOCATE(matrx)
263  ABI_DEALLOCATE(zhpev1)
264  ABI_DEALLOCATE(zhpev2)
265 
266 end subroutine interpolate_gkk