TABLE OF CONTENTS
- ABINIT/m_pptools
- m_pptools/print_fofr_cube
- m_pptools/print_fofr_ri
- m_pptools/print_fofr_xyzri
- m_pptools/printbxsf
- m_pptools/printvtk
- m_pptools/printxsf
- m_pptools/prmat
ABINIT/m_pptools [ Modules ]
NAME
m_pptools
FUNCTION
Helper functions used for post-processing.
COPYRIGHT
Copyright (C) 2002-2024 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