TABLE OF CONTENTS


ABINIT/m_pptools [ Modules ]

[ Top ] [ Modules ]

NAME

 m_pptools

FUNCTION

  Helper functions used for simple post-processing.

COPYRIGHT

 Copyright (C) 2002-2018 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 .

PARENTS

CHILDREN

SOURCE

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

PARENTS

      m_cut3d

CHILDREN

      wrap2_pmhalf

SOURCE

500 subroutine print_fofr_cube(nx,ny,nz,ldx,ldy,ldz,fofr,rprimd,natom,znucl_atom,xcart,unit)
501 
502 
503 !This section has been created automatically by the script Abilint (TD).
504 !Do not modify the following lines by hand.
505 #undef ABI_FUNC
506 #define ABI_FUNC 'print_fofr_cube'
507 !End of the abilint section
508 
509  implicit none
510 
511 !Arguments -----------------------------------------------
512 !scalars
513  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz,natom
514  integer,optional,intent(in) :: unit
515 !arrays
516  integer,intent(in) :: znucl_atom(natom)
517  real(dp),intent(in) :: rprimd(3,3),xcart(3,natom)
518  real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
519 
520 !Local variables-------------------------------
521 !scalars
522  integer,parameter :: cplx=2
523  integer :: ount,ix,iy,iz,iatom
524 !arrays
525 
526 ! *************************************************************************
527 
528  ount = std_out; if (PRESENT(unit)) ount = unit
529 
530 ! EXAMPLE FROM THE WEB
531 ! CPMD CUBE FILE.
532 ! OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z
533 ! 3    0.000000    0.000000    0.000000
534 ! 40    0.283459    0.000000    0.000000
535 ! 40    0.000000    0.283459    0.000000
536 ! 40    0.000000    0.000000    0.283459
537 ! 8    0.000000    5.570575    5.669178    5.593517
538 ! 1    0.000000    5.562867    5.669178    7.428055
539 ! 1    0.000000    7.340606    5.669178    5.111259
540 ! -0.25568E-04  0.59213E-05  0.81068E-05  0.10868E-04  0.11313E-04  0.35999E-05
541 
542  write(ount,'(a)') 'ABINIT generated cube file'
543  write(ount,'(a)') 'from cut3d tool'
544 
545  write(ount,'(i9,3(1x,f12.6))') natom,0.,0.,0.
546  write(ount,'(i9,3(1x,f12.6))') nx,(rprimd(iy,1)/nx, iy=1,3)
547  write(ount,'(i9,3(1x,f12.6))') ny,(rprimd(iy,2)/ny, iy=1,3)
548  write(ount,'(i9,3(1x,f12.6))') nz,(rprimd(iy,3)/nz, iy=1,3)
549 
550  do iatom=1,natom
551    write(ount,'(i9,4(3X,ES17.10))') znucl_atom(iatom),0.d0,xcart(1:3,iatom)
552  end do
553 
554 ! Note C ordering of the indexes
555  if (cplx==2) then
556    do ix=1,nx
557      do iy=1,ny
558        do iz=1,nz
559          write(ount,'(6(f12.6,2x))') sqrt(fofr(1,ix,iy,iz)**2 + fofr(2,ix,iy,iz)**2 )
560        end do
561      end do
562    end do
563  else
564    do ix=1,nx
565      do iy=1,ny
566        do iz=1,nz
567          write(ount,'(6(f12.6,2x))') fofr(1,ix,iy,iz)
568        end do
569      end do
570    end do
571  end if
572 
573 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

PARENTS

      m_cut3d

CHILDREN

      wrap2_pmhalf

SOURCE

303 subroutine print_fofr_ri(ri_mode,nx,ny,nz,ldx,ldy,ldz,fofr,unit)
304 
305 
306 !This section has been created automatically by the script Abilint (TD).
307 !Do not modify the following lines by hand.
308 #undef ABI_FUNC
309 #define ABI_FUNC 'print_fofr_ri'
310 !End of the abilint section
311 
312  implicit none
313 
314 !Arguments -----------------------------------------------
315 !scalars
316  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz
317  integer,optional,intent(in) :: unit
318  character(len=*),intent(in) :: ri_mode
319 !arrays
320  real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
321 
322 !Local variables-------------------------------
323 !scalars
324  integer :: ount,ix,iy,iz
325 !arrays
326 
327 ! *************************************************************************
328 
329  ount = std_out; if (PRESENT(unit)) ount = unit
330 
331  SELECT CASE (ri_mode)
332  CASE ("RI","ri")
333    do iz=1,nz
334      do iy=1,ny
335        do ix=1,nx
336          write(ount,'(2f20.16)') fofr(:,ix,iy,iz)
337        end do
338      end do
339    end do
340 
341  CASE ("R","r")
342    do iz=1,nz
343      do iy=1,ny
344        do ix=1,nx
345          write(ount,'(f20.16)') fofr(1,ix,iy,iz)
346        end do
347      end do
348    end do
349 
350  CASE ("I","i")
351    do iz=1,nz
352      do iy=1,ny
353        do ix=1,nx
354          write(ount,'(f20.16)') fofr(2,ix,iy,iz)
355        end do
356      end do
357    end do
358 
359  CASE DEFAULT
360    MSG_ERROR("Wrong ri_mode")
361  END SELECT
362 
363 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

