1! 2! Copyright (C) 2001-2020 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!---------------------------------------------------------------------------- 10SUBROUTINE setup() 11 !---------------------------------------------------------------------------- 12 !! This routine is called at the beginning of the calculation and: 13 ! 14 !! 1) determines various parameters of the calculation: 15 ! 16 !! * zv: charge of each atomic type; 17 !! * nelec: total number of electrons (if not given in input); 18 !! * nbnd: total number of bands (if not given in input); 19 !! * nbndx: max number of bands used in iterative diagonalization; 20 !! * tpiba: 2 pi / a (a = lattice parameter); 21 !! * tpiba2: square of tpiba; 22 !! * gcutm: cut-off in g space for charge/potentials; 23 !! * gcutms: cut-off in g space for smooth charge; 24 !! * ethr: convergence threshold for iterative diagonalization; 25 ! 26 !! 2) finds actual crystal symmetry: 27 ! 28 !! * s: symmetry matrices in the direct lattice vectors basis; 29 !! * nsym: number of crystal symmetry operations; 30 !! * nrot: number of lattice symmetry operations; 31 !! * ft: fractionary translations; 32 !! * irt: for each atom gives the corresponding symmetric; 33 !! * invsym: if true the system has inversion symmetry; 34 ! 35 !! 3) generates k-points corresponding to the actual crystal symmetry; 36 ! 37 !! 4) calculates various quantities used in magnetic, spin-orbit, PAW 38 !! electric-field, DFT+U(+V) calculations, and for parallelism. 39 ! 40 USE kinds, ONLY : DP 41 USE constants, ONLY : eps8, e2, fpi, pi, degspin 42 USE parameters, ONLY : npk 43 USE io_global, ONLY : stdout, ionode, ionode_id 44 USE io_files, ONLY : xmlfile 45 USE cell_base, ONLY : at, bg, alat, tpiba, tpiba2, ibrav 46 USE ions_base, ONLY : nat, tau, ntyp => nsp, ityp, zv 47 USE basis, ONLY : starting_pot, natomwfc 48 USE gvect, ONLY : gcutm, ecutrho 49 USE gvecw, ONLY : gcutw, ecutwfc 50 USE gvecs, ONLY : doublegrid, gcutms, dual 51 USE klist, ONLY : xk, wk, nks, nelec, degauss, lgauss, & 52 ltetra, lxkcry, nkstot, & 53 nelup, neldw, two_fermi_energies, & 54 tot_charge, tot_magnetization 55 USE ener, ONLY : ef, ef_up, ef_dw 56 USE electrons_base, ONLY : set_nelup_neldw 57 USE start_k, ONLY : nks_start, xk_start, wk_start, & 58 nk1, nk2, nk3, k1, k2, k3 59 USE ktetra, ONLY : tetra_type, opt_tetra_init, tetra_init 60 USE symm_base, ONLY : s, t_rev, irt, nrot, nsym, invsym, nosym, & 61 d1,d2,d3, time_reversal, sname, set_sym_bl, & 62 find_sym, inverse_s, no_t_rev, fft_fact, & 63 allfrac 64 USE wvfct, ONLY : nbnd, nbndx 65 USE control_flags, ONLY : tr2, ethr, lscf, lbfgs, lmd, david, lecrpa, & 66 isolve, niter, noinv, ts_vdw, & 67 lbands, use_para_diag, gamma_only, & 68 restart 69 USE cellmd, ONLY : calc 70 USE uspp_param, ONLY : upf, n_atom_wfc 71 USE uspp, ONLY : okvan 72 USE ldaU, ONLY : lda_plus_u, init_lda_plus_u 73 USE bp, ONLY : gdir, lberry, nppstr, lelfield, lorbm, nx_el,& 74 nppstr_3d,l3dstring, efield 75 USE fixed_occ, ONLY : f_inp, tfixed_occ, one_atom_occupations 76 USE mp_images, ONLY : intra_image_comm 77 USE mp_pools, ONLY : kunit 78 USE mp_bands, ONLY : intra_bgrp_comm, nyfft 79 USE mp, ONLY : mp_bcast 80 USE lsda_mod, ONLY : lsda, nspin, current_spin, isk, & 81 starting_magnetization 82 USE spin_orb, ONLY : lspinorb, domag 83 USE noncollin_module, ONLY : noncolin, npol, i_cons, m_loc, & 84 angle1, angle2, bfield, ux, nspin_lsda, & 85 nspin_gga, nspin_mag 86 USE qexsd_module, ONLY : qexsd_readschema 87 USE qexsd_copy, ONLY : qexsd_copy_efermi 88 USE qes_libs_module, ONLY : qes_reset 89 USE qes_types_module, ONLY : output_type 90 USE exx, ONLY : ecutfock 91 USE exx_base, ONLY : exx_grid_init, exx_mp_init, exx_div_check 92 USE funct, ONLY : dft_is_meta, dft_is_hybrid, dft_is_gradient 93 USE paw_variables, ONLY : okpaw 94 USE fcp_variables, ONLY : lfcpopt, lfcpdyn 95 USE extfield, ONLY : gate 96 USE additional_kpoints, ONLY : add_additional_kpoints 97 ! 98 IMPLICIT NONE 99 ! 100 INTEGER :: na, is, ierr, ibnd, ik, nrot_ 101 LOGICAL :: magnetic_sym, skip_equivalence=.FALSE. 102 REAL(DP) :: iocc, ionic_charge, one 103 ! 104 LOGICAL, EXTERNAL :: check_para_diag 105 ! 106 TYPE(output_type) :: output_obj 107 ! 108#if defined(__MPI) 109 LOGICAL :: lpara = .true. 110#else 111 LOGICAL :: lpara = .false. 112#endif 113 114 ! 115 ! ... okvan/okpaw = .TRUE. : at least one pseudopotential is US/PAW 116 ! 117 okvan = ANY( upf(1:ntyp)%tvanp ) 118 okpaw = ANY( upf(1:ntyp)%tpawp ) 119 ! 120 ! ... check for features not implemented with US-PP or PAW 121 ! 122 IF ( okvan .OR. okpaw ) THEN 123 IF ( dft_is_meta() ) CALL errore( 'setup', & 124 'Meta-GGA not implemented with USPP/PAW', 1 ) 125 IF ( noncolin .AND. lberry) CALL errore( 'setup', & 126 'Noncolinear Berry Phase/electric not implemented with USPP/PAW', 1 ) 127 IF (ts_vdw ) CALL errore ('setup',& 128 'Tkatchenko-Scheffler not implemented with USPP/PAW', 1 ) 129 IF ( lorbm ) CALL errore( 'setup', & 130 'Orbital Magnetization not implemented with USPP/PAW', 1 ) 131 END IF 132 133 IF ( dft_is_hybrid() ) THEN 134 IF ( lberry ) CALL errore( 'setup ', & 135 'hybrid XC not allowed in Berry-phase calculations',1 ) 136 IF ( lelfield ) CALL errore( 'setup ', & 137 'hybrid XC and electric fields untested',1 ) 138 IF ( allfrac ) CALL errore( 'setup ', & 139 'option use_all_frac incompatible with hybrid XC', 1 ) 140 IF (.NOT. lscf) CALL errore( 'setup ', & 141 'hybrid XC not allowed in non-scf calculations', 1 ) 142 IF ( ANY (upf(1:ntyp)%nlcc) ) CALL infomsg( 'setup ', 'BEWARE:' // & 143 & ' nonlinear core correction is not consistent with hybrid XC') 144 IF (okvan) THEN 145 IF (ecutfock /= 4.0_dp*ecutwfc) THEN 146 ecutfock = MIN(4.0_dp*ecutwfc,ecutrho) 147 CALL infomsg ('setup', & 148 'Warning: ecutfock not valid for US/PAW, ignored') 149 END IF 150 IF ( lmd .OR. lbfgs ) CALL errore & 151 ('setup','forces for hybrid functionals + US/PAW not implemented',1) 152 IF ( noncolin ) CALL errore & 153 ('setup','Noncolinear hybrid XC for USPP not implemented',1) 154 END IF 155 IF ( noncolin ) no_t_rev=.true. 156 END IF 157 ! 158 IF ( dft_is_meta() .AND. noncolin ) CALL errore( 'setup', & 159 'Non-collinear Meta-GGA not implemented', 1 ) 160 ! 161 ! ... Compute the ionic charge for each atom type and the total ionic charge 162 ! 163 zv(1:ntyp) = upf(1:ntyp)%zp 164 ! 165#if defined (__PGI) 166 ionic_charge = 0._DP 167 DO na = 1, nat 168 ionic_charge = ionic_charge + zv( ityp(na) ) 169 END DO 170#else 171 ionic_charge = SUM( zv(ityp(1:nat)) ) 172#endif 173 ! 174 ! ... set the number of electrons 175 ! 176 nelec = ionic_charge - tot_charge 177 ! 178 IF ( .NOT. lscf .OR. ( (lfcpopt .OR. lfcpdyn ) .AND. restart )) THEN 179 ! 180 ! ... in these cases, we need (or it is useful) to read the Fermi energy 181 ! 182 IF (ionode) CALL qexsd_readschema ( xmlfile(), ierr, output_obj ) 183 CALL mp_bcast(ierr, ionode_id, intra_image_comm) 184 IF ( ierr > 0 ) CALL errore ( 'setup', 'problem reading ef from file ' //& 185 & TRIM(xmlfile()), ierr ) 186 IF (ionode) CALL qexsd_copy_efermi ( output_obj%band_structure, & 187 nelec, ef, two_fermi_energies, ef_up, ef_dw ) 188 ! convert to Ry a.u. 189 ef = ef*e2 190 ef_up = ef_up*e2 191 ef_dw = ef_dw*e2 192 CALL mp_bcast(nelec, ionode_id, intra_image_comm) 193 CALL mp_bcast(ef, ionode_id, intra_image_comm) 194 CALL mp_bcast(two_fermi_energies, ionode_id, intra_image_comm) 195 CALL mp_bcast(ef_up, ionode_id, intra_image_comm) 196 CALL mp_bcast(ef_dw, ionode_id, intra_image_comm) 197 CALL qes_reset ( output_obj ) 198 ! 199 END IF 200 IF ( (lfcpopt .OR. lfcpdyn) .AND. restart ) THEN 201 tot_charge = ionic_charge - nelec 202 END IF 203 ! 204 ! ... magnetism-related quantities 205 ! 206 ! ... Set the domag variable to make a spin-orbit calculation with zero 207 ! ... magnetization 208 ! 209 IF ( noncolin ) THEN 210 domag = ANY ( ABS( starting_magnetization(1:ntyp) ) > 1.D-6 ) 211 ELSE 212 domag = .false. 213 END IF 214 ! 215 ! Set the different spin indices 216 ! 217 CALL set_spin_vars( lsda, noncolin, lspinorb, domag, & 218 npol, nspin, nspin_lsda, nspin_mag, nspin_gga, current_spin ) 219 ! 220 ! time reversal operation is set up to 0 by default 221 t_rev = 0 222 ! 223 ALLOCATE( m_loc( 3, nat ) ) 224 IF ( noncolin ) THEN 225 ! 226 ! gamma_only and noncollinear not allowed 227 ! 228 if (gamma_only) call errore('setup', & 229 'gamma_only and noncolin not allowed',1) 230 ! 231 DO na = 1, nat 232 ! 233 m_loc(1,na) = starting_magnetization(ityp(na)) * & 234 SIN( angle1(ityp(na)) ) * COS( angle2(ityp(na)) ) 235 m_loc(2,na) = starting_magnetization(ityp(na)) * & 236 SIN( angle1(ityp(na)) ) * SIN( angle2(ityp(na)) ) 237 m_loc(3,na) = starting_magnetization(ityp(na)) * & 238 COS( angle1(ityp(na)) ) 239 END DO 240 ! 241 ! initialize the quantization direction for gga 242 ! 243 ux=0.0_DP 244 if (dft_is_gradient()) call compute_ux(m_loc,ux,nat) 245 ! 246 ELSE 247 ! 248 IF (lspinorb) CALL errore( 'setup ', & 249 'spin orbit requires a non collinear calculation', 1 ) 250 ! 251 IF ( i_cons == 1) then 252 do na=1,nat 253 m_loc(1,na) = starting_magnetization(ityp(na)) 254 end do 255 end if 256 IF ( i_cons /= 0 .AND. nspin==1 ) & 257 CALL errore( 'setup', 'this i_cons requires a magnetic calculation ', 1 ) 258 IF ( i_cons /= 0 .AND. i_cons /= 1 ) & 259 CALL errore( 'setup', 'this i_cons requires a non colinear run', 1 ) 260 ! 261 END IF 262 ! 263 ! ... if this is not a spin-orbit calculation, all spin-orbit pseudopotentials 264 ! ... are transformed into standard pseudopotentials 265 ! 266 IF ( lspinorb ) THEN 267 IF ( ALL ( .NOT. upf(:)%has_so ) ) & 268 CALL infomsg ('setup','At least one non s.o. pseudo') 269 ELSE 270 CALL average_pp ( ntyp ) 271 END IF 272 ! 273 ! ... If the occupations are from input, check the consistency with the 274 ! ... number of electrons 275 ! 276 IF ( tfixed_occ ) THEN 277 ! 278 iocc = 0 279 ! 280 DO is = 1, nspin_lsda 281 ! 282#if defined (__PGI) 283 DO ibnd = 1, nbnd 284 iocc = iocc + f_inp(ibnd,is) 285 END DO 286#else 287 iocc = iocc + SUM( f_inp(1:nbnd,is) ) 288#endif 289 ! 290 DO ibnd = 1, nbnd 291 if (f_inp(ibnd,is) > 2.d0/nspin_lsda .or. f_inp(ibnd,is) < 0.d0) & 292 call errore('setup','wrong fixed occupations',is) 293 END DO 294 END DO 295 ! 296 IF ( ABS( iocc - nelec ) > 1D-5 ) & 297 CALL errore( 'setup', 'strange occupations: '//& 298 'number of electrons from occupations is wrong.', 1 ) 299 ! 300 END IF 301 ! 302 ! ... Check: if there is an odd number of electrons, the crystal is a metal 303 ! 304 IF ( lscf .AND. ABS( NINT( nelec / 2.D0 ) - nelec / 2.D0 ) > eps8 & 305 .AND. .NOT. lgauss .AND. .NOT. ltetra .AND. .NOT. tfixed_occ ) & 306 CALL infomsg( 'setup', 'the system is metallic, specify occupations' ) 307 ! 308 ! ... Check: spin-polarized calculations require either broadening or 309 ! fixed occupation 310 ! 311 IF ( lscf .AND. lsda & 312 .AND. .NOT. lgauss .AND. .NOT. ltetra & 313 .AND. .NOT. tfixed_occ .AND. .NOT. two_fermi_energies ) & 314 CALL errore( 'setup', 'spin-polarized system, specify occupations', 1 ) 315 ! 316 ! ... setting the number of up and down electrons 317 ! 318 call set_nelup_neldw ( tot_magnetization, nelec, nelup, neldw ) 319 ! 320 ! ... Set the number of occupied bands if not given in input 321 ! 322 IF ( nbnd == 0 ) THEN 323 ! 324 IF (nat==0) CALL errore('setup','free electrons: nbnd required in input',1) 325 ! 326 nbnd = MAX ( NINT( nelec / degspin ), NINT(nelup), NINT(neldw) ) 327 ! 328 IF ( lgauss .OR. ltetra ) THEN 329 ! 330 ! ... metallic case: add 20% more bands, with a minimum of 4 331 ! 332 nbnd = MAX( NINT( 1.2D0 * nelec / degspin ), & 333 NINT( 1.2D0 * nelup), NINT( 1.2d0 * neldw ), & 334 ( nbnd + 4 ) ) 335 ! 336 END IF 337 ! 338 ! ... In the case of noncollinear magnetism, bands are NOT 339 ! ... twofold degenerate : 340 ! 341 IF ( noncolin ) nbnd = INT( degspin ) * nbnd 342 ! 343 ELSE 344 ! 345 IF ( nbnd < NINT( nelec / degspin ) .AND. lscf ) & 346 CALL errore( 'setup', 'too few bands', 1 ) 347 ! 348 IF ( nbnd < NINT( nelup ) .AND. lscf ) & 349 CALL errore( 'setup', 'too few spin up bands', 1 ) 350 IF ( nbnd < NINT( neldw ) .AND. lscf ) & 351 CALL errore( 'setup', 'too few spin dw bands', 1 ) 352 ! 353 IF ( nbnd < NINT( nelec ) .AND. lscf .AND. noncolin ) & 354 CALL errore( 'setup', 'too few bands', 1 ) 355 ! 356 END IF 357 ! 358 ! ... Here we set the precision of the diagonalization for the first scf 359 ! ... iteration of for the first ionic step 360 ! ... for subsequent steps ethr is automatically updated in electrons 361 ! 362 IF ( nat==0 ) THEN 363 ! 364 ethr=1.0D-8 365 ! 366 ELSE IF ( .NOT. lscf ) THEN 367 ! 368 ! ... do not allow convergence threshold of scf and nscf to become too small 369 ! 370 IF ( ethr == 0.D0 ) ethr = MAX(1.D-13, 0.1D0 * MIN( 1.D-2, tr2 / nelec )) 371 ! 372 ELSE 373 ! 374 IF ( ethr == 0.D0 ) THEN 375 ! 376 IF ( starting_pot == 'file' ) THEN 377 ! 378 ! ... if you think that the starting potential is good 379 ! ... do not spoil it with a lousy first diagonalization : 380 ! ... set a strict ethr in the input file (diago_thr_init) 381 ! 382 ethr = 1.D-5 383 ! 384 ELSE 385 ! 386 ! ... starting atomic potential is probably far from scf 387 ! ... do not waste iterations in the first diagonalizations 388 ! 389 ethr = 1.0D-2 390 ! 391 END IF 392 ! 393 END IF 394 ! 395 END IF 396 ! 397 IF ( .NOT. lscf ) niter = 1 398 ! 399 ! ... set number of atomic wavefunctions 400 ! 401 natomwfc = n_atom_wfc( nat, ityp, noncolin ) 402 ! 403 ! ... set the max number of bands used in iterative diagonalization 404 ! 405 nbndx = nbnd 406 IF ( isolve == 0 ) nbndx = david * nbnd 407 ! 408 use_para_diag = check_para_diag( nbnd ) 409 ! 410 ! ... Set the units in real and reciprocal space 411 ! 412 tpiba = 2.D0 * pi / alat 413 tpiba2 = tpiba**2 414 ! 415 ! ... Compute the cut-off of the G vectors 416 ! 417 doublegrid = ( dual > 4.0_dp + eps8 ) 418 IF ( doublegrid .AND. ( .NOT.okvan .AND. .NOT.okpaw .AND. & 419 .NOT. ANY (upf(1:ntyp)%nlcc) ) ) & 420 CALL infomsg ( 'setup', 'no reason to have ecutrho>4*ecutwfc' ) 421 gcutm = dual * ecutwfc / tpiba2 422 gcutw = ecutwfc / tpiba2 423 ! 424 IF ( doublegrid ) THEN 425 ! 426 gcutms = 4.D0 * ecutwfc / tpiba2 427 ! 428 ELSE 429 ! 430 gcutms = gcutm 431 ! 432 END IF 433 ! 434 ! ... Test that atoms do not overlap 435 ! 436 call check_atoms ( nat, tau, bg ) 437 ! 438 ! ... generate transformation matrices for the crystal point group 439 ! ... First we generate all the symmetry matrices of the Bravais lattice 440 ! 441 call set_sym_bl ( ) 442 ! 443 ! ... If lecrpa is true, nosym must be set to true also 444 ! 445 IF ( lecrpa ) nosym = .TRUE. 446 IF ( lecrpa ) skip_equivalence=.TRUE. 447 ! 448 ! ... if nosym is true: do not reduce automatic k-point grids to the IBZ 449 ! ... using the symmetries of the lattice (only k <-> -k symmetry is used) 450 ! ... Does not change the number "nrot" of symmetries of the lattice 451 ! 452 IF ( nosym ) THEN 453 nrot_ = 1 454 ELSE 455 nrot_ = nrot 456 END IF 457 ! 458 ! ... time_reversal = use q=>-q symmetry for k-point generation 459 ! 460 magnetic_sym = noncolin .AND. domag 461 time_reversal = .NOT. noinv .AND. .NOT. magnetic_sym 462 ! 463 ! ... Automatic generation of k-points (if required) 464 ! 465 IF ( nks_start == 0 ) THEN 466 ! 467 IF (lelfield .OR. lorbm) THEN 468 ! 469 CALL kpoint_grid_efield (at,bg, npk, & 470 k1,k2,k3, nk1,nk2,nk3, nkstot, xk, wk, nspin) 471 nosym = .TRUE. 472 nrot_ = 1 473 ! 474 ELSE IF (lberry ) THEN 475 ! 476 CALL kp_strings( nppstr, gdir, nrot_, s, bg, npk, & 477 k1, k2, k3, nk1, nk2, nk3, nkstot, xk, wk ) 478 nosym = .TRUE. 479 nrot_ = 1 480 ! 481 ELSE 482 ! 483 CALL kpoint_grid ( nrot_,time_reversal, skip_equivalence, s, t_rev, bg,& 484 npk, k1,k2,k3, nk1,nk2,nk3, nkstot, xk, wk) 485 ! 486 END IF 487 ! 488 ELSE 489 nkstot = nks_start 490 xk(:,1:nkstot) = xk_start(:,1:nks_start) 491 wk(1:nkstot) = wk_start(1:nks_start) 492 ! 493 IF( lelfield) THEN 494 ! 495 IF(noncolin) THEN 496 allocate(nx_el(nkstot,3)) 497 ELSE 498 allocate(nx_el(nkstot*nspin,3)) 499 END IF 500 501 IF ( gdir<1 .OR. gdir>3 ) CALL errore('setup','invalid gdir value'& 502 &' (valid values: 1=x, 2=y, 3=z)',10) 503 DO ik=1,nkstot 504 nx_el(ik,gdir)=ik 505 END DO 506 ! sanity check (when nkstot==1 we /could/ just set nppstr=1): 507 IF(nppstr==0) CALL errore('setup', 'When lefield is true and kpoint are '& 508 &'specified manually you MUST set nppstr',1) 509 if(nspin==2) nx_el(nkstot+1:2*nkstot,:) = nx_el(1:nkstot,:) + nkstot 510 nppstr_3d(gdir)=nppstr 511 l3dstring=.false. 512 nosym = .TRUE. 513 nrot_ = 1 514 ! 515 END IF 516 END IF 517 ! 518 CALL add_additional_kpoints(nkstot, xk, wk) 519 ! 520 IF ( nat==0 ) THEN 521 ! 522 nsym=nrot 523 invsym=.true. 524 CALL inverse_s ( ) 525 ! 526 ELSE 527 ! 528 ! ... eliminate rotations that are not symmetry operations 529 ! 530 CALL find_sym ( nat, tau, ityp, magnetic_sym, m_loc, gate ) 531 ! 532 ! ... do not force FFT grid to be commensurate with fractional translations 533 ! 534 IF ( allfrac ) fft_fact(:) = 1 535 ! 536 END IF 537 ! 538 ! ... nosym: do not use any point-group symmetry (s(:,:,1) is the identity) 539 ! 540 IF ( nosym ) THEN 541 nsym = 1 542 invsym = .FALSE. 543 fft_fact(:) = 1 544 END IF 545 ! 546 IF ( nsym > 1 .AND. ibrav == 0 ) CALL infomsg('setup', & 547 'using ibrav=0 with symmetry is DISCOURAGED, use correct ibrav instead') 548 ! 549 ! ... Input k-points are assumed to be given in the IBZ of the Bravais 550 ! ... lattice, with the full point symmetry of the lattice. 551 ! ... If some symmetries of the lattice are missing in the crystal, 552 ! ... "irreducible_BZ" computes the missing k-points. 553 ! 554 IF ( .NOT. lbands ) THEN 555 CALL irreducible_BZ (nrot_, s, nsym, time_reversal, & 556 magnetic_sym, at, bg, npk, nkstot, xk, wk, t_rev) 557 ELSE 558 one = SUM (wk(1:nkstot)) 559 IF ( one > 0.0_dp ) wk(1:nkstot) = wk(1:nkstot) / one 560 END IF 561 ! 562 ! ... if dynamics is done the system should have no symmetries 563 ! ... (inversion symmetry alone is allowed) 564 ! 565 IF ( lmd .AND. ( nsym == 2 .AND. .NOT. invsym .OR. nsym > 2 ) & 566 .AND. .NOT. ( calc == 'mm' .OR. calc == 'nm' ) ) & 567 CALL infomsg( 'setup', 'Dynamics, you should have no symmetries' ) 568 ! 569 IF ( ltetra ) THEN 570 ! 571 ! ... Calculate quantities used in tetrahedra method 572 ! 573 IF (nks_start /= 0) CALL errore( 'setup ', 'tetrahedra need automatic k-point grid',1) 574 ! 575 IF (tetra_type == 0) then 576 CALL tetra_init( nsym, s, time_reversal, t_rev, at, bg, npk, k1,k2,k3, & 577 nk1, nk2, nk3, nkstot, xk ) 578 ELSE 579 CALL opt_tetra_init(nsym, s, time_reversal, t_rev, at, bg, npk, & 580 k1, k2, k3, nk1, nk2, nk3, nkstot, xk, 1) 581 END IF 582 ! 583 END IF 584 ! 585 IF ( lsda ) THEN 586 ! 587 ! ... LSDA case: two different spin polarizations, 588 ! ... each with its own kpoints 589 ! 590 if (nspin /= 2) call errore ('setup','nspin should be 2; check iosys',1) 591 ! 592 CALL set_kup_and_kdw( xk, wk, isk, nkstot, npk ) 593 ! 594 ELSE IF ( noncolin ) THEN 595 ! 596 ! ... noncolinear magnetism: potential and charge have dimension 4 (1+3) 597 ! 598 if (nspin /= 4) call errore ('setup','nspin should be 4; check iosys',1) 599 current_spin = 1 600 isk(:) = 1 601 ! 602 ELSE 603 ! 604 ! ... LDA case: the two spin polarizations are identical 605 ! 606 wk(1:nkstot) = wk(1:nkstot) * degspin 607 current_spin = 1 608 isk(:) = 1 609 ! 610 IF ( nspin /= 1 ) & 611 CALL errore( 'setup', 'nspin should be 1; check iosys', 1 ) 612 ! 613 END IF 614 ! 615 IF ( nkstot > npk ) CALL errore( 'setup', 'too many k points', nkstot ) 616 ! 617 ! ... distribute k-points (and their weights and spin indices) 618 ! 619 kunit = 1 620 CALL divide_et_impera ( nkstot, xk, wk, isk, nks ) 621 ! 622 IF ( dft_is_hybrid() ) THEN 623 IF ( nks == 0 ) CALL errore('setup','pools with no k-points not allowed for hybrid functionals',1) 624 CALL exx_grid_init() 625 CALL exx_mp_init() 626 CALL exx_div_check() 627 ENDIF 628 629 IF (one_atom_occupations) THEN 630 DO ik=1,nkstot 631 DO ibnd=natomwfc+1, nbnd 632 IF (f_inp(ibnd,ik)> 0.0_DP) CALL errore('setup', & 633 'no atomic wavefunction for some band',1) 634 ENDDO 635 ENDDO 636 ENDIF 637 ! 638 ! ... Set up Hubbard parameters for DFT+U(+V) calculation 639 ! 640 CALL init_lda_plus_u ( upf(1:ntyp)%psd, nspin, noncolin ) 641 ! 642 ! ... initialize d1 and d2 to rotate the spherical harmonics 643 ! 644 IF (lda_plus_u .or. okpaw .or. (okvan.and.dft_is_hybrid()) ) CALL d_matrix( d1, d2, d3 ) 645 ! 646 RETURN 647 ! 648END SUBROUTINE setup 649! 650!---------------------------------------------------------------------------- 651LOGICAL FUNCTION check_para_diag( nbnd ) 652 !----------------------------------------------------------------------------- 653 !! Some checks for parallel diagonalization. 654 ! 655 USE io_global, ONLY : stdout, ionode, ionode_id 656 USE mp_bands, ONLY : intra_bgrp_comm 657 USE mp_pools, ONLY : intra_pool_comm 658 659 IMPLICIT NONE 660 661 INCLUDE 'laxlib.fh' 662 663 INTEGER, INTENT(IN) :: nbnd 664 !! number of bands 665 ! 666 LOGICAL, SAVE :: first = .TRUE. 667 LOGICAL, SAVE :: saved_value = .FALSE. 668 INTEGER :: np_ortho(2), ortho_parent_comm 669 670#if defined(__MPI) 671 IF( .NOT. first ) THEN 672 check_para_diag = saved_value 673 RETURN 674 END IF 675 first = .FALSE. 676 ! 677 CALL laxlib_getval( np_ortho = np_ortho, ortho_parent_comm = ortho_parent_comm ) 678 ! 679 IF( np_ortho(1) > nbnd ) & 680 CALL errore ('check_para_diag', 'Too few bands for required ndiag',nbnd) 681 ! 682 check_para_diag = ( np_ortho( 1 ) > 1 .AND. np_ortho( 2 ) > 1 ) 683 saved_value = check_para_diag 684 ! 685 IF ( ionode ) THEN 686 ! 687 WRITE( stdout, '(/,5X,"Subspace diagonalization in iterative solution ",& 688 & "of the eigenvalue problem:")' ) 689 IF ( check_para_diag ) THEN 690 IF (ortho_parent_comm .EQ. intra_pool_comm) THEN 691 WRITE( stdout, '(5X,"one sub-group per k-point group (pool) will be used")' ) 692 ELSE IF (ortho_parent_comm .EQ. intra_bgrp_comm) THEN 693 WRITE( stdout, '(5X,"one sub-group per band group will be used")' ) 694 ELSE 695 CALL errore( 'setup','Unexpected sub-group communicator ', 1 ) 696 END IF 697#if defined(__ELPA) || defined(__ELPA_2015) || defined(__ELPA_2016) || defined(__ELPA_2017) || defined(__ELPA_2018) || defined(__ELPA_2019) 698 WRITE( stdout, '(5X,"ELPA distributed-memory algorithm ", & 699 & "(size of sub-group: ", I2, "*", I3, " procs)",/)') & 700 np_ortho(1), np_ortho(2) 701#elif defined(__SCALAPACK) 702 WRITE( stdout, '(5X,"scalapack distributed-memory algorithm ", & 703 & "(size of sub-group: ", I2, "*", I3, " procs)",/)') & 704 np_ortho(1), np_ortho(2) 705#else 706 WRITE( stdout, '(5X,"custom distributed-memory algorithm ", & 707 & "(size of sub-group: ", I2, "*", I3, " procs)",/)') & 708 np_ortho(1), np_ortho(2) 709#endif 710 ELSE 711 WRITE( stdout, '(5X,"a serial algorithm will be used",/)' ) 712 END IF 713 ! 714 END IF 715 ! 716#else 717 check_para_diag = .FALSE. 718#endif 719 RETURN 720END FUNCTION check_para_diag 721