LCOV - code coverage report
Current view: top level - src/98_main - fold2Bloch.F90 (source / functions) Coverage Total Hit
Test: coverage.info Lines: 81.3 % 139 113
Test Date: 2026-09-19 17:42:43 Functions: 100.0 % 2 2

            Line data    Source code
       1              : !!****p* ABINIT/fold2Bloch
       2              : !! NAME
       3              : !! fold2Bloch
       4              : !!
       5              : !! FUNCTION
       6              : !! Main routine for the unfolding of the wavefuntion.
       7              : !!
       8              : !! COPYRIGHT
       9              : !! Copyright (C) 2014-2026 ABINIT group (AB)
      10              : !! This file is distributed under the terms of the
      11              : !! GNU General Public License, see ~abinit/COPYING
      12              : !! or http://www.gnu.org/copyleft/gpl.txt .
      13              : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
      14              : !!
      15              : !! INPUTS
      16              : !!  (main program)
      17              : !!
      18              : !! OUTPUT
      19              : !!  (main program)
      20              : !!
      21              : !! NOTES
      22              : !! folds= Array of folds in X,Y, and Z directions
      23              : !!
      24              : !! SOURCE
      25              : 
      26              : #if defined HAVE_CONFIG_H
      27              : #include "config.h"
      28              : #endif
      29              : 
      30              : #include "abi_common.h"
      31              : 
      32            2 : program fold2Bloch
      33              : 
      34            2 :  use defs_basis
      35              :  use m_errors
      36              :  use m_abicore
      37              :  use m_wfk
      38              :  use m_xmpi
      39              :  use m_nctk
      40              :  use m_hdr
      41              :  use m_crystal
      42              :  use m_ebands
      43              :  use m_fold2block
      44              :  use netcdf
      45              : 
      46              :  use m_fstrings,       only : strcat
      47              :  use m_io_tools,       only : get_unit, iomode_from_fname, open_file, prompt
      48              : 
      49              :  implicit none
      50              : 
      51              : !Arguments --------------------------------------------------------------
      52              : 
      53              : !Local variables-------------------------------
      54              : !scalars
      55              : integer :: ikpt, iband,nspinor,nsppol,mband,nkpt,mcg,csppol, cspinor, nfold, iss, ii
      56              : integer :: comm, my_rank, nargs, iomode, ncid, ncerr, fform, kunf_varid, weights_varid, eigunf_varid
      57              : integer :: cg_b, count, outfile, outfile1, outfile2, lwcg, hicg, pos
      58              : character(fnlen) :: fname, outname,seedname
      59              : character(len=500) :: msg
      60            8 : type(wfk_t) :: wfk
      61          104 : type(crystal_t) :: cryst
      62            2 : type(ebands_t) :: ebands
      63              : !arrays
      64              : integer :: folds(3),fold_matrix(3,3)
      65            2 : integer, allocatable :: kg(:,:),nband(:), npwarr(:)
      66            4 : real(dp), allocatable :: cg(:,:), eig(:),kpts(:,:), weights(:),coefc(:,:), nkval(:,:)
      67              : 
      68              : !*************************************************************************
      69              : 
      70              : !0) Change communicator for I/O (mandatory!)
      71            2 :  call abi_io_redirect(new_io_comm=xmpi_world)
      72              : 
      73            2 :  call xmpi_init()
      74            2 :  comm = xmpi_world; my_rank = xmpi_comm_rank(xmpi_world)
      75              : 
      76              : !Initialize memory profiling if it is activated
      77              : !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
      78              : !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
      79              : #ifdef HAVE_MEM_PROFILING
      80              :  call abimem_init(0)
      81              : #endif
      82              : 
      83            2 :  if (xmpi_comm_size(comm) /= 1) then
      84            0 :    ABI_ERROR("fold2bloch not programmed for parallel execution.")
      85              :  end if
      86              : 
      87            2 :  nargs = command_argument_count()
      88              : 
      89            2 :  if (nargs == 0) then
      90            2 :    call prompt("Enter WFK file name:", fname)
      91            2 :    call prompt("Enter x y z integers giving the multiplicity:", folds)
      92              :  else
      93            0 :    call getargs(folds, fname) !Process command line arguments
      94              :  end if
      95              :  ! Use fold_matrix instead of folds(1:3) to prepare possible generalization.
      96            2 :  fold_matrix = 0
      97            8 :  do ii=1,3
      98            8 :    fold_matrix(ii,ii) = folds(ii)
      99              :  end do
     100              : 
     101              :  ! Test if the netcdf library supports MPI-IO
     102              :  !call nctk_test_mpiio()
     103              : 
     104            2 :  if (nctk_try_fort_or_ncfile(fname, msg) /= 0) then
     105            0 :    ABI_ERROR(msg)
     106              :  end if
     107              : 
     108            2 :  pos=INDEX(fname, "_")
     109            2 :  write(seedname,'(a)') fname(1:pos-1)
     110              : 
     111            2 :  write(std_out,*) '         '//achar(27)//'[97m ***********************' !print program header in pearl white
     112            2 :  write(std_out,*) '          ** Fold2Bloch V 1.1  **'
     113            2 :  write(std_out,*) '          **Build  Mar 16, 2015**'
     114            2 :  write(std_out,*) '          ***********************'//achar(27)//'[0m'
     115              : 
     116            2 :  ebands = wfk_read_ebands(fname, xmpi_comm_self)
     117            2 :  iomode = iomode_from_fname(fname)
     118            2 :  call wfk%open_read(fname, 0, iomode, get_unit(), comm)
     119              : 
     120            2 :  nkpt=wfk%hdr%nkpt
     121            2 :  ABI_MALLOC(npwarr,(nkpt))
     122            2 :  ABI_MALLOC(nband,(nkpt))
     123            2 :  ABI_MALLOC(kpts,(3,nkpt))
     124              : 
     125            2 :  nsppol=wfk%hdr%nsppol
     126            2 :  nspinor=wfk%hdr%nspinor
     127           88 :  npwarr=wfk%hdr%npwarr
     128          340 :  kpts=wfk%hdr%kptns
     129           88 :  nband=wfk%hdr%nband
     130           86 :  mband=maxval(nband)
     131           86 :  mcg=maxval(npwarr)*nspinor*mband
     132            8 :  nfold = product(folds)
     133              : 
     134            2 :  cryst = wfk%hdr%get_crystal()
     135              : 
     136            2 :  NCF_CHECK(nctk_open_create(ncid, strcat(seedname, "_FOLD2BLOCH.nc"), xmpi_comm_self))
     137            2 :  fform = fform_from_ext("FOLD2BLOCH.nc")
     138            2 :  NCF_CHECK(wfk%hdr%ncwrite(ncid, fform, nc_define=.True.))
     139            2 :  NCF_CHECK(cryst%ncwrite(ncid))
     140            2 :  NCF_CHECK(ebands%ncwrite(ncid))
     141              : 
     142              :  ncerr = nctk_def_dims(ncid, [ &
     143              :  nctkdim_t("nk_unfolded", nkpt * nfold), &
     144            6 :  nctkdim_t("nsppol_times_nspinor", wfk%hdr%nsppol * wfk%hdr%nspinor)], defmode=.True.)
     145            2 :  NCF_CHECK(ncerr)
     146              :  ncerr = nctk_def_arrays(ncid, [ &
     147              :  nctkarr_t("fold_matrix", "int", "number_of_reduced_dimensions, number_of_reduced_dimensions"), &
     148              :  nctkarr_t("reduced_coordinates_of_unfolded_kpoints", "dp", "number_of_reduced_dimensions, nk_unfolded"), &
     149              :  nctkarr_t("unfolded_eigenvalues", "dp", "max_number_of_states, nk_unfolded, number_of_spins"), &
     150              :  nctkarr_t("spectral_weights", "dp", "max_number_of_states, nk_unfolded, nsppol_times_nspinor") &
     151           10 :  ])
     152            2 :  NCF_CHECK(ncerr)
     153            2 :  NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_unfolded_kpoints", kunf_varid))
     154            2 :  NCF_CHECK(nf90_inq_varid(ncid, "unfolded_eigenvalues", eigunf_varid))
     155            2 :  NCF_CHECK(nf90_inq_varid(ncid, "spectral_weights", weights_varid))
     156            2 :  NCF_CHECK(nctk_set_datamode(ncid))
     157            2 :  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "fold_matrix"), fold_matrix))
     158              : 
     159            2 :  call cryst%free()
     160            2 :  call ebands%free()
     161              : 
     162            4 :  do csppol=1, nsppol
     163            2 :    if (nsppol==1) then !Determine spin polarization for output file
     164            2 :      outname=trim(seedname)//".f2b"
     165            0 :    elseif ((nsppol==2).and.(csppol==1)) then
     166            0 :      outname=trim(seedname)//"_UP.f2b"
     167            0 :      write(std_out,*) "     ===================="
     168            0 :      write(std_out,*) "     SPIN POLARIZATION UP"
     169            0 :      write(std_out,*) "     ===================="
     170            0 :    elseif ((nsppol==2).and.(csppol==2)) then
     171            0 :      outname=trim(seedname)//"_DOWN.f2b"
     172            0 :      write(std_out,*) "     ======================"
     173            0 :      write(std_out,*) "     SPIN POLARIZATION DOWN"
     174            0 :      write(std_out,*) "     ======================"
     175              :    end if
     176            2 :    if (nspinor==2) then
     177              :      !open output file
     178            0 :      if (open_file(trim(seedname)//"_SPOR_1.f2b", msg, newunit=outfile1, form="formatted", status="unknown") /= 0) then
     179            0 :        ABI_ERROR(msg)
     180              :      end if
     181            0 :      if (open_file(trim(seedname)//"_SPOR_2.f2b", msg, newunit=outfile2, form="formatted", status="unknown") /= 0) then
     182            0 :        ABI_ERROR(msg)
     183              :      end if
     184              :    else
     185            2 :      if (open_file(outname, msg, newunit=outfile1,form="formatted", status="unknown") /= 0) then
     186            0 :        ABI_ERROR(msg)
     187              :      end if
     188              :    end if
     189              : 
     190           86 :    do ikpt=1, nkpt !For each K point
     191           84 :      ABI_MALLOC(cg,(2,mcg))
     192           84 :      ABI_MALLOC(eig,((2*mband)**0*mband))
     193           84 :      ABI_MALLOC(kg,(3,npwarr(ikpt)))
     194           84 :      ABI_MALLOC(coefc,(2,nspinor*npwarr(ikpt)))
     195           84 :      ABI_MALLOC(weights, (nfold))
     196           84 :      ABI_MALLOC(nkval,(3, nfold))
     197           84 :      call progress(ikpt,nkpt,kpts(:,ikpt)) !Write progress information
     198              : 
     199              :      !Read a block of data
     200          252 :      call wfk%read_band_block([1, nband(ikpt)], ikpt, csppol, xmpio_single, kg_k=kg, cg_k=cg, eig_k=eig)
     201              : 
     202              :      !Determine unfolded K point states
     203           84 :      call newk(kpts(1,ikpt),kpts(2,ikpt),kpts(3,ikpt),folds(1),folds(2),folds(3),nkval)
     204           84 :      if (csppol == 1) then
     205          420 :        NCF_CHECK(nf90_put_var(ncid, kunf_varid, nkval, start=[1, 1 + (ikpt-1) * nfold], count=[3, nfold]))
     206              :      end if
     207              : 
     208           84 :      cg_b=1
     209          420 :      do iband=1, nband(ikpt) !Foe each Eigenvalue
     210       648264 :        coefc=cg(:,cg_b:(cg_b+nspinor*npwarr(ikpt)-1)) !Split coefficients per eigen value according to the number of "kg"
     211          672 :        do cspinor=1,nspinor
     212          336 :          if (cspinor==1) then
     213          336 :            outfile=outfile1
     214          336 :            lwcg=1
     215          336 :            hicg=npwarr(ikpt)
     216              :          else
     217              :           ! Move coefficient span to spinor 2
     218            0 :            outfile=outfile2
     219            0 :            lwcg=npwarr(ikpt)+1
     220            0 :            hicg=npwarr(ikpt)*nspinor
     221              :          end if
     222          336 :          call sortc(folds(1),folds(2),folds(3),kg,coefc(:,lwcg:hicg),npwarr(ikpt),weights)
     223              :          ! Write out results, format: new k states(x, y, and z), eigenvalue, weight
     224         2352 :          do count=1, nfold
     225         2352 :            write(outfile,50) nkval(1,count),nkval(2,count),nkval(3,count),eig(iband),weights(count)
     226              :            50 format(f11.6, f11.6, f11.6, f11.6, f11.6)
     227              :          end do
     228          336 :          iss = csppol; if (nspinor == 2) iss = cspinor
     229              :          ncerr = nf90_put_var(ncid, weights_varid, weights, start=[iband, 1 + (ikpt-1) * nfold, iss], &
     230         3360 :          stride=[mband, 1, 1], count=[1, nfold, 1])
     231          336 :          NCF_CHECK(ncerr)
     232          672 :          if (cspinor == 1) then
     233         2352 :            weights = eig(iband) ! Use weights as workspace array.
     234              :            ncerr = nf90_put_var(ncid, eigunf_varid, weights, start=[iband, 1 + (ikpt-1) * nfold, csppol], &
     235         3360 :            stride=[mband, 1, 1], count=[1, nfold, 1])
     236              :              !count=[1, nfold, 1])
     237          336 :            NCF_CHECK(ncerr)
     238              :          end if
     239              :        end do ! cspinor
     240          420 :        cg_b=cg_b+nspinor*npwarr(ikpt) !shift coefficient pointer for next eigenvalue
     241              :      end do ! iband
     242              : 
     243           84 :      ABI_FREE(cg)
     244           84 :      ABI_FREE(eig)
     245           84 :      ABI_FREE(kg)
     246           84 :      ABI_FREE(coefc)
     247           84 :      ABI_FREE(weights)
     248           86 :      ABI_FREE(nkval)
     249              :    end do
     250            4 :    if (nspinor==2) then
     251            0 :      close(outfile1) !close output file
     252            0 :      close(outfile2)
     253              :    else
     254            2 :      close(outfile1)
     255              :    end if
     256              :  end do
     257            2 :  call wfk%close()
     258              : 
     259            2 :  ABI_FREE(kpts)
     260            2 :  ABI_FREE(nband)
     261            2 :  ABI_FREE(npwarr)
     262              : 
     263              : ! Print summary
     264            2 :  write(std_out,*) '    '//achar(27)//'[97m Number of K points processed:', nkpt
     265            2 :  if (nsppol==2) then
     266            0 :    write(std_out,*) '     Data was written to: ', trim(seedname)//"_UP.f2b", " & ", trim(seedname)//"_DOWN.f2b"
     267              :  else
     268            2 :    if (nspinor==2) then
     269            0 :      write(std_out,*) '     Data was written to: ', trim(seedname)//"_SPOR_1.f2b", " & ", trim(seedname)//"_SPOR_2.f2b"
     270              :    else
     271            2 :      write(std_out,*) '     Data was written to: ', trim(seedname)//".f2b"
     272              :    end if
     273              :  end if
     274            2 :  write(std_out,*) '     Data format: KX, KY, KZ, Eigenvalue(Ha), Weight'//achar(27)//'[0m'
     275              : 
     276            2 :  NCF_CHECK(nf90_close(ncid))
     277              : 
     278              : !Write information on file about the memory before ending mpi module, if memory profiling is enabled
     279            2 :  call abinit_doctor("__fold2bloch")
     280              : 
     281            2 :  call xmpi_end()
     282              : 
     283            0 :  end program fold2Bloch
     284              : !!***
        

Generated by: LCOV version 2.3-1