1module m_ts_io 2 3 use precision, only : dp 4 use parallel, only : Node 5 6 implicit none 7 8 public :: ts_read_TSHS_opt 9 public :: ts_read_TSHS 10 public :: ts_write_TSHS 11 public :: fname_TSHS 12 public :: TSHS_version 13 public :: FC_index 14 15 private 16 17contains 18 19 subroutine ts_read_TSHS_opt(TSHS,DUMMY,na_u,no_u,no_s,nspin,n_nzs, & 20 xa, ucell, nsc, Qtot, Temp, Ef, & 21 Gamma,Gamma_TS,kscell,kdispl,onlyS,lasto, & 22 Bcast) 23 24 use m_os, only : file_exist 25#ifdef MPI 26 use mpi_siesta 27#endif 28 29! *********************** 30! * INPUT variables * 31! *********************** 32 character(len=*), intent(in) :: TSHS 33 34! *********************** 35! * OUTPUT variables * 36! *********************** 37 integer, optional :: DUMMY ! MUST NEVER BE PASSED 38 integer, intent(out), optional :: na_u, no_u, no_s, nspin, n_nzs, lasto(:), kscell(3,3), nsc(3) 39 real(dp), intent(out), optional :: xa(:,:), ucell(3,3), Qtot, Temp, Ef, kdispl(3) 40 logical, intent(out), optional :: Gamma, Gamma_TS, onlyS 41 logical, intent(in), optional :: Bcast 42 43! *********************** 44! * LOCAL variables * 45! *********************** 46 integer :: lna_u, lno_u, lno_s, lnspin, ln_nzs 47 integer :: uTSHS, version, tkscell(3,3), tnsc(3) 48 real(dp), allocatable :: txa(:,:) 49 real(dp) :: rtmp(3), tucell(3,3) 50 logical :: fGamma, bool(3), lonlyS 51#ifdef MPI 52 integer :: buffer_size, ipos 53 character(len=1), allocatable :: buffer(:) 54 integer :: MPIerror 55#endif 56 57 external :: io_assign, io_close 58 59 if ( present(DUMMY) ) call die('ts_read_TSHS_opt: Arguments has to be & 60 &named. Please correct sources.') 61 62 if ( .not. file_exist(TSHS, Bcast = Bcast) ) then 63 call die('ERROR: Could not read '//trim(TSHS)//'.') 64 end if 65 66#ifdef NCDF_4 67 ! If it is a netcdf file we read from that 68 ! instead 69 lna_u = len_trim(TSHS) 70 if ( TSHS(lna_u-1:lna_u) == 'nc' ) then 71 call ts_read_TSHS_opt_nc(TSHS,na_u=na_u,no_u=no_u,no_s=no_s, & 72 nspin=nspin,n_nzs=n_nzs,xa=xa, ucell=ucell, & 73 nsc=nsc, Qtot=Qtot, Temp=Temp, Ef=Ef, & 74 Gamma=Gamma,Gamma_TS=Gamma_TS, & 75 kscell=kscell,kdispl=kdispl,onlyS=onlyS,lasto=lasto, & 76 Bcast = Bcast) 77 return 78 end if 79#endif 80 81 if ( Node == 0 ) then 82 83 version = tshs_version(tshs) 84 85 select case ( version ) 86 case ( 0 , 1 ) 87 ! do nothing 88 case default 89 call die('Unsupported TSHS version file [0,1]') 90 end select 91 92 call io_assign(uTSHS) 93 open(file=trim(TSHS),unit=uTSHS,form='unformatted') 94 if ( version /= 0 ) then 95 read(uTSHS) version 96 end if 97 read(uTSHS) lna_u,lno_u,lno_s,lnspin,ln_nzs 98 allocate(txa(3,lna_u)) 99 if ( present(na_u) ) na_u = lna_u 100 if ( present(no_u) ) no_u = lno_u 101 if ( present(no_s) ) no_s = lno_s 102 if ( present(nspin) ) nspin = lnspin 103 if ( present(n_nzs) ) n_nzs = ln_nzs 104 if ( version == 0 ) then 105 read(uTSHS) txa 106 read(uTSHS) ! iza 107 read(uTSHS) tucell 108 tnsc = 0 109 else if ( version == 1 ) then 110 read(uTSHS) tnsc ! corrected nsc 111 read(uTSHS) tucell, txa 112 end if 113 if ( present(nsc) ) nsc = tnsc 114 if ( present(ucell) ) ucell = tucell 115 if ( present(xa) ) xa = txa 116 if ( version == 0 ) then 117 read(uTSHS) fGamma ! SIESTA_Gamma 118 if ( present(Gamma) ) Gamma = fGamma 119 read(uTSHS) lonlyS 120 if ( present(OnlyS) ) onlys = lonlyS 121 122 if ( present(Gamma_TS) ) then 123 read(uTSHS) Gamma_TS 124 else 125 read(uTSHS) ! Gamma_TS 126 end if 127 if ( present(kscell) ) then 128 read(uTSHS) kscell 129 else 130 read(uTSHS) ! ts_kscell_file 131 end if 132 if ( present(kdispl) ) then 133 read(uTSHS) kdispl 134 else 135 read(uTSHS) ! ts_kdispl_file 136 end if 137 else if ( version == 1 ) then 138 read(uTSHS) fGamma, bool(1), lonlyS 139 if ( present(Gamma) ) Gamma = fGamma 140 if ( present(Gamma_TS) ) Gamma_TS = bool(1) 141 if ( present(OnlyS) ) OnlyS = lonlyS 142 read(uTSHS) tkscell, rtmp 143 if ( present(kscell) ) kscell = tkscell 144 if ( present(kdispl) ) kdispl = rtmp 145 read(uTSHS) rtmp(1:3) 146 if ( present(Ef) ) Ef = rtmp(1) 147 if ( present(Qtot) ) Qtot = rtmp(2) 148 if ( present(Temp) ) Temp = rtmp(3) 149 end if 150 151 ! These quantities are not needed 152 read(uTSHS) ! istep, ia1 153 154 if ( present(lasto) ) then 155 if ( size(lasto) /= lna_u+1 ) call die('ts_read_TSHS: Wrong size of lasto') 156 read(uTSHS) lasto 157 else 158 read(uTSHS) ! lasto 159 end if 160 161 if ( version == 0 .and. .not. fGamma ) then 162 read(uTSHS) ! indxuo 163 end if 164 165 read(uTSHS) ! numh 166 167 if ( version == 0 ) then 168 read(uTSHS) rtmp(1:2) 169 if ( present(Qtot) ) Qtot = rtmp(1) 170 if ( present(Temp) ) Temp = rtmp(2) 171 172 if ( present(Ef) ) then 173 read(uTSHS) Ef 174 else 175 read(uTSHS) ! Ef 176 end if 177 end if 178 179 deallocate(txa) 180 call io_close(uTSHS) 181 182 end if 183 184#ifdef MPI 185 if ( present(Bcast) ) then 186 ! if we do not request broadcasting, then return... 187 if ( .not. Bcast ) return 188 end if 189 190 ! Broadcast na_u (for easy reference) 191 call MPI_Bcast(lna_u,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 192 193#ifdef MPI_OLD 194 if ( present(na_u) ) & 195 call MPI_Bcast(na_u,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 196 if ( present(no_u) ) & 197 call MPI_Bcast(no_u,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 198 if ( present(no_s) ) & 199 call MPI_Bcast(no_s,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 200 if ( present(nspin) ) & 201 call MPI_Bcast(nspin,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 202 if ( present(n_nzs) ) & 203 call MPI_Bcast(n_nzs,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 204 if ( present(nsc) ) & 205 call MPI_Bcast(nsc,3,MPI_Integer,0,MPI_Comm_World,MPIerror) 206 if ( present(xa) ) & 207 call MPI_Bcast(xa(1,1),3*lna_u,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 208 if ( present(ucell) ) & 209 call MPI_Bcast(ucell(1,1),9,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 210 if ( present(Gamma) ) & 211 call MPI_Bcast(Gamma,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 212 if ( present(OnlyS) ) & 213 call MPI_Bcast(OnlyS,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 214 if ( present(Gamma_TS) ) & 215 call MPI_Bcast(Gamma_TS,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 216 if ( present(kscell) ) & 217 call MPI_Bcast(kscell,9,MPI_Integer,0,MPI_Comm_World,MPIerror) 218 if ( present(kdispl) ) & 219 call MPI_Bcast(kdispl,3,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 220 if ( present(lasto) ) & 221 call MPI_Bcast(lasto(1),lna_u+1,MPI_Integer,0,MPI_Comm_World,MPIerror) 222 if ( present(Qtot) ) & 223 call MPI_Bcast(Qtot,1,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 224 if ( present(Temp) ) & 225 call MPI_Bcast(Temp,1,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 226 if ( present(Ef) ) & 227 call MPI_Bcast(Ef,1,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 228#else 229 230 ! this should be more than enough... 231 buffer_size = 8 * (lna_u * 6 + 103) 232 allocate(buffer(buffer_size)) 233 ! position of data in buffer... 234 ipos = 0 235 236 if ( Node == 0 ) then 237 if ( present(na_u) ) & ! 4 238 call MPI_Pack(na_u,1,MPI_Integer, & 239 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 240 if ( present(no_u) ) & ! 4 241 call MPI_Pack(no_u,1,MPI_Integer, & 242 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 243 if ( present(no_s) ) & ! 4 244 call MPI_Pack(no_s,1,MPI_Integer, & 245 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 246 if ( present(nspin) ) & ! 4 247 call MPI_Pack(nspin,1,MPI_Integer, & 248 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 249 if ( present(n_nzs) ) & ! 4 250 call MPI_Pack(n_nzs,1,MPI_Integer, & 251 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 252 if ( present(xa) ) & ! 8 * 3 * na_u 253 call MPI_Pack(xa(1,1),3*lna_u,MPI_Double_Precision, & 254 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 255 if ( present(ucell) ) & ! 8 * 3 * 3 256 call MPI_Pack(ucell(1,1),9,MPI_Double_Precision, & 257 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 258 if ( present(nsc) ) & ! 4 259 call MPI_Pack(nsc(1),3,MPI_Integer, & 260 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 261 if ( present(Gamma) ) & ! 4 262 call MPI_Pack(Gamma,1,MPI_Logical, & 263 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 264 if ( present(kscell) ) &! 4 * 3 * 3 265 call MPI_Pack(kscell(1,1),9,MPI_Integer, & 266 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 267 if ( present(kdispl) ) &! 8 * 3 268 call MPI_Pack(kdispl(1),3,MPI_Double_Precision, & 269 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 270 if ( present(OnlyS) ) & ! 4 271 call MPI_Pack(OnlyS,1,MPI_Logical, & 272 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 273 if ( present(Gamma_TS) ) &! 4 274 call MPI_Pack(Gamma_TS,1,MPI_Logical, & 275 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 276 if ( present(lasto) ) & ! 4 * (na_u+1) 277 call MPI_Pack(lasto(1),lna_u+1,MPI_Logical, & 278 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 279 if ( present(Qtot) ) & ! 8 280 call MPI_Pack(Qtot,1,MPI_Double_Precision, & 281 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 282 if ( present(Temp) ) & ! 8 283 call MPI_Pack(Temp,1,MPI_Double_Precision, & 284 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 285 if ( present(Ef) ) & ! 8 286 call MPI_Pack(Ef,1,MPI_Double_Precision, & 287 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 288 289 if ( ipos >= buffer_size .or. ipos < 0 .or. MPIerror /= MPI_Success ) then 290 call die('Error in estimating the buffer-size for the & 291 &TSHS reading. Please contact the developers') 292 end if 293 294 end if 295 296 call MPI_Bcast(buffer,buffer_size,MPI_Packed, & 297 0, MPI_Comm_World, MPIerror) 298 299 if ( Node /= 0 ) then 300 if ( present(na_u) ) & 301 call MPI_UnPack(buffer,buffer_size,ipos, & 302 na_u,1,MPI_Integer, & 303 MPI_Comm_World, MPIerror) 304 if ( present(no_u) ) & 305 call MPI_UnPack(buffer,buffer_size,ipos, & 306 no_u,1,MPI_Integer, & 307 MPI_Comm_World, MPIerror) 308 if ( present(no_s) ) & 309 call MPI_UnPack(buffer,buffer_size,ipos, & 310 no_s,1,MPI_Integer, & 311 MPI_Comm_World, MPIerror) 312 if ( present(nspin) ) & 313 call MPI_UnPack(buffer,buffer_size,ipos, & 314 nspin,1,MPI_Integer, & 315 MPI_Comm_World, MPIerror) 316 if ( present(n_nzs) ) & 317 call MPI_UnPack(buffer,buffer_size,ipos, & 318 n_nzs,1,MPI_Integer, & 319 MPI_Comm_World, MPIerror) 320 if ( present(xa) ) & 321 call MPI_UnPack(buffer,buffer_size,ipos, & 322 xa(1,1),3*lna_u,MPI_Double_Precision, & 323 MPI_Comm_World, MPIerror) 324 if ( present(ucell) ) & 325 call MPI_UnPack(buffer,buffer_size,ipos, & 326 ucell(1,1),9,MPI_Double_Precision, & 327 MPI_Comm_World, MPIerror) 328 if ( present(nsc) ) & 329 call MPI_UnPack(buffer,buffer_size,ipos, & 330 nsc(1),3,MPI_Integer, & 331 MPI_Comm_World, MPIerror) 332 if ( present(Gamma) ) & 333 call MPI_UnPack(buffer,buffer_size,ipos, & 334 Gamma,1,MPI_Logical, & 335 MPI_Comm_World, MPIerror) 336 if ( present(kscell) ) & 337 call MPI_UnPack(buffer,buffer_size,ipos, & 338 kscell(1,1),9,MPI_Integer, & 339 MPI_Comm_World, MPIerror) 340 if ( present(kdispl) ) & 341 call MPI_UnPack(buffer,buffer_size,ipos, & 342 kdispl(1),3,MPI_Double_Precision, & 343 MPI_Comm_World, MPIerror) 344 if ( present(OnlyS) ) & 345 call MPI_UnPack(buffer,buffer_size,ipos, & 346 OnlyS,1,MPI_Logical, & 347 MPI_Comm_World, MPIerror) 348 if ( present(Gamma_TS) ) & 349 call MPI_UnPack(buffer,buffer_size,ipos, & 350 Gamma_TS,1,MPI_Logical, & 351 MPI_Comm_World, MPIerror) 352 if ( present(lasto) ) & 353 call MPI_UnPack(buffer,buffer_size,ipos, & 354 lasto(1),lna_u+1,MPI_Logical, & 355 MPI_Comm_World, MPIerror) 356 if ( present(Qtot) ) & 357 call MPI_UnPack(buffer,buffer_size,ipos, & 358 Qtot,1,MPI_Double_Precision, & 359 MPI_Comm_World, MPIerror) 360 if ( present(Temp) ) & 361 call MPI_UnPack(buffer,buffer_size,ipos, & 362 Temp,1,MPI_Double_Precision, & 363 MPI_Comm_World, MPIerror) 364 if ( present(Ef) ) & 365 call MPI_UnPack(buffer,buffer_size,ipos, & 366 Ef,1,MPI_Double_Precision, & 367 MPI_Comm_World, MPIerror) 368 end if 369 370 deallocate(buffer) 371 372#endif 373 374#endif 375 376 end subroutine ts_read_TSHS_opt 377 378#ifdef NCDF_4 379 380 ! This is a routine for reading information from a 381 ! siesta NetCDF-4 format file 382 subroutine ts_read_TSHS_opt_nc(TSHS,na_u,no_u,no_s,nspin,n_nzs, & 383 xa, ucell, nsc, Qtot, Temp, Ef, & 384 Gamma,Gamma_TS,kscell,kdispl,onlyS,lasto, & 385 Bcast) 386 387 use netcdf_ncdf, ncdf_parallel => parallel 388#ifdef MPI 389 use mpi_siesta 390#endif 391! *********************** 392! * INPUT variables * 393! *********************** 394 character(len=*), intent(in) :: TSHS 395 396! *********************** 397! * OUTPUT variables * 398! *********************** 399 integer, intent(out), optional :: na_u, no_u, no_s, nspin, n_nzs, lasto(:), kscell(3,3), nsc(3) 400 real(dp), intent(out), optional :: xa(:,:), ucell(3,3), Qtot, Temp, Ef, kdispl(3) 401 logical, intent(out), optional :: Gamma, Gamma_TS, onlyS 402 logical, intent(in), optional :: Bcast 403 404! *********************** 405! * LOCAL variables * 406! *********************** 407 type(hNCDF) :: ncdf, grp 408 integer :: lna_u, tnsc(3) 409#ifdef MPI 410 integer :: buffer_size, ipos 411 character(len=1), allocatable :: buffer(:) 412 integer :: MPIerror 413#endif 414 415 if ( present(onlyS) ) onlyS = .false. 416 417 call ncdf_open(ncdf,trim(TSHS),mode=NF90_NOWRITE) 418 419 ! Now we read in the things 420 call ncdf_inq_dim(ncdf,'na_u',len=lna_u) 421 if ( present(na_u) ) na_u = lna_u 422 if ( present(no_u) ) & 423 call ncdf_inq_dim(ncdf,'no_u',len=no_u) 424 if ( present(no_s) ) & 425 call ncdf_inq_dim(ncdf,'no_s',len=no_s) 426 if ( present(nspin) ) & 427 call ncdf_inq_dim(ncdf,'spin',len=nspin) 428 if ( present(xa) ) & 429 call ncdf_get_var(ncdf,'xa',xa) 430 if ( present(ucell) ) & 431 call ncdf_get_var(ncdf,'cell',ucell) 432 call ncdf_get_var(ncdf,'nsc',tnsc) 433 if ( present(nsc) ) nsc = tnsc 434 if ( present(Gamma) ) then 435 Gamma = sum(tnsc) == 1 436 if ( present(Gamma_TS) ) Gamma_TS = Gamma 437 end if 438 if ( present(Ef) ) & 439 call ncdf_get_var(ncdf,'Ef',Ef) 440 if ( present(Qtot) ) & 441 call ncdf_get_var(ncdf,'Qtot',Qtot) 442 443 if ( Node == 0 ) then 444 if ( present(lasto) ) then 445 if ( size(lasto) /= lna_u+1 ) call die('ts_read_TSHS: Wrong size of lasto') 446 lasto(1) = 0 447 call ncdf_get_var(ncdf,'lasto',lasto(2:lna_u+1)) 448 end if 449 end if 450 451 if ( present(n_nzs) ) then 452 call ncdf_open_grp(ncdf,'SPARSE',grp) 453 call ncdf_get_var(grp,'nnzs',n_nzs) 454 end if 455 456 call ncdf_open_grp(ncdf,'SETTINGS',grp) 457 if ( present(Temp) ) & 458 call ncdf_get_var(grp,'ElectronicTemperature',Temp) 459 if ( present(kscell) ) & 460 call ncdf_get_var(grp,'BZ',kscell) 461 if ( present(kdispl) ) & 462 call ncdf_get_var(grp,'BZ_displ',kdispl) 463 464 call ncdf_close(ncdf) 465 466#ifdef MPI 467 if ( present(Bcast) ) then 468 ! if we do not request broadcasting, then return... 469 if ( .not. Bcast ) return 470 end if 471 472 ! Broadcast na_u (for easy reference) 473 call MPI_Bcast(lna_u,1,MPI_Integer,0,MPI_Comm_World,MPIerror) 474 475 ! this should be more than enough... 476 buffer_size = 8 * (lna_u * 6 + 103) 477 allocate(buffer(buffer_size)) 478 ! position of data in buffer... 479 ipos = 0 480 481 if ( Node == 0 ) then 482 if ( present(na_u) ) & ! 4 483 call MPI_Pack(na_u,1,MPI_Integer, & 484 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 485 if ( present(no_u) ) & ! 4 486 call MPI_Pack(no_u,1,MPI_Integer, & 487 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 488 if ( present(no_s) ) & ! 4 489 call MPI_Pack(no_s,1,MPI_Integer, & 490 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 491 if ( present(nspin) ) & ! 4 492 call MPI_Pack(nspin,1,MPI_Integer, & 493 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 494 if ( present(n_nzs) ) & ! 4 495 call MPI_Pack(n_nzs,1,MPI_Integer, & 496 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 497 if ( present(xa) ) & ! 8 * 3 * na_u 498 call MPI_Pack(xa(1,1),3*lna_u,MPI_Double_Precision, & 499 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 500 if ( present(ucell) ) & ! 8 * 3 * 3 501 call MPI_Pack(ucell(1,1),9,MPI_Double_Precision, & 502 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 503 if ( present(nsc) ) & ! 4 504 call MPI_Pack(nsc(1),3,MPI_Integer, & 505 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 506 if ( present(Gamma) ) & ! 4 507 call MPI_Pack(Gamma,1,MPI_Logical, & 508 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 509 if ( present(kscell) ) &! 4 * 3 * 3 510 call MPI_Pack(kscell(1,1),9,MPI_Integer, & 511 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 512 if ( present(kdispl) ) &! 8 * 3 513 call MPI_Pack(kdispl(1),3,MPI_Double_Precision, & 514 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 515 if ( present(Gamma_TS) ) &! 4 516 call MPI_Pack(Gamma_TS,1,MPI_Logical, & 517 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 518 if ( present(lasto) ) & ! 4 * (na_u+1) 519 call MPI_Pack(lasto(1),lna_u+1,MPI_Logical, & 520 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 521 if ( present(Qtot) ) & ! 8 522 call MPI_Pack(Qtot,1,MPI_Double_Precision, & 523 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 524 if ( present(Temp) ) & ! 8 525 call MPI_Pack(Temp,1,MPI_Double_Precision, & 526 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 527 if ( present(Ef) ) & ! 8 528 call MPI_Pack(Ef,1,MPI_Double_Precision, & 529 buffer,buffer_size, ipos, MPI_Comm_World, MPIerror) 530 531 if ( ipos >= buffer_size .or. ipos < 0 .or. MPIerror /= MPI_Success ) then 532 call die('Error in estimating the buffer-size for the & 533 &TSHS reading. Please contact the developers') 534 end if 535 536 end if 537 538 call MPI_Bcast(buffer,buffer_size,MPI_Packed, & 539 0, MPI_Comm_World, MPIerror) 540 541 if ( Node /= 0 ) then 542 if ( present(na_u) ) & 543 call MPI_UnPack(buffer,buffer_size,ipos, & 544 na_u,1,MPI_Integer, & 545 MPI_Comm_World, MPIerror) 546 if ( present(no_u) ) & 547 call MPI_UnPack(buffer,buffer_size,ipos, & 548 no_u,1,MPI_Integer, & 549 MPI_Comm_World, MPIerror) 550 if ( present(no_s) ) & 551 call MPI_UnPack(buffer,buffer_size,ipos, & 552 no_s,1,MPI_Integer, & 553 MPI_Comm_World, MPIerror) 554 if ( present(nspin) ) & 555 call MPI_UnPack(buffer,buffer_size,ipos, & 556 nspin,1,MPI_Integer, & 557 MPI_Comm_World, MPIerror) 558 if ( present(n_nzs) ) & 559 call MPI_UnPack(buffer,buffer_size,ipos, & 560 n_nzs,1,MPI_Integer, & 561 MPI_Comm_World, MPIerror) 562 if ( present(xa) ) & 563 call MPI_UnPack(buffer,buffer_size,ipos, & 564 xa(1,1),3*lna_u,MPI_Double_Precision, & 565 MPI_Comm_World, MPIerror) 566 if ( present(ucell) ) & 567 call MPI_UnPack(buffer,buffer_size,ipos, & 568 ucell(1,1),9,MPI_Double_Precision, & 569 MPI_Comm_World, MPIerror) 570 if ( present(nsc) ) & 571 call MPI_UnPack(buffer,buffer_size,ipos, & 572 nsc(1),3,MPI_Integer, & 573 MPI_Comm_World, MPIerror) 574 if ( present(Gamma) ) & 575 call MPI_UnPack(buffer,buffer_size,ipos, & 576 Gamma,1,MPI_Logical, & 577 MPI_Comm_World, MPIerror) 578 if ( present(kscell) ) & 579 call MPI_UnPack(buffer,buffer_size,ipos, & 580 kscell(1,1),9,MPI_Integer, & 581 MPI_Comm_World, MPIerror) 582 if ( present(kdispl) ) & 583 call MPI_UnPack(buffer,buffer_size,ipos, & 584 kdispl(1),3,MPI_Double_Precision, & 585 MPI_Comm_World, MPIerror) 586 if ( present(Gamma_TS) ) & 587 call MPI_UnPack(buffer,buffer_size,ipos, & 588 Gamma_TS,1,MPI_Logical, & 589 MPI_Comm_World, MPIerror) 590 if ( present(lasto) ) & 591 call MPI_UnPack(buffer,buffer_size,ipos, & 592 lasto(1),lna_u+1,MPI_Logical, & 593 MPI_Comm_World, MPIerror) 594 if ( present(Qtot) ) & 595 call MPI_UnPack(buffer,buffer_size,ipos, & 596 Qtot,1,MPI_Double_Precision, & 597 MPI_Comm_World, MPIerror) 598 if ( present(Temp) ) & 599 call MPI_UnPack(buffer,buffer_size,ipos, & 600 Temp,1,MPI_Double_Precision, & 601 MPI_Comm_World, MPIerror) 602 if ( present(Ef) ) & 603 call MPI_UnPack(buffer,buffer_size,ipos, & 604 Ef,1,MPI_Double_Precision, & 605 MPI_Comm_World, MPIerror) 606 end if 607 608 deallocate(buffer) 609 610#endif 611 612 end subroutine ts_read_TSHS_opt_nc 613#endif 614 615 subroutine ts_read_TSHS(filename, & 616 onlyS, Gamma, TSGamma, & 617 ucell, nsc, na_u, no_u, nspin, & 618 kscell, kdispl, & 619 xa, lasto, & 620 sp, H, S, isc_off, & 621 Ef, Qtot, Temp, & 622 istep, ia1, tag, & 623 Bcast) 624 625! ********************************************************************* 626! Saves the hamiltonian and overlap matrices, and other data required 627! to obtain the bands and density of states 628! Writen by J.Soler July 1997. 629! Note because of the new more compact method of storing H and S 630! this routine is NOT backwards compatible 631! Modified by M.Paulsson 2009 to: 632! 1: To include information of which FC step for phonon calculations 633! 2: To only save the overlap matrix if onlyS flag is set 634! (Used for e-ph coupling calculations) 635! 3: File format changed to unify Copenhagen/Barcelona Transiesta vers. 636! 4: Smaller files by writing arrays directly instead of element wise 637! *************************** INPUT ********************************** 638! logical Gamma : Is only gamma point used? 639! logical TSGamma : Is only TS gamma point used? 640! logical onlyS : Should only overlap matrix be saved? 641! ******************** INPUT or OUTPUT (depending on task) *********** 642! integer no_u : Number of basis orbitals per unit cell 643! integer no_s : Number of basis orbitals per supercell 644! integer Enspin : Spin polarization (1 or 2) 645! integer n_nzs : First dimension of listh, H, S and 646! integer numh(nuo) : Number of nonzero elements of each row 647! of hamiltonian matrix 648! integer listhptr(nuo) : Pointer to the start of each row (-1) 649! of hamiltonian matrix 650! integer listh(n_nzs) : Nonzero hamiltonian-matrix element column 651! indexes for each matrix row 652! real*8 H(n_nzs,Enspin) : Hamiltonian in sparse form 653! real*8 S(n_nzs) : Overlap in sparse form 654! real*8 qtot : Total number of electrons 655! real*8 temp : Electronic temperature for Fermi smearing 656! real*8 xij(3,n_nzs),isc_off : Vectors between orbital centers (sparse) 657! (not read/written if only gamma point) 658! TSS Begin 659! ********************* ADDED ARGUMENTS FOR TRANSIESTA **************** 660! integer fnlength : file name length 661! character(fnlength) fname : file nema for input or output 662! integer na_u : Number of atoms per unit cell 663! integer istep, ia1 : Force constant step and atom number 664! logical check_kcell : Do a check of the kscell and kdispl 665! TSS End 666! *************************** UNITS *********************************** 667! Units should be consistent between task='read' and 'write' 668! ********************************************************************* 669 670 use sys, only : die 671#ifdef MPI 672 use mpi_siesta 673#endif 674 use alloc, only : re_alloc 675 use geom_helper, only : iaorb, ucorb 676 use m_io_s 677 use m_os, only : file_exist 678 use class_Sparsity 679 use class_OrbitalDistribution 680 use class_dSpData1D 681 use class_dSpData2D 682 use m_sparse, only : list_col_correct, xij_offset, calc_nsc 683 684 implicit none 685 686! ********************** 687! * INPUT variables * 688! ********************** 689 character(len=*), intent(in) :: filename 690 logical, intent(out) :: onlyS, Gamma, TSGamma 691 real(dp), intent(out) :: ucell(3,3) 692 integer, intent(out) :: nsc(3), na_u, no_u, nspin 693 integer, intent(out) :: kscell(3,3) 694 real(dp), intent(out) :: kdispl(3) 695 real(dp), pointer :: xa(:,:) 696 integer, pointer :: lasto(:) ! (0:na_u) 697 type(Sparsity), intent(inout) :: sp 698 type(dSpData2D), intent(inout) :: H 699 type(dSpData1D), intent(inout) :: S 700 integer, pointer :: isc_off(:,:) 701 real(dp), intent(out) :: Ef, Qtot,Temp 702 ! These have to be set before entrance (makes it possible to read 703 ! in FCrun TSHS files...) 704 integer, intent(out) :: istep, ia1 705 character(len=*), intent(in), optional :: tag 706 ! If true it will broadcast every information within the code... 707 logical, intent(in), optional :: Bcast 708 709! ************************ 710! * LOCAL variables * 711! ************************ 712 type(OrbitalDistribution), pointer :: dit 713 type(dSpData2D) :: xij 714 integer :: iu, version, no_s, n_s, n_nzs 715 integer :: i, j, ind 716 integer, pointer :: ncol(:), l_ptr(:), l_col(:) 717 integer, allocatable :: indxuo(:) 718 real(dp), pointer :: lxij(:,:) 719 character(len=250) :: ltag 720 logical :: lBcast, exist 721#ifdef MPI 722 integer :: all_I(0:9) 723 integer :: MPIerror 724#endif 725 726 external :: io_assign, io_close 727 728 if ( .not. file_exist(filename, Bcast = Bcast ) ) then 729 call die('ERROR: Could not read '//trim(filename)//'.') 730 end if 731 732#ifdef NCDF_4 733 ! If it is a NetCDF file, we call the netCDF 734 ! routine 735 i = len_trim(filename) 736 if ( filename(i-1:i) == 'nc' ) then 737 call ts_read_TSHS_nc(filename, & 738 Gamma, TSGamma, & 739 ucell, nsc, na_u, no_u, nspin, & 740 kscell, kdispl, & 741 xa, lasto, & 742 sp, H, S, isc_off, & 743 Ef, Qtot, Temp, & 744 tag = tag, Bcast = Bcast) 745 746 OnlyS = .false. 747 Temp = 0._dp 748 istep = 0 749 ia1 = 0 750 751 return 752 753 end if 754#endif 755 756#ifdef TRANSIESTA_DEBUG 757 call write_debug( 'PRE ts_io_read' ) 758#endif 759 760 nullify(xa,lasto,ncol,l_ptr,l_col,isc_off) 761 762 ltag = trim(filename) 763 if ( present(tag) ) ltag = trim(tag) 764 765 ! Determine whether to broadcast afterwards 766 lBcast = .false. 767 if ( present(Bcast) ) lBcast = Bcast 768 769 ! Many variables are needed initialization due 770 ! to compiletime unawareness 771 nsc(:) = 0 772 773 exist = file_exist(filename, Bcast = Bcast ) 774 775 if ( Node == 0 ) then 776 777 if ( .not. exist ) then 778 call die('ERROR: Could not read '//trim(filename)//'.') 779 end if 780 781 ! Get file version 782 version = tshs_version(filename) 783 784 select case ( version ) 785 case ( 0 , 1 ) 786 ! do nothing 787 case default 788 call die('Unsupported TSHS version file [0,1]') 789 end select 790 791 ! Open file 792 call io_assign( iu ) 793 open( iu, file=filename, form='unformatted', status='old' ) 794 795 ! Read Dimensions Information 796 if ( version /= 0 ) then 797 read(iu) ! version 798 end if 799 read(iu) na_u, no_u, no_s, nspin, n_nzs 800 801 ! Read Geometry information 802 allocate(xa(3,na_u)) 803 call memory('A','D',3*na_u,'iohs') 804 call memory('A','I',na_u,'iohs') 805 if ( version == 0 ) then 806 read(iu) xa 807 read(iu) ! iza 808 read(iu) ucell 809 else if ( version == 1 ) then 810 read(iu) nsc ! original nsc, and the corrected nsc 811 read(iu) ucell, xa 812 end if 813 814 ! Read k-point sampling information 815 if ( version == 0 ) then 816 read(iu) Gamma 817 read(iu) onlyS 818 read(iu) TSGamma 819 read(iu) kscell 820 read(iu) kdispl 821 else if ( version == 1 ) then 822 read(iu) Gamma, TSGamma, onlyS 823 read(iu) kscell, kdispl 824 read(iu) Ef, Qtot, Temp 825 end if 826 read(iu) istep, ia1 827 828 allocate(lasto(0:na_u)) 829 call memory('A','I',1+na_u,'iohs') 830 read(iu) lasto 831 832 if ( version == 0 ) then 833 ! We check that the format of the supercell indices 834 ! are consistent with the general formalism 835 allocate(indxuo(no_s)) 836 read(iu) indxuo 837 do i = 1 , no_s 838 if ( indxuo(i) /= ucorb(i,no_u) ) & 839 call die('Error in indxuo, not recognized') 840 end do 841 deallocate(indxuo) 842 end if 843 844 end if 845 846#ifdef MPI 847 if ( lBcast ) then 848 ! Bcast initial sizes 849 if ( Node == 0 ) then 850 all_I(0:9) = (/version,na_u,no_u,nspin,n_nzs,nsc(1),nsc(2),nsc(3),istep,ia1/) 851 end if 852 call MPI_Bcast(all_I(0),10,MPI_Integer,0,MPI_Comm_World,MPIerror) 853 version = all_I(0) 854 na_u = all_I(1) 855 no_u = all_I(2) 856 nspin = all_I(3) 857 n_nzs = all_I(4) 858 nsc(1) = all_I(5) 859 nsc(2) = all_I(6) 860 nsc(3) = all_I(7) 861 istep = all_I(8) 862 ia1 = all_I(9) 863 call MPI_Bcast(ucell(1,1),9,MPI_Double_Precision,0, & 864 MPI_Comm_World,MPIerror) 865 call MPI_Bcast(Gamma,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 866 call MPI_Bcast(TSGamma,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 867 call MPI_Bcast(onlyS,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 868 call MPI_Bcast(kscell(1,1),9,MPI_Integer,0,MPI_Comm_World,MPIerror) 869 call MPI_Bcast(kdispl(1),3,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 870 871 if ( Node /= 0 ) then 872 allocate(xa(3,na_u)) 873 call memory('A','D',3*na_u,'iohs') 874 allocate(lasto(0:na_u)) 875 call memory('A','I',1+na_u,'iohs') 876 end if 877 call MPI_Bcast(xa(1,1),3*na_u,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 878 call MPI_Bcast(lasto(0),1+na_u,MPI_Integer,0, MPI_Comm_World,MPIerror) 879 880 end if 881#endif 882 883 if ( version == 0 ) then 884 885 nullify(ncol,l_col) 886 allocate(ncol(no_u),l_col(n_nzs)) 887 888 if ( Node == 0 ) then 889 890 read(iu) ncol 891 892 ! Read Electronic Structure Information 893 read(iu) Qtot,Temp 894 read(iu) Ef 895 896 j = 0 897 do i = 1 , no_u 898 read(iu) l_col(j+1:j+ncol(i)) 899 j = j + ncol(i) 900 end do 901 902 end if 903 904#ifdef MPI 905 if ( lBcast ) then 906 call MPI_Bcast(ncol,no_u,MPI_Integer,0, MPI_Comm_World,MPIerror) 907 call MPI_Bcast(l_col,n_nzs,MPI_Integer,0, MPI_Comm_World,MPIerror) 908 end if 909#endif 910 911 nullify(l_ptr) 912 allocate(l_ptr(no_u)) 913 914 l_ptr(1) = 0 915 do i = 2 , no_u 916 l_ptr(i) = l_ptr(i-1) + ncol(i-1) 917 end do 918 919 ! Create the sparsity pattern 920 call newSparsity(sp,no_u,no_u, & 921 n_nzs, ncol, l_ptr, l_col, trim(ltag)) 922 923 deallocate(ncol,l_ptr,l_col) 924 nullify(ncol,l_ptr,l_col) 925 926 else if ( version == 1 ) then 927 call io_read_Sp(iu,no_u,sp,tag=trim(ltag),Bcast=Bcast) 928 end if 929 930 ! Read in S 931 call io_read_d1D(iu,sp,S,trim(ltag)//': S',Bcast=Bcast) 932 933 ! Utilize the same distribution for the others 934 dit => dist(S) 935 936 if ( .not. onlyS ) then 937 ! Read in H 938 call io_read_d2D(iu,sp,H,nspin,trim(ltag)//': H',Bcast=Bcast,dit=dit) 939 end if 940 941 if ( .not. Gamma ) then 942 943 if ( version == 0 ) then 944 call io_read_d2D(iu,sp,xij,3,'xij',sparsity_dim=2, & 945 dit=dit, Bcast=Bcast) 946 947 call attach(sp,n_col=ncol) 948 lxij => val(xij) 949 950 ! A stupid "bug" in the very old TRANSIESTA code 951 ! was transposing the xij array before writing... sigh... 952 ind = 0 953 do i = 1 , no_u 954 lxij(:,ind+1:ind+ncol(i)) = transpose( & 955 reshape(lxij(:,ind+1:ind+ncol(i)), (/ncol(i),3/))) 956 ind = ind + ncol(i) 957 end do 958 959 ! We do not need (MUST NOT) do bcast the routines: 960 ! calc_nsc, list_col_correct, xij_offset 961 ! They are duplicated on all nodes, hence b-casting 962 ! will introduce wrong columns... 963 call calc_nsc(ucell,na_u,xa,lasto,xij,nsc) 964 965 ! Ensure that list_col is correctly formatted (not always 966 ! needed, but for consistency) 967 call list_col_correct(ucell,nsc,na_u,xa,lasto,xij) 968 969 call xij_offset(ucell,nsc,na_u,xa,lasto,xij,isc_off) 970 971 ! We do not need the xij array anymore... :) 972 call delete(xij) 973 974 else if ( version == 1 ) then 975 976 ! Number of supercells 977 n_s = product(nsc) 978 call re_alloc(isc_off,1,3,1,n_s) 979 980 if ( Node == 0 ) then 981 read(iu) isc_off 982 end if 983#ifdef MPI 984 if ( lBcast ) then 985 call MPI_Bcast(isc_off(1,1),3*n_s, MPI_Integer, 0, & 986 MPI_Comm_World,MPIerror) 987 end if 988#endif 989 990 end if 991 else 992 993 ! Number of supercells 994 n_s = product(nsc) 995 call re_alloc(isc_off,1,3,1,n_s) 996 isc_off(:,:) = 0 997 998 end if 999 1000 if ( Node == 0 ) then 1001 ! Close file 1002 call io_close( iu ) 1003 end if 1004 1005#ifdef MPI 1006 if ( lBcast ) then 1007 call MPI_Bcast(Qtot,1,MPI_Double_Precision,0, MPI_Comm_World,MPIerror) 1008 call MPI_Bcast(Temp,1,MPI_Double_Precision,0, MPI_Comm_World,MPIerror) 1009 call MPI_Bcast(Ef,1,MPI_Double_Precision,0, MPI_Comm_World,MPIerror) 1010 end if 1011#endif 1012 1013#ifdef TRANSIESTA_DEBUG 1014 call write_debug( 'POS ts_io_read' ) 1015#endif 1016 1017 end subroutine ts_read_TSHS 1018 1019#ifdef NCDF_4 1020 subroutine ts_read_TSHS_nc(filename, & 1021 Gamma, TSGamma, & 1022 ucell, nsc, na_u, no_u, nspin, & 1023 kscell, kdispl, & 1024 xa, lasto, & 1025 sp, H, S, isc_off, & 1026 Ef, Qtot, Temp, & 1027 tag, Bcast) 1028 1029 use m_ncdf_io, only : cdf_r_Sp, cdf_r_d1D, cdf_r_d2D 1030 1031 use netcdf_ncdf, ncdf_parallel => parallel 1032 use sys, only : die 1033#ifdef MPI 1034 use mpi_siesta 1035#endif 1036 use alloc, only : re_alloc 1037 use class_Sparsity 1038 use class_OrbitalDistribution 1039 use class_dSpData1D 1040 use class_dSpData2D 1041 1042 implicit none 1043 1044! ********************** 1045! * INPUT variables * 1046! ********************** 1047 character(len=*), intent(in) :: filename 1048 logical, intent(out) :: Gamma, TSGamma 1049 real(dp), intent(out) :: ucell(3,3) 1050 integer, intent(out) :: nsc(3), na_u, no_u, nspin 1051 integer, intent(out) :: kscell(3,3) 1052 real(dp), intent(out) :: kdispl(3) 1053 real(dp), pointer :: xa(:,:) 1054 integer, pointer :: lasto(:) ! (0:na_u) 1055 type(Sparsity), intent(inout) :: sp 1056 type(dSpData2D), intent(inout) :: H 1057 type(dSpData1D), intent(inout) :: S 1058 integer, pointer :: isc_off(:,:) 1059 real(dp), intent(out) :: Ef, Qtot,Temp 1060 character(len=*), intent(in), optional :: tag 1061 ! If true it will broadcast every information within the code... 1062 logical, intent(in), optional :: Bcast 1063 1064! ************************ 1065! * LOCAL variables * 1066! ************************ 1067 type(hNCDF) :: ncdf, grp 1068 type(OrbitalDistribution), pointer :: dit 1069 integer :: n_s 1070 character(len=250) :: ltag 1071 logical :: lBcast 1072#ifdef MPI 1073 integer :: all_I(6) 1074 integer :: MPIerror 1075#endif 1076 1077 nullify(xa,lasto,isc_off) 1078 1079 ltag = trim(filename) 1080 if ( present(tag) ) ltag = trim(tag) 1081 1082 ! Determine whether to broadcast afterwards 1083 lBcast = .false. 1084 if ( present(Bcast) ) lBcast = Bcast 1085 1086 call ncdf_open(ncdf,filename,mode=NF90_NOWRITE) 1087 1088 call ncdf_inq_dim(ncdf,'na_u',len=na_u) 1089 call ncdf_inq_dim(ncdf,'no_u',len=no_u) 1090 call ncdf_inq_dim(ncdf,'spin',len=nspin) 1091 1092 call ncdf_get_var(ncdf,'nsc',nsc) 1093 call ncdf_get_var(ncdf,'cell',ucell) 1094 call ncdf_get_var(ncdf,'Ef',Ef) 1095 call ncdf_get_var(ncdf,'Qtot',Qtot) 1096 1097 ! The Brillouin zone sampling for siesta 1098 call ncdf_open_grp(ncdf,'SETTINGS',grp) 1099 call ncdf_get_var(grp,'BZ',kscell) 1100 call ncdf_get_var(grp,'BZ_displ',kdispl) 1101 call ncdf_get_var(grp,'ElectronicTemperature',Temp) 1102 1103#ifdef MPI 1104 if ( lBcast ) then 1105 ! Bcast initial sizes 1106 if ( Node == 0 ) then 1107 all_I(1:6) = (/na_u,no_u,nspin,nsc(1),nsc(2),nsc(3)/) 1108 end if 1109 call MPI_Bcast(all_I(1),6,MPI_Integer,0,MPI_Comm_World,MPIerror) 1110 na_u = all_I(1) 1111 no_u = all_I(2) 1112 nspin = all_I(3) 1113 nsc(1) = all_I(4) 1114 nsc(2) = all_I(5) 1115 nsc(3) = all_I(6) 1116 call MPI_Bcast(ucell(1,1),9,MPI_Double_Precision,0, & 1117 MPI_Comm_World,MPIerror) 1118 call MPI_Bcast(Gamma,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 1119 call MPI_Bcast(TSGamma,1,MPI_Logical,0,MPI_Comm_World,MPIerror) 1120 call MPI_Bcast(kscell(1,1),9,MPI_Integer,0,MPI_Comm_World,MPIerror) 1121 call MPI_Bcast(kdispl(1),3,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 1122 call MPI_Bcast(Ef,1,MPI_Double_Precision,0, MPI_Comm_World,MPIerror) 1123 call MPI_Bcast(Qtot,1,MPI_Double_Precision,0, MPI_Comm_World,MPIerror) 1124 call MPI_Bcast(Temp,1,MPI_Double_Precision,0, MPI_Comm_World,MPIerror) 1125 1126 end if 1127#endif 1128 1129 if ( Node == 0 ) then 1130 ! Read Geometry information 1131 allocate(xa(3,na_u)) 1132 allocate(lasto(0:na_u)) 1133 call ncdf_get_var(ncdf,'xa',xa) 1134 call ncdf_get_var(ncdf,'lasto',lasto(1:na_u)) 1135 lasto(0) = 0 1136 end if 1137 1138 call ncdf_open_grp(ncdf,'SPARSE',grp) 1139 1140 ! Number of supercells 1141 n_s = product(nsc) 1142 call re_alloc(isc_off,1,3,1,n_s) 1143 call ncdf_get_var(grp,'isc_off',isc_off) 1144 1145 ! Calculate the Gamma and TSGamma 1146 Gamma = ( n_s == 1 ) 1147 TSGamma = sum(kscell) == 1 1148 1149#ifdef MPI 1150 if ( lBcast ) then 1151 if ( Node /= 0 ) then 1152 allocate(xa(3,na_u)) 1153 allocate(lasto(0:na_u)) 1154 end if 1155 call MPI_Bcast(xa(1,1),3*na_u,MPI_Double_Precision,0,MPI_Comm_World,MPIerror) 1156 call MPI_Bcast(lasto(0),1+na_u,MPI_Integer,0, MPI_Comm_World,MPIerror) 1157 call MPI_Bcast(isc_off(1,1),3*n_s, MPI_Integer, 0, & 1158 MPI_Comm_World,MPIerror) 1159 end if 1160#endif 1161 1162 call cdf_r_Sp(grp, no_u, sp, tag = trim(filename), Bcast = Bcast ) 1163 call cdf_r_d1D(grp, 'S', sp, S, tag = trim(filename)//': S', & 1164 Bcast = Bcast ) 1165 dit => dist(S) 1166 call cdf_r_d2D(grp, 'H', sp, H, nspin, tag = trim(filename)//': H', & 1167 Bcast = Bcast, dit = dit ) 1168 1169 call ncdf_close(ncdf) 1170 1171 end subroutine ts_read_TSHS_nc 1172#endif 1173 1174 subroutine ts_write_TSHS(filename, & 1175 onlyS, Gamma, TSGamma, & 1176 ucell, nsc, isc_off, na_u, no_s, nspin, & 1177 kscell, kdispl, & 1178 xa, lasto, & 1179 H2D, S1D, indxuo, & 1180 Ef, Qtot, Temp, & 1181 istep, ia1) 1182 1183! ********************************************************************* 1184! Saves the hamiltonian and overlap matrices, and other data required 1185! to obtain the bands and density of states 1186! Writen by J.Soler July 1997. 1187! Note because of the new more compact method of storing H and S 1188! this routine is NOT backwards compatible 1189! Modified by M.Paulsson 2009 to: 1190! 1: To include information of which FC step for phonon calculations 1191! 2: To only save the overlap matrix if onlyS flag is set 1192! (Used for e-ph coupling calculations) 1193! 3: File format changed to unify Copenhagen/Barcelona Transiesta vers. 1194! 4: Smaller files by writing arrays directly instead of element wise 1195! *************************** INPUT ********************************** 1196! logical Gamma : Is only gamma point used? 1197! logical TSGamma : Is only TS gamma point used? 1198! logical onlyS : Should only overlap matrix be saved? 1199! ******************** INPUT or OUTPUT (depending on task) *********** 1200! integer no_u : Number of basis orbitals per unit cell 1201! integer no_s : Number of basis orbitals per supercell 1202! integer Enspin : Spin polarization (1 or 2) 1203! integer indxuo(no_s) : Index of orbitals in supercell 1204! integer n_nzs : First dimension of l_col, H, S and 1205! second of xij 1206! integer ncol(nuo) : Number of nonzero elements of each row 1207! of hamiltonian matrix 1208! integer l_ptr(nuo) : Pointer to the start of each row (-1) 1209! of hamiltonian matrix 1210! integer l_col(n_nzs) : Nonzero hamiltonian-matrix element column 1211! indexes for each matrix row 1212! real*8 H(n_nzs,Enspin) : Hamiltonian in sparse form 1213! real*8 S(n_nzs) : Overlap in sparse form 1214! real*8 qtot : Total number of electrons 1215! real*8 temp : Electronic temperature for Fermi smearing 1216! real*8 xij(3,n_nzs) : Vectors between orbital centers (sparse) 1217! (not read/written if only gamma point) 1218! TSS Begin 1219! ********************* ADDED ARGUMENTS FOR TRANSIESTA **************** 1220! integer fnlength : file name length 1221! character(fnlength) fname : file nema for input or output 1222! integer na_u : Number of atoms per unit cell 1223! integer istep, ia1 : Force constant step and atom number 1224! logical check_kcell : Do a check of the kscell and kdispl 1225! TSS End 1226! ********************************************************************* 1227 1228 use class_Sparsity 1229 use class_OrbitalDistribution 1230 use class_dSpData1D 1231 use class_dSpData2D 1232 use geom_helper, only : ucorb 1233 use m_sparse, only : xij_offset 1234 1235 use m_io_s, only: io_write_Sp, io_write_d1D, io_write_d2D 1236#ifdef MPI 1237 use mpi_siesta 1238#endif 1239 1240! ********************** 1241! * INPUT variables * 1242! ********************** 1243 character(len=*), intent(in) :: filename 1244 logical, intent(in) :: onlyS 1245 logical, intent(in) :: Gamma, TSGamma 1246 real(dp), intent(in) :: ucell(3,3) 1247 integer, intent(in) :: kscell(3,3) 1248 real(dp), intent(in) :: kdispl(3) 1249 integer, intent(in) :: nsc(3), na_u, no_s, nspin 1250 integer, intent(in) :: isc_off(3,product(nsc)) 1251 real(dp), intent(in) :: xa(3,na_u) 1252 type(dSpData2D), intent(inout) :: H2D 1253 type(dSpData1D), intent(inout) :: S1D 1254 integer, intent(in) :: indxuo(no_s) 1255 integer, intent(in) :: lasto(0:na_u) 1256 real(dp), intent(in) :: Ef 1257 real(dp), intent(in) :: Qtot, Temp 1258 integer, intent(in) :: istep, ia1 1259 1260! ************************ 1261! * LOCAL variables * 1262! ************************ 1263 type(OrbitalDistribution), pointer :: dit 1264 type(Sparsity), pointer :: sp 1265 integer :: iu, no_l, no_u, n_nzs 1266 integer :: i, n_s 1267 integer :: n_nzsg 1268 integer, allocatable, target :: gncol(:) 1269 integer, pointer :: ncol(:), l_ptr(:), l_col(:) 1270#ifdef MPI 1271 integer :: MPIerror 1272#endif 1273 1274 external :: io_assign, io_close 1275 1276#ifdef TRANSIESTA_DEBUG 1277 call write_debug( 'PRE ts_io_write' ) 1278#endif 1279 1280 ! Gather sparse pattern 1281 dit => dist(H2D) 1282 sp => spar(H2D) 1283 call attach(sp,nrows=no_l,nrows_g=no_u,nnzs=n_nzs, & 1284 n_col=ncol,list_ptr=l_ptr,list_col=l_col) 1285 1286 ! This setting is (I am afraid) not constant 1287 ! with system. I had suspected that n_s ALWAYS would 1288 ! equal no_s / no_u, but this seems not to be the case. 1289 n_s = no_s / no_u 1290 if ( mod(no_s,no_u) /= 0 ) call die('Error in supercell orbitals, no_s') 1291 if ( n_s /= product(nsc) ) call die('Error in supercell orbitals, nsc') 1292 1293 ! Check that indxuo is constructed in a sane format 1294 do i = 1 , no_s 1295 if ( indxuo(i) /= ucorb(i,no_u) ) & 1296 call die('Error in indxuo, could not understand the format, & 1297 &please consult the developers.') 1298 end do 1299 1300#ifdef MPI 1301 ! get total number of non-zero elements 1302 call MPI_Reduce(n_nzs,n_nzsg,1,MPI_Integer, MPI_SUM, 0, & 1303 MPI_Comm_World,MPIerror) 1304#else 1305 n_nzsg = n_nzs 1306#endif 1307 1308 if ( Node == 0 ) then 1309 1310 ! Open file 1311 call io_assign( iu ) 1312 open( iu, file=filename, form='unformatted', status='unknown' ) 1313 1314 ! Write file version 1315 write(iu) 1 ! This is version ONE of the file format 1316 1317 ! Write Dimensions Information 1318 write(iu) na_u, no_u, no_s, nspin, n_nzsg 1319 write(iu) nsc 1320 1321 ! Write Geometry information 1322 write(iu) ucell, xa 1323 1324 ! Write k-point sampling information 1325 write(iu) Gamma, TSGamma, onlyS 1326 write(iu) kscell, kdispl 1327 1328 ! Write Electronic Structure Information 1329 write(iu) Ef, Qtot, Temp 1330 1331 ! The phonon-calculation 1332 write(iu) istep, ia1 1333 1334 write(iu) lasto 1335 1336 end if 1337 1338 allocate(gncol(no_u)) 1339 gncol(1) = -1 1340 1341 ! Write out sparsity pattern 1342 call io_write_Sp(iu,sp,dit,gncol=gncol) 1343 1344 ! Write overlap matrix 1345 call io_write_d1D(iu,S1D,gncol=gncol) 1346 1347 if ( .not. onlyS ) then 1348 1349 ! Write Hamiltonian 1350 call io_write_d2D(iu,H2D,gncol=gncol) 1351 1352 end if 1353 1354 deallocate(gncol) 1355 1356 if ( Node == 0 ) then 1357 1358 if ( .not. Gamma ) then 1359 write(iu) isc_off 1360 end if 1361 1362 ! Close file 1363 call io_close( iu ) 1364 1365 end if 1366 1367#ifdef TRANSIESTA_DEBUG 1368 call write_debug( 'POS ts_io_write' ) 1369#endif 1370 1371 end subroutine ts_write_TSHS 1372 1373 function fname_TSHS(slabel,istep,onlyS,ia1) result(fname) 1374 character(len=*), intent(in) :: slabel 1375 integer, intent(in), optional :: istep 1376 integer, intent(in), optional :: ia1 1377 logical, intent(in), optional :: onlyS 1378 character(len=255) :: fname 1379 integer :: fL 1380 logical :: lonlyS 1381 integer :: listep, lia1 1382 1383 ! Initialize... 1384 fname = ' ' 1385 1386 lonlyS = .false. 1387 if ( present(onlyS) ) lonlyS = onlyS 1388 listep = -1 1389 if ( present(istep) ) listep = istep 1390 lia1 = 0 1391 if ( present(ia1) ) lia1 = ia1 1392 1393 ! This is the criteria for not doing an FCrun 1394 ! There will never be an atom denoted by index 0 1395 if ( lia1 /= 0 .and. listep >= 0 ) then 1396 if ( listep == 0 ) then 1397 ! We need the extra 00000's to let python 1398 ! sort in a stringent way... :( 1399 write(fname,'(a,i5.5)') '.',0 1400 else 1401 ! Be sure to have the python sorting working by adding 1402 ! zeroes 1403 write(fname,'(a,i5.5,a,i1)') '.',lia1,'-',listep 1404 end if 1405 fname = trim(slabel)//trim(fname) 1406 else if ( listep >= 0 ) then 1407 ! Simply a consecutive number increase. 1408 write(fname,'(a,i0)') '.',listep 1409 fname = trim(slabel)//trim(fname) 1410 else 1411 fname = slabel 1412 end if 1413 1414 fL = len_trim(fname) 1415 if ( lonlyS ) then 1416 fname = trim(fname)//'.onlyS' 1417 else 1418 fname = trim(fname)//'.TSHS' 1419 end if 1420 1421 end function fname_TSHS 1422 1423 subroutine FC_index(istep,ia1,ostep,oa1) 1424 integer, intent(in) :: istep 1425 integer, intent(in) :: ia1 1426 integer, intent(out) :: ostep 1427 integer, intent(out) :: oa1 1428 1429 if ( istep == 0 ) then 1430 ostep = 0 1431 oa1 = ia1 1432 else 1433 ostep = mod(istep-1,6) + 1 1434 oa1 = (istep - mod(istep-1,6))/6 + ia1 1435 end if 1436 1437 end subroutine FC_index 1438 1439 1440 function TSHS_version(fname) result(version) 1441 character(len=*), intent(in) :: fname 1442 integer :: version 1443 integer :: iu 1444 integer :: na_u, no_u, no_s, nspin, n_nzs, err 1445 1446 external :: io_assign, io_close 1447 1448 ! Initialize 1449 version = -1 1450 if ( Node /= 0 ) return 1451 1452 ! Open file 1453 call io_assign( iu ) 1454 open( iu, file=fname, form='unformatted', status='unknown' ) 1455 1456 read(iu,iostat=err) na_u, no_u, no_s, nspin, n_nzs 1457 if ( err == 0 ) then 1458 ! we can successfully read 5 integers 1459 version = 0 1460 else 1461 backspace(iu) 1462 read(iu,iostat=err) version 1463 end if 1464 1465 call io_close(iu) 1466 1467 end function TSHS_version 1468 1469end module m_ts_io 1470