read_xml_cards.f90 73.7 KB
Newer Older
1 2 3
!
!
!-------------------------------------------------------------!
4
! This module handles the cards reading for xml input         !
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
!                                                             !
!   written by Simone Ziraldo (08/2010)                       !
!-------------------------------------------------------------!
!
! cards not yet implemented:
! KSOUT
! AUTOPILOT
! ATOMIC_FORCES
! PLOT_WANNIER
! WANNIER_AC
! DIPOLE
! ESR
!
! to implement these cards take inspiration from file read_cards.f90
!
MODULE read_xml_cards_module
  !
  !
23
  USE io_global, ONLY : xmlinputunit => qestdin
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
  USE iotk_module, ONLY : iotk_scan_begin, iotk_scan_end, iotk_scan_dat,&
       iotk_scan_dat_inside, iotk_scan_attr, iotk_attlenx
  USE read_xml_fields_module, ONLY : clean_str
  USE kinds, ONLY : DP
  !
  USE io_global, ONLY : stdout
  !
  USE input_parameters
  !
  !
  IMPLICIT NONE
  !
  SAVE
  !
  PRIVATE
  !
  PUBLIC :: card_xml_atomic_species, card_xml_atomic_list, card_xml_chain, card_xml_cell, &
       card_xml_kpoints, card_xml_occupations, card_xml_constraints, card_xml_climbing_images, &
42
       card_xml_plot_wannier, card_default, card_bcast
43 44 45 46 47 48 49
  !
  !
  !
CONTAINS
  !
  !
  !--------------------------------------------------------------------------!
50
  !   This subroutine sets all the cards default values; as an input         !
51 52 53 54 55 56 57
  !   takes the card name that you want to set                               !
  !--------------------------------------------------------------------------!
  SUBROUTINE card_default( card )
    !
    !
    USE autopilot, ONLY : init_autopilot
    !
58 59
    USE read_namelists_module, ONLY : sm_not_set
    !
60 61 62 63 64 65 66 67 68 69 70 71 72
    !
    IMPLICIT NONE
    !
    !
    CHARACTER( len = * ),INTENT( IN ) :: card
    !
    !
    SELECT CASE ( trim(card) )
       !
    CASE ('INIT_AUTOPILOT')
       CALL init_autopilot()
       !
    CASE ('ATOMIC_LIST')
73 74 75 76 77 78
       !
       ! ... nothing to initialize
       ! ... because we don't have nat
       !
    CASE ('CHAIN' )
       !
79
       ! ... nothing to initialize
80
       ! ... because we don't have nat
81 82 83 84 85 86 87
       !
    CASE ('CELL')
       trd_ht = .false.
       rd_ht = 0.0_DP
       !
    CASE ('ATOMIC_SPECIES')
       atom_mass = 0.0_DP
88
       hubbard_u = 0.0_DP
89
       hubbard_j0 = 0.0_DP
90
       hubbard_alpha = 0.0_DP
91
       hubbard_beta = 0.0_DP
92 93 94 95 96 97 98 99 100
       starting_magnetization = sm_not_set
       starting_ns_eigenvalue = -1.0_DP
       angle1 = 0.0_DP
       angle2 = 0.0_DP
       ion_radius = 0.5_DP
       nhgrp = 0
       fnhscl = -1.0_DP
       tranp = .false.
       amprp = 0.0_DP
101 102 103
       !
    CASE ('K_POINTS')
       k_points = 'gamma'
104
       tk_inp   = .false.
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
       nkstot   = 1
       nk1      = 0
       nk2      = 0
       nk3      = 0
       k1       = 0
       k2       = 0
       k3       = 0
       !
    CASE ('OCCUPATIONS')
       tf_inp = .FALSE.
       !
    CASE ('CONSTRAINTS')
       nconstr_inp    = 0
       constr_tol_inp = 1.E-6_DP
       !
    CASE ('CLIMBING_IMAGES')
       ! ... nothing to initialize
       !
123 124 125 126
    CASE ('PLOT_WANNIER')
       !
       !       wannier_index =
       !
