electrons_base.f90 16.4 KB
Newer Older
1
!
2
! Copyright (C) 2002-2009 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 electrons_base
!------------------------------------------------------------------------------!

12
      USE kinds, ONLY: DP
13 14 15 16
!
      IMPLICIT NONE
      SAVE

cavazzon's avatar
cavazzon committed
17 18
      INTEGER :: nbnd       = 0    !  number electronic bands, each band contains
                                   !  two spin states
19
      INTEGER :: nbndx      = 0    !  array dimension nbndx >= nbnd
cavazzon's avatar
cavazzon committed
20
      INTEGER :: nspin      = 0    !  nspin = number of spins (1=no spin, 2=LSDA)
21 22 23 24
      INTEGER :: nel(2)     = 0    !  number of electrons (up, down)
      INTEGER :: nelt       = 0    !  total number of electrons ( up + down )
      INTEGER :: nupdwn(2)  = 0    !  number of states with spin up (1) and down (2)
      INTEGER :: iupdwn(2)  = 0    !  first state with spin (1) and down (2)
cavazzon's avatar
cavazzon committed
25
      INTEGER :: nudx       = 0    !  max (nupdw(1),nupdw(2))
cavazzon's avatar
cavazzon committed
26
      INTEGER :: nbsp       = 0    !  total number of electronic states 
27
                                   !  (nupdwn(1)+nupdwn(2))
cavazzon's avatar
cavazzon committed
28
      INTEGER :: nbspx      = 0    !  array dimension nbspx >= nbsp
29 30 31 32 33 34 35 36
      !
      INTEGER :: nupdwn_bgrp(2)  = 0    !  number of states with spin up (1) and down (2) in this band group
      INTEGER :: iupdwn_bgrp(2)  = 0    !  first state with spin (1) and down (2) in this band group
      INTEGER :: nudx_bgrp       = 0    !  max (nupdw_bgrp(1),nupdw_bgrp(2)) in this band group
      INTEGER :: nbsp_bgrp       = 0    !  total number of electronic states 
                                        !  (nupdwn_bgrp(1)+nupdwn_bgrp(2)) in this band group
      INTEGER :: nbspx_bgrp      = 0    !  array dimension nbspx_bgrp >= nbsp_bgrp local to the band group
      INTEGER :: i2gupdwn_bgrp(2)= 0    !  global index of the first local band
cavazzon's avatar
cavazzon committed
37 38

      LOGICAL :: telectrons_base_initval = .FALSE.
39 40
      LOGICAL :: keep_occ = .FALSE.  ! if .true. when reading restart file keep 
                                     ! the occupations calculated in initval
cavazzon's avatar
cavazzon committed
41

42
      REAL(DP), ALLOCATABLE :: f(:)   ! occupation numbers ( at gamma )
43
      REAL(DP) :: qbac = 0.0_DP       ! background neutralizing charge
44
      INTEGER, ALLOCATABLE :: ispin(:) ! spin of each state
45 46 47

      REAL(DP), ALLOCATABLE :: f_bgrp(:)     ! occupation numbers ( at gamma )
      INTEGER, ALLOCATABLE  :: ispin_bgrp(:) ! spin of each state
48
      INTEGER, ALLOCATABLE :: ibgrp_g2l(:)    ! local index of the i-th global band index
49
!
cavazzon's avatar
cavazzon committed
50 51 52 53
!------------------------------------------------------------------------------!
  CONTAINS
!------------------------------------------------------------------------------!

cavazzon's avatar
cavazzon committed
54

55
    SUBROUTINE electrons_base_initval( zv_ , na_ , nsp_ , nbnd_ , nspin_ , &
56
          occupations_ , f_inp, tot_charge_, tot_magnetization_ )
cavazzon's avatar
cavazzon committed
57

58 59
      USE constants,         ONLY   : eps8
      USE io_global,         ONLY   : stdout
cavazzon's avatar
cavazzon committed
60

