TABLE OF CONTENTS


ABINIT/m_nesting [ Modules ]

[ Top ] [ Modules ]

NAME

  m_nesting

FUNCTION

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (MG, MJV)
 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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 MODULE m_nesting
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_kptrank
32  use m_sort
33 
34  use m_numeric_tools,  only : wrap2_zero_one, interpol3d
35  use m_io_tools,       only : open_file
36  use m_bz_mesh,        only : make_path
37  use m_pptools,        only : printxsf
38 
39  implicit none
40 
41  private
42 
43  public :: bfactor
44  public :: mknesting
45  public :: outnesting

m_nesting/bfactor [ Functions ]

[ Top ] [ m_nesting ] [ Functions ]

NAME

 bfactor

FUNCTION

 Calculate the nesting factor

INPUTS

  nkptfull = number of k-points in full grid
  kptfull(3,nkptfull) = k-point grid
  nqpt = number of qpoints
  qpt(3,nqpt) = q-point grid (must be a subgrid of the k grid),
                the nesting factor will be calculated for each q point in this array
  nkpt = eventually reduced number of k-points
  weight(nband,nkpt) =  integration weights for each k-point and band (NOT NORMALISED!!!)
  nband = number of bands

OUTPUT

  nestfactor(nqpt) = array containing the nesting factor values

NOTES

 Inspired to nmsq_gam_sumfs and mkqptequiv
  TODO : better use of symmetries to reduce the computational effort
 Must be called with kpt = full grid! Reduction by symmetry is not possible for q-dependent quantities (or not easy :)

PARENTS

      m_nesting,outelph

CHILDREN

      make_path,printxsf,wrap2_zero_one

SOURCE

 86 subroutine bfactor(nkptfull,kptfull,nqpt,qpt,kptrank_t,nkpt,weight,nband,nestfactor)
 87 
 88 
 89 !This section has been created automatically by the script Abilint (TD).
 90 !Do not modify the following lines by hand.
 91 #undef ABI_FUNC
 92 #define ABI_FUNC 'bfactor'
 93 !End of the abilint section
 94 
 95  implicit none
 96 
 97 !Arguments ------------------------------------
 98 !scalars
 99  integer,intent(in) :: nband,nkptfull,nqpt,nkpt
100 !arrays
101  real(dp),intent(in) :: kptfull(3,nkptfull),qpt(3,nqpt),weight(nband,nkpt)
102  real(dp),intent(out) :: nestfactor(nqpt)
103  type(kptrank_type), intent(in) :: kptrank_t
104 
105 !Local variables-------------------------------
106 !scalars
107  integer :: ib1,ib2,ikplusq_irr,ikpt
108  integer :: irank_kpt,ikpt_irr,iqpt,symrank_kpt
109  real(dp) :: w1,w2
110  !character(len=500) :: message
111 !arrays
112  real(dp) :: kptpq(3)
113 
114 ! *************************************************************************
115 
116  nestfactor(:)=zero
117 
118  do iqpt=1,nqpt
119    do ikpt=1,nkptfull
120      call get_rank_1kpt (kptfull(:,ikpt),irank_kpt,kptrank_t)
121      ikpt_irr = kptrank_t%invrank(irank_kpt)
122 
123      kptpq(:) = kptfull(:,ikpt) + qpt(:,iqpt)
124      call get_rank_1kpt (kptpq,symrank_kpt,kptrank_t)
125 
126      ikplusq_irr = kptrank_t%invrank(symrank_kpt)
127      if (ikplusq_irr == -1) then
128        MSG_ERROR('It looks like no kpoint equiv to k+q!')
129      end if
130 
131      do ib1=1,nband
132        w1 = weight(ib1,ikpt_irr) !weight for distance from the Fermi surface
133        if (w1 < tol6 ) cycle
134        do ib2=1,nband
135          w2 = weight(ib2,ikplusq_irr) !weight for distance from the Fermi surface
136          if (w1 < tol6 ) cycle
137          nestfactor(iqpt) = nestfactor(iqpt) + w1*w2
138        end do !ib2
139      end do !ib1
140 
141    end do !ikpt
142  end do !iqpt
143 
144 !need prefactor of (1/nkptfull) for normalisation of integration
145  nestfactor(:) = (one/nkptfull) * nestfactor(:)
146 
147 end subroutine bfactor