127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
    CASE ('KSOUT')
       ! ... not yet implemented in xml reading
       CALL allocate_input_iprnks( 0, nspin )
       nprnks  = 0
       !
    CASE ('DIPOLE')
       ! ... not yet implemented in xml reading
       tdipole_card = .FALSE.
    CASE ('ESR')
       ! ... not yet implemented in xml reading
       iesr_inp = 1
       !
    CASE ('ION_VELOCITIES')
       ! ... not yet implemented in xml reading
       tavel = .false.
       !
    CASE DEFAULT
       CALL errore ( 'card_default', 'You want to initialize a card that does &
            &not exist or is not yet implemented ( '//trim(card)//' card)', 1 )
       !
    END SELECT
    !
    !
  END SUBROUTINE card_default
  !
  !
  !
  !
  !---------------------------------------------------------------------------!
  !    This subroutine broadcasts the varibles defined in the various cards;  !
  !    the input string is the name of the card that you want to broadcast    !
  !---------------------------------------------------------------------------!
  SUBROUTINE card_bcast( card )
    !
    !
    USE io_global, ONLY : ionode, ionode_id                                                           
    !
    USE mp,        ONLY : mp_bcast
165
    USE mp_world,  ONLY : world_comm
166 167 168 169 170 171 172 173 174 175 176 177
    !
    IMPLICIT NONE
    !
    !
    CHARACTER( len = * ),INTENT( IN ) :: card
    INTEGER :: nspin0
    !
    !
    SELECT CASE ( trim(card) )
       !
       !
    CASE ( 'CELL' )
178 179 180 181 182 183 184 185 186 187 188
       CALL mp_bcast( ibrav, ionode_id, world_comm )
       CALL mp_bcast( celldm, ionode_id, world_comm )
       CALL mp_bcast( A, ionode_id, world_comm )
       CALL mp_bcast( B, ionode_id, world_comm )
       CALL mp_bcast( C, ionode_id, world_comm )
       CALL mp_bcast( cosAB, ionode_id, world_comm )
       CALL mp_bcast( cosAC, ionode_id, world_comm )
       CALL mp_bcast( cosBC, ionode_id, world_comm )
       CALL mp_bcast( cell_units, ionode_id, world_comm )
       CALL mp_bcast( rd_ht, ionode_id, world_comm )
       CALL mp_bcast( trd_ht, ionode_id, world_comm )
189 190
       !
    CASE ( 'ATOMIC_SPECIES' )
191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208
       CALL mp_bcast( ntyp, ionode_id, world_comm )
       CALL mp_bcast( atom_mass, ionode_id, world_comm )
       CALL mp_bcast( atom_pfile, ionode_id, world_comm )
       CALL mp_bcast( atom_label, ionode_id, world_comm )
       CALL mp_bcast( taspc, ionode_id, world_comm )
       CALL mp_bcast( hubbard_u, ionode_id, world_comm )
       CALL mp_bcast( hubbard_j0, ionode_id, world_comm )
       CALL mp_bcast( hubbard_alpha, ionode_id, world_comm )
       CALL mp_bcast( hubbard_beta, ionode_id, world_comm )
       CALL mp_bcast( starting_magnetization, ionode_id, world_comm )
       CALL mp_bcast( starting_ns_eigenvalue, ionode_id, world_comm )
       CALL mp_bcast( angle1, ionode_id, world_comm )
       CALL mp_bcast( angle2, ionode_id, world_comm )
       CALL mp_bcast( ion_radius, ionode_id, world_comm )
       CALL mp_bcast( nhgrp, ionode_id, world_comm )
       CALL mp_bcast( fnhscl, ionode_id, world_comm )
       CALL mp_bcast( tranp, ionode_id, world_comm )
       CALL mp_bcast( amprp, ionode_id, world_comm )
209 210
       !
    CASE ( 'ATOMIC_LIST' )
211 212 213
       CALL mp_bcast( atomic_positions, ionode_id, world_comm )
       CALL mp_bcast( nat, ionode_id, world_comm )
!       CALL mp_bcast( num_of_images, ionode_id, world_comm )
214 215 216 217
       ! ... ionode has already done it inside card_xml_atomic_list
       IF (.not.ionode) THEN
          CALL allocate_input_ions( ntyp, nat )
       END IF
218
!       CALL mp_bcast( pos, ionode_id )
219 220 221 222 223 224 225
       CALL mp_bcast( if_pos, ionode_id, world_comm )
       CALL mp_bcast( na_inp, ionode_id, world_comm )
       CALL mp_bcast( sp_pos, ionode_id, world_comm )
       CALL mp_bcast( rd_pos, ionode_id, world_comm )
       CALL mp_bcast( sp_vel, ionode_id, world_comm )
       CALL mp_bcast( rd_vel, ionode_id, world_comm )
       CALL mp_bcast( tapos, ionode_id, world_comm )
226
       !
227
    CASE ( 'CONSTRAINTS' )
228 229
       CALL mp_bcast( nconstr_inp, ionode_id, world_comm )
       CALL mp_bcast( constr_tol_inp, ionode_id, world_comm )
230
       IF ( .not.ionode ) CALL allocate_input_constr()
231 232 233 234
       CALL mp_bcast( constr_type_inp, ionode_id, world_comm )
       CALL mp_bcast( constr_target_inp, ionode_id, world_comm )
       CALL mp_bcast( constr_target_set, ionode_id, world_comm )
       CALL mp_bcast( constr_inp, ionode_id, world_comm )
235 236
       !
    CASE ( 'K_POINTS' )
237 238 239 240 241 242 243 244
       CALL mp_bcast( k_points, ionode_id, world_comm )
       CALL mp_bcast( nkstot, ionode_id, world_comm )
       CALL mp_bcast( nk1, ionode_id, world_comm )
       CALL mp_bcast( nk2, ionode_id, world_comm )
       CALL mp_bcast( nk3, ionode_id, world_comm )
       CALL mp_bcast( k1, ionode_id, world_comm )
       CALL mp_bcast( k2, ionode_id, world_comm )
       CALL mp_bcast( k3, ionode_id, world_comm )
245
       IF ( .not.ionode ) ALLOCATE( xk(3,MAX(1,nkstot)), wk(MAX(1,nkstot)) ) 
246 247
       CALL mp_bcast( xk, ionode_id, world_comm )
       CALL mp_bcast( wk, ionode_id, world_comm )
248 249 250 251 252 253 254
       !
    CASE ( 'OCCUPATIONS' )
       IF ( .not.ionode ) THEN
          nspin0 = nspin
          if ( nspin == 4 ) nspin0 = 1
          ALLOCATE( f_inp (nbnd, nspin0 ) )
       END IF
255
       CALL mp_bcast( f_inp, ionode_id, world_comm )
256
       !
257
    CASE ( 'PLOT_WANNIER' )
258
       CALL mp_bcast( wannier_index, ionode_id, world_comm )
259
       !
260 261 262 263 264 265 266 267 268 269 270 271 272
    CASE DEFAULT
       CALL errore ( 'card_bcast', 'You want to broadcast a card that does &
            &not exist or is not yet implemented', 1 )
       !
       !
    END SELECT
    !
    !
    !
  END SUBROUTINE card_bcast
  !
  !
  !-------------------------------------------------------------------------!
273
  ! Hereafter there are the reading of the xml cards                        !
274 275 276 277 278 279 280 281 282 283 284 285
  ! For more information see the Help file                                  !
  !-------------------------------------------------------------------------!
  !                                                                         !
  !                                                                         !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !                                                                         !
  !  CELL  (compulsory)                                                     !
  !                                                                         !
  !   specify the cell of your calculation                                  !
  !                                                                         !
  ! Syntax:                                                                 !
  !                                                                         !
286
  !   <cell type="type" sym="sym">                                          !
287
  !       depends on the type                                               !
288
  !   </cell>                                                               !
289 290 291 292 293 294 295
  !                                                                         !
  !  sym can be cubic or exagonal                                           !
  !                                                                         !
  !  if:                                                                    !
  !                                                                         !
  !  1) type is qecell, inside CELL node there is:                          !
  !                                                                         !
296 297
  !      <qecell ibrav="ibrav" alat="celldm(1)">                            !
  !         <real rank="1" n1="6">                                          !
298
  !            celldm(2) celldm(3) celldm(4) celldm(5) celldm(6)            !
299 300
  !         </real>                                                         !
  !      </qecell>                                                          !
301 302 303
  !                                                                         !
  !  2) type is abc, inside CELL node there is:                             !
  !                                                                         !
304
  !      <abc ibrav="ibrav">                                                !
305
  !          A B C cosAB cosAC cosBC                                        !
306
  !      </abc>                                                             !
307 308 309
  !                                                                         !
  !  3) type is matrix, inside there will be:                               !
  !                                                                         !
310
  !      <matrix units="units" alat="alat">                                   !
311
  !        <real rank="2" n1="3" n2="3">                                    !
312 313 314
  !          HT(1,1) HT(1,2) HT(1,3)                                        !
  !          HT(2,1) HT(2,2) HT(2,3)                                        !
  !          HT(3,1) HT(3,2) HT(3,3)                                        !
