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-2018 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

PARENTS

CHILDREN

      abi_io_redirect,abimem_init,abinit_doctor,destroy_mpi_enreg,flush_unit
      hdr_echo,hdr_fort_read,hdr_fort_write,hdr_free,herald,initmpi_seq
      wfk_close,wfk_open_read,wfk_read_eigk,wrtout,xmpi_end,xmpi_init

SOURCE

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