PARENTS

      m_cut3d

CHILDREN

      wrap2_pmhalf

SOURCE

398 subroutine print_fofr_xyzri(ri_mode,nx,ny,nz,ldx,ldy,ldz,fofr,rprimd,conv_fact,unit)
399 
400 
401 !This section has been created automatically by the script Abilint (TD).
402 !Do not modify the following lines by hand.
403 #undef ABI_FUNC
404 #define ABI_FUNC 'print_fofr_xyzri'
405 !End of the abilint section
406 
407  implicit none
408 
409 !Arguments -----------------------------------------------
410 !scalars
411  integer,intent(in) :: nx,ny,nz,ldx,ldy,ldz
412  integer,optional,intent(in) :: unit
413  real(dp),optional,intent(in) :: conv_fact
414  character(len=*),intent(in) :: ri_mode
415 !arrays
416  real(dp),intent(in) :: rprimd(3,3)
417  real(dp),intent(in) :: fofr(2,ldx,ldy,ldz)
418 
419 !Local variables-------------------------------
420 !scalars
421  integer :: ount,ix,iy,iz
422  real(dp) :: xnow,ynow,znow,my_cfact
423 
424 ! *************************************************************************
425 
426  ount = std_out; if (PRESENT(unit)) ount = unit
427  my_cfact = one; if (PRESENT(conv_fact)) my_cfact = conv_fact
428 
429  SELECT CASE (ri_mode)
430  CASE ("RI","ri")
431    do iz=1,nz
432      do iy=1,ny
433        do ix=1,nx
434          xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
435          ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
436          znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
437          write(ount,'(3f16.10,2f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(:,ix,iy,iz)
438        end do
439      end do
440    end do
441 
442  CASE ("R","r")
443    do iz=1,nz
444      do iy=1,ny
445        do ix=1,nx
446          xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
447          ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
448          znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
449          write(ount,'(3f16.10,f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(1,ix,iy,iz)
450        end do
451      end do
452    end do
453 
454  CASE ("I","i")
455    do iz=1,nz
456      do iy=1,ny
457        do ix=1,nx
458          xnow = rprimd(1,1)*(ix-1)/nx + rprimd(1,2)*(iy-1)/ny + rprimd(1,3)*(iz-1)/nz
459          ynow = rprimd(2,1)*(ix-1)/nx + rprimd(2,2)*(iy-1)/ny + rprimd(2,3)*(iz-1)/nz
460          znow = rprimd(3,1)*(ix-1)/nx + rprimd(3,2)*(iy-1)/ny + rprimd(3,3)*(iz-1)/nz
461          write(ount,'(3f16.10,f20.16)') my_cfact*xnow, my_cfact*ynow, my_cfact*znow,fofr(2,ix,iy,iz)
462        end do
463      end do
464    end do
465 
466  CASE DEFAULT
467    MSG_ERROR("Wrong ri_mode")
468  END SELECT
469 
470 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.

PARENTS

      m_ebands,m_ifc

CHILDREN

      wrap2_pmhalf

SOURCE

613 subroutine printbxsf(eigen,ewind,fermie,gprimd,kptrlatt,mband,&
614 & nkptirred,kptirred,nsym,use_afm,symrec,symafm,use_tr,nsppol,shiftk,nshiftk,fname,ierr)
615 
616 
617 !This section has been created automatically by the script Abilint (TD).
618 !Do not modify the following lines by hand.
619 #undef ABI_FUNC
620 #define ABI_FUNC 'printbxsf'
621 !End of the abilint section
622 
623  implicit none
624 
625 !Arguments ------------------------------------
626 !scalars
627  integer,intent(in) :: mband,nkptirred,nshiftk,nsppol,nsym
628  integer,intent(out) :: ierr
629  real(dp),intent(in) :: ewind,fermie
630  logical,intent(in) :: use_afm,use_tr
631  character(len=*),intent(in) :: fname
632 !arrays
633  integer,intent(in) :: kptrlatt(3,3),symafm(nsym),symrec(3,3,nsym)
634  real(dp),intent(in) :: eigen(mband,nkptirred,nsppol),gprimd(3,3)
635  real(dp),intent(in) :: kptirred(3,nkptirred),shiftk(3,nshiftk)
636 
637 !Local variables-------------------------------
638 !scalars
639  integer :: iband,ik1,ik2,ik3,ikgrid,ikpt,indx
640  integer :: isppol,isym,maxband,minband,nk1,nk2,nk3,nkptfull,ubxsf,timrev
641  integer :: symkptrank, nsymfm, isymfm
642  real(dp) :: ene
643  character(len=500) :: msg
644  type(kptrank_type) :: kptrank_t
645 !arrays
646  integer,allocatable :: fulltoirred(:),symrecfm(:,:,:)
647  real(dp) :: kptgrid(3),gmet(3,3)
648 
649 ! *************************************************************************
650 
651  ierr = 0
652 
653 ! Error if klatt is no simple orthogonal lattice (in red space)
654 ! for generalization to MP grids, need new version of XCrysDen
655 
656  if ( kptrlatt(1,2)/=0 .or. kptrlatt(1,3)/=0 .or. kptrlatt(2,1)/=0 .or. &
657 & kptrlatt(2,3)/=0 .or. kptrlatt(3,1)/=0 .or. kptrlatt(3,2)/=0 ) then
658    write(msg,'(3a)')&
659 &   'kptrlatt should be diagonal, for the FS calculation ',ch10,&
660 &   'Action: use an orthogonal k-grid for the GS calculation '
661    MSG_COMMENT(msg)
662    ierr = ierr + 1
663  end if
664 
665 ! Error if there are not at least 2 kpts in each direction:
666 ! kptrank will fail for the intermediate points below
667  if ( abs(kptrlatt(1,1))<2 .or. abs(kptrlatt(2,2))<2 .or. abs(kptrlatt(3,3))<2) then
668    write(msg,'(3a)')&
669 &   'You need at least 2 points in each direction in k space to output BXSF files ',ch10,&
670 &   'Action: use an augmented k-grid for the GS calculation (at least 2x2x2) '
671    MSG_COMMENT(msg)
672    ierr = ierr + 1
673  end if
674 
675  if (ANY(ABS(shiftk) > tol10)) then
676    write(msg,'(3a)')&
677 &   'Origin of the k-grid should be (0,0,0) for the FS calculation ',ch10,&
678 &   'Action: use a non-shifted k-grid for the GS calculation. Returning '
679    MSG_COMMENT(msg)
680    ierr = ierr + 1
681  end if
682 
683  if (ierr /= 0) return
684 
685  ! Compute reciprocal space metric.
686  gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
687 
688  if (use_afm) then
689    nsymfm = 0
690    do isym = 1, nsym
691      if (symafm(isym) == 1) nsymfm = nsymfm+1
692    end do
693    ABI_MALLOC(symrecfm,(3,3,nsymfm))
694    isymfm = 0
695    do isym = 1, nsym
696      if (symafm(isym) == 1) then
697        isymfm = isymfm + 1
698        symrecfm(:,:,isymfm) = symrec(:,:,isym)
699      end if
700    end do
701  else
702    nsymfm = nsym
703    ABI_MALLOC(symrecfm,(3,3,nsymfm))
704    symrecfm = symrec
705  end if
706 
707 !Xcrysden uses aperiodical data-grid
708  nk1 = kptrlatt(1,1)
709  nk2 = kptrlatt(2,2)
710  nk3 = kptrlatt(3,3)
711  nkptfull=(nk1+1)*(nk2+1)*(nk3+1)
712 
713  ABI_MALLOC(fulltoirred,(nkptfull))
714  timrev=0; if (use_tr) timrev=1
715 
716  call mkkptrank (kptirred,nkptirred,kptrank_t, nsym=nsymfm, symrec=symrecfm, time_reversal=use_tr)
717 
718 !Xcrysden employs the C-ordering for the Fermi Surface.
719  ikgrid=0
720  do ik1=0,nk1
721    do ik2=0,nk2
722      do ik3=0,nk3
723 
724        ikgrid=ikgrid+1
725        kptgrid(1)=DBLE(ik1)/kptrlatt(1,1)
726        kptgrid(2)=DBLE(ik2)/kptrlatt(2,2)
727        kptgrid(3)=DBLE(ik3)/kptrlatt(3,3)
728 
729        ! Find correspondence between the Xcrysden grid and the IBZ ===
730        call get_rank_1kpt (kptgrid, symkptrank, kptrank_t)
731        fulltoirred(ikgrid) = kptrank_t%invrank(symkptrank)
732 
733        if (fulltoirred(ikgrid) < 1) then
734          write(msg,'(a,3es16.8,2a,i0,2a)')&
735 &         'kpt = ',kptgrid,ch10,' with rank ', symkptrank, ch10,&
736 &         'has no symmetric among the k-points used in the GS calculation '
737          ierr=ierr + 1
738          MSG_WARNING(msg)
739        end if
740 
741      end do !ik1
742    end do !ik2
743  end do !ik3
744 
745  call destroy_kptrank(kptrank_t)
746 
747  if (ierr/=0) then
748    ABI_FREE(fulltoirred)
749    MSG_ERROR("Bug")
750    RETURN
751  end if
752 
753  if (abs(ewind) < tol12 ) then ! Keep all bands.
754    minband=1
755    maxband=mband
756  else ! Select a subset of bands.
757    minband = mband
758    maxband = 0
759    ene=abs(ewind)
760    do isppol=1,nsppol
761      do iband=1,mband
762        if(minval(eigen(iband,:,isppol))-fermie < -ene) then
763          minband = iband
764        end if
765      end do
766      do iband=mband,1,-1
767        if (maxval(eigen(iband,:,isppol))-fermie > ene) then
768          maxband = iband
769        end if
770      end do
771    end do ! isppol
772 
773  end if ! abs(energy_window)
774 
775  ! Dump the results on file
776  if (open_file(fname,msg,newunit=ubxsf,status='unknown',form='formatted') /=0) then
777    MSG_WARNING(msg)
778    ierr=ierr +1; RETURN
779  end if
780 
781 ! Write header
782  write(ubxsf,*)' BEGIN_INFO'
783  write(ubxsf,*)'   #'
784  write(ubxsf,*)'   # this is a Band-XCRYSDEN-Structure-File for Visualization of Fermi Surface'
785  write(ubxsf,*)'   # generated by the ABINIT package'
786  write(ubxsf,*)'   #'
787  write(ubxsf,*)'   #  bands between ',minband,' and ',maxband
788  write(ubxsf,*)'   #'
789  if (nsppol == 2 ) then
790    write(ubxsf,*)'   # NOTE: the first band is relative to spin-up electrons,'
791    write(ubxsf,*)'   # the second band to spin-down and so on .. '
792    write(ubxsf,*)'   #'
793  end if
794  write(ubxsf,*)'   # Launch as: xcrysden --bxsf '
795  write(ubxsf,*)'   #'
796  write(ubxsf,'(a,es16.8)')'   Fermi Energy: ',fermie
797  write(ubxsf,*)' END_INFO'
798  write(ubxsf,*)' '
799  write(ubxsf,*)' BEGIN_BLOCK_BANDGRID_3D'
800  write(ubxsf,*)' band_energies'
801  write(ubxsf,*)' BEGIN_BANDGRID_3D'
802 
803  write(ubxsf,*)' ',(maxband-minband+1)*nsppol
804  write(ubxsf,*)' ',nk1+1,nk2+1,nk3+1
805  write(ubxsf,*)' ',shiftk(:,1)
806 !NOTE : Angstrom units are used in the BXSF format
807  write(ubxsf,*)' ',gprimd(:,1)/Bohr_Ang
808  write(ubxsf,*)' ',gprimd(:,2)/Bohr_Ang
809  write(ubxsf,*)' ',gprimd(:,3)/Bohr_Ang
810 
811 !print out data for all relevant bands and full kpt grid (redundant, yes)
812 !for each kpt in full zone, find equivalent irred kpt and print eigenval
813  indx=0
814  do iband=minband,maxband
815    do isppol=1,nsppol
816      write(ubxsf,*)' BAND: ',indx+minband
817      write(ubxsf,'(7(es16.8))')(eigen(iband,fulltoirred(ikpt),isppol),ikpt=1,nkptfull)
818      indx=indx+1
819    end do
820  end do
821 
822  write(ubxsf,*)'  END_BANDGRID_3D'
823  write(ubxsf,*)' END_BLOCK_BANDGRID_3D'
824  close (ubxsf)
825 
826  ABI_FREE(fulltoirred)
827  ABI_FREE(symrecfm)
828 
829 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.

PARENTS

      elphon

CHILDREN

      wrap2_pmhalf

SOURCE

 867 subroutine printvtk(eigen,v_surf,ewind,fermie,gprimd,kptrlatt,mband,&
 868 & nkptirred,kptirred,nsym,use_afm,symrec,symafm,use_tr,nsppol,shiftk,nshiftk,fname,ierr)
 869 
 870 
 871 !This section has been created automatically by the script Abilint (TD).
 872 !Do not modify the following lines by hand.
 873 #undef ABI_FUNC
 874 #define ABI_FUNC 'printvtk'
 875 !End of the abilint section
 876 
 877  implicit none
 878 
 879 !Arguments ------------------------------------
 880 !scalars
 881  integer,intent(in) :: mband,nkptirred,nshiftk,nsppol,nsym
 882  integer,intent(out) :: ierr
 883  real(dp),intent(in) :: ewind,fermie
 884  logical,intent(in) :: use_afm,use_tr
 885  character(len=*),intent(in) :: fname
 886 !arrays
 887  integer,intent(in) :: kptrlatt(3,3),symafm(nsym),symrec(3,3,nsym)
 888  real(dp),intent(in) :: eigen(mband,nkptirred,nsppol),gprimd(3,3)
 889  real(dp),intent(in) :: v_surf(mband,kptrlatt(1,1)+1,kptrlatt(2,2)+1,kptrlatt(3,3)+1,3,nsppol)
 890  real(dp),intent(in) :: kptirred(3,nkptirred),shiftk(3,nshiftk)
 891 
 892 !Local variables-------------------------------
 893 !scalars
 894  integer :: iband,ikgrid,ikpt1,indx
 895  integer :: ikpt,jkpt,kkpt,ikpt_fine, ik1, ik2, ik3
 896  integer :: isppol,isym,itim,maxband,minband,nk1,nk2,nk3,nkptfull,uvtk,timrev
 897  real(dp) :: ene,res,ss,timsign
 898  logical :: found
 899  character(len=500) :: msg, format_str
 900 !arrays
 901  integer,allocatable :: fulltoirred(:)
 902  real(dp) :: kconv(3),kpt(3),kptgrid(3),kptsym(3)
 903 
 904 ! *************************************************************************
 905 
 906  ierr=0
 907 
 908 !Error if klatt is no simple orthogonal lattice (in red space)
 909 !for generalization to MP grids, need new version of XCrysDen
 910 
 911  if ( kptrlatt(1,2)/=0 .or. kptrlatt(1,3)/=0 .or. kptrlatt(2,1)/=0 .or. &
 912 & kptrlatt(2,3)/=0 .or. kptrlatt(3,1)/=0 .or. kptrlatt(3,2)/=0 ) then
 913    write(msg,'(3a)')&
 914 &   'kptrlatt should be diagonal, for the FS calculation ',ch10,&
 915 &   'action: use an orthogonal k-grid for the GS calculation '
 916    MSG_COMMENT(msg)
 917    ierr=ierr+1
 918  end if
 919 
 920  if (ANY(ABS(shiftk(:,:))>tol10)) then
 921    write(msg,'(3a)')&
 922 &   'Origin of the k-grid should be (0,0,0) for the FS calculation ',ch10,&
 923 &   'Action: use a non-shifted k-grid for the GS calculation. Returning '
 924    MSG_COMMENT(msg)
 925    ierr=ierr+1
 926  end if
 927 
 928  if (ierr/=0) RETURN
 929 
 930 !Xcrysden uses aperiodical data-grid
 931  nk1 = kptrlatt(1,1)
 932  nk2 = kptrlatt(2,2)
 933  nk3 = kptrlatt(3,3)
 934  nkptfull=(nk1+1)*(nk2+1)*(nk3+1)
 935 
 936  ABI_ALLOCATE(fulltoirred,(nkptfull))
 937  timrev=0; if (use_tr) timrev=1
 938 
 939 !Xcrysden employs the C-ordering for the Fermi Surface.
 940  ierr = 0
 941  ikgrid=0
 942  do ik1=0,nk1
 943    do ik2=0,nk2
 944      do ik3=0,nk3
 945 
 946        ikgrid=ikgrid+1
 947        kptgrid(1)=DBLE(ik1)/kptrlatt(1,1)
 948        kptgrid(2)=DBLE(ik2)/kptrlatt(2,2)
 949        kptgrid(3)=DBLE(ik3)/kptrlatt(3,3)
 950        call wrap2_pmhalf(kptgrid(1),kpt(1),res)
 951        call wrap2_pmhalf(kptgrid(2),kpt(2),res)
 952        call wrap2_pmhalf(kptgrid(3),kpt(3),res)
 953 
 954 !      === Find correspondence between the Xcrysden grid and the IBZ ===
 955 !      * If AFM case, use only Ferromagetic symmetries.
 956        found=.FALSE.
 957        irred: do ikpt1=1,nkptirred
 958          do itim=0,timrev
 959            do isym=1,nsym
 960              if (use_afm.and.symafm(isym)==-1) CYCLE
 961              timsign = one-two*itim
 962              kptsym(:) = timsign*(symrec(:,1,isym)*kptirred(1,ikpt1) + &
 963 &             symrec(:,2,isym)*kptirred(2,ikpt1) + &
 964 &             symrec(:,3,isym)*kptirred(3,ikpt1))
 965              call wrap2_pmhalf(kptsym(1),kconv(1),res)
 966              call wrap2_pmhalf(kptsym(2),kconv(2),res)
 967              call wrap2_pmhalf(kptsym(3),kconv(3),res)
 968 !            * is kconv equivalent to kpt?
 969              ss= (kpt(1)-kconv(1))**2 + (kpt(2)-kconv(2))**2 + (kpt(3)-kconv(3))**2
 970              if (ss < tol6) then
 971                found=.TRUE.
 972                fulltoirred(ikgrid)=ikpt1
 973                exit irred
 974              end if
 975 
 976            end do !itim
 977          end do !isym
 978        end do irred
 979 
 980        if (.not.found) then
 981          write(msg,'(a,3es16.8,2a)')&
 982 &         ' kpt = ',kpt,ch10,' has no symmetric among the irred k-points used in the GS calculation '
 983          ierr=ierr+1
 984          MSG_ERROR(msg)
 985        end if
 986 
 987      end do !ik1
 988    end do !ik2
 989  end do !ik3
 990 
 991 
 992  if (ierr/=0) then
 993    ABI_DEALLOCATE(fulltoirred)
 994    RETURN
 995  end if
 996 
 997  if (abs(ewind) < tol12 ) then ! Keep all bands.
 998    minband=1
 999    maxband=mband
1000  else ! Select a subset of bands.
1001    minband = mband
1002    maxband = 0
1003    ene=abs(ewind)
1004    do isppol=1,nsppol
1005      do iband=1,mband
1006        if(minval(eigen(iband,:,isppol))-fermie < -ene) then
1007          minband = iband
1008        end if
1009      end do
1010      do iband=mband,1,-1
1011        if (maxval(eigen(iband,:,isppol))-fermie > ene) then
1012          maxband = iband
1013        end if
1014      end do
1015    end do ! isppol
1016 
1017  end if ! abs(energy_window)
1018 
1019 !=== Dump the results on file ===
1020  if (open_file(fname,msg,newunit=uvtk,status='unknown',form='formatted') /= 0) then
1021    ABI_DEALLOCATE(fulltoirred)
1022    MSG_WARNING(msg)
1023    ierr=ierr +1; RETURN
1024  end if
1025 
1026 !write header
1027  write(uvtk,"(a)") '# vtk DataFile Version 2.0'
1028  write(uvtk,"(a)") 'Eigen values for the Fermi surface'
1029  write(uvtk,"(a)") 'ASCII'
1030  write(uvtk,*) ''
1031  write(uvtk,"(a)") 'DATASET STRUCTURED_GRID'
1032  write(uvtk,"(a,3i6)") 'DIMENSIONS', nk1+1,nk2+1,nk3+1
1033  write(uvtk,"(a,i6,a)") 'POINTS',nkptfull,' float'
1034 
1035  do ik3 = 0, nk3
1036    do ik2 = 0, nk2
1037      do ik1 = 0, nk1
1038        write(uvtk,'(3es16.8)') dble(ik1)/nk1*gprimd(1,1)+ &
1039 &       dble(ik2)/nk2*gprimd(1,2)+ &
1040 &       dble(ik3)/nk3*gprimd(1,3), &
1041 &       dble(ik1)/nk1*gprimd(2,1)+ &
1042 &       dble(ik2)/nk2*gprimd(2,2)+ &
1043 &       dble(ik3)/nk3*gprimd(2,3), &
1044 &       dble(ik1)/nk1*gprimd(3,1)+ &
1045 &       dble(ik2)/nk2*gprimd(3,2)+ &
1046 &       dble(ik3)/nk3*gprimd(3,3)
1047      end do
1048    end do
1049  end do
1050 
1051 !print out data for all relevant bands and full kpt grid (redundant, yes)
1052 !for each kpt in full zone, find equivalent irred kpt and print eigenval
1053  write(uvtk,*) ''
1054  write(uvtk,"(a,i6)") 'POINT_DATA',nkptfull
1055  indx=0
1056  do iband=minband,maxband
1057    do isppol=1,nsppol
1058      if (minband+indx < 10) then
1059        format_str="(a14,i1,1X,a)"
1060      else
1061        format_str="(a14,i2,1X,a)"
1062      end if
1063      write(uvtk,format_str) 'SCALARS eigval', minband+indx, 'float 1'
1064      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1065      write(uvtk,*) ' '
1066      do kkpt = nk3/2+1, nk3+nk3/2+1
1067        do jkpt = nk2/2+1, nk2+nk2/2+1
1068          do ikpt = nk1/2+1, nk1+nk1/2+1
1069            ik1 = ikpt
1070            ik2 = jkpt
1071            ik3 = kkpt
1072            if (ikpt > nk1+1) ik1 = ikpt - nk1
1073            if (jkpt > nk2+1) ik2 = jkpt - nk2
1074            if (kkpt > nk3+1) ik3 = kkpt - nk3
1075 !          get the index with zyx order
1076            ikpt_fine = (ik1-1)*(nk2+1)*(nk3+1) + (ik2-1)*(nk3+1) + ik3
1077            write(uvtk,'(es16.8)') eigen(iband,fulltoirred(ikpt_fine),isppol)
1078          end do
1079        end do
1080      end do
1081      indx=indx+1
1082    end do
1083  end do
1084 
1085  write(uvtk,*) ''
1086  indx=0
1087  do iband=minband,maxband
1088    do isppol=1,nsppol
1089      if (minband+indx < 10) then
1090        format_str="(a10,i1,1X,a)"
1091      else
1092        format_str="(a10,i2,1X,a)"
1093      end if
1094      write(uvtk,format_str) 'SCALARS ve', minband+indx, 'float'
1095      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1096      write(uvtk,*) ' '
1097      do kkpt = nk3/2+1, nk3+nk3/2+1
1098        do jkpt = nk2/2+1, nk2+nk2/2+1
1099          do ikpt = nk1/2+1, nk1+nk1/2+1
1100            ik1 = ikpt
1101            ik2 = jkpt
1102            ik3 = kkpt
1103            if (ikpt > nk1+1) ik1 = ikpt - nk1
1104            if (jkpt > nk2+1) ik2 = jkpt - nk2
1105            if (kkpt > nk3+1) ik3 = kkpt - nk3
1106 !          write(uvtk,'(3i6,3es16.8)') ik1,ik2,ik3,v_surf(iband,ik1,ik2,ik3,1,isppol), &
1107 !          &                                                 v_surf(iband,ik1,ik2,ik3,2,isppol), &
1108 !          &                                                 v_surf(iband,ik1,ik2,ik3,3,isppol)
1109            write(uvtk,'(es16.8)') sqrt(v_surf(iband,ik1,ik2,ik3,1,isppol)*v_surf(iband,ik1,ik2,ik3,1,isppol)+ &
1110 &           v_surf(iband,ik1,ik2,ik3,2,isppol)*v_surf(iband,ik1,ik2,ik3,2,isppol)+ &
1111 &           v_surf(iband,ik1,ik2,ik3,3,isppol)*v_surf(iband,ik1,ik2,ik3,3,isppol))
1112          end do
1113        end do
1114      end do
1115      write(uvtk,format_str) 'SCALARS vx', minband+indx, 'float'
1116      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1117      write(uvtk,*) ' '
1118      do kkpt = nk3/2+1, nk3+nk3/2+1
1119        do jkpt = nk2/2+1, nk2+nk2/2+1
1120          do ikpt = nk1/2+1, nk1+nk1/2+1
1121            ik1 = ikpt
1122            ik2 = jkpt
1123            ik3 = kkpt
1124            if (ikpt > nk1+1) ik1 = ikpt - nk1
1125            if (jkpt > nk2+1) ik2 = jkpt - nk2
1126            if (kkpt > nk3+1) ik3 = kkpt - nk3
1127            write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,1,isppol)
1128          end do
1129        end do
1130      end do
1131      write(uvtk,format_str) 'SCALARS vy', minband+indx, 'float'
1132      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1133      write(uvtk,*) ' '
1134      do kkpt = nk3/2+1, nk3+nk3/2+1
1135        do jkpt = nk2/2+1, nk2+nk2/2+1
1136          do ikpt = nk1/2+1, nk1+nk1/2+1
1137            ik1 = ikpt
1138            ik2 = jkpt
1139            ik3 = kkpt
1140            if (ikpt > nk1+1) ik1 = ikpt - nk1
1141            if (jkpt > nk2+1) ik2 = jkpt - nk2
1142            if (kkpt > nk3+1) ik3 = kkpt - nk3
1143            write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,2,isppol)
1144          end do
1145        end do
1146      end do
1147      write(uvtk,format_str) 'SCALARS vz', minband+indx, 'float'
1148      write(uvtk,"(a)") 'LOOKUP_TABLE default'
1149      write(uvtk,*) ' '
1150      do kkpt = nk3/2+1, nk3+nk3/2+1
1151        do jkpt = nk2/2+1, nk2+nk2/2+1
1152          do ikpt = nk1/2+1, nk1+nk1/2+1
1153            ik1 = ikpt
1154            ik2 = jkpt
1155            ik3 = kkpt
1156            if (ikpt > nk1+1) ik1 = ikpt - nk1
1157            if (jkpt > nk2+1) ik2 = jkpt - nk2
1158            if (kkpt > nk3+1) ik3 = kkpt - nk3
1159            write(uvtk,'(es16.8)') v_surf(iband,ik1,ik2,ik3,3,isppol)
1160          end do
1161        end do
1162      end do
1163      indx=indx+1
1164    end do
1165  end do
1166 
1167  close (uvtk)
1168  ABI_DEALLOCATE(fulltoirred)
1169 
1170 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

