TABLE OF CONTENTS


ABINIT/mrggkk [ Programs ]

[ Top ] [ Programs ]

NAME

 mrggkk

FUNCTION

 This program merges a GS file and several 1WF or GKK files for
 different q-vectors and perturbations.

COPYRIGHT

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

INPUTS

  (main routine)

OUTPUT

  (main routine)

NOTES

 GKK file structure is composed of header records and eigenvalue arrays,
 in binary or ascii:
   GS header = hdr
   GS eigenvalues = eigen
   number of perturbations = ntot
   for each perturbation
      1WF header = hdr1
      1st order eigenvalues = eigen1

SOURCE

 34 #if defined HAVE_CONFIG_H
 35 #include "config.h"
 36 #endif
 37 
 38 #include "abi_common.h"
 39 
 40 program mrggkk
 41 
 42  use defs_basis
 43  use m_build_info
 44  use m_abicore
 45  use m_xmpi
 46  use m_errors
 47  use m_wfk
 48  use m_nctk
 49 #ifdef HAVE_NETCDF
 50  use netcdf
 51 #endif
 52  use m_hdr
 53 
 54  use m_specialmsg,      only : specialmsg_getcount, herald
 55  use m_fstrings,        only : endswith, sjoin
 56  use m_io_tools,        only : flush_unit, open_file, file_exists
 57  use m_mpinfo,          only : destroy_mpi_enreg, initmpi_seq
 58 
 59  implicit none
 60 
 61 !Arguments ------------------------------------
 62 
 63 !Local variables-------------------------------
 64 !scalars
 65  integer,parameter :: unit1wf=22,unitgkk=24,unitgs=21,unitout=23,formeig0=0,formeig1=1
 66  integer :: binascii,fform,headform,i1wf,igkk,ik_ibz,ios,spin,mband
 67  integer :: n1wf,ngkk,ntot,ntotgkk,comm,iomode
 68  integer :: iband,jband,nband_k,rdwrout,ierr,ipos !base,
 69  real(dp) :: tolgkk=tol6
 70  character(len=1),parameter :: comment="#"
 71  character(len=24) :: codename
 72  character(len=500) :: message
 73  character(len=fnlen) :: file1wf,filegkk,filegs,outfile
 74  type(hdr_type) :: hdr,hdr1
 75  type(wfk_t) :: GS_wfk,PH_wfk
 76 !arrays
 77  real(dp),allocatable :: eig_k(:)
 78 
 79 ! *************************************************************************
 80 
 81 !Change communicator for I/O (mandatory!)
 82  call abi_io_redirect(new_io_comm=xmpi_world)
 83 
 84 !Initialize MPI
 85  call xmpi_init()
 86  comm = xmpi_world
 87 
 88 !Initialize memory profiling if it is activated
 89 !if a full abimem.mocc report is desired, set the argument of abimem_init to "2" instead of "0"
 90 !note that abimem.mocc files can easily be multiple GB in size so don't use this option normally
 91 #ifdef HAVE_MEM_PROFILING
 92  call abimem_init(0)
 93 #endif
 94 
 95  codename='MRGGKK'//repeat(' ',18)
 96 
 97 !write greating,read the file names, etc.
 98  call herald(codename,abinit_version,std_out)
 99 
