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