1! 2! Copyright (C) 2002-2009 Quantum ESPRESSO group 3! This file is distributed under the terms of the 4! GNU General Public License. See the file `License' 5! in the root directory of the present distribution, 6! or http://www.gnu.org/copyleft/gpl.txt . 7! 8!------------------------------------------------------------------------------! 9 MODULE electrons_base 10!------------------------------------------------------------------------------! 11 12 USE kinds, ONLY: DP 13! 14 IMPLICIT NONE 15 SAVE 16 17 INTEGER :: nbnd = 0 ! number electronic bands, each band contains 18 ! two spin states 19 INTEGER :: nbndx = 0 ! array dimension nbndx >= nbnd 20 INTEGER :: nspin = 0 ! nspin = number of spins (1=no spin, 2=LSDA) 21 INTEGER :: nel(2) = 0 ! number of electrons (up, down) 22 INTEGER :: nelt = 0 ! total number of electrons ( up + down ) 23 INTEGER :: nupdwn(2) = 0 ! number of states with spin up (1) and down (2) 24 INTEGER :: iupdwn(2) = 0 ! first state with spin (1) and down (2) 25 INTEGER :: nudx = 0 ! max (nupdw(1),nupdw(2)) 26 INTEGER :: nbsp = 0 ! total number of electronic states 27 ! (nupdwn(1)+nupdwn(2)) 28 INTEGER :: nbspx = 0 ! array dimension nbspx >= nbsp 29 ! 30 INTEGER :: nupdwn_bgrp(2) = 0 ! number of states with spin up (1) and down (2) in this band group 31 INTEGER :: iupdwn_bgrp(2) = 0 ! first state with spin (1) and down (2) in this band group 32 INTEGER :: nudx_bgrp = 0 ! max (nupdw_bgrp(1),nupdw_bgrp(2)) in this band group 33 INTEGER :: nbsp_bgrp = 0 ! total number of electronic states 34 ! (nupdwn_bgrp(1)+nupdwn_bgrp(2)) in this band group 35 INTEGER :: nbspx_bgrp = 0 ! array dimension nbspx_bgrp >= nbsp_bgrp local to the band group 36 INTEGER :: i2gupdwn_bgrp(2)= 0 ! global index of the first local band 37 38 LOGICAL :: telectrons_base_initval = .FALSE. 39 LOGICAL :: keep_occ = .FALSE. ! if .true. when reading restart file keep 40 ! the occupations calculated in initval 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 REAL(DP), ALLOCATABLE :: f_bgrp(:) ! occupation numbers ( at gamma ) 47 INTEGER, ALLOCATABLE :: ispin_bgrp(:) ! spin of each state 48 INTEGER, ALLOCATABLE :: ibgrp_g2l(:) ! local index of the i-th global band index 49! 50!------------------------------------------------------------------------------! 51 CONTAINS 52!------------------------------------------------------------------------------! 53 54 55 SUBROUTINE electrons_base_initval( zv_ , na_ , nsp_ , nbnd_ , nspin_ , & 56 occupations_ , f_inp, tot_charge_, tot_magnetization_ ) 57 58 USE constants, ONLY : eps8 59 USE io_global, ONLY : stdout 60 61 REAL(DP), INTENT(IN) :: zv_ (:), tot_charge_ 62 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_ 66 CHARACTER(LEN=*), INTENT(IN) :: occupations_ 67 68 REAL(DP) :: nelec, nelup, neldw, ocp, fsum 69 INTEGER :: iss, i, in 70 71 nspin = nspin_ 72 ! 73 ! ... set nelec 74 ! 75 nelec = 0.0_DP 76 DO i = 1, nsp_ 77 nelec = nelec + na_ ( i ) * zv_ ( i ) 78 END DO 79 nelec = nelec - tot_charge_ 80 ! 81 ! ... set nelup/neldw 82 ! 83 nelup = 0._dp 84 neldw = 0._dp 85 call set_nelup_neldw (tot_magnetization_, nelec, nelup, neldw ) 86 87 IF( ABS( nelec - ( nelup + neldw ) ) > eps8 ) THEN 88 CALL errore(' electrons_base_initval ',' inconsistent n. of electrons ', 2 ) 89 END IF 90 ! 91 ! Compute the number of bands 92 ! 93 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 END IF 98 99 100 IF( nelec < 1 ) THEN 101 CALL errore(' electrons_base_initval ',' nelec less than 1 ', 1 ) 102 END IF 103 ! 104 IF( ABS( NINT( nelec ) - nelec ) > eps8 ) THEN 105 CALL errore(' electrons_base_initval ',' nelec must be integer', 2 ) 106 END IF 107 ! 108 IF( nbnd < 1 ) & 109 CALL errore(' electrons_base_initval ',' nbnd out of range ', 1 ) 110 ! 111 112 IF ( nspin /= 1 .AND. nspin /= 2 ) THEN 113 WRITE( stdout, * ) 'nspin = ', nspin 114 CALL errore( ' electrons_base_initval ', ' nspin out of range ', 1 ) 115 END IF 116 117 IF( MOD( nbnd, 2 ) == 0 ) THEN 118 nbspx = nbnd * nspin 119 ELSE 120 nbspx = ( nbnd + 1 ) * nspin 121 END IF 122 123 ALLOCATE( f ( nbspx ) ) 124 ALLOCATE( ispin ( nbspx ) ) 125 f = 0.0_DP 126 ispin = 0 127 128 iupdwn ( 1 ) = 1 129 nel = 0 130 131 SELECT CASE ( TRIM(occupations_) ) 132 CASE ('bogus') 133 ! 134 ! bogus to ensure \sum_i f_i = Nelec (nelec is integer) 135 ! 136 f ( : ) = nelec / nbspx 137 nel (1) = nint( nelec ) 138 nupdwn (1) = nbspx 139 if ( nspin == 2 ) then 140 ! 141 ! bogus to ensure Nelec = Nup + Ndw 142 ! 143 nel (1) = ( nint(nelec) + 1 ) / 2 144 nel (2) = nint(nelec) / 2 145 nupdwn (1)=nbnd 146 nupdwn (2)=nbnd 147 iupdwn (2)=nbnd+1 148 end if 149 ! 150 keep_occ = .true. 151 ! 152 CASE ('from_input') 153 ! 154 ! occupancies have been read from input 155 ! 156 ! count electrons 157 ! 158 IF( nspin == 1 ) THEN 159 nelec = SUM( f_inp( :, 1 ) ) 160 nelup = nelec / 2.0_DP 161 neldw = nelec / 2.0_DP 162 ELSE 163 nelup = SUM ( f_inp ( :, 1 ) ) 164 neldw = SUM ( f_inp ( :, 2 ) ) 165 nelec = nelup + neldw 166 END IF 167 ! 168 ! consistency check 169 ! 170 IF( nspin == 1 ) THEN 171 IF( f_inp( 1, 1 ) <= 0.0_DP ) & 172 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 173 ELSE 174 IF( f_inp( 1, 1 ) < 0.0_DP ) & 175 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 176 IF( f_inp( 1, 2 ) < 0.0_DP ) & 177 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 178 IF( ( f_inp( 1, 1 ) + f_inp( 1, 2 ) ) == 0.0_DP ) & 179 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 180 END IF 181 DO i = 2, nbnd 182 IF( nspin == 1 ) THEN 183 IF( f_inp( i, 1 ) > 0.0_DP .AND. f_inp( i-1, 1 ) <= 0.0_DP ) & 184 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 185 ELSE 186 IF( f_inp( i, 1 ) > 0.0_DP .AND. f_inp( i-1, 1 ) <= 0.0_DP ) & 187 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 188 IF( f_inp( i, 2 ) > 0.0_DP .AND. f_inp( i-1, 2 ) <= 0.0_DP ) & 189 CALL errore(' electrons_base_initval ',' Zero or negative occupation are not allowed ', 1 ) 190 END IF 191 END DO 192 ! 193 ! count bands 194 ! 195 nupdwn (1) = 0 196 nupdwn (2) = 0 197 DO i = 1, nbnd 198 IF( nspin == 1 ) THEN 199 IF( f_inp( i, 1 ) > 0.0_DP ) nupdwn (1) = nupdwn (1) + 1 200 ELSE 201 IF( f_inp( i, 1 ) > 0.0_DP ) nupdwn (1) = nupdwn (1) + 1 202 IF( f_inp( i, 2 ) > 0.0_DP ) nupdwn (2) = nupdwn (2) + 1 203 END IF 204 END DO 205 ! 206 if( nspin == 1 ) then 207 nel (1) = nint( nelec ) 208 iupdwn (1) = 1 209 else 210 nel (1) = nint(nelup) 211 nel (2) = nint(neldw) 212 iupdwn (1) = 1 213 iupdwn (2) = nupdwn (1) + 1 214 end if 215 ! 216 DO iss = 1, nspin 217 DO in = iupdwn ( iss ), iupdwn ( iss ) - 1 + nupdwn ( iss ) 218 f( in ) = f_inp( in - iupdwn ( iss ) + 1, iss ) 219 END DO 220 END DO 221 ! 222 CASE ('fixed') 223 224 if( nspin == 1 ) then 225 nel(1) = nint(nelec) 226 nupdwn(1) = nbnd 227 iupdwn(1) = 1 228 else 229 IF ( nelup + neldw /= nelec ) THEN 230 CALL errore(' electrons_base_initval ',' wrong # of up and down spin', 1 ) 231 END IF 232 nel(1) = nint(nelup) 233 nel(2) = nint(neldw) 234 nupdwn(1) = nint(nelup) 235 nupdwn(2) = nint(neldw) 236 iupdwn(1) = 1 237 iupdwn(2) = nupdwn(1) + 1 238 end if 239 240! if( (nspin == 1) .and. MOD( nint(nelec), 2 ) /= 0 ) & 241! CALL errore(' electrons_base_initval ', & 242! ' must use nspin=2 for odd number of electrons', 1 ) 243 244 ! ocp = 2 for spinless systems, ocp = 1 for spin-polarized systems 245 ocp = 2.0_DP / nspin 246 ! 247 ! default filling: attribute ocp electrons to each states 248 ! until the good number of electrons is reached 249 do iss = 1, nspin 250 fsum = 0.0_DP 251 do in = iupdwn ( iss ), iupdwn ( iss ) - 1 + nupdwn ( iss ) 252 if ( fsum + ocp < nel ( iss ) + 0.0001_DP ) then 253 f (in) = ocp 254 else 255 f (in) = max( nel ( iss ) - fsum, 0.0_DP ) 256 end if 257 fsum = fsum + f(in) 258 end do 259 end do 260 ! 261 CASE ('ensemble','ensemble-dft','edft') 262 263 if ( nspin == 1 ) then 264 ! 265 f ( : ) = nelec / nbnd 266 nel (1) = nint(nelec) 267 nupdwn (1) = nbnd 268 ! 269 else 270 ! 271 if (nelup.ne.0) then 272 if ((nelup+neldw).ne.nelec) then 273 CALL errore(' electrons_base_initval ',' nelup+neldw .ne. nelec', 1 ) 274 end if 275 nel (1) = nelup 276 nel (2) = neldw 277 else 278 nel (1) = ( nint(nelec) + 1 ) / 2 279 nel (2) = nint(nelec) / 2 280 end if 281 ! 282 nupdwn (1) = nbnd 283 nupdwn (2) = nbnd 284 iupdwn (2) = nbnd+1 285 ! 286 do iss = 1, nspin 287 do i = iupdwn ( iss ), iupdwn ( iss ) - 1 + nupdwn ( iss ) 288 f (i) = nel (iss) / DBLE (nupdwn (iss)) 289 end do 290 end do 291 ! 292 end if 293 294 CASE DEFAULT 295 CALL errore(' electrons_base_initval ',' occupation method not implemented', 1 ) 296 END SELECT 297 298 299 do iss = 1, nspin 300 do in = iupdwn(iss), iupdwn(iss) - 1 + nupdwn(iss) 301 ispin(in) = iss 302 end do 303 end do 304 305 nbndx = nupdwn (1) 306 nudx = nupdwn (1) 307 nbsp = nupdwn (1) + nupdwn (2) 308 309 IF ( nspin == 1 ) THEN 310 nelt = nel(1) 311 ELSE 312 nelt = nel(1) + nel(2) 313 END IF 314 315 IF( nupdwn(1) < nupdwn(2) ) & 316 CALL errore(' electrons_base_initval ',' nupdwn(1) should be greater or equal nupdwn(2) ', 1 ) 317 318 IF( nbnd < nupdwn(1) ) & 319 CALL errore(' electrons_base_initval ',' inconsistent nbnd, should be .GE. than nupdwn(1) ', 1 ) 320 321 IF( nbspx < ( nupdwn(1) * nspin ) ) & 322 CALL errore(' electrons_base_initval ',' inconsistent nbspx, should be .GE. than nspin * nupdwn(1) ', 1 ) 323 324 IF( ( 2 * nbnd ) < nelt ) & 325 CALL errore(' electrons_base_initval ',' too few states ', 1 ) 326 327 IF( nbsp < INT( nelec * nspin / 2.0_DP ) ) & 328 CALL errore(' electrons_base_initval ',' too many electrons ', 1 ) 329 330 telectrons_base_initval = .TRUE. 331 332 RETURN 333 334 END SUBROUTINE electrons_base_initval 335 336!---------------------------------------------------------------------------- 337! 338 subroutine set_nelup_neldw ( tot_magnetization_, nelec_, nelup_, neldw_ ) 339 ! 340 USE kinds, ONLY : DP 341 USE constants, ONLY : eps8 342 ! 343 REAL (KIND=DP), intent(IN) :: tot_magnetization_ 344 REAL (KIND=DP), intent(IN) :: nelec_ 345 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 IF ( integer_charge) THEN 353 nelup_ = INT( nelec_ + 1 ) / 2 354 neldw_ = nelec_ - nelup_ 355 ELSE 356 nelup_ = nelec_ / 2 357 neldw_ = nelup_ 358 END IF 359 ELSE 360 ! tot_magnetization specified in input 361 ! 362 if ( (tot_magnetization_ > 0) .and. (nspin==1) ) & 363 CALL errore(' set_nelup_neldw ', & 364 'tot_magnetization is inconsistent with nspin=1 ', 2 ) 365 integer_magnetization = ( ABS( tot_magnetization_ - & 366 NINT(tot_magnetization_) ) < eps8 ) 367 IF ( integer_charge .AND. integer_magnetization) THEN 368 ! 369 ! odd tot_magnetization requires an odd number of electrons 370 ! even tot_magnetization requires an even number of electrons 371 ! 372 if ( ((MOD(NINT(tot_magnetization_),2) == 0) .and. & 373 (MOD(NINT(nelec_),2)==1)) .or. & 374 ((MOD(NINT(tot_magnetization_),2) == 1) .and. & 375 (MOD(NINT(nelec_),2)==0)) ) & 376 CALL infomsg(' set_nelup_neldw ', & 377 'BEWARE: non-integer number of up and down electrons!' ) 378 ! 379 ! ... setting nelup/neldw 380 ! 381 nelup_ = ( INT(nelec_) + tot_magnetization_ ) / 2 382 neldw_ = ( INT(nelec_) - tot_magnetization_ ) / 2 383 ELSE 384 ! 385 nelup_ = ( nelec_ + tot_magnetization_ ) / 2 386 neldw_ = ( nelec_ - tot_magnetization_ ) / 2 387 END IF 388 END IF 389 390 return 391 end subroutine set_nelup_neldw 392 393!---------------------------------------------------------------------------- 394 395 396 SUBROUTINE deallocate_elct() 397 IF( ALLOCATED( f ) ) DEALLOCATE( f ) 398 IF( ALLOCATED( ispin ) ) DEALLOCATE( ispin ) 399 IF( ALLOCATED( f_bgrp ) ) DEALLOCATE( f_bgrp ) 400 IF( ALLOCATED( ispin_bgrp ) ) DEALLOCATE( ispin_bgrp ) 401 IF( ALLOCATED( ibgrp_g2l ) ) DEALLOCATE( ibgrp_g2l ) 402 telectrons_base_initval = .FALSE. 403 RETURN 404 END SUBROUTINE deallocate_elct 405 406!---------------------------------------------------------------------------- 407 408 SUBROUTINE distribute_bands( nbgrp, my_bgrp_id ) 409 INTEGER, INTENT(IN) :: nbgrp, my_bgrp_id 410 INTEGER, EXTERNAL :: ldim_block, gind_block 411 INTEGER :: iss, n1, n2, m1, m2, ilocal, iglobal 412 ! 413 IF( .NOT. telectrons_base_initval ) & 414 CALL errore( ' distribute_bands ', ' electrons_base_initval not yet called ', 1 ) 415 416 nupdwn_bgrp = nupdwn 417 iupdwn_bgrp = iupdwn 418 nudx_bgrp = nudx 419 nbsp_bgrp = nbsp 420 nbspx_bgrp = nbspx 421 i2gupdwn_bgrp= 1 422 423 DO iss = 1, nspin 424 nupdwn_bgrp( iss ) = ldim_block( nupdwn( iss ), nbgrp, my_bgrp_id ) 425 i2gupdwn_bgrp( iss ) = gind_block( 1, nupdwn( iss ), nbgrp, my_bgrp_id ) 426 END DO 427 ! 428 iupdwn_bgrp(1) = 1 429 IF( nspin > 1 ) THEN 430 iupdwn_bgrp(2) = iupdwn_bgrp(1) + nupdwn_bgrp( 1 ) 431 END IF 432 nudx_bgrp = nupdwn_bgrp( 1 ) 433 nbsp_bgrp = nupdwn_bgrp( 1 ) + nupdwn_bgrp ( 2 ) 434 nbspx_bgrp = nbsp_bgrp 435 IF( MOD( nbspx_bgrp, 2 ) /= 0 ) nbspx_bgrp = nbspx_bgrp + 1 436 437 ALLOCATE( f_bgrp ( nbspx_bgrp ) ) 438 ALLOCATE( ispin_bgrp ( nbspx_bgrp ) ) 439 ALLOCATE( ibgrp_g2l ( nbspx ) ) 440 f_bgrp = 0.0 441 ispin_bgrp = 0 442 ibgrp_g2l = 0 443 ! 444 DO iss = 1, nspin 445 n1 = iupdwn_bgrp(iss) 446 n2 = n1 + nupdwn_bgrp(iss) - 1 447 m1 = iupdwn(iss)+i2gupdwn_bgrp(iss) - 1 448 m2 = m1 + nupdwn_bgrp(iss) - 1 449 f_bgrp(n1:n2) = f(m1:m2) 450 ispin_bgrp(n1:n2) = ispin(m1:m2) 451 ilocal = n1 452 DO iglobal = m1, m2 453 ibgrp_g2l( iglobal ) = ilocal 454 ilocal = ilocal + 1 455 END DO 456 END DO 457 458 RETURN 459 460 END SUBROUTINE distribute_bands 461 462 463!------------------------------------------------------------------------------! 464 END MODULE electrons_base 465!------------------------------------------------------------------------------! 466