TABLE OF CONTENTS


ABINIT/m_pptools [ Modules ]

[ Top ] [ Modules ]

NAME

 m_pptools

FUNCTION

  Helper functions used for post-processing.

COPYRIGHT

 Copyright (C) 2002-2022 ABINIT group (MG, ZL, MJV, BXu)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .
 For the initials of contributors, see ~abinit/doc/developers/contributors.txt .

SOURCE

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 MODULE m_pptools
24 
25  use defs_basis
26  use m_errors
27  use m_abicore
28  use m_krank
29 
30  use m_io_tools,        only : open_file
31  use m_fstrings,        only : sjoin, itoa
32  use m_numeric_tools,   only : wrap2_pmhalf
33 
34  implicit none
35 
36  private
37 
38  public :: prmat                ! print real(dp) matrices in an attractive format.
39  public :: printxsf             ! Write a generic array in the XSF format (XCrysden format)
40  public :: print_fofr_ri        ! Print the [real, imaginary] part of an array
41  public :: print_fofr_xyzri     ! Print the Cartesian coordinates and the [real,imaginary] part of an array
42  public :: print_fofr_cube      ! Print ||fofr|| in CUBE format.
43  public :: printbxsf            ! Print band structure energies in XCrysDen format.
44  public :: printvtk             ! Print band structure energies and velocities in VTK format.
45 
46 CONTAINS  !===========================================================

