TABLE OF CONTENTS


ABINIT/fold2Bloch [ Programs ]

[ Top ] [ Programs ]

NAME

 fold2Bloch

FUNCTION

 Main routine for the unfolding of the wavefuntion.

COPYRIGHT

 Copyright (C) 2014-2018 ABINIT group (AB)
 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 .

INPUTS

  (main program)

OUTPUT

  (main program)

NOTES

 folds= Array of folds in X,Y, and Z directions

PARENTS

CHILDREN

      abi_io_redirect,abimem_init,abinit_doctor,crystal_free,crystal_from_hdr
      ebands_free,getargs,newk,progress,prompt,sortc,wfk_close,wfk_open_read
      wfk_read_band_block,xmpi_end,xmpi_init

SOURCE

 34 #if defined HAVE_CONFIG_H
 35 #include "config.h"
 36 #endif
 37 
 38 #include "abi_common.h"
 39 
 40 program fold2Bloch
 41 
 42  use defs_basis
 43  use m_errors
 44  use m_abicore
 45  use m_wfk
 46  use m_xmpi
 47  use m_nctk
 48  use m_hdr
 49  use m_crystal
 50  use m_crystal_io
 51  use m_ebands
 52  use m_fold2block
 53 #ifdef HAVE_NETCDF
 54  use netcdf
 55 #endif
 56 
 57  use m_fstrings,       only : strcat
 58  use m_io_tools,       only : get_unit, iomode_from_fname, open_file, prompt
 59  use defs_datatypes,   only : ebands_t
 60 
 61 !This section has been created automatically by the script Abilint (TD).
 62 !Do not modify the following lines by hand.
 63 #undef ABI_FUNC
 64 #define ABI_FUNC 'fold2Bloch'
 65 !End of the abilint section
 66 
 67 implicit none
 68 
 69 !Arguments --------------------------------------------------------------
 70 
 71 !Local variables-------------------------------
 72 !scalars
 73 integer :: ikpt, iband,nspinor,nsppol,mband,nkpt,mcg,csppol, cspinor, nfold, iss, ii
 74 integer :: comm, my_rank, nargs, iomode, ncid, ncerr, fform, timrev, kunf_varid, weights_varid, eigunf_varid
 75 integer :: cg_b, count, outfile, outfile1, outfile2, lwcg, hicg, pos
 76 character(fnlen) :: fname, outname,seedname
 77 character(len=500) :: msg
 78 type(wfk_t) :: wfk
 79 type(crystal_t) :: cryst
 80 type(ebands_t) :: ebands
 81 !arrays
 82 integer :: folds(3),fold_matrix(3,3)
 83 integer, allocatable :: kg(:,:),nband(:), npwarr(:)
 84 real(dp), allocatable :: cg(:,:), eig(:),kpts(:,:), weights(:),coefc(:,:), nkval(:,:)
 85 
 86 !*************************************************************************
 87 
 88 !0) Change communicator for I/O (mandatory!)
 89  call abi_io_redirect(new_io_comm=xmpi_world)
 90 
 91  call xmpi_init()
 92  comm = xmpi_world; my_rank = xmpi_comm_rank(xmpi_world)
 93 
 94 !Initialize memory profiling if it is activated
 95 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
 96 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
 97 #ifdef HAVE_MEM_PROFILING
 98  call abimem_init(0)
 99 #endif