61 62
      REAL(DP),         INTENT(IN) :: zv_ (:), tot_charge_
      REAL(DP),         INTENT(IN) :: f_inp(:,:)
63
      REAL(DP),         INTENT(IN) :: tot_magnetization_
64
      INTEGER,          INTENT(IN) :: na_ (:) , nsp_
65
      INTEGER,          INTENT(IN) :: nbnd_ , nspin_
cavazzon's avatar
cavazzon committed
66
      CHARACTER(LEN=*), INTENT(IN) :: occupations_
67 68 69

      REAL(DP)                     :: nelec, nelup, neldw, ocp, fsum
      INTEGER                      :: iss, i, in
cavazzon's avatar
cavazzon committed
70

71
      nspin = nspin_
72 73 74
      !
      ! ... set nelec
      !
75 76 77 78 79
      nelec = 0.0_DP
      DO i = 1, nsp_
         nelec = nelec + na_ ( i ) * zv_ ( i )
      END DO 
      nelec = nelec - tot_charge_
80 81 82
      !
      ! ... set nelup/neldw
      !
83 84
      nelup = 0._dp
      neldw = 0._dp
85
      call set_nelup_neldw (tot_magnetization_, nelec, nelup, neldw )
cavazzon's avatar
cavazzon committed
86

87 88 89 90 91 92 93
      IF( ABS( nelec - ( nelup + neldw ) ) > eps8 ) THEN
         CALL errore(' electrons_base_initval ',' inconsistent n. of electrons ', 2 )
      END IF
      !
      !   Compute the number of bands
      !
      IF( nbnd_ /= 0 ) THEN
94
        nbnd  = nbnd_                          ! nbnd is given from input
95
      ELSE
96
        nbnd  = NINT( MAX( nelup, neldw ) )    ! take the maximum between up and down states
97 98 99
      END IF


cavazzon's avatar
cavazzon committed
100 101 102
      IF( nelec < 1 ) THEN
         CALL errore(' electrons_base_initval ',' nelec less than 1 ', 1 )
      END IF
103 104
      !
      IF( ABS( NINT( nelec ) - nelec ) > eps8 ) THEN
cavazzon's avatar
cavazzon committed
105 106
         CALL errore(' electrons_base_initval ',' nelec must be integer', 2 )
      END IF
107
      !
cavazzon's avatar
cavazzon committed
108 109
      IF( nbnd < 1 ) &
        CALL errore(' electrons_base_initval ',' nbnd out of range ', 1 )
110
      !
cavazzon's avatar
cavazzon committed
111

112 113 114 115
      IF ( nspin /= 1 .AND. nspin /= 2 ) THEN
        WRITE( stdout, * ) 'nspin = ', nspin
        CALL errore( ' electrons_base_initval ', ' nspin out of range ', 1 )
      END IF
cavazzon's avatar
cavazzon committed
116

117 118 119 120 121
      IF( MOD( nbnd, 2 ) == 0 ) THEN
         nbspx = nbnd * nspin
      ELSE
         nbspx = ( nbnd + 1 ) * nspin
      END IF
cavazzon's avatar
cavazzon committed
122 123

      ALLOCATE( f     ( nbspx ) )
124
      ALLOCATE( ispin ( nbspx ) )
125
      f     = 0.0_DP
126
      ispin = 0
cavazzon's avatar
cavazzon committed
127 128 129 130 131 132 133 134 135

      iupdwn ( 1 ) = 1
      nel          = 0

      SELECT CASE ( TRIM(occupations_) )
      CASE ('bogus')
         !
         ! bogus to ensure \sum_i f_i = Nelec  (nelec is integer)
         !
136
         f ( : )    = nelec / nbspx
cavazzon's avatar
cavazzon committed
137
         nel (1)    = nint( nelec )
138
         nupdwn (1) = nbspx