100  write(message,'(17a)')&
101 & ' Files file format: ',ch10,ch10,&
102 & '  Name of the output file',ch10,&
103 & '  Integer flag: 0 --> binary output,   1 --> ascii formatted output',ch10,&
104 & '  Name of the groud state wavefunction file WF',ch10,&
105 & '  Number of 1WF, of GKK files, and number of 1WF files in all the GKK files',ch10,&
106 & '  Names of the 1WF files...',ch10,&
107 & '  Names of the GKK files...',ch10,ch10,&
108 & ' Enter name of output file: '
109  call wrtout(std_out,message,'COLL')
110 
111 !get file with filenames and number of 1wf files
112  read(*,'(a)') outfile
113  ipos=INDEX(outfile,comment)
114  if (ipos/=0) outfile=outfile(:ipos-1)
115 
116  read(*,*) binascii
117 
118  read(*,'(a)') filegs
119  ipos=INDEX(filegs,comment)
120  if (ipos/=0) filegs=filegs(:ipos-1)
121 
122  read(*,*) n1wf,ngkk,ntotgkk
123 
124  write(message,'(7a,i4,2a,i4,2a,i4,a)')&
125 & ' Output                     = ',trim(outfile),ch10,&
126 & ' Ground State file          = ',trim(filegs),ch10,&
127 & ' Number of 1WF files        = ',n1wf,ch10,&
128 & ' Number of GKK files        = ',ngkk,ch10,&
129 & ' Total Number of 1WF in GKK = ',ntotgkk,ch10
130  call wrtout(std_out,message,'COLL')
131 
132  iomode = IO_MODE_FORTRAN
133 #ifdef HAVE_MPI_IO
134  iomode = IO_MODE_MPI
135 #endif
136  !iomode = IO_MODE_FORTRAN
137 
138  ! Trick needed so that we can run the automatic tests with both Fortran and netcdf
139  if (.not. file_exists(filegs) .and. file_exists(nctk_ncify(filegs))) then
140    write(message, "(3a)")"- File: ",trim(filegs)," does not exist but found netcdf file with similar name."
141    call wrtout(std_out,message,'COLL')
142    filegs = nctk_ncify(filegs)
143    iomode = IO_MODE_ETSF
144  end if
145 
146  !output without rewinding the file
147  if (binascii == 0) then
148    ! open output file
149    ios = open_file(outfile,message,unit=unitout,form='unformatted')
150    rdwrout = 6
151  else if (binascii == 1) then
152    !  rdwrout=4 ! use for screen output and change writes of eigen to (*,*)
153    !  MJV 27/5/2008 removed 'new' constraint on gkk files: presume competent user!
154    ios = open_file(outfile,message,unit=unitout,form='formatted')
155    rdwrout = 4
156  else if (binascii == 2) then
157    !  this is for simple "short" output of the matrices, without headers or imaginary part
158    ios = open_file(outfile,message,unit=unitout,form='formatted')
159    rdwrout = 4
160  else
161    ABI_ERROR(' binascii must be between 0 and 2')
162  end if
163 
164  ABI_CHECK(ios==0,message)
165  rewind (unitout)
166 
167 !-------------------------------------------------------
168 !now read and write information for GS file
169 !-------------------------------------------------------
170 
171 !open GS wf file
172  call wrtout(std_out,' normal input for GS file',"COLL")
173 
174  call wfk_open_read(GS_wfk,filegs,formeig0,iomode,unitgs,comm)
175 
176 !Copy header of GS file to output.
177  if (binascii /= 2) then
178    if (rdwrout == 4) then
179      call GS_wfk%Hdr%echo(GS_wfk%fform, rdwrout)
180    else
181      call GS_wfk%Hdr%fort_write(unitout, GS_wfk%fform, ierr)
182      ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
183    end if
184  end if
185 
186  call wrtout(std_out,' header echoed to output file',"COLL")
187 
188  ABI_MALLOC(eig_k,(GS_wfk%mband))
189 
190 !Retrieve GS eigenvalues from GS wf file and echo to output
191  do spin=1,GS_wfk%nsppol
192    do ik_ibz=1,GS_wfk%nkpt
193      nband_k = GS_wfk%nband(ik_ibz,spin)
194 
195      call GS_wfk%read_eigk(ik_ibz,spin,xmpio_single,eig_k)
196      if (binascii==0) then
197        write(unitout) eig_k(1:nband_k)
198      else
199        write(unitout,*) eig_k(1:nband_k)
200      end if
201    end do
202  end do
203 
204  ABI_FREE(eig_k)
205 
206 !Close GS wf file
207  call GS_wfk%close()
208 
209  ntot = n1wf + ntotgkk
210  if (binascii==0) then
211    write (unitout) ntot
212  else
213    write (unitout,*) ntot
214  end if
215 
216 !-------------------------------------------------------
217 !now read and write information for 1WF files
218 !-------------------------------------------------------
219  do i1wf=1,n1wf
220 !  for each 1wf file, get name...
221    read(*,'(a)') file1wf
222    ipos=INDEX(file1wf,comment)
223    if (ipos/=0) file1wf=file1wf(:ipos-1)
224 
225    ! Trick needed so that we can run the automatic tests with both Fortran and netcdf
226    if (.not. file_exists(file1wf) .and. file_exists(nctk_ncify(file1wf))) then
227      write(message, "(3a)")"- File: ",trim(file1wf)," does not exist but found netcdf file with similar name."
228      call wrtout(std_out,message,'COLL')
229      file1wf = nctk_ncify(file1wf)
230      iomode = IO_MODE_ETSF
231    end if
232 
233 !  open 1wf file
234    call wrtout(std_out,' normal input for 1WF file ',"COLL")
235 
236    call wfk_open_read(PH_wfk,file1wf,formeig1,iomode,unit1wf,comm,Hdr_out=hdr1)
237 
238 !  copy header of 1WF file to output
239    if (binascii /= 2) then
240      if (rdwrout == 4) then
241        call hdr1%echo(PH_wfk%fform, rdwrout)
242      else
243        call hdr1%fort_write(unitout, PH_wfk%fform, ierr)
244        ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
245      end if
246    else
247      write (unitout,'(a,3E20.10)') "qpt ", hdr1%qptn
248      write (unitout,'(a,I6)') "pertnum ", hdr1%pertcase
249    end if
250 
251 !  retrieve 1WF <psi_k+q | H | psi_k> from 1wf file and echo to output
252    mband = maxval(hdr1%nband)
253    headform=hdr1%headform
254    ABI_MALLOC(eig_k,(2*mband*mband))
255 
256    ABI_CHECK(ALL(PH_wfk%nband == PH_wfk%nband(1,1)),"nband must be constant")
257 
258    do spin=1,hdr1%nsppol
259      do ik_ibz=1,hdr1%nkpt
260 !      write(std_out,*) 'spin,ik_ibz = ', spin,ik_ibz
261        nband_k = PH_wfk%nband(ik_ibz,spin)
262 
263        call PH_wfk%read_eigk(ik_ibz,spin,xmpio_single,eig_k)
264 
265        !base = 0
266        !do jband=1,nband_k
267        !  base = 2*(jband-1)*nband_k
268        !  do iband=1,2*nband_k
269        !    write(777,*) iband,jband,eig_k(base+iband)
270        !  end do
271        !end do
272 
273        if (binascii==0) then
274          write(unitout) eig_k(1:2*nband_k**2)
275        else if (binascii==1) then
276          write(unitout,*) eig_k(1:2*nband_k**2)
277        else if (binascii==2) then
278          do iband=1,nband_k
279            do jband=1,nband_k
280              if (abs(eig_k(2*nband_k*(iband-1)+2*(jband-1)+1))>tolgkk) then
281                write(unitout,'(E18.7, 2x)', ADVANCE='NO') eig_k(2*nband_k*(iband-1)+2*(jband-1)+1)
282              else
283                write(unitout,'(I18, 2x)', ADVANCE='NO') 0
284              end if
285              if (abs(eig_k(2*nband_k*(iband-1)+2*(jband-1)+2))>tolgkk) then
286                write(unitout,'(E18.7, 2x)', ADVANCE='NO') eig_k(2*nband_k*(iband-1)+2*(jband-1)+2)
287              else
288                write(unitout,'(I18, 2x)', ADVANCE='NO') 0
289              end if
290            end do
291            write(unitout,*)
292          end do
293          write(unitout,*)
294        end if
295 !
296      end do
297      if (binascii==2) write(unitout,'(2a)') ch10, ch10
298    end do
299 
300    ABI_FREE(eig_k)
301 
302    ! clean header to deallocate everything
303    call hdr1%free()
304    call PH_wfk%close()
305  end do
306 
307 !-------------------------------------------------------
308 !now read and write information for small GKK files
309 !-------------------------------------------------------
310  do igkk=1,ngkk
311 !
312 !  for each gkk file, get name...
313    read(*,'(a)') filegkk
314    ipos=INDEX(filegkk,comment)
315    if (ipos/=0) filegkk=filegkk(:ipos-1)
316 
317 !  open gkk file
318    call wrtout(std_out,' normal input for GKK file',"COLL")
319 
320    if (open_file(filegkk,message,unit=unitgkk,form='unformatted',status='old') /= 0) then
321      ABI_ERROR(message)
322    end if
323    rewind (unitgkk)
324 
325 !  read in header of GS file and eigenvalues
326 !  could force a comparison of header with global header above for consistency
327    call hdr_fort_read(hdr, unitgkk, fform)
328    ABI_CHECK(fform /= 0, sjoin("fform == 0 while reading:", filegkk))
329 
330    mband = maxval(hdr%nband)
331    ABI_MALLOC(eig_k,(mband))
332    call wrtout(std_out,'mrggkk : try to reread GS eigenvalues','COLL')
333 
334    do spin=1,hdr%nsppol
335      do ik_ibz=1,hdr%nkpt
336        nband_k = hdr%nband(ik_ibz + (spin-1)* hdr%nkpt)
337        read (unitgkk,IOSTAT=ierr) eig_k(1:nband_k)
338        ABI_CHECK(ierr==0,'error reading eigen from gkk file')
339      end do
340    end do
341    ABI_FREE(eig_k)
342 
343    read(unitgkk,IOSTAT=ierr) n1wf
344    ABI_CHECK(ierr==0,'error reading n1wf record')
345 
346    ABI_MALLOC(eig_k,(2*mband*mband))
347    do i1wf=1,n1wf
348 !    read in header of 1WF file
349      call hdr_fort_read(hdr1, unitgkk, fform)
350      if (fform == 0) then
351        write(message,'(a,i0,a)')' 1WF header number ',i1wf,' was mis-read. fform == 0'
352        ABI_ERROR(message)
353      end if
354 
355 !    copy header of 1WF file to output
356      if (binascii /= 2) then
357        if (rdwrout == 4) then
358          call hdr1%echo(fform, rdwrout)
359        else
360          call hdr1%fort_write(unitout, fform, ierr)
361          ABI_CHECK(ierr == 0 , "hdr_fort_write returned ierr != 0")
362        end if
363      else
364        write (unitout,'(a,3E20.10)') "qpt ", hdr1%qptn
365        write (unitout,'(a,I6)') "pertnum ", hdr1%pertcase
366      end if
367 
368 !    retrieve 1WF <psi_k+q | H | psi_k> from gkk file and echo to output
369      do spin=1,hdr1%nsppol
370        do ik_ibz=1,hdr1%nkpt
371          nband_k = hdr%nband(ik_ibz + (spin-1)* hdr1%nkpt)
372          read (unitgkk,IOSTAT=ierr) eig_k(1:2*nband_k**2)
373          if (ierr /= 0) write (std_out,*) 'error reading eigen2 from gkk file',spin,ik_ibz
374 
375          if (binascii==0) then
376            write (unitout) eig_k(1:2*nband_k**2)
377          else if (binascii==1) then
378            write (unitout,*) eig_k(1:2*nband_k**2)
379          else if (binascii==2) then
380            do iband=1,nband_k
381              do jband=1,nband_k
382                if (abs(eig_k(2*nband_k*(iband-1)+2*(jband-1)+1))>tolgkk) then
383                  write(unitout,'(E18.7, 2x)', ADVANCE='NO') eig_k(2*nband_k*(iband-1)+2*(jband-1)+1)
384                else
385                  write(unitout,'(I18, 2x)', ADVANCE='NO') 0
386                end if
387                if (abs(eig_k(2*nband_k*(iband-1)+2*(jband-1)+2))>tolgkk) then
388                  write(unitout,'(E18.7, 2x)', ADVANCE='NO') eig_k(2*nband_k*(iband-1)+2*(jband-1)+2)
389                else
390                  write(unitout,'(I18, 2x)', ADVANCE='NO') 0
391                end if
392              end do
393              write(unitout,*)
394            end do
395            write(unitout,*)
396          end if
397 !
398        end do
399        if (binascii==2) write(unitout,'(2a)') ch10, ch10
400      end do
401      call hdr1%free()
402    end do !  end loop over 1wf segments in small gkk file
403 
404    ABI_FREE(eig_k)
405 
406    close (unitgkk)
407    call hdr%free()
408  end do !end loop over small gkk files
409 
410  close(unitout)
411 
412  write(message,'(2a)')ch10,' Done'
413  call wrtout(std_out,message,'COLL')
414 
415  call flush_unit(std_out)
416 
417  call abinit_doctor("__mrggkk")
418 
419  call xmpi_end()
420 
421  end program mrggkk