Branch data 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 : : !!***
|