m_pptools/print_fofr_cube [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

  print_fofr_cube

FUNCTION

  Print array fofr in the cube file format

INPUTS

  nx,ny,nz,ldx,ldy,ldz = Logical and physical dimensions of the array.
  fofr(2,ldx,ldy,ldz) = Input data
  rprimd(3,3)=Lattive vectors in Bohr
  [unit] = Fortran unit number. Default: std_out

OUTPUT

  Only writing

SOURCE

433 subroutine print_fofr_cube(nx,ny,nz,ldx,ldy,ldz,fofr,rprimd,natom,znucl_atom,xcart,unit)
434 
435 
436 !Arguments -----------------------------------------------
437 !scalars
438  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,natom
439  integer,optional,intent(in) :: unit
440 !arrays
441  integer,intent(in) :: znucl_atom(natom)
442  real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
443  real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
444 
445 !Local variables-------------------------------
446 !scalars
447  integer,parameter :: cplx=2
448  integer :: ount,ix,iy,iz,iatom
449 !arrays
450 
451 ! *************************************************************************
452 
453  ount = std_out; if (PRESENT(unit)) ount = unit
454 
455 ! EXAMPLE FROM THE WEB
456 ! CPMD CUBE FILE.
457 ! OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
458 ! 3    0.000000    0.000000    0.000000
459 ! 40    0.283459    0.000000    0.000000
460 ! 40    0.000000    0.283459    0.000000
461 ! 40    0.000000    0.000000    0.283459
462 ! 8    0.000000    5.570575    5.669178    5.593517
463 ! 1    0.000000    5.562867    5.669178    7.428055
464 ! 1    0.000000    7.340606    5.669178    5.111259
465 ! -0.25568E-04  0.59213E-05  0.81068E-05  0.10868E-04  0.11313E-04  0.35999E-05
466 
467  write(ount,'(a)') 'ABINIT generated cube file'
468  write(ount,'(a)') 'from cut3d tool'
469 
470  write(ount,'(i9,3(1x,f12.6))') natom,0.,0.,0.
471  write(ount,'(i9,3(1x,f12.6))') nx,(rprimd(iy,1)/nx, iy=1,3)
472  write(ount,'(i9,3(1x,f12.6))') ny,(rprimd(iy,2)/ny, iy=1,3)
473  write(ount,'(i9,3(1x,f12.6))') nz,(rprimd(iy,3)/nz, iy=1,3)
474 
475  do iatom=1,natom
476    write(ount,'(i9,4(3X,ES17.10))') znucl_atom(iatom),0.d0,xcart(1:3,iatom)
477  end do
478 
479 ! Note C ordering of the indexes
480  if (cplx==2) then
481    do ix=1,nx
482      do iy=1,ny
483        do iz=1,nz
484          write(ount,'(6(f12.6,2x))') sqrt(fofr(1,ix,iy,iz)**2 + fofr(2,ix,iy,iz)**2 )
485        end do
486      end do
487    end do
488  else
489    do ix=1,nx
490      do iy=1,ny
491        do iz=1,nz
492          write(ount,'(6(f12.6,2x))') fofr(1,ix,iy,iz)
493        end do
494      end do
495    end do
496  end if
497 
498 end subroutine print_fofr_cube

m_pptools/print_fofr_ri [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

  print_fofr_ri

FUNCTION

  Print the [real,imaginary] part of fofr on unit unit

INPUTS

  ri_mode =
    "RI" if both real and imag part are wanted
    "R"  for real part
    "I"  for imaginary part
  nx,ny,nz,ldx,ldy,ldz = Logical and physical dimensions of the array.
  fofr(2,ldx,ldy,ldz) = Input data
  [unit] = Fortran unit number. Default: std_out

OUTPUT

  Only writing

SOURCE

264 subroutine print_fofr_ri(ri_mode,nx,ny,nz,ldx,ldy,ldz,fofr,unit)
265 
266 
267 !Arguments -----------------------------------------------
268 !scalars
269  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz
270  integer,optional,intent(in) :: unit
271  character(len=*),intent(in) :: ri_mode
272 !arrays
273  real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
274 
275 !Local variables-------------------------------
276 !scalars
277  integer :: ount,ix,iy,iz
278 !arrays
279 
280 ! *************************************************************************
281 
282  ount = std_out; if (PRESENT(unit)) ount = unit
283 
284  SELECT CASE (ri_mode)
285  CASE ("RI","ri")
286    do iz=1,nz
287      do iy=1,ny
288        do ix=1,nx
289          write(ount,'(2f20.16)') fofr(:,ix,iy,iz)
290        end do
291      end do
292    end do
293 
294  CASE ("R","r")
295    do iz=1,nz
296      do iy=1,ny
297        do ix=1,nx
298          write(ount,'(f20.16)') fofr(1,ix,iy,iz)
299        end do
300      end do
301    end do
302 
303  CASE ("I","i")
304    do iz=1,nz
305      do iy=1,ny
306        do ix=1,nx
307          write(ount,'(f20.16)') fofr(2,ix,iy,iz)
308        end do
309      end do
310    end do
311 
312  CASE DEFAULT
313    ABI_ERROR("Wrong ri_mode")
314  END SELECT
315 
316 end subroutine print_fofr_ri

m_pptools/print_fofr_xyzri [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

  print_fofr_xyzri

FUNCTION

  Print the Cartesian coordinates and the [real,imaginary] part of fofr on unit unit

INPUTS

  ri_mode =
    "RI" if both real and imag part are wanted
    "R"  for real part
    "I"  for imaginary part
  nx,ny,nz,ldx,ldy,ldz = Logical and physical dimensions of the array.
  fofr(2,ldx,ldy,ldz) = Input data
  rprimd(3,3)=Lattive vectors in Bohr
  [conv_fact] = Conversion factor for rprimd (rprimd is multiplied by conv_fact). Default is one
  [unit] = Fortran unit number. Default: std_out

OUTPUT

  Only writing

SOURCE

345 subroutine print_fofr_xyzri(ri_mode,nx,ny,nz,ldx,ldy,ldz,fofr,rprimd,conv_fact,unit)
346 
347 
348 !Arguments -----------------------------------------------
349 !scalars
350  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz
351  integer,optional,intent(in) :: unit
352  real(dp),optional,intent(in) :: conv_fact
353  character(len=*),intent(in) :: ri_mode
354 !arrays
355  real(dp),intent(in) :: rprimd(3,3)
356  real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
357 
358 !Local variables-------------------------------
359 !scalars
360  integer :: ount,ix,iy,iz
361  real(dp) :: xnow,ynow,znow,my_cfact
362 
363 ! *************************************************************************
364 
365  ount = std_out; if (PRESENT(unit)) ount = unit
366  my_cfact = one; if (PRESENT(conv_fact)) my_cfact = conv_fact
367 
368  SELECT CASE (ri_mode)
369  CASE ("RI","ri")
370    do iz=1,nz
371      do iy=1,ny
372        do ix=1,nx
373          xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
374          ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
375          znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
376          write(ount,'(3f16.10,2f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(:,ix,iy,iz)
377        end do
378      end do
379    end do
380 
381  CASE ("R","r")
382    do iz=1,nz
383      do iy=1,ny
384        do ix=1,nx
385          xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
386          ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
387          znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
388          write(ount,'(3f16.10,f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(1,ix,iy,iz)
389        end do
390      end do
391    end do
392 
393  CASE ("I","i")
394    do iz=1,nz
395      do iy=1,ny
396        do ix=1,nx
397          xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
398          ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
399          znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
400          write(ount,'(3f16.10,f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(2,ix,iy,iz)
401        end do
402      end do
403    end do
404 
405  CASE DEFAULT
406    ABI_ERROR("Wrong ri_mode")
407  END SELECT
408 
409 end subroutine print_fofr_xyzri

m_pptools/printbxsf [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

 printbxsf

FUNCTION

  Print band structure energies in XCrysDen format.

INPUTS

  eigen(mband,nkpt,nsppol) = eigenvalues in hartree
  ewind = energy window around the fermi level.
          if ewind /= 0 ==> a band is considered in the plot of FSurf
                            only if it is inside [ ef-ewind, ef+ewind ] for some k point
          if ewind == 0 ==> all bands will be keept in the _BXSF file
  fermie = Fermi energy (Hartree)
  gprimd(3,3) = dimensional primitive translations for reciprocal space (bohr^-1)
  kptrlatt(3,3) = reciprocal of lattice vectors for full kpoint grid
  mband = maximum number of bands
  nsppol = 1 for unpolarized, 2 for spin-polarized
  shiftk(3,nshiftk) =shift vector for k point grid
  fname = filename for the fortran file
  symafm(nsym)=(Anti)ferromagnetic symmetries.
  use_afm=.TRUE. if (anti)ferromagnetic symmetries are used.

OUTPUT

  ierr=Status error.
  BXSF file.

SOURCE

532 subroutine printbxsf(eigen,ewind,fermie,gprimd,kptrlatt,mband,&
533 & nkptirred,kptirred,nsym,use_afm,symrec,symafm,use_tr,nsppol,shiftk,nshiftk,fname,ierr)
534 
535 
536 !Arguments ------------------------------------
537 !scalars
538  integer,intent(in) :: mband,nkptirred,nshiftk,nsppol,nsym
539  integer,intent(out) :: ierr
540  real(dp),intent(in) :: ewind,fermie
541  logical,intent(in) :: use_afm,use_tr
542  character(len=*),intent(in) :: fname
543 !arrays
544  integer,intent(in) :: kptrlatt(3,3),symafm(nsym),symrec(3,3,nsym)
545  real(dp),intent(in) :: eigen(mband,nkptirred,nsppol),gprimd(3,3)
546  real(dp),intent(in) :: kptirred(3,nkptirred),shiftk(3,nshiftk)
547 
548 !Local variables-------------------------------
549 !scalars
550  integer,parameter :: enough = 50
551  integer :: iband,ik1,ik2,ik3,ikgrid,ikpt,indx
552  integer :: isppol,isym,maxband,minband,nk1,nk2,nk3,nkptfull,ubxsf,timrev
553  integer :: symkptrank, nsymfm, isymfm
554  real(dp) :: ene
555  character(len=500) :: msg
556  type(krank_t) :: krank
557 !arrays
558  integer,allocatable :: fulltoirred(:),symrecfm(:,:,:)
559  real(dp) :: kptgrid(3),gmet(3,3)
560 
561 ! *************************************************************************
562 
563  ierr = 0
564 
565  ! Error if klatt is no simple orthogonal lattice (in red space)
566  ! for generalization to MP grids, need new version of XCrysDen
567  if (kptrlatt(1,2)/=0 .or. kptrlatt(1,3)/=0 .or. kptrlatt(2,1)/=0 .or. &
568      kptrlatt(2,3)/=0 .or. kptrlatt(3,1)/=0 .or. kptrlatt(3,2)/=0 ) then
569    write(msg,'(3a)')&
570     'kptrlatt should be diagonal, for the FS calculation ',ch10,&
571     'Action: use an orthogonal k-grid for the GS calculation '
572    ABI_COMMENT(msg)
573    ierr = ierr + 1
574  end if
575 
576  ! Error if there are not at least 2 kpts in each direction:
577  ! kptrank will fail for the intermediate points below
578  if (abs(kptrlatt(1,1)) < 2 .or. abs(kptrlatt(2,2)) < 2 .or. abs(kptrlatt(3,3)) < 2) then
579    write(msg,'(3a)')&
580     'You need at least 2 points in each direction in k space to output BXSF files ',ch10,&
581     'Action: use an augmented k-grid for the GS calculation (at least 2x2x2) '
582    ABI_COMMENT(msg)
583    ierr = ierr + 1
584  end if
585 
586  if (ANY(ABS(shiftk) > tol10)) then
587    write(msg,'(3a)')&
588     'Origin of the k-grid should be (0,0,0) for the FS calculation ',ch10,&
589     'Action: use a non-shifted k-grid for the GS calculation. Returning '
590    ABI_COMMENT(msg)
591    ierr = ierr + 1
592  end if
593 
594  if (ierr /= 0) return
595 
596  ! Compute reciprocal space metric.
597  gmet = MATMUL(TRANSPOSE(gprimd), gprimd)
598 
599  if (use_afm) then
600    nsymfm = 0
601    do isym = 1, nsym
602      if (symafm(isym) == 1) nsymfm = nsymfm+1
603    end do
604    ABI_MALLOC(symrecfm, (3,3,nsymfm))
605    isymfm = 0
606    do isym = 1, nsym
607      if (symafm(isym) == 1) then
608        isymfm = isymfm + 1
609        symrecfm(:,:,isymfm) = symrec(:,:,isym)
610      end if
611    end do
612  else
613    nsymfm = nsym
614    ABI_MALLOC(symrecfm, (3,3,nsymfm))
615    symrecfm = symrec
616  end if
617 
618  ! Xcrysden uses aperiodic data-grid (images are included in the grid)
619  nk1 = kptrlatt(1,1); nk2 = kptrlatt(2,2); nk3 = kptrlatt(3,3)
620  nkptfull = (nk1+1) * (nk2+1) * (nk3+1)
621 
622  ABI_MALLOC(fulltoirred, (nkptfull))
623  timrev = 0; if (use_tr) timrev=1
624 
625  krank = krank_new(nkptirred, kptirred, nsym=nsymfm, symrec=symrecfm, time_reversal=use_tr)
626 
627  ! Xcrysden employs the C-ordering for the Fermi Surface (x-y-z)
628  ikgrid=0
629  do ik1=0,nk1
630    do ik2=0,nk2
631      do ik3=0,nk3
632 
633        ikgrid = ikgrid+1
634        kptgrid(1) = DBLE(ik1)/kptrlatt(1,1)
635        kptgrid(2) = DBLE(ik2)/kptrlatt(2,2)
636        kptgrid(3) = DBLE(ik3)/kptrlatt(3,3)
637 
638        ! Find correspondence between the Xcrysden grid and the IBZ
639        symkptrank = krank%get_rank(kptgrid)
640        fulltoirred(ikgrid) = krank%invrank(symkptrank)
641 
642        if (fulltoirred(ikgrid) < 1) then
643          if (ierr <= enough) then
644            write(msg,'(a,3es16.8,2a,i0,2a)')&
645             'kpt = ',kptgrid,ch10,' with rank ', symkptrank, ch10,&
646             'has no symmetric among the k-points used in the GS calculation '
647            ABI_WARNING(msg)
648          end if
649          ierr = ierr + 1
650        end if
651 
652      end do !ik1
653    end do !ik2
654  end do !ik3
655 
656  call krank%free()
657 
658  ABI_CHECK(ierr == 0, "See above warnings")
659 
660  if (abs(ewind) < tol12 ) then
661    ! Keep all bands.
662    minband=1
663    maxband=mband
664  else
665    ! Select a subset of bands.
666    minband = mband
667    maxband = 0
668    ene=abs(ewind)
669    do isppol=1,nsppol
670      do iband=1,mband
671        if(minval(eigen(iband,:,isppol))-fermie < -ene) minband = iband
672      end do
673      do iband=mband,1,-1
674        if (maxval(eigen(iband,:,isppol))-fermie > ene) maxband = iband
675      end do
676    end do ! isppol
677 
678  end if ! abs(energy_window)
679 
680  ! Dump results to file
681  if (open_file(fname,msg, newunit=ubxsf, status='unknown', action="write", form='formatted') /= 0 ) then
682    ABI_WARNING(msg)
683    ierr=ierr +1; RETURN
684  end if
685 
686  ! Write header
687  write(ubxsf,*)' BEGIN_INFO'
688  write(ubxsf,*)'   #'
689  write(ubxsf,*)'   # this is a Band-XCRYSDEN-Structure-File for Visualization of Fermi Surface'
690  write(ubxsf,*)'   # generated by the ABINIT package'
691  write(ubxsf,*)'   #'
692  write(ubxsf,*)'   #  bands between ',minband,' and ',maxband
693  write(ubxsf,*)'   #'
694  if (nsppol == 2 ) then
695    write(ubxsf,*)'   # NOTE: the first band is relative to spin-up electrons,'
696    write(ubxsf,*)'   # the second band to spin-down and so on .. '
697    write(ubxsf,*)'   #'
698  end if
699  write(ubxsf,*)'   # Launch as: xcrysden --bxsf '
700  write(ubxsf,*)'   #'
701  write(ubxsf,'(a,es16.8)')'   Fermi Energy: ',fermie
702  write(ubxsf,*)' END_INFO'
703  write(ubxsf,*)' '
704  write(ubxsf,*)' BEGIN_BLOCK_BANDGRID_3D'
705  write(ubxsf,*)' band_energies'
706  write(ubxsf,*)' BEGIN_BANDGRID_3D'
707 
708  write(ubxsf,*)' ',(maxband-minband+1)*nsppol
709  write(ubxsf,*)' ',nk1+1,nk2+1,nk3+1
710  write(ubxsf,*)' ',shiftk(:,1)
711  ! Angstrom units are used in the BXSF format
712  write(ubxsf,*)' ',gprimd(:,1)/Bohr_Ang
713  write(ubxsf,*)' ',gprimd(:,2)/Bohr_Ang
714  write(ubxsf,*)' ',gprimd(:,3)/Bohr_Ang
715 
716  ! print out data for all relevant bands and full kpt grid (redundant, yes)
717  ! for each kpt in full zone, find equivalent irred kpt and print eigenval
718  indx = 0
719  do iband=minband,maxband
720    do isppol=1,nsppol
721      write(ubxsf,*)' BAND: ',indx+minband
722      write(ubxsf,'(7(es16.8))')(eigen(iband,fulltoirred(ikpt),isppol),ikpt=1,nkptfull)
723      indx=indx+1
724    end do
725  end do
726 
727  write(ubxsf,*)'  END_BANDGRID_3D'
728  write(ubxsf,*)' END_BLOCK_BANDGRID_3D'
729  close(ubxsf)
730 
731  ABI_FREE(fulltoirred)
732  ABI_FREE(symrecfm)
733 
734 end subroutine printbxsf

m_pptools/printvtk [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

 printvtk

FUNCTION

  Print band structure energies and velocities in VTK format.

INPUTS

  eigen(mband,nkpt,nsppol) = eigenvalues in hartree
  ewind = energy window around the fermi level.
          if ewind /= 0 ==> a band is considered in the plot of FSurf
                            only if it is inside [ ef-ewind, ef+ewind ] for some k point
          if ewind == 0 ==> all bands will be keept in the _BXSF file
  fermie = Fermi energy (Hartree)
  gprimd(3,3) = dimensional primitive translations for reciprocal space (bohr^-1)
  kptrlatt(3,3) = reciprocal of lattice vectors for full kpoint grid
  mband = maximum number of bands
  nsppol = 1 for unpolarized, 2 for spin-polarized
  shiftk(3,nshiftk) =shift vector for k point grid
  fname = filename for the fortran file
  symafm(nsym)=(Anti)ferromagnetic symmetries.
  use_afm=.TRUE. if (anti)ferromagnetic symmetries are used.

OUTPUT

  ierr=Status error.
  BXSF file.

SOURCE

 766 subroutine printvtk(eigen,v_surf,ewind,fermie,gprimd,kptrlatt,mband,&
 767 & nkptirred,kptirred,nsym,use_afm,symrec,symafm,use_tr,nsppol,shiftk,nshiftk,fname,ierr)
 768 
 769 
 770 !Arguments ------------------------------------
 771 !scalars
 772  integer,intent(in) :: mband,nkptirred,nshiftk,nsppol,nsym
 773  integer,intent(out) :: ierr
 774  real(dp),intent(in) :: ewind,fermie
 775  logical,intent(in) :: use_afm,use_tr
 776  character(len=*),intent(in) :: fname
 777 !arrays
 778  integer,intent(in) :: kptrlatt(3,3),symafm(nsym),symrec(3,3,nsym)
 779  real(dp),intent(in) :: eigen(mband,nkptirred,nsppol),gprimd(3,3)
 780  real(dp),intent(in) :: v_surf(mband,kptrlatt(1,1)+1,kptrlatt(2,2)+1,kptrlatt(3,3)+1,3,nsppol)
 781  real(dp),intent(in) :: kptirred(3,nkptirred),shiftk(3,nshiftk)
 782 
 783 !Local variables-------------------------------
 784 !scalars
 785  integer :: iband,ikgrid,ikpt1,indx
 786  integer :: ikpt,jkpt,kkpt,ikpt_fine, ik1, ik2, ik3
 787  integer :: isppol,isym,itim,maxband,minband,nk1,nk2,nk3,nkptfull,uvtk,timrev
 788  real(dp) :: ene,res,ss,timsign
 789  logical :: found
 790  character(len=500) :: msg, format_str
 791 !arrays
 792  integer,allocatable :: fulltoirred(:)
 793  real(dp) :: kconv(3),kpt(3),kptgrid(3),kptsym(3)
 794 
 795 ! *************************************************************************
 796 
 797  ierr=0
 798 
 799 !Error if klatt is no simple orthogonal lattice (in red space)
 800 !for generalization to MP grids, need new version of XCrysDen
 801 
 802  if (kptrlatt(1,2)/=0 .or. kptrlatt(1,3)/=0 .or. kptrlatt(2,1)/=0 .or. &
 803      kptrlatt(2,3)/=0 .or. kptrlatt(3,1)/=0 .or. kptrlatt(3,2)/=0 ) then
 804    write(msg,'(3a)')&
 805     'kptrlatt should be diagonal, for the FS calculation ',ch10,&
 806     'action: use an orthogonal k-grid for the GS calculation '
 807    ABI_COMMENT(msg)
 808    ierr=ierr+1
 809  end if
 810 
 811  if (ANY(ABS(shiftk(:,:))>tol10)) then
 812    write(msg,'(3a)')&
 813     'Origin of the k-grid should be (0,0,0) for the FS calculation ',ch10,&
 814     'Action: use a non-shifted k-grid for the GS calculation. Returning '
 815    ABI_COMMENT(msg)
 816    ierr=ierr+1
 817  end if
 818 
 819  if (ierr/=0) RETURN
 820 
 821  ! Xcrysden uses aperiodical data-grid
 822  nk1 = kptrlatt(1,1)
 823  nk2 = kptrlatt(2,2)
 824  nk3 = kptrlatt(3,3)
 825  nkptfull=(nk1+1)*(nk2+1)*(nk3+1)
 826 
 827  ABI_MALLOC(fulltoirred,(nkptfull))
 828  timrev=0; if (use_tr) timrev=1
 829 
 830  !Xcrysden employs C-ordering for the Fermi Surface.
 831  ierr = 0
 832  ikgrid=0
 833  do ik1=0,nk1
 834    do ik2=0,nk2
 835      do ik3=0,nk3
 836 
 837        ikgrid=ikgrid+1
 838        kptgrid(1)=DBLE(ik1)/kptrlatt(1,1)
 839        kptgrid(2)=DBLE(ik2)/kptrlatt(2,2)
 840        kptgrid(3)=DBLE(ik3)/kptrlatt(3,3)
 841        call wrap2_pmhalf(kptgrid(1),kpt(1),res)
 842        call wrap2_pmhalf(kptgrid(2),kpt(2),res)
 843        call wrap2_pmhalf(kptgrid(3),kpt(3),res)
 844 
 845        ! === Find correspondence between the Xcrysden grid and the IBZ ===
 846        ! If AFM case, use only Ferromagetic symmetries.
 847        found=.FALSE.
 848        irred: do ikpt1=1,nkptirred
 849          do itim=0,timrev
 850            do isym=1,nsym
 851              if (use_afm.and.symafm(isym)==-1) CYCLE
 852              timsign = one-two*itim
 853              kptsym(:) = timsign*(symrec(:,1,isym)*kptirred(1,ikpt1) + &
 854                                   symrec(:,2,isym)*kptirred(2,ikpt1) + &
 855                                   symrec(:,3,isym)*kptirred(3,ikpt1))
 856              call wrap2_pmhalf(kptsym(1),kconv(1),res)
 857              call wrap2_pmhalf(kptsym(2),kconv(2),res)
 858              call wrap2_pmhalf(kptsym(3),kconv(3),res)
 859              ! is kconv equivalent to kpt?
 860              ss= (kpt(1)-kconv(1))**2 + (kpt(2)-kconv(2))**2 + (kpt(3)-kconv(3))**2
 861              if (ss < tol6) then
 862                found=.TRUE.
 863                fulltoirred(ikgrid)=ikpt1
 864                exit irred
 865              end if
 866 
 867            end do !itim
 868          end do !isym
 869        end do irred
 870 
 871        if (.not.found) then
 872          write(msg,'(a,3es16.8,2a)')&
 873           ' kpt = ',kpt,ch10,' has no symmetric among the irred k-points used in the GS calculation '
 874          ierr=ierr+1
 875          ABI_ERROR(msg)
 876        end if
 877 
 878      end do !ik1
 879    end do !ik2
 880  end do !ik3
 881 
 882 
 883  if (ierr/=0) then
 884    ABI_FREE(fulltoirred)
 885    RETURN
 886  end if
 887 
 888  if (abs(ewind) < tol12 ) then
 889    ! Keep all bands.
 890    minband=1
 891    maxband=mband
 892  else
 893    ! Select a subset of bands.
 894    minband = mband
 895    maxband = 0
 896    ene=abs(ewind)
 897    do isppol=1,nsppol
 898      do iband=1,mband
 899        if(minval(eigen(iband,:,isppol))-fermie < -ene) then
 900          minband = iband
 901        end if
 902      end do
 903      do iband=mband,1,-1
 904        if (maxval(eigen(iband,:,isppol))-fermie > ene) then
 905          maxband = iband
 906        end if
 907      end do
 908    end do ! isppol
 909 
 910  end if ! abs(energy_window)
 911 
 912  ! Dump the results on file ===
 913  if (open_file(fname,msg,newunit=uvtk,status='unknown',form='formatted') /= 0) then
 914    ABI_FREE(fulltoirred)
 915    ABI_WARNING(msg)
 916    ierr=ierr +1; RETURN
 917  end if
 918 
 919  ! write header
 920  write(uvtk,"(a)") '# vtk DataFile Version 2.0'
 921  write(uvtk,"(a)") 'Eigen values for the Fermi surface'
 922  write(uvtk,"(a)") 'ASCII'
 923  write(uvtk,*) ''
 924  write(uvtk,"(a)") 'DATASET STRUCTURED_GRID'
 925  write(uvtk,"(a,3i6)") 'DIMENSIONS', nk1+1,nk2+1,nk3+1
 926  write(uvtk,"(a,i6,a)") 'POINTS',nkptfull,' float'
 927 
 928  do ik3 = 0, nk3
 929    do ik2 = 0, nk2
 930      do ik1 = 0, nk1
 931        write(uvtk,'(3es16.8)') dble(ik1)/nk1*gprimd(1,1)+ &
 932                                dble(ik2)/nk2*gprimd(1,2)+ &
 933                                dble(ik3)/nk3*gprimd(1,3), &
 934                                dble(ik1)/nk1*gprimd(2,1)+ &
 935                                dble(ik2)/nk2*gprimd(2,2)+ &
 936                                dble(ik3)/nk3*gprimd(2,3), &
 937                                dble(ik1)/nk1*gprimd(3,1)+ &
 938                                dble(ik2)/nk2*gprimd(3,2)+ &
 939                                dble(ik3)/nk3*gprimd(3,3)
 940      end do
 941    end do
 942  end do
 943 
 944 !print out data for all relevant bands and full kpt grid (redundant, yes)
 945 !for each kpt in full zone, find equivalent irred kpt and print eigenval
 946  write(uvtk,*) ''
 947  write(uvtk,"(a,i6)") 'POINT_DATA',nkptfull
 948  indx=0
 949  do iband=minband,maxband
 950    do isppol=1,nsppol
 951      if (minband+indx < 10) then
 952        format_str="(a14,i1,1X,a)"
 953      else
 954        format_str="(a14,i2,1X,a)"
 955      end if
 956      write(uvtk,format_str) 'SCALARS eigval', minband+indx, 'float 1'
 957      write(uvtk,"(a)") 'LOOKUP_TABLE default'
 958      write(uvtk,*) ' '
 959      do kkpt = nk3/2+1, nk3+nk3/2+1
 960        do jkpt = nk2/2+1, nk2+nk2/2+1
 961          do ikpt = nk1/2+1, nk1+nk1/2+1
 962            ik1 = ikpt
 963            ik2 = jkpt
 964            ik3 = kkpt
 965            if (ikpt > nk1+1) ik1 = ikpt - nk1
 966            if (jkpt > nk2+1) ik2 = jkpt - nk2
 967            if (kkpt > nk3+1) ik3 = kkpt - nk3
 968 !          get the index with zyx order
 969            ikpt_fine = (ik1-1)*(nk2+1)*(nk3+1) + (ik2-1)*(nk3+1) + ik3
 970            write(uvtk,'(es16.8)') eigen(iband,fulltoirred(ikpt_fine),isppol)
 971          end do
 972        end do
 973      end do
 974      indx=indx+1
 975    end do
 976  end do
 977 
 978  write(uvtk,*) ''
 979  indx=0
 980  do iband=minband,maxband
 981    do isppol=1,nsppol
 982      if (minband+indx < 10) then
 983        format_str="(a10,i1,1X,a)"
 984      else
 985        format_str="(a10,i2,1X,a)"
 986      end if
 987      write(uvtk,format_str) 'SCALARS ve', minband+indx, 'float'
 988      write(uvtk,"(a)") 'LOOKUP_TABLE default'
 989      write(uvtk,*) ' '
 990      do kkpt = nk3/2+1, nk3+nk3/2+1
 991        do jkpt = nk2/2+1, nk2+nk2/2+1
 992          do ikpt = nk1/2+1, nk1+nk1/2+1
 993            ik1 = ikpt
 994            ik2 = jkpt
 995            ik3 = kkpt
 996            if (ikpt > nk1+1) ik1 = ikpt - nk1
 997            if (jkpt > nk2+1) ik2 = jkpt - nk2
 998            if (kkpt > nk3+1) ik3 = kkpt - nk3
 999 !          write(uvtk,'(3i6,3es16.8)') ik1,ik2,ik3,v_surf(iband,ik1,ik2,ik3,1,isppol), &
1000 !          &                                                 v_surf(iband,ik1,ik2,ik3,2,isppol), &
1001 !          &                                                 v_surf(iband,ik1,ik2,ik3,3,isppol)
1002            write(uvtk,'(es16.8)') sqrt(v_surf(iband,ik1,ik2,ik3,1,isppol)*v_surf(iband,ik1,ik2,ik3,1,isppol)+ &
1003 &           v_surf(iband,ik1,ik2,ik3,2,isppol)*v_surf(iband,ik1,ik2,ik3,2,isppol)+ &
1004 &           v_surf(iband,ik1,ik2,ik3,3,isppol)*v_surf(iband,ik1,ik2,ik3,3,isppol))
1005          end do
1006        end do
1007      end do
1008      write(uvtk,format_str) 'SCALARS vx', minband+indx, 'float'
1009      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1010      write(uvtk,*) ' '
1011      do kkpt = nk3/2+1, nk3+nk3/2+1
1012        do jkpt = nk2/2+1, nk2+nk2/2+1
1013          do ikpt = nk1/2+1, nk1+nk1/2+1
1014            ik1 = ikpt
1015            ik2 = jkpt
1016            ik3 = kkpt
1017            if (ikpt > nk1+1) ik1 = ikpt - nk1
1018            if (jkpt > nk2+1) ik2 = jkpt - nk2
1019            if (kkpt > nk3+1) ik3 = kkpt - nk3
1020            write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,1,isppol)
1021          end do
1022        end do
1023      end do
1024      write(uvtk,format_str) 'SCALARS vy', minband+indx, 'float'
1025      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1026      write(uvtk,*) ' '
1027      do kkpt = nk3/2+1, nk3+nk3/2+1
1028        do jkpt = nk2/2+1, nk2+nk2/2+1
1029          do ikpt = nk1/2+1, nk1+nk1/2+1
1030            ik1 = ikpt
1031            ik2 = jkpt
1032            ik3 = kkpt
1033            if (ikpt > nk1+1) ik1 = ikpt - nk1
1034            if (jkpt > nk2+1) ik2 = jkpt - nk2
1035            if (kkpt > nk3+1) ik3 = kkpt - nk3
1036            write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,2,isppol)
1037          end do
1038        end do
1039      end do
1040      write(uvtk,format_str) 'SCALARS vz', minband+indx, 'float'
1041      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1042      write(uvtk,*) ' '
1043      do kkpt = nk3/2+1, nk3+nk3/2+1
1044        do jkpt = nk2/2+1, nk2+nk2/2+1
1045          do ikpt = nk1/2+1, nk1+nk1/2+1
1046            ik1 = ikpt
1047            ik2 = jkpt
1048            ik3 = kkpt
1049            if (ikpt > nk1+1) ik1 = ikpt - nk1
1050            if (jkpt > nk2+1) ik2 = jkpt - nk2
1051            if (kkpt > nk3+1) ik3 = kkpt - nk3
1052            write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,3,isppol)
1053          end do
1054        end do
1055      end do
1056      indx=indx+1
1057    end do
1058  end do
1059 
1060  close (uvtk)
1061  ABI_FREE(fulltoirred)
1062 
1063 end subroutine printvtk

m_pptools/printxsf [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

 printxsf

FUNCTION

 Write a generic array in the XSF format (XCrysden format)

INPUTS

 basis(3,3) = basis vectors of the direct real lattice or of the reciprocal lattice (fortran convention)
              (Bohr units if realrecip=0, Bohr^-1 if realrecip=1, see below)
 realrecip = 0  for a plot in real space
             1  for a plot in reciprocal space
 nunit   = unit number of the output file (already open by the caller, not closed here!)
 n1=grid size along x
 n2=grid size along y
 n3=grid size along z
 origin(3) = origin of the grid
 datagrid(n1*n2*n3) = datagrid values stored using the fortran convention

OUTPUT

 Only write

SOURCE

140 subroutine printxsf(n1,n2,n3,datagrid,basis,origin,natom,ntypat,typat,xcart,znucl,nunit,realrecip)
141 
142 
143 !Arguments ------------------------------------
144 !scalars
145  integer,intent(in) :: n1,n2,n3,nunit,realrecip
146  integer,intent(in) :: natom,ntypat
147 !arrays
148  integer,intent(in) :: typat(natom)
149  real(dp),intent(in) :: basis(3,3),datagrid(n1*n2*n3),origin(3)
150  real(dp),intent(in) :: xcart(3,natom), znucl(ntypat)
151 
152 !Local variables-------------------------------
153 !scalars
154  integer :: ix,iy,iz,nslice,nsym,iatom
155  real(dp) :: fact
156 !arrays
157  real(dp) :: tau(3,natom)
158 
159 ! *************************************************************************
160 
161  DBG_ENTER("COLL")
162 
163  if (all(realrecip/= [0,1])) then
164    ABI_BUG(sjoin('The argument realrecip should be 0 or 1, received:', itoa(realrecip)))
165  end if
166 
167 !conversion between ABINIT default units and XCrysden units
168  fact=Bohr_Ang; if (realrecip ==1) fact=one/fact  !since we are in reciprocal space
169 
170 !TODO insert crystalline structure and dummy atoms in case of reciprocal space
171 !need to convert basis too
172 
173  write(nunit,'(1X,A)')  'DIM-GROUP'
174  write(nunit,*) '3  1'
175  write(nunit,'(1X,A)') 'PRIMVEC'
176  do iy = 1,3
177    write(nunit,'(3(ES17.10,2X))') (Bohr_Ang*basis(ix,iy), ix=1,3)
178  end do
179 !
180 !generate translated coordinates to fit origin shift
181 !
182  do iatom = 1,natom
183    tau (:,iatom) = xcart(:,iatom) - origin(:)
184  end do
185 
186  write(nunit,'(1X,A)') 'PRIMCOORD'
187  write(nunit,*) natom, ' 1'
188  do iatom = 1,natom
189    write(nunit,'(i9,3(3X,ES17.10))') NINT(znucl(typat(iatom))), &  ! WARNING alchemy not supported by XCrysden
190      Bohr_Ang*tau(1,iatom), &
191      Bohr_Ang*tau(2,iatom), &
192      Bohr_Ang*tau(3,iatom)
193  end do
194  write(nunit,'(1X,A)') 'ATOMS'
195  do iatom = 1,natom
196    write(nunit,'(i9,3(3X,ES17.10))') NINT(znucl(typat(iatom))), & ! WARNING alchemy not supported by XCrysden
197      Bohr_Ang*tau(1,iatom), &
198      Bohr_Ang*tau(2,iatom), &
199      Bohr_Ang*tau(3,iatom)
200  end do
201 
202  write(nunit,'(a)')' BEGIN_BLOCK_DATAGRID3D'
203  write(nunit,'(a)')' datagrid'
204  write(nunit,'(a)')' DATAGRID_3D_DENSITY'
205 !NOTE: XCrysden uses aperiodical data grid
206  write(nunit,*)n1+1,n2+1,n3+1
207  write(nunit,*)origin
208  write(nunit,*)basis(:,1)*fact
209  write(nunit,*)basis(:,2)*fact
210  write(nunit,*)basis(:,3)*fact
211 
212  nslice=1
213  do iz=1,n3
214    do iy=1,n2
215      write(nunit,'(8es16.8)') datagrid(1+n1*(nslice-1):n1+n1*(nslice-1)),datagrid(1+n1*(nslice-1))
216      nslice=nslice+1
217    end do
218    nsym=nslice-n2
219    write (nunit,'(8es16.8)') datagrid(1+n1*(nsym-1):n1+n1*(nsym-1)),datagrid(1+n1*(nsym-1))
220  end do
221 
222 !Now write upper plane
223  nslice=1
224  do iy=1,n2
225    write (nunit,'(8es16.8)') datagrid(1+n1*(nslice-1):n1+n1*(nslice-1)),datagrid(1+n1*(nslice-1))
226    nslice=nslice+1
227  end do
228 
229  nsym=nslice-n2
230  write (nunit,'(8es16.8)') datagrid(1+n1*(nsym-1):n1+n1*(nsym-1)),datagrid(1+n1*(nsym-1))
231 
232  write (nunit,'(a)')' END_DATAGRID_3D'
233  write (nunit,'(a)')' END_BLOCK_DATAGRID3D'
234 
235  DBG_EXIT("COLL")
236 
237 end subroutine printxsf

m_pptools/prmat [ Functions ]

[ Top ] [ m_pptools ] [ Functions ]

NAME

 prmat

FUNCTION

 This subroutine prints real*8 matrices in an attractive format.

INPUTS

  mat(mi,nj)= matrix to be printed
  mi        = no rows of mat
  ni        = no rows to print
  nj        = no colums of mat
  unitm     = unit to print to, if not provided std_out is chosen

OUTPUT

  (only writing)

SOURCE

 68 subroutine prmat (mat, ni, nj, mi, unitm)
 69 
 70 
 71 !Arguments ------------------------------------
 72 !scalars
 73  integer,intent(in)           :: mi,ni,nj
 74  integer,intent(in), optional :: unitm
 75 !arrays
 76  real(dp),intent(in)          :: mat(mi,nj)
 77 
 78 !Local variables-------------------------------
 79 !scalars
 80  character(len=1000)    :: message
 81  integer,parameter      :: nline=10
 82  integer                :: ii,jj,jstart,jstop,unitn
 83 
 84 ! *************************************************************************
 85 
 86  if (present(unitm)) then ! standard printing to std_out
 87    unitn=unitm
 88  else
 89    unitn=std_out
 90  end if
 91 
 92  do  jstart = 1, nj, nline
 93    jstop = min(nj, jstart+nline-1)
 94    write(message, '(3x,10(i4,8x))' ) (jj,jj=jstart,jstop)
 95    call wrtout(unitn,message,'COLL')
 96  end do
 97 
 98  do ii = 1,ni
 99    do jstart= 1, nj, nline
100      jstop = min(nj, jstart+nline-1)
101      if (jstart==1) then
102        write(message, '(i3,1p,10e12.4)' ) ii, (mat(ii,jj),jj=jstart,jstop)
103        call wrtout(unitn,message,'COLL')
104      else
105        write(message, '(3x,1p,10e12.4)' )    (mat(ii,jj),jj=jstart,jstop)
106        call wrtout(unitn,message,'COLL')
107      end if
108    end do
109  end do
110 
111 end subroutine prmat