cavazzon's avatar
cavazzon committed
139 140 141 142 143 144 145 146 147 148 149
         if ( nspin == 2 ) then
            !
            ! bogus to ensure Nelec = Nup + Ndw
            !
            nel (1) = ( nint(nelec) + 1 ) / 2
            nel (2) =   nint(nelec)       / 2
            nupdwn (1)=nbnd
            nupdwn (2)=nbnd
            iupdwn (2)=nbnd+1
         end if
         !
150 151
         keep_occ = .true.
         !
cavazzon's avatar
cavazzon committed
152 153 154 155
      CASE ('from_input')
         !
         ! occupancies have been read from input
         !
156 157 158 159
         ! count electrons
         !
         IF( nspin == 1 ) THEN
            nelec = SUM( f_inp( :, 1 ) )
160 161
            nelup = nelec / 2.0_DP
            neldw = nelec / 2.0_DP
162 163 164 165 166 167
         ELSE
            nelup = SUM ( f_inp ( :, 1 ) )
            neldw = SUM ( f_inp ( :, 2 ) )
            nelec = nelup + neldw 
         END IF
         !
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192
         ! consistency check
         !
         IF( nspin == 1 ) THEN
           IF( f_inp( 1, 1 ) <= 0.0_DP )  &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
         ELSE
           IF( f_inp( 1, 1 ) < 0.0_DP )  &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
           IF( f_inp( 1, 2 ) < 0.0_DP )  &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
           IF( ( f_inp( 1, 1 ) + f_inp( 1, 2 ) ) == 0.0_DP )  &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
         END IF
         DO i = 2, nbnd
           IF( nspin == 1 ) THEN
             IF( f_inp( i, 1 ) > 0.0_DP .AND. f_inp( i-1, 1 ) <= 0.0_DP )  &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
           ELSE
             IF( f_inp( i, 1 ) > 0.0_DP .AND. f_inp( i-1, 1 ) <= 0.0_DP ) &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
             IF( f_inp( i, 2 ) > 0.0_DP .AND. f_inp( i-1, 2 ) <= 0.0_DP ) &
               CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 )
           END IF
         END DO
         !
193 194 195 196 197 198
         ! count bands
         !
         nupdwn (1) = 0
         nupdwn (2) = 0
         DO i = 1, nbnd
           IF( nspin == 1 ) THEN
199
             IF( f_inp( i, 1 ) > 0.0_DP ) nupdwn (1) = nupdwn (1) + 1
200
           ELSE
201 202
             IF( f_inp( i, 1 ) > 0.0_DP ) nupdwn (1) = nupdwn (1) + 1
             IF( f_inp( i, 2 ) > 0.0_DP ) nupdwn (2) = nupdwn (2) + 1
203 204 205
           END IF
         END DO
         !
cavazzon's avatar
cavazzon committed
206
         if( nspin == 1 ) then
207 208
           nel (1)    = nint( nelec )
           iupdwn (1) = 1
cavazzon's avatar
cavazzon committed
209
         else
210 211 212 213
           nel (1)    = nint(nelup)
           nel (2)    = nint(neldw)
           iupdwn (1) = 1
           iupdwn (2) = nupdwn (1) + 1
cavazzon's avatar
cavazzon committed
214 215
         end if
         !
216 217 218 219 220 221
         DO iss = 1, nspin
           DO in = iupdwn ( iss ), iupdwn ( iss ) - 1 + nupdwn ( iss )
             f( in ) = f_inp( in - iupdwn ( iss ) + 1, iss )
           END DO
         END DO
         !
cavazzon's avatar
cavazzon committed
222 223 224
      CASE ('fixed')

         if( nspin == 1 ) then
225 226 227
            nel(1)    = nint(nelec)
            nupdwn(1) = nbnd
            iupdwn(1) = 1
cavazzon's avatar
cavazzon committed
228 229 230 231
         else
            IF ( nelup + neldw /= nelec  ) THEN
               CALL errore(' electrons_base_initval ',' wrong # of up and down spin', 1 )
            END IF
