xml_io_base.f90 41.9 KB
Newer Older
1
!
2
! Copyright (C) 2005-2014 Quantum ESPRESSO group
3 4 5 6 7 8 9 10 11
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!----------------------------------------------------------------------------
MODULE xml_io_base
  !----------------------------------------------------------------------------
  !
12
  ! ... this module contains some common subroutines used to read and write
13
  ! ... in XML format the data produced by Quantum ESPRESSO package
14 15
  !
  ! ... written by Carlo Sbraccia (2005)
16
  ! ... modified by Andrea Ferretti (2006-08)
17 18 19 20
  !
  USE iotk_module
  !
  USE kinds,     ONLY : DP
21
  USE io_files,  ONLY : tmp_dir, prefix, iunpun, xmlpun
22
  USE io_global, ONLY : ionode, ionode_id, stdout
23
  USE mp,        ONLY : mp_bcast
24
  USE parser,    ONLY : version_compare
25 26
  !
  IMPLICIT NONE
27
  PRIVATE
28
  !
29
  CHARACTER(iotk_attlenx)  :: attr
giannozz's avatar
giannozz committed
30
  LOGICAL,       SAVE      :: rho_binary = .TRUE.
31
  !
32
  PUBLIC :: rho_binary
33
  PUBLIC :: attr
34
  !
35 36
  PUBLIC :: read_wfc, write_wfc, read_rho, write_rho, &
       save_print_counter, read_print_counter
37
  PUBLIC :: create_directory, check_file_exst, restart_dir
38
  !
39 40 41 42 43 44
  CONTAINS
    !
    !------------------------------------------------------------------------
    SUBROUTINE create_directory( dirname )
      !------------------------------------------------------------------------
      !
45
      USE wrappers,  ONLY : f_mkdir_safe
46
      USE mp,        ONLY : mp_barrier
47
      USE mp_images, ONLY : me_image, intra_image_comm
48
      USE io_files,  ONLY : check_writable
49
      !
50
      CHARACTER(LEN=*), INTENT(IN) :: dirname
51
      !
52
      INTEGER                    :: ierr
53
      !
54
      CHARACTER(LEN=6), EXTERNAL :: int_to_char
55
      !
56
      IF ( ionode ) ierr = f_mkdir_safe( TRIM( dirname ) )
57
      CALL mp_bcast ( ierr, ionode_id, intra_image_comm )
58
      !