m_nesting/mknesting [ Functions ]

[ Top ] [ m_nesting ] [ Functions ]

NAME

 mknesting

FUNCTION

  Calculate the nesting factor over the dense k-grid, interpolate the values along a given q path
  and write the data on file in the X-Y format or in the XCrysden format (XSF)

INPUTS

  nkpt = number of k points
  kpt(3,nkpt) = k points
  nkx, nky, nkz = number of k-point along each direction
  nband = number of bands to be considered in the calculation
  weight(nband,nkpt) =  integration weights for each k-point and band
  nqpath = number of points requested along the trajectory
  qpath_vertices = vertices of the reciprocal space trajectory
  base_name = prefix of the output file
  gprimd(3,3) dimensional reciprocal lattice vectors
  gmet = metric in reciprocal space
  prtnest = flags governing the format of the output file

OUTPUT

   Write data to file.

PARENTS

      elphon,m_ebands

CHILDREN

      make_path,printxsf,wrap2_zero_one

SOURCE

184 subroutine mknesting(nkpt,kpt,kptrlatt,nband,weight,nqpath,&
185 & qpath_vertices,nqptfull,qptfull,base_name,gprimd,gmet,prtnest,qptrlatt,&
186 & nsym,symrec) ! optional
187 
188 
189 !This section has been created automatically by the script Abilint (TD).
190 !Do not modify the following lines by hand.
191 #undef ABI_FUNC
192 #define ABI_FUNC 'mknesting'
193 !End of the abilint section
194 
195  implicit none
196 
197 !Arguments ------------------------------------
198 !scalars
199  integer,intent(in) :: nband,nkpt,nqpath,prtnest
200  integer, intent(in) :: nqptfull
201  integer, intent(in), optional :: nsym
202  character(len=*),intent(in) :: base_name
203 !arrays
204  integer,intent(in) :: kptrlatt(3,3)
205  integer,intent(in),optional :: symrec(3,3,*)
206  real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt)
207  real(dp),intent(in) :: qptfull(3,nqptfull)
208  real(dp),intent(in) :: gmet(3,3)
209  real(dp),intent(in) :: qpath_vertices(3,nqpath)
210  real(dp),intent(in) :: weight(nband,nkpt)
211  integer,intent(in)  :: qptrlatt(3,3)
212 
213 !Local variables-------------------------------
214 !scalars
215  integer :: ikpt,jkpt
216  integer :: ik1, ik2, ik3, nkptfull
217  character(len=500) :: message
218  type(kptrank_type) :: kptrank_t
219 !arrays
220  integer,allocatable :: tmprank(:),ktable(:)
221  character(len=fnlen) :: tmpname
222  real(dp),allocatable :: nestfactor(:),nestordered(:)
223  real(dp), allocatable :: kptfull(:,:)
224 
225 ! *************************************************************************
226 
227  if (kptrlatt(1,2) /= 0 .or. kptrlatt(1,3) /= 0 .or. kptrlatt(2,1) /= 0 .or. &
228 &    kptrlatt(2,3) /= 0 .or. kptrlatt(3,1) /= 0 .or. kptrlatt(3,2) /= 0 ) then
229    write (message,'(4a)')&
230 &   'kptrlatt should be diagonal in order to calculate the nesting factor,',ch10,&
231 &   'skipping the nesting factor calculation ',ch10
232    MSG_WARNING(message)
233    return
234  end if
235 
236  if (prtnest /= 1 .and. prtnest /= 2) then
237    MSG_BUG('prtnest should be 1 or 2')
238  end if
239 
240  !write(message,'(a,9(i0,1x))')' mknesting : kptrlatt = ',kptrlatt
241  !call wrtout(std_out,message,'COLL')
242 
243  nkptfull = kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3)
244  ABI_MALLOC(nestordered,(nkptfull))
245  nestordered(:)=zero
246  ABI_MALLOC(kptfull,(3,nkptfull))
247 
248  ikpt = 0
249  do ik3 = 0, kptrlatt(3,3)-1
250    do ik2 = 0, kptrlatt(2,2)-1
251      do ik1 = 0, kptrlatt(1,1)-1
252        ikpt = ikpt+1
253        kptfull(:,ikpt) = (/dble(ik1)/dble(kptrlatt(1,1)), dble(ik2)/dble(kptrlatt(2,2)),dble(ik3)/dble(kptrlatt(3,3))/)
254      end do
255    end do
256  end do
257 
258 !NOTE: input weights are not normalised, the normalisation factor in introduced in bfactor
259 !new version now puts kptfull in correct order before bfactor, so no need to re-order...
260  if (present(symrec)) then
261    ABI_CHECK(present(nsym), "error - provide nsym and symrec arguments together")
262    call mkkptrank (kpt,nkpt,kptrank_t, nsym, symrec)
263  else
264    call mkkptrank (kpt,nkpt,kptrank_t)
265  end if
266 
267  call bfactor(nkptfull,kptfull,nkptfull,kptfull,kptrank_t,nkpt,weight,nband,nestordered)
268 
269 !================================================================================================
270 !use linear interpolation to plot the bfactor along the given q-path
271 !1) order the kpoints of the grid putting them in increasing x, then y, then z (FORTRAN convention)
272 !2) make table from input kpts to ordered kpts
273 !3) perform interpolation
274 !================================================================================================
275 
276  call outnesting(base_name,gmet,gprimd,kptrlatt,nestordered,nkptfull,nqpath,prtnest,qpath_vertices)
277  ABI_FREE(nestordered)
278 !
279 !now do the same, but for the nesting factor over the phonon qpoints only
280 !
281  ABI_MALLOC(nestfactor,(nqptfull))
282  call bfactor(nkptfull,kptfull,nqptfull,qptfull,kptrank_t,nkpt,weight,nband,nestfactor)
283 
284  call destroy_kptrank (kptrank_t)
285  ABI_FREE(kptfull)
286 
287  call mkkptrank (qptfull,nqptfull,kptrank_t)
288 
289  ABI_MALLOC(ktable,(nqptfull))
290  do ikpt=1,nqptfull
291    ktable(ikpt) = ikpt
292  end do
293 
294  ABI_MALLOC(tmprank, (nqptfull))
295  tmprank = kptrank_t%rank
296  call sort_int(nqptfull, tmprank, ktable)
297  ABI_FREE(tmprank)
298  call destroy_kptrank (kptrank_t)
299 
300 !fill the datagrid for the nesting factor using the Fortran convention and the conventional unit cell
301 !NOTE: the Fortran convention is a must if we want to plot the data
302 !in the BXSF format, useful for the linear interpolation since we use interpol3d.F90
303  ABI_MALLOC(nestordered,(nqptfull))
304  nestordered(:)=zero
305  do jkpt=1,nqptfull
306    ikpt = ktable(jkpt)
307    nestordered(ikpt)=nestfactor(jkpt)
308  end do
309  ABI_FREE(nestfactor)
310  ABI_FREE(ktable)
311 
312  tmpname = trim(base_name)//"kplusq"
313  call outnesting(tmpname,gmet,gprimd,qptrlatt,nestordered,nqptfull,nqpath,prtnest,qpath_vertices)
314 
315  ABI_FREE(nestordered)
316 
317 end subroutine mknesting