232 233 234 235 236 237
            nel(1)    = nint(nelup)
            nel(2)    = nint(neldw)
            nupdwn(1) = nint(nelup)
            nupdwn(2) = nint(neldw)
            iupdwn(1) = 1
            iupdwn(2) = nupdwn(1) + 1
cavazzon's avatar
cavazzon committed
238 239
         end if

240 241 242
!         if( (nspin == 1) .and. MOD( nint(nelec), 2 ) /= 0 ) &
!              CALL errore(' electrons_base_initval ', &
!              ' must use nspin=2 for odd number of electrons', 1 )
243
         
cavazzon's avatar
cavazzon committed
244
         ! ocp = 2 for spinless systems, ocp = 1 for spin-polarized systems
245
         ocp = 2.0_DP / nspin
246
         !
cavazzon's avatar
cavazzon committed
247 248 249
         ! default filling: attribute ocp electrons to each states
         !                  until the good number of electrons is reached
         do iss = 1, nspin
250
            fsum = 0.0_DP
cavazzon's avatar
cavazzon committed
251
            do in = iupdwn ( iss ), iupdwn ( iss ) - 1 + nupdwn ( iss )
252
               if ( fsum + ocp < nel ( iss ) + 0.0001_DP ) then
cavazzon's avatar
cavazzon committed
253 254
                  f (in) = ocp
               else
255
                  f (in) = max( nel ( iss ) - fsum, 0.0_DP )
cavazzon's avatar
cavazzon committed
256
               end if
257
                fsum = fsum + f(in)
cavazzon's avatar
cavazzon committed
258 259
            end do
         end do
260
         !
cavazzon's avatar
cavazzon committed
261 262 263
      CASE ('ensemble','ensemble-dft','edft')

          if ( nspin == 1 ) then
264 265 266 267 268
            !
            f ( : )    = nelec / nbnd
            nel (1)    = nint(nelec)
            nupdwn (1) = nbnd
            !
cavazzon's avatar
cavazzon committed
269
          else
270
            !
cavazzon's avatar
cavazzon committed
271 272 273 274 275 276 277 278 279 280
            if (nelup.ne.0) then
              if ((nelup+neldw).ne.nelec) then
                 CALL errore(' electrons_base_initval ',' nelup+neldw .ne. nelec', 1 )
              end if
              nel (1) = nelup
              nel (2) = neldw
            else
              nel (1) = ( nint(nelec) + 1 ) / 2
              nel (2) =   nint(nelec)       / 2
            end if
281
            !
cavazzon's avatar
cavazzon committed
282 283 284
            nupdwn (1) = nbnd
            nupdwn (2) = nbnd
            iupdwn (2) = nbnd+1
285
            !
cavazzon's avatar
cavazzon committed
286 287
            do iss = 1, nspin
             do i = iupdwn ( iss ), iupdwn ( iss ) - 1 + nupdwn ( iss )
288
                f (i) =  nel (iss) / DBLE (nupdwn (iss))
cavazzon's avatar
cavazzon committed
289 290
             end do
            end do
291
            !
cavazzon's avatar
cavazzon committed
292 293 294 295 296 297 298 299 300
          end if

      CASE DEFAULT
         CALL errore(' electrons_base_initval ',' occupation method not implemented', 1 )
      END SELECT


      do iss = 1, nspin
         do in = iupdwn(iss), iupdwn(iss) - 1 + nupdwn(iss)
301
            ispin(in) = iss
cavazzon's avatar
cavazzon committed
302 303 304
         end do
      end do

305 306
      nbndx = nupdwn (1)
      nudx  = nupdwn (1)
307
      nbsp  = nupdwn (1) + nupdwn (2)
cavazzon's avatar
cavazzon committed
308 309 310

      IF ( nspin == 1 ) THEN 
        nelt = nel(1)
cavazzon's avatar
cavazzon committed
311
      ELSE
cavazzon's avatar
cavazzon committed
312
        nelt = nel(1) + nel(2)
cavazzon's avatar
cavazzon committed
313
      END IF
cavazzon's avatar
cavazzon committed
314