315 316
  !        </real>                                                          !
  !      </matrix>                                                          !
317 318 319 320 321 322
  !                                                                         !
  !                                                                         !
  !      Where:                                                             !
  !      HT(i,j) (real)  cell dimensions ( in a.u. ),                       !
  !                      note the relation with lattice vectors:            !
  !                      HT(1,:) = A1, HT(2,:) = A2, HT(3,:) = A3           !
323
  !      units            can be bohr (default) or alat (in this case you   !
324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340
  !                      have to specify alat)                              !
  !                                                                         !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !
  SUBROUTINE card_xml_cell ( )
    !
    IMPLICIT NONE
    !
    !
    CHARACTER( LEN = iotk_attlenx ) :: attr, attr2
    CHARACTER( LEN = 20 ) :: option,option2
    INTEGER :: i, j, ierr
    LOGICAL :: found
    REAL( kind = DP ), DIMENSION(6) :: vect_tmp
    !
    !
    !
341 342
    CALL iotk_scan_begin( xmlinputunit, 'cell', attr = attr, found = found, ierr = ierr )
    IF ( ierr /= 0 ) CALL errore( 'read_xml_cell', 'error scanning begin of cell &
343 344 345 346 347 348
         &card', ABS( ierr ) )
    !
    IF ( found ) THEN
       !
       CALL iotk_scan_attr( attr, 'type', option, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning type &
349
            &attribute of cell node', abs(ierr) )
350 351 352 353
       !
       !
       IF ( trim(option) == 'qecell' ) THEN
          !
354
          CALL iotk_scan_begin( xmlinputunit, 'qecell', attr2, ierr = ierr )
355
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning begin &
356
               &of qecell node', abs(ierr) )
357 358 359
          !
          CALL iotk_scan_attr( attr2, 'ibrav', ibrav, ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading ibrav &
360
               &attribute of qecell node', abs(ierr) )
361 362 363
          !
          CALL iotk_scan_attr(attr2, 'alat', celldm(1), ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading alat &
364
               &attribute of qecell node', abs(ierr) )
365 366 367
          !
          CALL iotk_scan_dat_inside( xmlinputunit, celldm(2:6), ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading data inside &
368
               &qecell node', abs(ierr) )
369
          !
370
          CALL iotk_scan_end( xmlinputunit, 'qecell', ierr = ierr )
371
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning end of &
372
               &qecell node', abs(ierr) )
373 374 375
          !
       ELSE IF ( trim(option) == 'abc' ) THEN
          !
376
          CALL iotk_scan_begin(xmlinputunit, 'abc', attr2, ierr = ierr )
377
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning begin &
378
               &of abc node', abs(ierr) )
379 380 381
          !
          CALL iotk_scan_attr( attr2, 'ibrav', ibrav, ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading ibrav &
382
               &attribute of abc node', abs(ierr) )
383 384 385
          !
          CALL iotk_scan_dat_inside( xmlinputunit, vect_tmp, ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading data inside &
386
               &abc node', abs(ierr) )
387 388 389 390 391 392 393 394
          !
          A = vect_tmp(1) 
          B = vect_tmp(2)
          C = vect_tmp(3)
          cosAB = vect_tmp(4)
          cosAC = vect_tmp(5)
          cosBC = vect_tmp(6)
          !
395
          CALL iotk_scan_end(xmlinputunit,'abc', ierr = ierr)
396
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning end of &
397
               &abc node', abs(ierr) )
398 399 400 401 402
          !
       ELSE IF (trim(option)=='matrix') THEN
          !
          ibrav = 0
          !
403
          CALL iotk_scan_begin( xmlinputunit, 'matrix', attr2, ierr = ierr )
404
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning begin &
405
               &of matrix node', abs(ierr) )
406 407 408
          !
          CALL iotk_scan_attr( attr2, 'units', option2, found = found, ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading units attribute &
409
               &of matrix node', abs(ierr) )
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437
          !
          IF (found) THEN
             IF ( trim(option2) == 'alat' ) THEN
                !
                cell_units = 'alat'
                !
                CALL iotk_scan_attr(attr2, 'alat', celldm(1), ierr = ierr )
                IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading alat&
                     &attribute of MATRIX node', abs(ierr) )
                !
             ELSE IF ( trim(option2) == 'bohr' ) THEN
                !
                cell_units = 'bohr'
                !
             ELSE
                !
                CALL errore( 'card_xml_cell', 'invalid units attribute', abs(ierr) )
                !
             END IF
          ELSE
             !
             cell_units = 'bohr'
             !
          END IF
          !
          !
          CALL iotk_scan_dat_inside( xmlinputunit, rd_ht, ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error reading data inside &
438
               &matrix node', abs(ierr) )
439 440 441 442
          !
          rd_ht = transpose( rd_ht )
          trd_ht = .TRUE.
          !
443
          CALL iotk_scan_end( xmlinputunit, 'matrix', ierr = ierr )
444
          IF ( ierr /= 0 ) CALL errore( 'card_xml_cell', 'error scanning end of &
445
               &matrix node', abs(ierr) )
446 447
          !
       ELSE
448
          CALL errore( 'card_xml_cell', 'type '//trim(option)//' in cell node does not exist', 1 )
449 450
       END IF
       !
451 452
       CALL iotk_scan_end( xmlinputunit, 'cell', ierr = ierr)
       IF ( ierr /= 0 ) CALL errore( 'read_xml_pw', 'error scanning end of cell &
453 454 455
            &card', ABS( ierr ) )
    ELSE
       !
456
       CALL errore( 'read_xml_pw', 'cell card not found', 1 )
457 458 459 460 461 462 463 464 465 466 467 468 469 470
       !
    END IF
    !
    !
    RETURN
    !
  END SUBROUTINE card_xml_cell
  !
  !
  !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !                                                                         !
  ! ATOMIC_SPECIES  (compulsory)                                            !
  !                                                                         !
471
  !   set the atomic species and their pseudopotential files                !
472 473 474
  !                                                                         !
  ! Syntax:                                                                 !
  !                                                                         !
475 476
  !    <atomic_species ntyp="ntyp">                                         !
  !                                                                         !
477 478 479 480 481 482 483
  !       <specie name="label(i)">                                          !
  !                                                                         !
  !           <property name="mass">                                        !
  !              <real>                                                     !
  !                 mass(i)                                                 !
  !              </real>                                                    !
  !           </property>                                                   !
484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557
  !                                                                         !
  !           <property name="pseudofile">                                  !
  !              <string>                                                   !
  !                 psfile(i)                                               !
  !              </string>                                                  !
  !           </property>                                                   !
  ![ optional                                                               !
  !           <property name="starting_magnetization">                      !
  !              <real>                                                     !
  !                 starting_magnetization(i)                               !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="hubbard_alpha">                               !
  !              <real>                                                     !
  !                 hubbard_alpha(i)                                        !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="hubbard_u">                                   !
  !              <real>                                                     !
  !                 hubbard_alpha(i)                                        !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="starting_ns_eigenvalue" ispin="" ns="">       !
  !              <real>                                                     !
  !                 starting_ns_eigenvalue(ns , ispin, i )                  !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="angle1">                                      !
  !              <real>                                                     !
  !                 angle1(i)                                               !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="angle2">                                      !
  !              <real>                                                     !
  !                 angle2(i)                                               !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="ion_radius">                                  !
  !              <real>                                                     !
  !                 ion_radius(i)                                           !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="nhgrp">                                       !
  !              <integer>                                                  !
  !                 nhgrp(i)                                                !
  !              </integer>                                                 !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="fnhscl">                                      !
  !              <real>                                                     !
  !                 fnhscl(i)                                               !
  !              </real>                                                    !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="tranp">                                       !
  !              <logical>                                                  !
  !                 tranp(i)                                                !
  !              </logical>                                                 !
  !           </property>                                                   !
  !                                                                         !
  !           <property name="amprp">                                       !
  !              <real>                                                     !
  !                 amprp(i)                                                !
  !              </real>                                                    !
  !           </property>                                                   !
  !]                                                                        !
  !       </specie>                                                         !
558 559
  !       ....                                                              !
  !       ....                                                              !
560
  !    </atomic_species>                                                    !
561 562 563
  !                                                                         !
  ! Where:                                                                  !
  !                                                                         !
564 565
  !      only the pseudofile property is compulsory, the others are optional!
  !                                                                         !
566 567 568
  !      label(i)  ( character(len=4) )  label of the atomic species        !
  !      mass(i)   ( real )              atomic mass                        !
  !                                      ( in u.m.a, carbon mass is 12.0 )  !
569
  !      psfile(i) ( character(len=80) ) pseudopotential filename           !
570 571 572 573 574 575 576 577
  !                                                                         !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !
  SUBROUTINE card_xml_atomic_species( )
    !
    IMPLICIT NONE
    !
    !
578
    INTEGER            :: is, ip, ierr, direction
579 580 581
    CHARACTER( LEN = 4 )   :: lb_pos
    CHARACTER( LEN = 256 ) :: psfile
    CHARACTER( LEN = iotk_attlenx ) :: attr, attr2
582
    LOGICAL :: found, psfile_found
583 584 585
    !
    !
    !
586 587
    CALL iotk_scan_begin( xmlinputunit, 'atomic_species', attr = attr, found = found, ierr = ierr )
    IF ( ierr /= 0 ) CALL errore( 'read_xml_pw', 'error scanning begin of atomic_species &
588 589 590 591 592 593
         &card', ABS( ierr ) )
    !
    IF ( found ) THEN
       !
       CALL iotk_scan_attr( attr, 'ntyp', ntyp, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
594
            &reading ntyp attribute inside atomic_species node', abs( ierr ) )
595 596 597 598 599 600 601
       !
       IF( ntyp < 0 .OR. ntyp > nsx ) &
            CALL errore( 'card_xml_atomic_species', &
            ' ntyp is too large', MAX( ntyp, 1) )
       !
       DO is = 1, ntyp
          !
602
          CALL iotk_scan_begin( xmlinputunit, 'specie', attr = attr2, ierr = ierr )
603
          IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
604
               &scanning specie node', abs( ierr ) )
605 606 607
          !
          CALL iotk_scan_attr( attr2, 'name', lb_pos, ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
608
               &reading name attribute of specie node', abs( ierr ) )
609
          !
610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638
          psfile_found = .false.
          !
          DO
             CALL iotk_scan_begin( xmlinputunit, 'property', attr = attr2, &
                  direction = direction, ierr = ierr )
             IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
                  &scanning begin property node', abs( ierr ) )
             !
             IF (direction == -1) EXIT
             !
             CALL read_property( attr2 )
             !
             !
             CALL iotk_scan_end( xmlinputunit, 'property', ierr = ierr )
             IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
                  &scanning end of property node', abs( ierr ) )

          END DO
          !
          CALL iotk_scan_end( xmlinputunit, 'property', ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
               &scanning end of property node', abs( ierr ) )
          !
          CALL iotk_scan_end( xmlinputunit, 'specie', ierr = ierr )
          IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_species',  'error &
               &scanning end of specie node', abs( ierr ) )
          !
          IF  (.not. psfile_found ) CALL errore( 'card_xml_atomic_species', &
               'no pseudofile found', abs( is ) )
639 640 641 642 643
          !
          atom_pfile(is) = trim( psfile )
          lb_pos         = adjustl( lb_pos )
          atom_label(is) = trim( lb_pos )
          !
644

645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660
          !
          DO ip = 1, is - 1
             !
             IF ( atom_label(ip) == atom_label(is) ) THEN
                CALL errore( ' card_xml_atomic_species ',  &
                     ' two occurrences of the same atomic label',  is )
             ENDIF
          ENDDO
          !
       ENDDO
       !
       ! ... this variable is necessary to mantain compatibility.
       ! ... With new xml input the compulsory of atomic_species is already given
       !
       taspc = .true.
       !
661 662 663
       CALL iotk_scan_end( xmlinputunit, 'atomic_species', ierr = ierr )
       IF (ierr/=0)  CALL errore( 'card_xml_atomic_species', 'error scanning end of &
            &atomic_species node', ABS( ierr ) )
664 665 666
       !
    ELSE
       !
667
       CALL errore( 'read_xml_pw', 'atomic_species  card not found', 1 )
668 669 670 671 672
       !
    ENDIF
    !
    RETURN
    !
673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688
  CONTAINS
    !
    SUBROUTINE read_property ( attr_in)
      !
      IMPLICIT NONE
      !
      CHARACTER( len = * ), INTENT( in ) :: attr_in
      INTEGER :: index1, index2
      CHARACTER( len = 50 ) :: prop_name
      !
      CALL iotk_scan_attr( attr_in, 'name', prop_name, ierr = ierr )
      IF (ierr/=0)  CALL errore( 'card_xml_atomic_species', 'error reading name &
           &attribute of property node', ABS( is ) )

      SELECT CASE ( trim(prop_name) )
         !
689 690 691
      CASE ( 'mass' )
         CALL iotk_scan_dat_inside( xmlinputunit, atom_mass(is) , ierr = ierr)
         !
692 693 694 695 696 697 698 699 700 701 702 703 704
      CASE ( 'pseudofile' )
         CALL iotk_scan_dat_inside( xmlinputunit, psfile, ierr = ierr)
         psfile = clean_str( psfile )
         psfile_found = .true.
         !
      CASE ( 'starting_magnetization' )
         CALL iotk_scan_dat_inside( xmlinputunit, starting_magnetization( is ),&
              ierr = ierr)
         !
      CASE ( 'hubbard_alpha' )
         CALL iotk_scan_dat_inside( xmlinputunit, hubbard_alpha( is ),&
              ierr = ierr)
         !
705 706 707 708
      CASE ( 'hubbard_beta' )
         CALL iotk_scan_dat_inside( xmlinputunit, hubbard_beta( is ),&
              ierr = ierr)
         !
709 710 711 712
      CASE ( 'hubbard_u' )
         CALL iotk_scan_dat_inside( xmlinputunit, hubbard_u( is ),&
              ierr = ierr)
         !
713 714 715 716
      CASE ( 'hubbard_j0' )
         CALL iotk_scan_dat_inside( xmlinputunit, hubbard_j0( is ),&
              ierr = ierr)
         !
717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773
      CASE ( 'starting_ns_eigenvalue' )
         !
         CALL iotk_scan_attr( attr_in, 'ns', index1, ierr = ierr )
         IF (ierr/=0)  CALL errore( 'card_xml_atomic_species', 'error reading ns &
              &attribute of property node', ABS( is ) )
         !
         CALL iotk_scan_attr( attr_in, 'ispin', index2, ierr = ierr )
         IF (ierr/=0)  CALL errore( 'card_xml_atomic_species', 'error reading ispin &
              &attribute of property node', ABS( is ) )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, &
              starting_ns_eigenvalue( index1, index2, is), ierr = ierr)
         !
      CASE ( 'angle1' )
         CALL iotk_scan_dat_inside( xmlinputunit, angle1( is ),&
              ierr = ierr)
         !
      CASE ( 'angle2' )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, angle2( is ),&
              ierr = ierr)
         !
      CASE ( 'ion_radius' )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, ion_radius( is ),&
              ierr = ierr)
         !
      CASE ( 'nhgrp' )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, nhgrp( is ),&
              ierr = ierr)
         !
      CASE ( 'fnhscl' )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, fnhscl( is ),&
              ierr = ierr)
         !
      CASE ( 'tranp' )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, tranp( is ),&
              ierr = ierr)
         !
      CASE ( 'amprp' )
         !
         CALL iotk_scan_dat_inside( xmlinputunit, amprp( is ),&
              ierr = ierr)
         !
      CASE DEFAULT
         CALL errore( 'card_xml_atomic_species', 'property '&
              //trim(prop_name)//' not known', abs( is ) )
      END SELECT
      !
      !
      IF ( ierr /= 0 )  CALL errore( 'card_xml_atomic_species', 'error reading ' &
           //trim(prop_name)//' data', abs( is ) )
      !
    END SUBROUTINE read_property
774 775 776 777 778 779 780 781 782 783 784 785 786 787
    !
  END SUBROUTINE card_xml_atomic_species
  !
  !
  !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !                                                                         !
  !                                                                         !
  ! ATOMIC_LIST (compulsory for PW)                                         !
  !                                                                         !
  !   set the atomic positions                                              !
  !                                                                         !
  ! Syntax:                                                                 !
  !                                                                         !
788 789 790 791
  !  <atomic_list units="units_option" nat="natom">                         !
  !     <atom name="label(1)">                                              !
  !        <position ifx="mbl(1,1)" ify="mbl(2,1)" ifz="mbl(3,1)">          !
  !           <real rank="1" n1="3">                                        !
792
  !              tau(1,1)  tau(2,1)  tau(3,1)                               !
793 794 795
  !           </real>                                                       !
  !        </position>                                                      !
  !     </atom>                                                             !
796
  !     ...                                                                 !
797
  !  </atomic_list>                                                         !
798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822
  !                                                                         !
  ! Where:                                                                  !
  !                                                                         !
  !   units_option == crystal   position are given in scaled units          !
  !   units_option == bohr      position are given in Bohr                  !
  !   units_option == angstrom  position are given in Angstrom              !
  !   units_option == alat      position are given in units of alat         !
  !                                                                         !
  !   label(k) ( character(len=4) )  atomic type                            !
  !   tau(:,k) ( real )              coordinates  of the k-th atom          !
  !   mbl(:,k) ( integer )           mbl(i,k) > 0 the i-th coord. of the    !
  !                                  k-th atom is allowed to be moved       !
  !                                                                         !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !
  SUBROUTINE card_xml_atomic_list( )
    !
    IMPLICIT NONE
    !
    !
    CHARACTER( len = iotk_attlenx ) :: attr
    INTEGER :: ierr, is
    LOGICAL :: found
    !
    !
823
    CALL iotk_scan_begin( xmlinputunit, 'atomic_list', attr, ierr = ierr )
824
    IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_list', 'error scanning begin &
825
         &of atomic_list node', abs(ierr) )
826 827 828
    !
    CALL iotk_scan_attr( attr, 'units', atomic_positions, found = found, ierr = ierr )
    IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_list', 'error reading units &
829
         &attribute of atomic_list node', abs(ierr) )
830 831 832 833 834 835 836 837 838
    !
    IF ( found ) THEN
       IF ( (trim( atomic_positions ) == 'crystal') .or. &
            (trim( atomic_positions ) == 'bohr') .or. &
            (trim( atomic_positions ) == 'angstrom').or. &
            (trim( atomic_positions ) == 'alat') ) THEN
          atomic_positions = trim( atomic_positions )
       ELSE
          CALL errore( 'car_xml_atom_lists',  &
839
               'error in units attribute of atomic_list node, unknow '&
840 841 842 843 844 845 846 847 848
               & //trim(atomic_positions)//' units', 1 )
       ENDIF
    ELSE
       ! ... default value
       atomic_positions = 'alat'
    ENDIF
    !
    CALL iotk_scan_attr( attr, 'nat', nat, ierr = ierr )
    IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_list', 'error reading nat attribute &
849
         &of atomic_list node', abs(ierr) )
850 851 852 853 854 855 856 857 858 859 860
    !
    IF ( nat < 1 ) THEN
       CALL errore( 'card_xml_atomic_list',  'nat out of range',  nat )
    END IF
    !
    ! ... allocation of needed arrays
    CALL allocate_input_ions( ntyp, nat )
    !
    if_pos = 1
    sp_pos = 0
    rd_pos = 0.0_DP
861 862
    sp_vel = 0
    rd_vel = 0.0_DP
863 864 865
    na_inp = 0
    !
    !
866
    CALL read_image( 1, rd_pos, rd_vel )
867
    !
868
    CALL iotk_scan_end( xmlinputunit, 'atomic_list', ierr = ierr )
869
    IF ( ierr /= 0 ) CALL errore( 'card_xml_atomic_list', 'error scanning end of &
870
         &atomic_list node', abs( ierr ) )
871 872 873 874 875 876 877 878 879 880 881 882 883
    !
    !
    tapos = .true.
    !
    RETURN
    !
    !
  END SUBROUTINE card_xml_atomic_list
  !
  !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-!
  !                                                                            !
  !                                                                            !
884
  ! CHAIN  (used in neb and smd calculation) OBSOLETE, NOT IMPLEMENTED         !
885 886 887 888 889
  !                                                                            !
  !   set the atomic positions for a chian                                     !
  !                                                                            !
  ! Syntax:                                                                    !
  !                                                                            !
890 891 892 893 894
  !  <chain num_of_images="">                                                  !
  !     <atomic_list units="units_option" nat="natom" num="1">                 !
  !        <atom name="label(1)" ifx="mbl(1,1)" ify="mbl(2,1)" ifz="mbl(3,1)"> !
  !          <position>                                                        !
  !             <real rank="1" n1="3">                                         !
895
  !                 tau(1,1)  tau(2,1)  tau(3,1)                               !
896 897 898
  !             </real>                                                        !
  !          </position>
  !        </atom>                                                             !
899
  !        ...                                                                 !
900 901
  !     </atomic_list>                                                         !
  !     <atomic_list num="2">                                                  !
902
  !        ...                                                                 !
903
  !     </atomic_list>                                                         !
904
  !     ...                                                                    !
905
  !  </chain>                                                                  !
906 907 908 909
  !                                                                            !
  !                                                                            !
  ! Where:                                                                     !
  !                                                                            !
910 911
  ! notation of atomic_list node is the same of the atomic_list cards.         !
  ! the difference is that inside the chain card you put more atomic_list node !
912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933
  ! with the attribute num that indicates the number of the image              !
  !                                                                            !
  !                                                                            !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-!
  !
  SUBROUTINE card_xml_chain( )
    !
    IMPLICIT NONE
    !
    !
    CHARACTER( LEN = iotk_attlenx ) :: attr
    LOGICAL :: found,end_of_chain
    INTEGER :: ierr
    REAL (DP), DIMENSION( :, :), ALLOCATABLE :: tmp_image
    !
    !
    end_of_chain = .false.

    RETURN
    !
  END SUBROUTINE card_xml_chain
  !
934
  ! ... Subroutine that reads a single image inside chain node
935
  !
936
  SUBROUTINE read_image( image, image_pos, image_vel )
937 938 939 940
    !
    IMPLICIT NONE
    !
    INTEGER, INTENT( in ) :: image
941 942
    REAL( DP ), INTENT( inout ), DIMENSION( 3, nat ) :: image_pos
    REAL( DP ), INTENT( inout ), DIMENSION( 3, nat ), OPTIONAL :: image_vel
943 944 945 946 947
    !
    !
    INTEGER :: ia, idx, ierr, is, direction
    CHARACTER( len = iotk_attlenx ) :: attr
    CHARACTER( len = 4 ) :: lb_pos
948
    LOGICAL :: found_vel, read_vel
949 950 951 952 953 954
    REAL( DP ) :: default
    !
    default = 1.0_DP
    !
    ia = 0
    !
955 956 957
    read_vel = .true.
    IF (present(image_vel)) read_vel = .true.
    !
958 959
    DO
       !
960
       CALL iotk_scan_begin( xmlinputunit, 'atom', attr = attr, &
961 962
            direction = direction, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'read_image', 'error scanning begin of &
963
            &atom node', abs(ierr) )
964 965 966
       !
       IF (direction == -1) THEN
          IF (ia < nat) CALL errore( 'read_image', &
967
               'less atoms than axpected in atomic_list', image )
968 969 970 971 972 973
          EXIT
       END IF
       !
       ia = ia + 1
       !
       IF ( ia > nat) CALL errore( 'read_image', &
974
            'more atoms than axpected in atomic_list', image )
975 976 977 978
       !
       ! ... compulsory name attribute
       CALL iotk_scan_attr( attr, 'name', lb_pos, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'read_image', 'error reading &
979 980 981 982 983
            &name attribute of atom node', abs(ierr) )
       !
       CALL iotk_scan_dat( xmlinputunit,'position', image_pos( 1:3, ia ), attr = attr, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'read_image', 'error reading position data of &
            &atom node', abs(ierr) )
984
       !
985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002
       IF (read_vel) THEN
          CALL iotk_scan_begin( xmlinputunit, 'velocity', &
               found = found_vel, ierr = ierr)
          IF ( ierr /= 0 ) CALL errore( 'read_al_image', 'error scanning begin of &
               &velocity node', abs(ierr) )
          !
          IF (found_vel) THEN
             !
             CALL iotk_scan_dat_inside( xmlinputunit, image_vel( 1:3, ia ), ierr = ierr )
             IF ( ierr /= 0 ) CALL errore( 'read_al_image', 'error reading &
                  &velocity', abs(ierr) )
             !
             CALL iotk_scan_end( xmlinputunit, 'velocity', ierr = ierr)
             IF ( ierr /= 0 ) CALL errore( 'read_al_image', 'error scanning end of &
                  &velocity node', abs(ierr) )
             !
          ENDIF
       ENDIF
1003
       !
1004 1005

       CALL iotk_scan_end( xmlinputunit, 'atom', ierr = ierr )
1006
       IF ( ierr /= 0 ) CALL errore( 'read_image', 'error scanning end of &
1007
            &atom node', abs(ierr) )
1008 1009 1010 1011 1012 1013
       !
       !
       IF ( image  == 1 ) THEN
          !
          CALL iotk_scan_attr( attr, 'ifx', if_pos(1,ia), default = 1, ierr=ierr )
          IF ( ierr /= 0) CALL errore( 'read_image', &
1014
               'error reading ifx attribute of atom node', image )
1015 1016 1017
          !
          CALL iotk_scan_attr( attr, 'ify', if_pos(2,ia), default = 1, ierr = ierr )
          IF ( ierr /= 0) CALL errore( 'read_image', &
1018
               'error reading ify attribute of atom node', image )
1019 1020 1021
          !
          CALL iotk_scan_attr( attr, 'ifz', if_pos(3,ia), default = 1, ierr = ierr )
          IF ( ierr /= 0) CALL errore( 'read_image', &
1022
               'error reading ifz attribute of atom node', image )
1023 1024 1025 1026 1027 1028 1029 1030
          !
          lb_pos = adjustl( lb_pos )
          !
          match_label_path: DO is = 1, ntyp
             !
             IF ( trim( lb_pos ) == trim( atom_label(is) ) ) THEN
                !
                sp_pos( ia ) = is
1031
                IF (found_vel .and. read_vel) sp_vel( ia) = is 
1032 1033 1034 1035 1036 1037 1038 1039
                !
                EXIT match_label_path
                !
             ENDIF
             !
          ENDDO match_label_path
          !
          IF ( ( sp_pos( ia ) < 1 ) .or. ( sp_pos( ia ) > ntyp ) ) CALL errore( &
1040
               'read_image', 'wrong name in atomic_list node', ia )
1041 1042 1043 1044 1045 1046 1047 1048 1049
          !
          is = sp_pos( ia )
          !
          na_inp( is ) = na_inp( is ) + 1
          !
       ENDIF
       !
    ENDDO
    !
1050
    CALL iotk_scan_end( xmlinputunit, 'atom', ierr = ierr )
1051
    IF ( ierr /= 0 ) CALL errore( 'read_image', 'error scanning end of &
1052
         &atom node', abs(ierr) )
1053
    !
1054 1055 1056 1057 1058 1059 1060 1061
    IF ( image == 1) THEN
       DO is = 1, ntyp
          IF( na_inp( is ) < 1 ) &
               CALL errore( 'read_image', 'no atom found in atomic_list for '&
               //trim(atom_label(is))//' specie', is )
       ENDDO
    ENDIF
   !
1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075
    RETURN
    !
  END SUBROUTINE read_image
  !
  !
  !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !                                                                         !
  !   K_POINTS                                                              !
  !                                                                         !
  !   use the specified set of k points                                     !
  !                                                                         !
  ! Syntax:                                                                 !
  !                                                                         !
1076
  ! <k_points type="mesh_option">                                           !
1077 1078
  !                                                                         !
  ! if mesh_option = tpiba, crystal, tpiba_b or crystal_b :                 !
1079 1080
  !      <mesh npoints="n">                                                 !
  !         <real rank="2" n1="4" n2="n">                                   !
1081 1082 1083 1084
  !                                                                         !
  !            xk(1,1) xk(2,1) xk(3,1) wk(1)                                !
  !            ...     ...     ...     ...                                  !
  !            xk(1,n) xk(2,n) xk(3,n) wk(n)                                !
1085 1086
  !         </real>                                                         !
  !      </mesh>                                                            !
1087 1088
  !                                                                         !
  ! else if mesh_option = automatic                                         !
1089 1090
  !      <mesh>                                                             !
  !         <real rank="1" n1="6">                                          !
1091
  !             nk1 nk2 nk3 k1 k2 k3                                        !
1092 1093
  !         </real>                                                         !
  !      </mesh>                                                            !
1094
  !                                                                         !
1095
  ! </k_points>                                                             !
1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139
  !                                                                         !
  !                                                                         !
  ! Where:                                                                  !
  !                                                                         !
  !   mesh_option == automatic  k points mesh is generated automatically    !
  !                             with Monkhorst-Pack algorithm               !
  !   mesh_option == crystal    k points mesh is given in stdin in scaled   !
  !                             units                                       !
  !   mesh_option == tpiba      k points mesh is given in stdin in units    !
  !                             of ( 2 PI / alat )                          !
  !   mesh_option == gamma      only gamma point is used ( default in       !
  !                             CPMD simulation )                           !
  !   mesh_option == tpiba_b    as tpiba but the weights gives the          !
  !                             number of points between this point         !
  !                             and the next                                !
  !   mesh_option == crystal_b  as crystal but the weights gives the        !
  !                             number of points between this point and     !
  !                             the next                                    !
  !                                                                         !
  !   n       ( integer )  number of k points                               !
  !   xk(:,i) ( real )     coordinates of i-th k point                      !
  !   wk(i)   ( real )     weights of i-th k point                          !
  !                                                                         !
  !                                                                         !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !
  SUBROUTINE card_xml_kpoints( attr )
    !
    IMPLICIT NONE
    !
    CHARACTER( len = * ), INTENT( in ) :: attr
    !
    LOGICAL :: kband = .FALSE.
    CHARACTER( len = 20 ) :: type
    CHARACTER( len = iotk_attlenx ) :: attr2
    INTEGER :: i,j, nkaux, ierr
    INTEGER, DIMENSION( 6 ) :: tmp
    INTEGER, DIMENSION( : ), ALLOCATABLE :: wkaux
    REAL( DP ), DIMENSION( : , : ), ALLOCATABLE :: points_tmp, xkaux
    REAL( DP ) :: delta
    !
    !
    CALL iotk_scan_attr(attr, 'type', type, ierr = ierr )
    IF ( ierr /= 0 ) CALL errore( 'card_xml_kpoints', 'error reading type attribute &
1140
         &of k_points node', abs( ierr ) )
1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178
    !
    SELECT CASE ( trim( type ) )
       !
    CASE ('automatic')
       !automatic generation of k-points
       k_points = 'automatic'
       !
    CASE ('crystal')
       !  input k-points are in crystal (reciprocal lattice) axis
       k_points = 'crystal'
       !
    CASE ('crystal_b')
       k_points = 'crystal'
       kband=.true.
       !
    CASE ('tpiba')
       !  input k-points are in 2pi/a units
       k_points = 'tpiba'
       !
    CASE ('tpiba_b')
       k_points = 'tpiba'
       kband=.true.
       !
    CASE ('gamma')
       !  Only Gamma (k=0) is used
       k_points = 'gamma'
       !
    CASE DEFAULT
       !  by default, input k-points are in 2pi/a units
       k_points = 'tpiba'
       !
    END SELECT
    !
    IF ( k_points == 'automatic' ) THEN
       !
       ! ... automatic generation of k-points
       !
       nkstot = 0
1179
       ALLOCATE ( xk(3,1), wk(1) )
1180 1181
       CALL iotk_scan_dat( xmlinputunit, 'mesh', tmp, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'card_xml_kpoints', 'error reading data inside mesh &
1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204
            &node', abs( ierr ) )
       !
       nk1 = tmp( 1 )
       nk2 = tmp( 2 )
       nk3 = tmp( 3 )
       k1  = tmp( 4 )
       k2  = tmp( 5 )
       k3  = tmp( 6 )
       !
       ! ... some checks
       !
       IF ( k1 < 0 .or. k1 > 1 .or. &
               k2 < 0 .or. k2 > 1 .or. &
               k3 < 0 .or. k3 > 1 ) CALL errore &
               ('card_xml_kpoints', 'invalid offsets: must be 0 or 1', 1)
       !
       IF ( nk1 <= 0 .or. nk2 <= 0 .or. nk3 <= 0 ) CALL errore &
            ('card_xml_kpoints', 'invalid values for nk1, nk2, nk3', 1)
       !
    ELSE IF ( ( k_points == 'tpiba' ) .OR. ( k_points == 'crystal' ) ) THEN
       !
       ! ... input k-points are in 2pi/a units
       !
1205 1206
       CALL iotk_scan_begin( xmlinputunit, 'mesh', attr2, ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'card_xml_kpoints', 'error scanning begin of mesh &
1207 1208 1209
            &node', abs( ierr ) )
       !
       CALL iotk_scan_attr( attr2, 'npoints', nkstot, ierr = ierr )
1210
       IF ( ierr /= 0 ) CALL errore( 'card_xml_kpoints', 'error reading attribute npoints of mesh &
1211 1212 1213
            &node', abs( ierr ) )
       !
       !
1214 1215
       !IF ( nkstot > size( xk, 2 )  ) CALL errore &
       !     ('card_xml_kpoints', 'too many k-points', nkstot)
1216 1217 1218 1219
       !
       allocate( points_tmp(4,nkstot) )
       !
       CALL iotk_scan_dat_inside( xmlinputunit, points_tmp, ierr = ierr )
1220
       IF ( ierr /= 0 ) CALL errore( 'card_xml_kpoints', 'error reading data inside mesh &
1221 1222
            &node', abs( ierr ) )
       !
1223 1224
       ALLOCATE ( xk(3,nkstot), wk(nkstot) )
       !
1225 1226 1227 1228 1229
       xk( :, 1:nkstot ) = points_tmp( 1:3, : )
       wk( 1:nkstot ) = points_tmp( 4, : )
       !
       deallocate( points_tmp )
       !
1230 1231
       CALL iotk_scan_end( xmlinputunit, 'mesh', ierr = ierr )
       IF ( ierr /= 0 ) CALL errore( 'card_xml_kpoints', 'error scanning end of mesh &
1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273
            &node', abs( ierr ) )
       !
       !
       IF ( kband ) THEN
          !
          nkaux=nkstot
          !
          allocate( xkaux( 3, nkstot ) )
          allocate( wkaux( nkstot ) )
          !
          xkaux( :, 1:nkstot ) = xk( :, 1:nkstot )
          wkaux( 1:nkstot ) = nint( wk(1:nkstot) )
          nkstot = 0
          !
          DO i = 1, nkaux-1
             !
             delta = 1.0_DP/wkaux(i)
             !
             DO j=0, wkaux(i)-1
                !
                nkstot=nkstot+1
                IF ( nkstot > SIZE (xk,2)  ) CALL errore &
                     ('card_xml_kpoints', 'too many k-points',nkstot)
                !
                xk( :, nkstot ) = xkaux( :, i ) + delta*j*( xkaux(:,i+1) - xkaux(:,i) ) 
                wk(nkstot)=1.0_DP
                !
             ENDDO
             !
          ENDDO
          !
          nkstot = nkstot + 1
          xk( :, nkstot ) = xkaux( :, nkaux )
          wk( nkstot ) = 1.0_DP
          !
          deallocate(xkaux)
          deallocate(wkaux)
       ENDIF
       !
    ELSE IF ( k_points == 'gamma' ) THEN
       !
       nkstot = 1
1274
       ALLOCATE ( xk(3,1), wk(1) )
1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296
       xk(:, 1) = 0.0_DP
       wk(1) = 1.0_DP
       !
    ENDIF
    !
    tk_inp = .TRUE.
    !
    RETURN
    !
    !
  END SUBROUTINE card_xml_kpoints
  !
  !
  !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !                                                                         !
  ! OCCUPATIONS (optional)                                                  !
  !                                                                         !
  !   use the specified occupation numbers for electronic states.           !
  !                                                                         !
  ! Syntax (nspin == 1) or (nspin == 4):                                    !
  !                                                                         !
1297 1298
  !   <occupations>                                                         !
  !      <real rank="1" n1="nbnd">                                          !
1299 1300 1301 1302
  !         f(1)                                                            !
  !         ....                                                            !
  !         ....                                                            !
  !         f(nbnd)                                                         !
1303 1304
  !      </real>                                                            !
  !   </occupations>                                                        !
1305 1306 1307
  !                                                                         !
  ! Syntax (nspin == 2):                                                    !
  !                                                                         !
1308 1309
  !   <occupations>                                                         !
  !         <real rank="2" n1="nbnd" n2="2">                                !
1310 1311
  !            u(1) ... u(nbnd)                                             !
  !            d(1) ... d(nbnd)                                             !
1312 1313
  !         </real>                                                         !
  !   </occupations>                                                        !
1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348
  !                                                                         !
  ! Where:                                                                  !
  !                                                                         !
  !      f(:) (real)  these are the occupation numbers                      !
  !                   for LDA electronic states.                            !
  !                                                                         !
  !      u(:) (real)  these are the occupation numbers                      !
  !                   for LSD spin == 1 electronic states                   !
  !      d(:) (real)  these are the occupation numbers                      !
  !                   for LSD spin == 2 electronic states                   !
  !                                                                         !
  !                                                                         !
  !_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_!
  !
  SUBROUTINE card_xml_occupations( )
    !
    !
    IMPLICIT NONE
    ! 
    INTEGER :: nspin0, ierr
    REAL( DP ), ALLOCATABLE :: tmp_data(:)
    !
    !
    nspin0 = nspin
    IF (nspin == 4) nspin0 = 1
    !
    IF (nbnd==0) CALL errore( 'card_xml_occupation', 'nbdn is not defined ', 1 )
    !
    allocate ( f_inp ( nbnd, nspin0 ) )
    !
    IF ( nspin0 == 2 ) THEN
       !
       CALL iotk_scan_dat_inside( xmlinputunit, f_inp, ierr = ierr )
       !
       IF ( ierr /= 0 ) CALL errore( 'card_xml_occupations