100 
101  if (xmpi_comm_size(comm) /= 1) then
102    MSG_ERROR("fold2bloch not programmed for parallel execution.")
103  end if
104 
105  nargs = command_argument_count()
106 
107  if (nargs == 0) then
108    call prompt("Enter WFK file name:", fname)
109    call prompt("Enter x y z integers giving the multiplicity:", folds)
110  else
111    call getargs(folds, fname) !Process command line arguments
112  end if
113  ! Use fold_matrix instead of folds(1:3) to prepare possible generalization.
114  fold_matrix = 0
115  do ii=1,3
116    fold_matrix(ii,ii) = folds(ii)
117  end do
118 
119  ! Test if the netcdf library supports MPI-IO
120  !call nctk_test_mpiio()
121 
122  if (nctk_try_fort_or_ncfile(fname, msg) /= 0) then
123    MSG_ERROR(msg)
124  end if
125 
126  pos=INDEX(fname, "_")
127  write(seedname,'(a)') fname(1:pos-1)
128 
129  write(std_out,*) '         '//achar(27)//'[97m ***********************' !print program header in pearl white
130  write(std_out,*) '          ** Fold2Bloch V 1.1  **'
131  write(std_out,*) '          **Build  Mar 16, 2015**'
132  write(std_out,*) '          ***********************'//achar(27)//'[0m'
133 
134  ebands = wfk_read_ebands(fname, xmpi_comm_self)
135  iomode = iomode_from_fname(fname)
136  call wfk_open_read(wfk,fname,0,iomode,get_unit(),comm)
137 
138  nkpt=wfk%hdr%nkpt
139  ABI_ALLOCATE(npwarr,(nkpt))
140  ABI_ALLOCATE(nband,(nkpt))
141  ABI_ALLOCATE(kpts,(3,nkpt))
142 
143  nsppol=wfk%hdr%nsppol
144  nspinor=wfk%hdr%nspinor
145  npwarr=wfk%hdr%npwarr
146  kpts=wfk%hdr%kptns
147  nband=wfk%hdr%nband
148  mband=maxval(nband)
149  mcg=maxval(npwarr)*nspinor*mband
150  nfold = product(folds)
151 
152 #ifdef HAVE_NETCDF
153  timrev = 2; if (any(wfk%hdr%kptopt == [3, 4])) timrev = 1
154  call crystal_from_hdr(cryst, wfk%hdr, timrev)
155 
156  NCF_CHECK(nctk_open_create(ncid, strcat(seedname, "_FOLD2BLOCH.nc"), xmpi_comm_self))
157  fform = fform_from_ext("FOLD2BLOCH.nc")
158  NCF_CHECK(hdr_ncwrite(wfk%hdr, ncid, fform, nc_define=.True.))
159  NCF_CHECK(crystal_ncwrite(cryst, ncid))
160  NCF_CHECK(ebands_ncwrite(ebands, ncid))
161 
162  ncerr = nctk_def_dims(ncid, [ &
163  nctkdim_t("nk_unfolded", nkpt * nfold), &
164  nctkdim_t("nsppol_times_nspinor", wfk%hdr%nsppol * wfk%hdr%nspinor)], defmode=.True.)
165  NCF_CHECK(ncerr)
166  ncerr = nctk_def_arrays(ncid, [ &
167  nctkarr_t("fold_matrix", "int", "number_of_reduced_dimensions, number_of_reduced_dimensions"), &
168  nctkarr_t("reduced_coordinates_of_unfolded_kpoints", "dp", "number_of_reduced_dimensions, nk_unfolded"), &
169  nctkarr_t("unfolded_eigenvalues", "dp", "max_number_of_states, nk_unfolded, number_of_spins"), &
170  nctkarr_t("spectral_weights", "dp", "max_number_of_states, nk_unfolded, nsppol_times_nspinor") &
171  ])
172  NCF_CHECK(ncerr)
173  NCF_CHECK(nf90_inq_varid(ncid, "reduced_coordinates_of_unfolded_kpoints", kunf_varid))
174  NCF_CHECK(nf90_inq_varid(ncid, "unfolded_eigenvalues", eigunf_varid))
175  NCF_CHECK(nf90_inq_varid(ncid, "spectral_weights", weights_varid))
176  NCF_CHECK(nctk_set_datamode(ncid))
177  NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "fold_matrix"), fold_matrix))
178  call crystal_free(cryst)
179 #endif
180 
181  call ebands_free(ebands)
182 
183  do csppol=1, nsppol
184    if (nsppol==1) then !Determine spin polarization for output file
185      outname=trim(seedname)//".f2b"
186    elseif ((nsppol==2).and.(csppol==1)) then
187      outname=trim(seedname)//"_UP.f2b"
188      write(std_out,*) "     ===================="
189      write(std_out,*) "     SPIN POLARIZATION UP"
190      write(std_out,*) "     ===================="
191    elseif ((nsppol==2).and.(csppol==2)) then
192      outname=trim(seedname)//"_DOWN.f2b"
193      write(std_out,*) "     ======================"
194      write(std_out,*) "     SPIN POLARIZATION DOWN"
195      write(std_out,*) "     ======================"
196    end if
197    if (nspinor==2) then
198      !open output file
199      if (open_file(trim(seedname)//"_SPOR_1.f2b", msg, newunit=outfile1, form="formatted", status="unknown") /= 0) then
200        MSG_ERROR(msg)
201      end if
202      if (open_file(trim(seedname)//"_SPOR_2.f2b", msg, newunit=outfile2, form="formatted", status="unknown") /= 0) then
203        MSG_ERROR(msg)
204      end if
205    else
206      if (open_file(outname, msg, newunit=outfile1,form="formatted", status="unknown") /= 0) then
207        MSG_ERROR(msg)
208      end if
209    end if
210 
211    do ikpt=1, nkpt !For each K point
212      ABI_ALLOCATE(cg,(2,mcg))
213      ABI_ALLOCATE(eig,((2*mband)**0*mband))
214      ABI_ALLOCATE(kg,(3,npwarr(ikpt)))
215      ABI_ALLOCATE(coefc,(2,nspinor*npwarr(ikpt)))
216      ABI_ALLOCATE(weights, (nfold))
217      ABI_ALLOCATE(nkval,(3, nfold))
218      call progress(ikpt,nkpt,kpts(:,ikpt)) !Write progress information
219 
220      !Read a block of data
221      call wfk_read_band_block(wfk, [1, nband(ikpt)], ikpt, csppol, xmpio_single, kg_k=kg, cg_k=cg, eig_k=eig)
222 
223      !Determine unfolded K point states
224      call newk(kpts(1,ikpt),kpts(2,ikpt),kpts(3,ikpt),folds(1),folds(2),folds(3),nkval)
225 #ifdef HAVE_NETCDF
226      if (csppol == 1) then
227        NCF_CHECK(nf90_put_var(ncid, kunf_varid, nkval, start=[1, 1 + (ikpt-1) * nfold], count=[3, nfold]))
228      end if
229 #endif
230 
231      cg_b=1
232      do iband=1, nband(ikpt) !Foe each Eigenvalue
233        coefc=cg(:,cg_b:(cg_b+nspinor*npwarr(ikpt)-1)) !Split coefficients per eigen value according to the number of "kg"
234        do cspinor=1,nspinor
235          if (cspinor==1) then
236            outfile=outfile1
237            lwcg=1
238            hicg=npwarr(ikpt)
239          else
240           ! Move coefficient span to spinor 2
241            outfile=outfile2
242            lwcg=npwarr(ikpt)+1
243            hicg=npwarr(ikpt)*nspinor
244          end if
245          call sortc(folds(1),folds(2),folds(3),kg,coefc(:,lwcg:hicg),npwarr(ikpt),weights)
246          ! Write out results, format: new k states(x, y, and z), eigenvalue, weight
247          do count=1, nfold
248            write(outfile,50) nkval(1,count),nkval(2,count),nkval(3,count),eig(iband),weights(count)
249            50 format(f11.6, f11.6, f11.6, f11.6, f11.6)
250          end do
251 #ifdef HAVE_NETCDF
252          iss = csppol; if (nspinor == 2) iss = cspinor
253          ncerr = nf90_put_var(ncid, weights_varid, weights, start=[iband, 1 + (ikpt-1) * nfold, iss], &
254          stride=[mband, 1, 1], count=[1, nfold, 1])
255          NCF_CHECK(ncerr)
256          if (cspinor == 1) then
257            weights = eig(iband) ! Use weights as workspace array.
258            ncerr = nf90_put_var(ncid, eigunf_varid, weights, start=[iband, 1 + (ikpt-1) * nfold, csppol], &
259            stride=[mband, 1, 1], count=[1, nfold, 1])
260              !count=[1, nfold, 1])
261            NCF_CHECK(ncerr)
262          end if
263 #endif
264        end do ! cspinor
265        cg_b=cg_b+nspinor*npwarr(ikpt) !shift coefficient pointer for next eigenvalue
266      end do ! iband
267 
268      ABI_DEALLOCATE(cg)
269      ABI_DEALLOCATE(eig)
270      ABI_DEALLOCATE(kg)
271      ABI_DEALLOCATE(coefc)
272      ABI_DEALLOCATE(weights)
273      ABI_DEALLOCATE(nkval)
274    end do
275    if (nspinor==2) then
276      close(outfile1) !close output file
277      close(outfile2)
278    else
279      close(outfile1)
280    end if
281  end do
282  call wfk_close(wfk)
283 
284  ABI_FREE(kpts)
285  ABI_FREE(nband)
286  ABI_FREE(npwarr)
287 
288 ! Print summary
289  write(std_out,*) '    '//achar(27)//'[97m Number of K points processed:', nkpt
290  if (nsppol==2) then
291    write(std_out,*) '     Data was written to: ', trim(seedname)//"_UP.f2b", " & ", trim(seedname)//"_DOWN.f2b"
292  else
293    if (nspinor==2) then
294      write(std_out,*) '     Data was written to: ', trim(seedname)//"_SPOR_1.f2b", " & ", trim(seedname)//"_SPOR_2.f2b"
295    else
296      write(std_out,*) '     Data was written to: ', trim(seedname)//".f2b"
297    end if
298  end if
299  write(std_out,*) '     Data format: KX, KY, KZ, Eigenvalue(Ha), Weight'//achar(27)//'[0m'
300 
301 #ifdef HAVE_NETCDF
302  NCF_CHECK(nf90_close(ncid))
303 #endif
304 
305 !Write information on file about the memory before ending mpi module, if memory profiling is enabled
306  call abinit_doctor("__fold2bloch")
307 
308  call xmpi_end()
309 
310  end program fold2Bloch