! This file is copied and modified from QUANTUM ESPRESSO ! Kun Cao, Henry Lambert, Feliciano Giustino ! ! Copyright (C) 2008 Quantum ESPRESSO group ! 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 . ! !---------------------------------------------------------------------------- SUBROUTINE setup_nscf ( newgrid, xq ) !---------------------------------------------------------------------------- ! ! ... This routine initializes variables for the non-scf calculations at k ! ... and k+q required by the linear response calculation at finite q. ! ... In particular: finds the symmetry group of the crystal that leaves ! ... the phonon q-vector (xq) or the single atomic displacement (modenum) ! ... unchanged; determines the k- and k+q points in the irreducible BZ ! ... Needed on input (read from data file): ! ... "nsym" crystal symmetries s, ftau, t_rev, "nrot" lattice symetries "s" ! ... "nkstot" k-points in the irreducible BZ wrt lattice symmetry ! ... Produced on output: ! ... symmetries ordered with the "nsymq" phonon symmetries first ! ... "nkstot" k- and k+q-points in the IBZ calculated for the phonon sym.) ! ... Misc. data needed for running the non-scf calculation ! USE kinds, ONLY : DP USE constants, ONLY : eps8 USE parameters, ONLY : npk USE io_global, ONLY : stdout USE constants, ONLY : pi, degspin USE cell_base, ONLY : at, bg, alat, tpiba, tpiba2, ibrav, omega USE ions_base, ONLY : nat, tau, ntyp => nsp, ityp, zv USE force_mod, ONLY : force USE basis, ONLY : natomwfc USE klist, ONLY : xk, wk, nks, nelec, degauss, lgauss, & nkstot, qnorm USE lsda_mod, ONLY : lsda, nspin, current_spin, isk, & starting_magnetization USE symm_base, ONLY : s, t_rev, irt, ftau, nrot, nsym, & time_reversal, sname, d1, d2, d3, invsym, & copy_sym, s_axis_to_cart USE wvfct, ONLY : nbnd, nbndx, npwx USE control_flags, ONLY : ethr, isolve, david, max_cg_iter, & noinv, modenum, use_para_diag, io_level USE mp_global, ONLY : kunit USE spin_orb, ONLY : domag USE noncollin_module, ONLY : noncolin, npol USE start_k, ONLY : nks_start, xk_start, wk_start, & nk1, nk2, nk3, k1, k2, k3 USE paw_variables, ONLY : okpaw USE modes, ONLY : nsymq, invsymq, minus_q USE uspp_param, ONLY : n_atom_wfc USE control_ph, ONLY : symoff, man_kpoints, reduce_io USE wavefunctions_module, ONLY: evc0 USE opt_tetra_mod, ONLY : tetra_type, opt_tetra_init USE ktetra, ONLY : ltetra, tetra, ntetra ! IMPLICIT NONE ! REAL (DP), INTENT(IN) :: xq(3) LOGICAL, INTENT (IN) :: newgrid ! REAL (DP), ALLOCATABLE :: rtau (:,:,:) LOGICAL :: magnetic_sym, sym(48) LOGICAL :: skip_equivalence ! IF ( .NOT. ALLOCATED( force ) ) ALLOCATE( force( 3, nat ) ) ! ! ... threshold for diagonalization ethr - should be good for all cases ! ethr= 1.0D-9 / nelec ! ! ... variables for iterative diagonalization ! ... Davdson: isolve=0, david=4 ; CG: isolve=1, david=1 isolve = 0 david = 4 nbndx = david*nbnd max_cg_iter=20 natomwfc = n_atom_wfc( nat, ityp, noncolin ) ! #ifdef __MPI IF ( use_para_diag ) CALL check_para_diag( nbnd ) #else use_para_diag = .FALSE. #endif IF(man_kpoints .or. symoff) THEN time_reversal = .false. magnetic_sym = .true. nsymq = 1 minus_q = .false. invsymq = .false. nrot = 1 write(stdout,*)'nsymq, nsym, minus_q',nsymq,nsym, minus_q t_rev = 0 nsym = 1 ELSE ! ! ... Symmetry and k-point section ! ! ... time_reversal = use q=>-q symmetry for k-point generation ! magnetic_sym = noncolin .AND. domag time_reversal = .NOT. noinv .AND. .NOT. magnetic_sym ! ! ... smallg_q flags in symmetry operations of the crystal ! ... that are not symmetry operations of the small group of q ! CALL set_small_group_of_q(nsymq,invsymq,minus_q) write(stdout,*)'nsymq,nsym,minus_q',nsymq,nsym, minus_q ! small group of q nsym=nsymq invsym = invsymq END IF !CALL mp_bcast (nsymq, root, world_comm) ! ! ... Input k-points are assumed to be given in the IBZ of the Bravais ! ... lattice, with the full point symmetry of the lattice. ! !print*, nks_start, newgrid ! if( nks_start > 0 .AND. .NOT. newgrid ) then ! ! In this case I keep the same points of the Charge density ! calculations ! ! nkstot = nks_start ! xk(:,1:nkstot) = xk_start(:,1:nkstot) ! wk(1:nkstot) = wk_start(1:nkstot) ! else ! ! In this case I generate a new set of k-points ! ! In the case of electron-phonon matrix element with ! wannier functions the k-points should not be reduced ! ! t_rev(:) = 0 IF(.NOT. man_kpoints)THEN skip_equivalence=.false. CALL kpoint_grid ( nrot, time_reversal, skip_equivalence, s, t_rev, & bg, nk1*nk2*nk3, k1,k2,k3, nk1,nk2,nk3, nkstot, xk, wk) ! endif ! ! ... If some symmetries of the lattice are missing in the crystal, ! ... "irreducible_BZ" computes the missing k-points. ! CALL irreducible_BZ (nrot, s, nsymq, minus_q, magnetic_sym, & at, bg, npk, nkstot, xk, wk, t_rev) ENDIF ! ... add k+q to the list of k ! CALL set_kplusq( xk, wk, xq, nkstot, npk ) IF ( ABS( xq(1) ) < eps8 .AND. ABS( xq(2) ) < eps8 .AND. & ABS( xq(3) ) < eps8 ) THEN ! kunit = 1 ! ELSE ! kunit = 2 ! ENDIF IF(ltetra .AND. (tetra_type /= 0)) THEN ntetra = 6 * nk1 * nk2 * nk3 IF(ALLOCATED(tetra)) DEALLOCATE(tetra) ALLOCATE( tetra( 20, ntetra ) ) CALL opt_tetra_init(nsymq, s, time_reversal .AND. minus_q, t_rev, at, bg, npk, k1, k2, k3, & & nk1, nk2, nk3, nkstot, xk, tetra, kunit) END IF ! IF ( lsda ) THEN ! ! ... LSDA case: two different spin polarizations, ! ... each with its own kpoints ! if (nspin /= 2) call errore ('setup_nscf','nspin should be 2; check iosys',1) ! CALL set_kup_and_kdw( xk, wk, isk, nkstot, npk ) ! ELSE IF ( noncolin ) THEN ! ! ... noncolinear magnetism: potential and charge have dimension 4 (1+3) ! if (nspin /= 4) call errore ('setup_nscf','nspin should be 4; check iosys',1) current_spin = 1 ! ELSE ! ! ... LDA case: the two spin polarizations are identical ! wk(1:nkstot) = wk(1:nkstot) * degspin current_spin = 1 ! IF ( nspin /= 1 ) & CALL errore( 'setup_nscf', 'nspin should be 1; check iosys', 1 ) ! END IF ! IF ( nkstot > npk ) CALL errore( 'setup_nscf', 'too many k points', nkstot ) ! ! ...notice: qnorm is used by allocate_nlpot to determine ! the correct size of the interpolation table "qrad" ! qnorm = sqrt(xq(1)**2 + xq(2)**2 + xq(3)**2) ! #ifdef __MPI ! ! ... set the granularity for k-point distribution ! ! ... distribute k-points (and their weights and spin indices) ! CALL divide_et_impera( xk, wk, isk, lsda, nkstot, nks ) ! #else ! nks = nkstot ! #endif IF(reduce_io)THEN !allocate (evc0 ( npwx*npol , nbnd, nks)) io_level = -5 END IF !KC: allocate the space for storing wavefunctions in memory ! RETURN ! END SUBROUTINE setup_nscf