315
      IF( nupdwn(1) < nupdwn(2) ) &
giannozz's avatar
giannozz committed
316
        CALL errore(' electrons_base_initval ',' nupdwn(1) should be greater or equal nupdwn(2) ', 1 )
317 318 319 320 321 322

      IF( nbnd < nupdwn(1) ) &
        CALL errore(' electrons_base_initval ',' inconsistent nbnd, should be .GE. than  nupdwn(1) ', 1 )

      IF( nbspx < ( nupdwn(1) * nspin ) ) &
        CALL errore(' electrons_base_initval ',' inconsistent nbspx, should be .GE. than  nspin * nupdwn(1) ', 1 )
cavazzon's avatar
cavazzon committed
323 324 325 326

      IF( ( 2 * nbnd ) < nelt ) &
        CALL errore(' electrons_base_initval ',' too few states ',  1  )

327
      IF( nbsp < INT( nelec * nspin / 2.0_DP ) ) &
328 329
        CALL errore(' electrons_base_initval ',' too many electrons ', 1 )

cavazzon's avatar
cavazzon committed
330 331
      telectrons_base_initval = .TRUE.

cavazzon's avatar
cavazzon committed
332
      RETURN
cavazzon's avatar
cavazzon committed
333

334
    END SUBROUTINE electrons_base_initval
cavazzon's avatar
cavazzon committed
335

336
!----------------------------------------------------------------------------
337
!
338
    subroutine set_nelup_neldw ( tot_magnetization_, nelec_, nelup_, neldw_ )
339
      !
340 341
      USE kinds,     ONLY : DP
      USE constants, ONLY : eps8
342
      !
343 344 345
      REAL (KIND=DP), intent(IN)  :: tot_magnetization_
      REAL (KIND=DP), intent(IN)  :: nelec_
      REAL (KIND=DP), intent(OUT) :: nelup_, neldw_
346
      LOGICAL :: integer_charge, integer_magnetization
347
      !
348
      integer_charge = ( ABS (nelec_ - NINT(nelec_)) < eps8 )
349
      !
350
      IF ( tot_magnetization_ < 0 ) THEN
351
         ! default when tot_magnetization is unspecified
352 353 354 355 356
         IF ( integer_charge) THEN
            nelup_ = INT( nelec_ + 1 ) / 2
            neldw_ = nelec_ - nelup_
         ELSE
            nelup_ = nelec_ / 2
357
            neldw_ = nelup_
358
         END IF
359
      ELSE
360 361 362 363 364
         ! tot_magnetization specified in input
         !
         if ( (tot_magnetization_ > 0) .and. (nspin==1) ) &
                 CALL errore(' set_nelup_neldw  ', &
                 'tot_magnetization is inconsistent with nspin=1 ', 2 )
365 366 367
         integer_magnetization = ( ABS( tot_magnetization_ - &
                                   NINT(tot_magnetization_) ) < eps8 )
         IF ( integer_charge .AND. integer_magnetization) THEN
368 369 370
            !
            ! odd  tot_magnetization requires an odd  number of electrons
            ! even tot_magnetization requires an even number of electrons
371
            !
372 373 374 375
            if ( ((MOD(NINT(tot_magnetization_),2) == 0) .and. &
                  (MOD(NINT(nelec_),2)==1))               .or. &
                 ((MOD(NINT(tot_magnetization_),2) == 1) .and. &
                  (MOD(NINT(nelec_),2)==0))      ) &
376 377
              CALL infomsg(' set_nelup_neldw ',                          &
             'BEWARE: non-integer number of up and down electrons!' )
378 379 380 381 382 383 384 385 386 387
            !
            ! ... setting nelup/neldw
            !
            nelup_ = ( INT(nelec_) + tot_magnetization_ ) / 2
            neldw_ = ( INT(nelec_) - tot_magnetization_ ) / 2
         ELSE
            !
            nelup_ = ( nelec_ + tot_magnetization_ ) / 2
            neldw_ = ( nelec_ - tot_magnetization_ ) / 2
         END IF