PARENTS

      denfgr,exc_plot,m_kxc,m_nesting,m_wfd,spin_current

CHILDREN

      wrap2_pmhalf

SOURCE

165 subroutine printxsf(n1,n2,n3,datagrid,basis,origin,natom,ntypat,typat,xcart,znucl,nunit,realrecip)
166 
167 
168 !This section has been created automatically by the script Abilint (TD).
169 !Do not modify the following lines by hand.
170 #undef ABI_FUNC
171 #define ABI_FUNC 'printxsf'
172 !End of the abilint section
173 
174  implicit none
175 
176 !Arguments ------------------------------------
177 !scalars
178  integer,intent(in) :: n1,n2,n3,nunit,realrecip
179  integer,intent(in) :: natom,ntypat
180 !arrays
181  integer,intent(in) :: typat(natom)
182  real(dp),intent(in) :: basis(3,3),datagrid(n1*n2*n3),origin(3)
183  real(dp),intent(in) :: xcart(3,natom), znucl(ntypat)
184 
185 !Local variables-------------------------------
186 !scalars
187  integer :: ix,iy,iz,nslice,nsym,iatom
188  real(dp) :: fact
189 !arrays
190  real(dp) :: tau(3,natom)
191 
192 ! *************************************************************************
193 
194  DBG_ENTER("COLL")
195 
196  if (all(realrecip/= [0,1])) then
197    MSG_BUG(sjoin('The argument realrecip should be 0 or 1, received:', itoa(realrecip)))
198  end if
199 
200 !conversion between ABINIT default units and XCrysden units
201  fact=Bohr_Ang; if (realrecip ==1) fact=one/fact  !since we are in reciprocal space
202 
203 !TODO insert crystalline structure and dummy atoms in case of reciprocal space
204 !need to convert basis too
205 
206  write(nunit,'(1X,A)')  'DIM-GROUP'
207  write(nunit,*) '3  1'
208  write(nunit,'(1X,A)') 'PRIMVEC'
209  do iy = 1,3
210    write(nunit,'(3(ES17.10,2X))') (Bohr_Ang*basis(ix,iy), ix=1,3)
211  end do
212 !
213 !generate translated coordinates to fit origin shift
214 !
215  do iatom = 1,natom
216    tau (:,iatom) = xcart(:,iatom) - origin(:)
217  end do
218 
219  write(nunit,'(1X,A)') 'PRIMCOORD'
220  write(nunit,*) natom, ' 1'
221  do iatom = 1,natom
222    write(nunit,'(i9,3(3X,ES17.10))') NINT(znucl(typat(iatom))), &  ! WARNING alchemy not supported by XCrysden
223 &  Bohr_Ang*tau(1,iatom), &
224 &  Bohr_Ang*tau(2,iatom), &
225 &  Bohr_Ang*tau(3,iatom)
226  end do
227  write(nunit,'(1X,A)') 'ATOMS'
228  do iatom = 1,natom
229    write(nunit,'(i9,3(3X,ES17.10))') NINT(znucl(typat(iatom))), & ! WARNING alchemy not supported by XCrysden
230 &  Bohr_Ang*tau(1,iatom), &
231 &  Bohr_Ang*tau(2,iatom), &
232 &  Bohr_Ang*tau(3,iatom)
233  end do
234 
235  write(nunit,'(a)')' BEGIN_BLOCK_DATAGRID3D'
236  write(nunit,'(a)')' datagrid'
237  write(nunit,'(a)')' DATAGRID_3D_DENSITY'
238 !NOTE: XCrysden uses aperiodical data grid
239  write(nunit,*)n1+1,n2+1,n3+1
240  write(nunit,*)origin
241  write(nunit,*)basis(:,1)*fact
242  write(nunit,*)basis(:,2)*fact
243  write(nunit,*)basis(:,3)*fact
244 
245  nslice=1
246  do iz=1,n3
247    do iy=1,n2
248      write(nunit,'(8es16.8)') datagrid(1+n1*(nslice-1):n1+n1*(nslice-1)),datagrid(1+n1*(nslice-1))
249      nslice=nslice+1
250    end do
251    nsym=nslice-n2
252    write (nunit,'(8es16.8)') datagrid(1+n1*(nsym-1):n1+n1*(nsym-1)),datagrid(1+n1*(nsym-1))
253  end do
254 
255 !Now write upper plane
256  nslice=1
257  do iy=1,n2
258    write (nunit,'(8es16.8)') datagrid(1+n1*(nslice-1):n1+n1*(nslice-1)),datagrid(1+n1*(nslice-1))
259    nslice=nslice+1
260  end do
261 
262  nsym=nslice-n2
263  write (nunit,'(8es16.8)') datagrid(1+n1*(nsym-1):n1+n1*(nsym-1)),datagrid(1+n1*(nsym-1))
264 
265  write (nunit,'(a)')' END_DATAGRID_3D'
266  write (nunit,'(a)')' END_BLOCK_DATAGRID3D'
267 
268  DBG_EXIT("COLL")
269 
270 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)