m_nesting/outnesting [ Functions ]

[ Top ] [ m_nesting ] [ Functions ]

NAME

 outnesting

FUNCTION

  Write ou the nesting factors calculated in mknesting
  Data on file in the X-Y format (prtnest 1) or
  in the XCrysden format (XSF)   (prtnest 2)

INPUTS

  base_name = prefix of the output file
  gmet = metric in reciprocal space
  gprimd(3,3) dimensional reciprocal lattice vectors
  kptrlatt(3,3) basis vectors for k-grid
  nestordered = nesting function on full grid, points ordered in x, then y, then z
  nkpt = number of k points
  nqpath = number of points requested along the trajectory
  prtnest = flags governing the format of the output file
  qpath_vertices = vertices of the reciprocal space trajectory

OUTPUT

  only write to file

PARENTS

      m_nesting

CHILDREN

      make_path,printxsf,wrap2_zero_one

SOURCE

353 subroutine outnesting(base_name,gmet,gprimd,kptrlatt,nestordered,nkpt,nqpath,prtnest,qpath_vertices)
354 
355 
356 !This section has been created automatically by the script Abilint (TD).
357 !Do not modify the following lines by hand.
358 #undef ABI_FUNC
359 #define ABI_FUNC 'outnesting'
360 !End of the abilint section
361 
362  implicit none
363 
364 !Arguments ------------------------------------
365  integer,intent(in) :: nqpath,prtnest,nkpt
366  character(len=*),intent(in) :: base_name
367 !arrays
368  integer,intent(in) :: kptrlatt(3,3)
369  real(dp),intent(in) :: gprimd(3,3)
370  real(dp),intent(in) :: gmet(3,3)
371  real(dp),intent(in) :: qpath_vertices(3,nqpath)
372  real(dp),intent(in) :: nestordered(nkpt)
373 
374 !Local variables-------------------------------
375 !scalars
376  integer :: unit_nest,nkx,nky,nkz
377  integer :: indx,ii,ipoint,npt_tot,realrecip
378  character(len=fnlen) :: fname
379  character(len=500) :: message
380  real(dp) :: res, kval
381 !arrays
382  integer :: ndiv(nqpath-1)
383  real(dp),allocatable :: finepath(:,:)
384  real(dp) :: tmpkpt(3)
385  real(dp) :: origin(3),qpt(3)
386 ! dummy variables for call to printxsf
387  integer :: natom, ntypat, typat(1)
388  real(dp) :: xcart (3,1), znucl(1)
389 
390 ! *************************************************************************
391 
392 !===================================================================
393 !Definition of the q path along which ph linwid will be interpolated
394 !===================================================================
395  call make_path(nqpath,qpath_vertices,gmet,'G',20,ndiv,npt_tot,finepath)
396 
397  nkx=kptrlatt(1,1)
398  nky=kptrlatt(2,2)
399  nkz=kptrlatt(3,3)
400 
401  if (nkpt /= nkx*nky*nkz) then
402    write(message,'(a,9(i0,1x),2x,i0)')' Wrong input value for kptrlatt  ',kptrlatt, nkpt
403    MSG_BUG(message)
404  end if
405 
406  ! open output file and write header
407  if (open_file(base_name,message,newunit=unit_nest,status="unknown",form="formatted",action="write") /= 0) then
408     MSG_ERROR(message)
409  end if
410 
411  write (unit_nest,'(a)')'#'
412  write (unit_nest,'(a)')'# ABINIT package : Nesting factor file'
413  write (unit_nest,'(a)')'#'
414  write (unit_nest,'(a,i10,a)')'# Nesting factor calculated on ',npt_tot,' Q-points'
415  write (unit_nest,'(a)')'# Description of the Q-path :'
416  write (unit_nest,'(a,i10)')'# Number of line segments = ',nqpath-1
417  write (unit_nest,'(a)')'# Vertices of the Q-path and corresponding index = '
418  indx=1
419  do ii=1,nqpath
420    write (unit_nest,'(a,3(E16.6,1x),i8)')'#  ',qpath_vertices(:,ii),indx
421    if(ii<nqpath) indx=indx+ndiv(ii)
422  end do
423  write (unit_nest,'(a)')'#'
424 
425 !Get qpoint along the q-path from finepath and interpolate the nesting factor
426  indx=1
427 
428  do ipoint=1, npt_tot
429    qpt(:) = finepath(:,ipoint)
430    call wrap2_zero_one(qpt(1),tmpkpt(1),res)
431    call wrap2_zero_one(qpt(2),tmpkpt(2),res)
432    call wrap2_zero_one(qpt(3),tmpkpt(3),res)
433 
434    kval = interpol3d(tmpkpt,nkx,nky,nkz,nestordered)
435 
436    write(unit_nest,'(i5,18e16.5)')indx,kval
437    indx = indx+1
438  end do
439 
440  close (unit_nest)
441  ABI_FREE(finepath)
442 
443  if (prtnest==2) then !write also the nest factor in the XSF format
444    fname=trim(base_name) // '_NEST_XSF'
445 
446    if (open_file(fname,message,newunit=unit_nest,status="unknown",form="formatted",action="write") /= 0) then
447       MSG_ERROR(message)
448    end if
449 
450    origin(:)=zero
451    realrecip=1 !reciprocal space
452    natom = 1
453    ntypat = 1
454    typat = (/1/)
455    xcart = reshape ((/zero, zero, zero/), (/3,1/))
456    znucl = (/one/)
457    call printxsf(nkx,nky,nkz,nestordered,gprimd,origin,natom, ntypat, typat, xcart, znucl, unit_nest,realrecip)
458 
459    close (unit_nest)
460  end if
461 
462 end subroutine outnesting