TABLE OF CONTENTS
ABINIT/mrggkk [ 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