PARENTS

      chiscwrt,ioniondist,newkpt,pawuj_det,pawuj_red,pawuj_utils,shellstruct

CHILDREN

      wrap2_pmhalf

SOURCE

 79 subroutine prmat (mat, ni, nj, mi, unitm)
 80 
 81 
 82 !This section has been created automatically by the script Abilint (TD).
 83 !Do not modify the following lines by hand.
 84 #undef ABI_FUNC
 85 #define ABI_FUNC 'prmat'
 86 !End of the abilint section
 87 
 88  implicit none
 89 
 90 !Arguments ------------------------------------
 91 !scalars
 92  integer,intent(in)           :: mi,ni,nj
 93  integer,intent(in), optional :: unitm
 94 !arrays
 95  real(dp),intent(in)          :: mat(mi,nj)
 96 
 97 !Local variables-------------------------------
 98 !scalars
 99  character(len=1000)    :: message
100  integer,parameter      :: nline=10
101  integer                :: ii,jj,jstart,jstop,unitn
102 
103 ! *************************************************************************
104 
105  if (present(unitm)) then ! standard printing to std_out
106    unitn=unitm
107  else
108    unitn=std_out
109  end if
110 
111  do  jstart = 1, nj, nline
112    jstop = min(nj, jstart+nline-1)
113    write(message, '(3x,10(i4,8x))' ) (jj,jj=jstart,jstop)
114    call wrtout(unitn,message,'COLL')
115  end do
116 
117  do ii = 1,ni
118    do jstart= 1, nj, nline
119      jstop = min(nj, jstart+nline-1)
120      if (jstart==1) then
121        write(message, '(i3,1p,10e12.4)' ) ii, (mat(ii,jj),jj=jstart,jstop)
122        call wrtout(unitn,message,'COLL')
123      else
124        write(message, '(3x,1p,10e12.4)' )    (mat(ii,jj),jj=jstart,jstop)
125        call wrtout(unitn,message,'COLL')
126      end if
127    end do
128  end do
129 
130 end subroutine prmat