388 389 390 391 392 393
      END IF

      return
    end subroutine set_nelup_neldw

!----------------------------------------------------------------------------
cavazzon's avatar
cavazzon committed
394 395


cavazzon's avatar
cavazzon committed
396 397
    SUBROUTINE deallocate_elct()
      IF( ALLOCATED( f ) ) DEALLOCATE( f )
398
      IF( ALLOCATED( ispin ) ) DEALLOCATE( ispin )
399 400
      IF( ALLOCATED( f_bgrp ) ) DEALLOCATE( f_bgrp )
      IF( ALLOCATED( ispin_bgrp ) ) DEALLOCATE( ispin_bgrp )
401
      IF( ALLOCATED( ibgrp_g2l ) ) DEALLOCATE( ibgrp_g2l )
cavazzon's avatar
cavazzon committed
402
      telectrons_base_initval = .FALSE.
cavazzon's avatar
cavazzon committed
403
      RETURN
404
    END SUBROUTINE deallocate_elct
cavazzon's avatar
cavazzon committed
405

406 407 408 409 410
!----------------------------------------------------------------------------

    SUBROUTINE distribute_bands( nbgrp, my_bgrp_id )
      INTEGER, INTENT(IN) :: nbgrp, my_bgrp_id
      INTEGER, EXTERNAL :: ldim_block, gind_block
411
      INTEGER :: iss, n1, n2, m1, m2, ilocal, iglobal
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 438
      !
      IF( .NOT. telectrons_base_initval ) &
        CALL errore( ' distribute_bands ', ' electrons_base_initval not yet called ', 1 )

      nupdwn_bgrp  = nupdwn
      iupdwn_bgrp  = iupdwn
      nudx_bgrp    = nudx
      nbsp_bgrp    = nbsp
      nbspx_bgrp   = nbspx
      i2gupdwn_bgrp= 1

      DO iss = 1, nspin
         nupdwn_bgrp( iss )  = ldim_block( nupdwn( iss ), nbgrp, my_bgrp_id )
         i2gupdwn_bgrp( iss ) = gind_block( 1, nupdwn( iss ), nbgrp, my_bgrp_id )
      END DO
      !
      iupdwn_bgrp(1) = 1
      IF( nspin > 1 ) THEN
         iupdwn_bgrp(2) = iupdwn_bgrp(1) + nupdwn_bgrp( 1 )
      END IF
      nudx_bgrp = nupdwn_bgrp( 1 )
      nbsp_bgrp = nupdwn_bgrp( 1 ) + nupdwn_bgrp ( 2 )
      nbspx_bgrp = nbsp_bgrp
      IF( MOD( nbspx_bgrp, 2 ) /= 0 ) nbspx_bgrp = nbspx_bgrp + 1

      ALLOCATE( f_bgrp     ( nbspx_bgrp ) )
      ALLOCATE( ispin_bgrp ( nbspx_bgrp ) )
439
      ALLOCATE( ibgrp_g2l ( nbspx ) )
440 441
      f_bgrp = 0.0
      ispin_bgrp = 0
442
      ibgrp_g2l = 0
443 444 445 446 447 448 449 450
      !
      DO iss = 1, nspin
         n1 = iupdwn_bgrp(iss)
         n2 = n1 + nupdwn_bgrp(iss) - 1
         m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1
         m2 = m1 + nupdwn_bgrp(iss) - 1
         f_bgrp(n1:n2) = f(m1:m2)
         ispin_bgrp(n1:n2) = ispin(m1:m2)
451 452 453 454 455
         ilocal = n1
         DO iglobal = m1, m2
            ibgrp_g2l( iglobal ) = ilocal
            ilocal = ilocal + 1
         END DO
456 457 458 459 460 461
      END DO
      
      RETURN

    END SUBROUTINE distribute_bands

462 463 464 465

!------------------------------------------------------------------------------!
  END MODULE electrons_base
!------------------------------------------------------------------------------!