59
      CALL errore( 'create_directory', &
60
           'unable to create directory ' // TRIM( dirname ), ierr )
61
      !
62
      ! ... syncronize all jobs (not sure it is really useful)
63
      !
64
      CALL mp_barrier( intra_image_comm )
65
      !
66 67
      ! ... check whether the scratch directory is writable
      !
68
      IF ( ionode ) ierr = check_writable ( dirname, me_image )
69
      CALL mp_bcast( ierr, ionode_id, intra_image_comm )
70
      !
71
      CALL errore( 'create_directory:', &
72
                   TRIM( dirname ) // ' non existent or non writable', ierr )
73
      !
74 75
      RETURN
      !
76 77 78
    END SUBROUTINE create_directory
    !
    !------------------------------------------------------------------------
79
    FUNCTION restart_dir( outdir, runit )
80 81
      !------------------------------------------------------------------------
      !
82
      ! KNK_nimage
83
      ! USE mp_images, ONLY:  my_image_id
84
      CHARACTER(LEN=256)           :: restart_dir
85
      CHARACTER(LEN=*), INTENT(IN) :: outdir
86 87
      INTEGER,          INTENT(IN) :: runit
      !
88 89
      CHARACTER(LEN=256)         :: dirname
      INTEGER                    :: strlen
90
      CHARACTER(LEN=6), EXTERNAL :: int_to_char
91 92 93
      !
      ! ... main restart directory
      !
94
      dirname = TRIM( prefix ) // '_' // TRIM( int_to_char( runit ) )// '.save/'
95
      !
96
      IF ( LEN( outdir ) > 1 ) THEN
97
         !
98
         strlen = INDEX( outdir, ' ' ) - 1
99
         !
100
         dirname = outdir(1:strlen) // '/' // dirname
101 102 103 104 105
         !
      END IF
      !
      restart_dir = TRIM( dirname )
      !
106 107
      RETURN
      !
108 109
    END FUNCTION restart_dir
    !
110
    !------------------------------------------------------------------------
111
    FUNCTION check_restartfile( outdir, ndr )
112 113 114
      !------------------------------------------------------------------------
      !
      USE io_global, ONLY : ionode, ionode_id
115
      USE mp_images, ONLY : intra_image_comm
116 117 118 119 120
      !
      IMPLICIT NONE
      !
      LOGICAL                      :: check_restartfile
      INTEGER,          INTENT(IN) :: ndr
121
      CHARACTER(LEN=*), INTENT(IN) :: outdir
122 123 124 125
      CHARACTER(LEN=256)           :: filename
      LOGICAL                      :: lval
      !
      !
126
      filename = restart_dir( outdir, ndr )
127 128 129
      !
      IF ( ionode ) THEN
         !
130
         filename = TRIM( filename ) // '/' // TRIM( xmlpun )
131 132 133 134 135
         !
         INQUIRE( FILE = TRIM( filename ), EXIST = lval )
         !
      END IF
      !
136
      CALL mp_bcast( lval, ionode_id, intra_image_comm )
137 138 139 140 141 142 143 144
      !
      check_restartfile = lval
      !
      RETURN
      !
    END FUNCTION check_restartfile
    !
    !------------------------------------------------------------------------
145 146 147 148
    FUNCTION check_file_exst( filename )
      !------------------------------------------------------------------------
      !
      USE io_global, ONLY : ionode, ionode_id
149
      USE mp_images, ONLY : intra_image_comm
150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170
      !
      IMPLICIT NONE
      !
      LOGICAL          :: check_file_exst
      CHARACTER(LEN=*) :: filename
      !
      LOGICAL :: lexists
      !
      IF ( ionode ) THEN 
         !
         INQUIRE( FILE = TRIM( filename ), EXIST = lexists )
         !
      ENDIF
      !
      CALL mp_bcast ( lexists, ionode_id, intra_image_comm )
      !
      check_file_exst = lexists
      RETURN
      !
    END FUNCTION check_file_exst
    !
171 172
    !
    !------------------------------------------------------------------------
173
    SUBROUTINE save_print_counter( iter, outdir, wunit )
174 175
      !------------------------------------------------------------------------
      !
176
      ! ... a counter indicating the last successful printout iteration is saved
177 178
      !
      USE io_global, ONLY : ionode, ionode_id
179
      USE mp_images, ONLY : intra_image_comm
180 181 182 183 184
      USE mp,        ONLY : mp_bcast
      !
      IMPLICIT NONE
      !
      INTEGER,          INTENT(IN) :: iter
185
      CHARACTER(LEN=*), INTENT(IN) :: outdir
186 187 188 189 190 191
      INTEGER,          INTENT(IN) :: wunit
      !
      INTEGER            :: ierr
      CHARACTER(LEN=256) :: filename, dirname
      !
      !
192
      dirname = restart_dir( outdir, wunit )
193
      !
194 195
      CALL create_directory( TRIM( dirname ) )
      !
196 197
      IF ( ionode ) THEN
         !
198
         filename = TRIM( dirname ) // 'print_counter.xml'
199 200 201 202
         !
         CALL iotk_open_write( iunpun, FILE = filename, &
                             & ROOT = "PRINT_COUNTER",  IERR = ierr )
         !
203
      END IF
204
      !
205
      CALL mp_bcast( ierr, ionode_id, intra_image_comm )
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224
      !
      CALL errore( 'save_print_counter', &
                   'cannot open restart file for writing', ierr )
      !
      IF ( ionode ) THEN
         !
         CALL iotk_write_begin( iunpun, "LAST_SUCCESSFUL_PRINTOUT" )
         CALL iotk_write_dat(   iunpun, "STEP", iter )
         CALL iotk_write_end(   iunpun, "LAST_SUCCESSFUL_PRINTOUT" )
         !
         CALL iotk_close_write( iunpun )
         !
      END IF
      !
      RETURN
      !
    END SUBROUTINE save_print_counter
    !
    !------------------------------------------------------------------------
225
    SUBROUTINE read_print_counter( nprint_nfi, outdir, runit )
226 227
      !------------------------------------------------------------------------
      !
228
      ! ... the counter indicating the last successful printout iteration 
229 230 231
      ! ... is read here
      !
      USE io_global, ONLY : ionode, ionode_id
232
      USE mp_images, ONLY : intra_image_comm
233 234 235 236 237
      USE mp,        ONLY : mp_bcast
      !
      IMPLICIT NONE
      !
      INTEGER,          INTENT(OUT) :: nprint_nfi
238
      CHARACTER(LEN=*), INTENT(IN)  :: outdir
239 240 241 242 243 244
      INTEGER,          INTENT(IN)  :: runit
      !
      INTEGER            :: ierr
      CHARACTER(LEN=256) :: filename, dirname
      !
      !
245
      dirname = restart_dir( outdir, runit )
246 247 248
      !
      IF ( ionode ) THEN
         !
249
         filename = TRIM( dirname ) // 'print_counter.xml'
250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268
         !
         CALL iotk_open_read( iunpun, FILE = filename, IERR = ierr )
         !
         IF ( ierr > 0 ) THEN
            !
            nprint_nfi = -1
            !
         ELSE
            !
            CALL iotk_scan_begin( iunpun, "LAST_SUCCESSFUL_PRINTOUT" )
            CALL iotk_scan_dat(   iunpun, "STEP", nprint_nfi )
            CALL iotk_scan_end(   iunpun, "LAST_SUCCESSFUL_PRINTOUT" )
            !
            CALL iotk_close_read( iunpun )
            !
         END IF
         !
      END IF
      !
269
      CALL mp_bcast( nprint_nfi, ionode_id, intra_image_comm )
270 271 272 273 274 275
      !
      RETURN
      !
    END SUBROUTINE read_print_counter   
    !
    !------------------------------------------------------------------------
276
    SUBROUTINE set_kpoints_vars( ik, nk, kunit, ngwl, igl, &
277 278
                                 ngroup, ikt, iks, ike, igwx, ipmask, ipsour, &
                                 ionode, root_in_group, intra_group_comm, inter_group_comm, parent_group_comm )
279 280 281 282 283
      !------------------------------------------------------------------------
      !
      ! ... set working variables for k-point index (ikt) and 
      ! ... k-points number (nkt)
      !
284
      USE mp,         ONLY : mp_sum, mp_get, mp_max, mp_rank, mp_size
285 286 287 288 289
      !
      IMPLICIT NONE
      !
      INTEGER, INTENT(IN)  :: ik, nk, kunit
      INTEGER, INTENT(IN)  :: ngwl, igl(:)
290
      INTEGER, INTENT(OUT) :: ngroup
291 292
      INTEGER, INTENT(OUT) :: ikt, iks, ike, igwx
      INTEGER, INTENT(OUT) :: ipmask(:), ipsour
293 294
      LOGICAL, INTENT(IN)  :: ionode
      INTEGER, INTENT(IN)  :: root_in_group, intra_group_comm, inter_group_comm, parent_group_comm
295 296 297
      !
      INTEGER :: ierr, i
      INTEGER :: nkl, nkr, nkbl, nkt
298
      INTEGER :: nproc_parent, nproc_group, my_group_id, me_in_group, me_in_parent, io_in_parent
299
      !
300 301 302 303 304 305 306 307 308 309 310
      nproc_parent = mp_size( parent_group_comm )
      nproc_group  = mp_size( intra_group_comm )
      my_group_id  = mp_rank( inter_group_comm )
      me_in_group  = mp_rank( intra_group_comm )
      me_in_parent = mp_rank( parent_group_comm )
      !
      ! find the ID (io_in_parent) of the io PE ( where ionode == .true. )
      !
      io_in_parent = 0
      IF( ionode ) io_in_parent = me_in_parent
      CALL mp_sum( io_in_parent, parent_group_comm )
311 312 313 314 315 316
      !
      ikt = ik
      nkt = nk
      !
      ! ... find out the number of pools
      !
317
      ngroup = nproc_parent / nproc_group 
318 319 320 321 322 323 324
      !
      ! ... find out number of k points blocks
      !
      nkbl = nkt / kunit  
      !
      ! ... k points per pool
      !
325
      nkl = kunit * ( nkbl / ngroup )
326 327 328
      !
      ! ... find out the reminder
      !
329
      nkr = ( nkt - nkl * ngroup ) / kunit
330 331 332
      !
      ! ... Assign the reminder to the first nkr pools
      !
333
      IF ( my_group_id < nkr ) nkl = nkl + kunit
334 335 336
      !
      ! ... find out the index of the first k point in this pool
      !
337
      iks = nkl * my_group_id + 1
338
      !
339
      IF ( my_group_id >= nkr ) iks = iks + nkr * kunit
340 341 342 343 344 345
      !
      ! ... find out the index of the last k point in this pool
      !
      ike = iks + nkl - 1
      !
      ipmask = 0
346
      ipsour = io_in_parent
347 348 349 350
      !
      ! ... find out the index of the processor which collect the data 
      ! ... in the pool of ik
      !
351
      IF ( ngroup > 1 ) THEN
352 353 354
         !
         IF ( ( ikt >= iks ) .AND. ( ikt <= ike ) ) THEN
            !
355
            IF ( me_in_group == root_in_group ) ipmask( me_in_parent + 1 ) = 1
356 357 358
            !
         END IF
         !
359
         ! ... Collect the mask for all proc in the image
360
         !
361
         CALL mp_sum( ipmask, parent_group_comm )
362
         !
363
         DO i = 1, nproc_parent
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389
            !
            IF( ipmask(i) == 1 ) ipsour = ( i - 1 )
            !
         END DO
         !
      END IF
      !
      igwx = 0
      ierr = 0
      !
      IF ( ( ikt >= iks ) .AND. ( ikt <= ike ) ) THEN
         !
         IF ( ngwl > SIZE( igl ) ) THEN
            !
            ierr = 1
            !
         ELSE
            !
            igwx = MAXVAL( igl(1:ngwl) )
            !
         END IF
         !
      END IF
      !
      ! ... get the maximum index within the pool
      !
390
      CALL mp_max( igwx, intra_group_comm )
391 392 393
      !
      ! ... now notify all procs if an error has been found 
      !
394
      CALL mp_max( ierr, parent_group_comm )
395
      !
giannozz's avatar
giannozz committed
396
      CALL errore( 'set_kpoint_vars ', 'wrong size ngl', ierr )
397
      !
398 399
      IF ( ipsour /= io_in_parent ) &
         CALL mp_get( igwx, igwx, me_in_parent, io_in_parent, ipsour, 1, parent_group_comm )
400 401 402 403 404
      !
      RETURN
      !
    END SUBROUTINE set_kpoints_vars
    !
405 406
    !
    !------------------------------------------------------------------------
407
    SUBROUTINE write_rho( dirname, rho, nspin, extension )
408 409 410
      !------------------------------------------------------------------------
      !
      ! ... this routine writes the charge-density in xml format into the
411
      ! ... $dirname directory - $dirname must exist and end with '/'
412 413 414 415 416 417 418 419 420
      !
      USE fft_base, ONLY : dfftp
      USE io_global,ONLY : ionode
      USE mp_bands, ONLY : intra_bgrp_comm, inter_bgrp_comm
      !
      IMPLICIT NONE
      !
      INTEGER,          INTENT(IN)           :: nspin
      REAL(DP),         INTENT(IN)           :: rho(dfftp%nnr,nspin)
421
      CHARACTER(LEN=*), INTENT(IN)           :: dirname
422 423
      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: extension
      !
424
      CHARACTER(LEN=256)    :: file_base
425
      CHARACTER(LEN=6)      :: ext
426 427 428 429 430 431
      REAL(DP), ALLOCATABLE :: rhoaux(:)
      !
      !
      ext = ' '
      IF ( PRESENT( extension ) ) ext = '.' // TRIM( extension )
      !
432
      file_base = TRIM( dirname ) // 'charge-density' // TRIM( ext )
433 434 435
      !
      IF ( nspin == 1 ) THEN
         !
436
         CALL write_rho_xml( file_base, rho(:,1), dfftp, ionode, inter_bgrp_comm )
437 438 439 440 441 442 443
         !
      ELSE IF ( nspin == 2 ) THEN
         !
         ALLOCATE( rhoaux( dfftp%nnr ) )
         !
         rhoaux(:) = rho(:,1) + rho(:,2)
         !
444
         CALL write_rho_xml( file_base, rhoaux, dfftp, ionode, inter_bgrp_comm )
445
         !
446
         file_base = TRIM( dirname ) // 'spin-polarization' // TRIM( ext )
447 448 449
         !
         rhoaux(:) = rho(:,1) - rho(:,2)
         !
450
         CALL write_rho_xml( file_base, rhoaux,  dfftp, ionode, inter_bgrp_comm )
451 452 453 454 455
         !
         DEALLOCATE( rhoaux )
         !
      ELSE IF ( nspin == 4 ) THEN
         !
456
         CALL write_rho_xml( file_base, rho(:,1), dfftp, ionode, inter_bgrp_comm )
457
         !
458
         file_base = TRIM( dirname ) // 'magnetization.x' // TRIM( ext )
459
         !
460
         CALL write_rho_xml( file_base, rho(:,2), dfftp, ionode, inter_bgrp_comm )
461
         !
462
         file_base = TRIM( dirname ) // 'magnetization.y' // TRIM( ext )
463
         !
464
         CALL write_rho_xml( file_base, rho(:,3), dfftp, ionode, inter_bgrp_comm )
465
         !
466
         file_base = TRIM( dirname ) // 'magnetization.z' // TRIM( ext )
467
         !
468
         CALL write_rho_xml( file_base, rho(:,4), dfftp, ionode, inter_bgrp_comm )
469 470 471 472 473 474 475 476
         !
      END IF
      !
      RETURN
      !
    END SUBROUTINE write_rho
    !
    !------------------------------------------------------------------------
477
    SUBROUTINE read_rho( dirname, rho, nspin, extension)
478 479 480 481 482 483 484 485 486 487 488
      !------------------------------------------------------------------------
      !
      ! ... this routine reads the charge-density in xml format from the
      ! ... files saved into the '.save' directory
      !
      USE fft_base,  ONLY : dfftp
      USE io_global, ONLY : ionode
      !
      IMPLICIT NONE
      !
      INTEGER,          INTENT(IN)           :: nspin
489
      CHARACTER(LEN=*), INTENT(IN)           :: dirname
490
      CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: extension
491
      REAL(DP),         INTENT(OUT)          :: rho(dfftp%nnr,nspin)
492
      !
493 494
      CHARACTER(LEN=256)  :: file_base
      CHARACTER(LEN=6)    :: ext
495 496
      REAL(DP), ALLOCATABLE :: rhoaux(:)
      !
497 498
      ext = ' '
      IF ( PRESENT( extension ) ) ext = '.' // TRIM( extension )
499
      !
500
      file_base = TRIM( dirname ) // 'charge-density' // TRIM( ext )
501
      CALL read_rho_xml ( file_base, dfftp, rho(:,1) ) 
502 503 504 505 506 507 508
      !
      IF ( nspin == 2 ) THEN
         !
         rho(:,2) = rho(:,1)
         !
         ALLOCATE( rhoaux( dfftp%nnr ) )
         !
509
         file_base = TRIM( dirname ) // 'spin-polarization' // TRIM( ext )
510
         CALL read_rho_xml ( file_base, dfftp, rhoaux ) 
511 512 513 514 515 516 517 518
         !
         rho(:,1) = 0.5D0*( rho(:,1) + rhoaux(:) )
         rho(:,2) = 0.5D0*( rho(:,2) - rhoaux(:) )
         !
         DEALLOCATE( rhoaux )
         !
      ELSE IF ( nspin == 4 ) THEN
         !
519
         file_base = TRIM( dirname ) // 'magnetization.x' // TRIM( ext )
520
         CALL read_rho_xml ( file_base, dfftp, rho(:,2) ) 
521
         !
522
         file_base = TRIM( dirname ) // 'magnetization.y' // TRIM( ext )
523
         CALL read_rho_xml ( file_base, dfftp, rho(:,3) ) 
524
         !
525
         file_base = TRIM( dirname ) // 'magnetization.z' // TRIM( ext )
526
         CALL read_rho_xml ( file_base, dfftp, rho(:,4) ) 
527 528 529 530 531 532
         !
      END IF
      !
      RETURN
      !
    END SUBROUTINE read_rho
533
    !------------------------------------------------------------------------
534
    SUBROUTINE write_rho_xml( rho_file_base, rho, fft_desc, ionode, inter_group_comm )
535
      !------------------------------------------------------------------------
536
      !
537 538 539 540
      ! ... Writes charge density rho, one plane at a time.
      ! ... If ipp and npp are specified, planes are collected one by one from
      ! ... all processors, avoiding an overall collect of the charge density
      ! ... on a single proc.
541
      !
542
      USE mp,        ONLY : mp_get, mp_sum, mp_rank, mp_size
543
#if defined __HDF5
544
      USE hdf5_qe,  ONLY  : write_rho_hdf5, h5fclose_f, &
545
                            prepare_for_writing_final, add_attributes_hdf5, rho_hdf5_write  
546
#endif
547
      USE fft_types
548 549 550
      !
      IMPLICIT NONE
      !
551
      CHARACTER(LEN=*),  INTENT(IN) :: rho_file_base
552
      REAL(DP),          INTENT(IN) :: rho(:)
553 554
      TYPE(fft_type_descriptor),INTENT(IN) :: fft_desc
      INTEGER,           INTENT(IN) :: inter_group_comm
555
      LOGICAL,           INTENT(IN) :: ionode
556
      !
557 558
      INTEGER               :: nr1,nr2,nr3, nr1x, nr2x,nr3x
      INTEGER               :: rhounit, ierr, i, j, jj, k, kk, ldr, ip
559
      CHARACTER(LEN=256)    :: rho_file
560
      CHARACTER(LEN=256)    :: rho_file_hdf5
561
      CHARACTER(LEN=10)     :: rho_extension
562 563
      REAL(DP), ALLOCATABLE :: rho_plane(:)
      INTEGER,  ALLOCATABLE :: kowner(:)
564 565 566
      INTEGER               :: my_group_id, me_group, me_group2, me_group3, &
                                            nproc_group, nproc_group2, nproc_group3, &
                                            io_group_id, io_group2, io_group3
567
      INTEGER,  EXTERNAL    :: find_free_unit
568
      !
569

570
      my_group_id = mp_rank( inter_group_comm )
571 572 573 574 575 576

      me_group = fft_desc%mype ; me_group2 = fft_desc%mype2 ; me_group3 = fft_desc%mype3
      nproc_group = fft_desc%nproc ; nproc_group2 = fft_desc%nproc2 ; nproc_group3 = fft_desc%nproc3
      !
      nr1  = fft_desc%nr1  ; nr2  = fft_desc%nr2  ; nr3  = fft_desc%nr3
      nr1x = fft_desc%nr1x ; nr2x = fft_desc%nr2x ; nr3x = fft_desc%nr3x
577
      !
578 579 580 581
      rho_extension = '.dat'
      IF ( .NOT. rho_binary ) rho_extension = '.xml'
      !
      rho_file = TRIM( rho_file_base ) // TRIM( rho_extension )
582
      rhounit = find_free_unit ()
583
      !
584
      IF ( ionode ) THEN 
585
#if defined  __HDF5
586 587 588 589 590
         rho_file_hdf5 = TRIM( rho_file_base ) // '.hdf5'
         CALL prepare_for_writing_final(rho_hdf5_write, 0 ,rho_file_hdf5)
         CALL add_attributes_hdf5(rho_hdf5_write,nr1,"nr1")
         CALL add_attributes_hdf5(rho_hdf5_write,nr2,"nr2")
         CALL add_attributes_hdf5(rho_hdf5_write,nr3,"nr3")
591
#else
592
         CALL iotk_open_write( rhounit, FILE = rho_file,  BINARY = rho_binary, IERR = ierr )
593
         CALL errore( 'write_rho_xml', 'cannot open ' // TRIM( rho_file ) // ' file for writing', ierr )
594
#endif
595
      END IF 
596
      !
597
#if !defined __HDF5
598 599 600 601
      IF ( ionode ) THEN
         !
         CALL iotk_write_begin( rhounit, "CHARGE-DENSITY" )
         !
602 603 604
         CALL iotk_write_attr( attr, "nr1", nr1, FIRST = .TRUE. )
         CALL iotk_write_attr( attr, "nr2", nr2 )
         CALL iotk_write_attr( attr, "nr3", nr3 )
605 606 607 608
         !
         CALL iotk_write_empty( rhounit, "INFO", attr )
         !
      END IF
609
#endif
610
      !
611
      ALLOCATE( rho_plane( nr1*nr2 ) )
612 613
      ALLOCATE( kowner( nr3 ) )
      !
614
      ! ... find the index of the group (pool) that will write rho
615
      !
616
      io_group_id = 0
617
      !
618
      IF ( ionode ) io_group_id = my_group_id
619
      !
620 621
      CALL mp_sum( io_group_id, fft_desc%comm )
      CALL mp_sum( io_group_id, inter_group_comm ) ! io_group_id is the (pool) group that contains the ionode
622
      !
623
      ! ... find the index of the ionode within Y and Z  groups
624
      !
625 626 627 628
      io_group2 = 0 ; IF ( ionode ) io_group2 = me_group2
      CALL mp_sum( io_group2, fft_desc%comm )  ! io_group2 is the group index of the ionode in the Y group (nproc2)
      io_group3 = 0 ; IF ( ionode ) io_group3 = me_group3
      CALL mp_sum( io_group3, fft_desc%comm )  ! io_group3 is the group index of the ionode in the Z group (nproc3)
629
      !
630 631
      ! ... find out the owner of each "z" plane
      !
632
      DO ip = 1, nproc_group3
633
         !
634
         kowner( (fft_desc%i0r3p(ip)+1):(fft_desc%i0r3p(ip)+fft_desc%nr3p(ip)) ) = ip - 1
635
         !
636
      END DO
637
      !
638
      ldr = nr1x*fft_desc%my_nr2p
639
      !
640
      IF ( ( my_group_id == io_group_id ) ) THEN ! only the group of ionode collects and writes the data
641
         !
642
         DO k = 1, nr3
643
            !
644
            !  Only one subgroup write the charge density
645
            ! 
646 647
            rho_plane = 0.d0
            IF( ( kowner(k) == me_group3 ) ) THEN
648
               !
649 650 651
               kk = k - fft_desc%my_i0r3p
               ! 
               DO jj = 1, fft_desc%my_nr2p
652
                  !
653 654 655 656 657 658
                  j = jj + fft_desc%my_i0r2p
                  DO i = 1, nr1
                     !
                     rho_plane(i+(j-1)*nr1) = rho(i+(jj-1)*nr1x+(kk-1)*ldr)
                     !
                  END DO
659
                  !
660
               END DO
661
               call mp_sum(rho_plane, fft_desc%comm2 ) ! collect the data over the Y group (nproc2)
662
               !
663
            END IF
664
            !
665 666 667 668 669
            ! if this processor is in the same comm3 group as ionode (me_group2==io_group2) 
            IF ( kowner(k) /= io_group3 .and. me_group2==io_group2) & 
               CALL mp_get( rho_plane, rho_plane, me_group3, io_group3, kowner(k), k, fft_desc%comm3 )
            !
            IF ( ionode ) THEN
670
#if defined __HDF5
671
            CALL write_rho_hdf5(rho_hdf5_write,k,rho_plane)
672
#else
673
               CALL iotk_write_dat( rhounit, "z" // iotk_index( k ), rho_plane )
674
#endif
675 676 677
            ENDIF
            !
         END DO
678
         !
679
      END IF
680
      !
681 682
      DEALLOCATE( rho_plane )
      DEALLOCATE( kowner )
683
      !
684
      IF ( ionode ) THEN
685 686 687 688
#if defined __HDF5
         CALL h5fclose_f(rho_hdf5_write%file_id,ierr)
#else
   !
689 690 691
         CALL iotk_write_end( rhounit, "CHARGE-DENSITY" )
         !
         CALL iotk_close_write( rhounit )
692
#endif       
693 694
         !
      END IF
695
      !
696
      RETURN
697
      !
698 699
    END SUBROUTINE write_rho_xml
    !
700
    !------------------------------------------------------------------------
701
    SUBROUTINE read_rho_xml( rho_file_base, fft_desc, rho )
702 703
      !------------------------------------------------------------------------
      !
704 705
      ! ... Reads charge density rho, one plane at a time, to avoid 
      ! ... collecting the entire charge density on a single processor
706
      !
707
      USE io_global, ONLY : ionode, ionode_id
708
      USE mp_images, ONLY : intra_image_comm
709
      USE mp,        ONLY : mp_put, mp_sum, mp_rank, mp_size
710
#if defined __HDF5
711 712
      USE hdf5_qe,   ONLY : read_rho_hdf5, read_attributes_hdf5, &
           prepare_for_reading_final, h5fclose_f, rho_hdf5_write, hdf5_type
713
#endif
714
      USE fft_types
715 716 717
      !
      IMPLICIT NONE
      !
718
      CHARACTER(LEN=*),  INTENT(IN)  :: rho_file_base
719
      TYPE(fft_type_descriptor),INTENT(IN) :: fft_desc
720 721
      REAL(DP),          INTENT(OUT) :: rho(:)
      !
722 723 724 725
      INTEGER               :: rhounit, ierr, i, j, jj, k, kk, ldr, ip
      INTEGER               :: nr( 3 ), nr1_, nr2_, nr3_, nr1, nr2, nr3, nr1x, nr2x, nr3x
      INTEGER               :: me_group, me_group2, me_group3, &
                               nproc_group, nproc_group2, nproc_group3
726
      CHARACTER(LEN=256)    :: rho_file
727
      CHARACTER(LEN=256)    :: rho_file_hdf5
728 729
      REAL(DP), ALLOCATABLE :: rho_plane(:)
      INTEGER,  ALLOCATABLE :: kowner(:)
730
      LOGICAL               :: exst
731
      INTEGER,  EXTERNAL    :: find_free_unit
732 733 734
#if defined(__HDF5)
      TYPE(hdf5_type),ALLOCATABLE   :: h5desc
#endif
735
      !
736 737 738 739 740
      me_group = fft_desc%mype ; me_group2 = fft_desc%mype2 ; me_group3 = fft_desc%mype3
      nproc_group = fft_desc%nproc ; nproc_group2 = fft_desc%nproc2 ; nproc_group3 = fft_desc%nproc3
      !
      nr1  = fft_desc%nr1  ; nr2  = fft_desc%nr2  ; nr3  = fft_desc%nr3
      nr1x = fft_desc%nr1x ; nr2x = fft_desc%nr2x ; nr3x = fft_desc%nr3x
741
      !
742 743 744 745 746
#if defined(__HDF5)
      rho_file_hdf5 = TRIM( rho_file_base ) // '.hdf5'
      exst = check_file_exst(TRIM(rho_file_hdf5))
      IF ( .NOT. exst ) CALL errore ('read_rho_xml', 'searching for '// TRIM(rho_file_hdf5),10)
#else 
747
      rhounit = find_free_unit ( )
748 749 750 751
      rho_file = TRIM( rho_file_base ) // ".dat"
      exst = check_file_exst( TRIM(rho_file) ) 
      !
      IF ( .NOT. exst ) CALL errore('read_rho_xml', 'searching for '//TRIM(rho_file), 10)
752
#endif
753
      !
754
      IF ( ionode ) THEN
755 756 757 758 759 760 761
#if defined (__HDF5)
         ALLOCATE ( h5desc)
         CALL prepare_for_reading_final(h5desc, 0 ,rho_file_hdf5)
         CALL read_attributes_hdf5(h5desc, nr1_,"nr1")
         CALL read_attributes_hdf5(h5desc, nr2_,"nr2")
         CALL read_attributes_hdf5(h5desc, nr3_,"nr3")
         nr = [nr1_,nr2_,nr3_]
762
#else
763 764
         CALL iotk_open_read( rhounit, FILE = rho_file, IERR = ierr )
         CALL errore( 'read_rho_xml', 'cannot open ' // TRIM( rho_file ) // ' file for reading', ierr )
765 766 767 768
         CALL iotk_scan_begin( rhounit, "CHARGE-DENSITY" )
         !
         CALL iotk_scan_empty( rhounit, "INFO", attr )
         !
769 770 771
         CALL iotk_scan_attr( attr, "nr1", nr(1) )
         CALL iotk_scan_attr( attr, "nr2", nr(2) )
         CALL iotk_scan_attr( attr, "nr3", nr(3) )
772
#endif
773
         !
774 775 776
         IF ( nr1 /= nr(1) .OR. nr2 /= nr(2) .OR. nr3 /= nr(3) ) &
            CALL errore( 'read_rho_xml', 'dimensions do not match', 1 )
         !
777 778
      END IF
      !
779
      ALLOCATE( rho_plane( nr1*nr2 ) )
780 781
      ALLOCATE( kowner( nr3 ) )
      !
782
      DO ip = 1, nproc_group3
783
         !
784
         kowner( (fft_desc%i0r3p(ip)+1):(fft_desc%i0r3p(ip)+fft_desc%nr3p(ip)) ) = ip - 1
785
         !
786
      END DO
787
      !
788
      ldr = nr1x*fft_desc%my_nr2p
789
      !
790
      ! ... explicit initialization to zero is needed because the physical
791
      ! ... dimensions of rho may exceed the true size of the FFT grid 
792
      !
793
      rho(:) = 0.0_DP
794
      !
795 796
      DO k = 1, nr3
         !
797
         ! ... only ionode reads the charge planes
798
         !
799 800
         IF ( ionode ) THEN
#if defined __HDF5
801
            CALL  read_rho_hdf5(h5desc , k,rho_plane)
802
#else
803
            CALL iotk_scan_dat( rhounit, "z" // iotk_index( k ), rho_plane )
804 805
#endif
         ENDIF
806
         !
807
         ! ... planes are sent to the destination processor (all processors in this image)
808
         !
809
         CALL mp_bcast( rho_plane, ionode_id, intra_image_comm )
810
         !
811
         IF( kowner(k) == me_group3 ) THEN
812
            !
813 814 815
            kk = k - fft_desc%my_i0r3p
            DO jj = 1, fft_desc%my_nr2p
               j = jj + fft_desc%my_i0r2p
816
               DO i = 1, nr1
817
                  rho(i+(jj-1)*nr1x+(kk-1)*ldr) = rho_plane(i+(j-1)*nr1)
818 819 820 821 822 823 824 825 826 827 828 829
               END DO
            END DO
            !
         END IF
         !
      END DO
      !
      DEALLOCATE( rho_plane )
      DEALLOCATE( kowner )
      !
      IF ( ionode ) THEN
         !
830
#if defined __HDF5
831 832
         CALL h5fclose_f(h5desc%file_id,ierr)
         DEALLOCATE ( h5desc)
833 834
#else
   !
835 836 837
         CALL iotk_scan_end( rhounit, "CHARGE-DENSITY" )
         !
         CALL iotk_close_read( rhounit )
838
#endif    
839 840 841 842 843 844
      END IF
      !
      RETURN
      !
    END SUBROUTINE read_rho_xml
    !
845
    !------------------------------------------------------------------------
846 847 848
    ! ... methods to write and read wavefunctions
    !
    !------------------------------------------------------------------------
849 850 851 852
    SUBROUTINE write_wfc( iuni, ik, nk, kunit, ispin, nspin, wf0, ngw,   &
                          gamma_only, nbnd, igl, ngwl, filename, scalef, &
                          ionode, root_in_group, intra_group_comm,       &
                          inter_group_comm, parent_group_comm )
853 854 855
      !------------------------------------------------------------------------
      !
      USE mp_wave,    ONLY : mergewf
856
      USE mp,         ONLY : mp_get, mp_size, mp_rank, mp_sum
857
      USE control_flags,     ONLY : lwfnscf, lwfpbe0nscf  ! Lingzhu Kong
858 859 860 861 862 863 864 865 866
#if defined  __HDF5
      !USE hdf5_qe,    ONLY : evc_hdf5, read_data_hdf5, write_data_hdf5, &
      !                        evc_hdf5_write,  &
      !                       setup_file_property_hdf5, &
      !                       write_final_data, prepare_for_writing_final, &
      USE hdf5_qe                
      USE mp_global,    ONLY : inter_pool_comm, world_comm
      USE HDF5
#endif
867 868 869 870 871 872