1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculation of band structures 8!> \par History 9!> 2015.06 created [JGH] 10!> \author JGH 11! ************************************************************************************************** 12MODULE qs_band_structure 13 USE cell_types, ONLY: cell_type,& 14 get_cell 15 USE cp_blacs_env, ONLY: cp_blacs_env_retain,& 16 cp_blacs_env_type 17 USE cp_control_types, ONLY: dft_control_type 18 USE cp_files, ONLY: close_file,& 19 open_file 20 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type 21 USE cp_fm_types, ONLY: cp_fm_type 22 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit 23 USE cp_para_env, ONLY: cp_para_env_retain 24 USE cp_para_types, ONLY: cp_para_env_type 25 USE dbcsr_api, ONLY: dbcsr_p_type 26 USE input_section_types, ONLY: section_vals_get,& 27 section_vals_get_subs_vals,& 28 section_vals_type,& 29 section_vals_val_get 30 USE kinds, ONLY: default_string_length,& 31 dp 32 USE kpoint_methods, ONLY: kpoint_env_initialize,& 33 kpoint_init_cell_index,& 34 kpoint_initialize_mos 35 USE kpoint_types, ONLY: get_kpoint_info,& 36 kpoint_create,& 37 kpoint_env_type,& 38 kpoint_release,& 39 kpoint_sym_create,& 40 kpoint_type 41 USE machine, ONLY: m_walltime 42 USE mathconstants, ONLY: twopi 43 USE message_passing, ONLY: mp_sum 44 USE physcon, ONLY: angstrom,& 45 evolt 46 USE qs_environment_types, ONLY: get_qs_env,& 47 qs_env_release,& 48 qs_environment_type 49 USE qs_gamma2kp, ONLY: create_kp_from_gamma 50 USE qs_matrix_pools, ONLY: mpools_get 51 USE qs_mo_types, ONLY: get_mo_set,& 52 init_mo_set,& 53 mo_set_p_type 54 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 55 USE qs_scf_diagonalization, ONLY: do_general_diag_kp 56 USE qs_scf_types, ONLY: qs_scf_env_type 57 USE scf_control_types, ONLY: scf_control_type 58 USE string_utilities, ONLY: uppercase 59#include "./base/base_uses.f90" 60 61 IMPLICIT NONE 62 63 PRIVATE 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_band_structure' 66 67 PUBLIC :: calculate_band_structure, calculate_kp_orbitals 68 69! ************************************************************************************************** 70 71CONTAINS 72 73! ************************************************************************************************** 74!> \brief Main routine for band structure calculation 75!> \param qs_env ... 76!> \author JGH 77! ************************************************************************************************** 78 SUBROUTINE calculate_band_structure(qs_env) 79 TYPE(qs_environment_type), POINTER :: qs_env 80 81 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_band_structure', & 82 routineP = moduleN//':'//routineN 83 84 LOGICAL :: do_kpoints, explicit 85 TYPE(qs_environment_type), POINTER :: qs_env_kp 86 TYPE(section_vals_type), POINTER :: bs_input 87 88 bs_input => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%BAND_STRUCTURE") 89 CALL section_vals_get(bs_input, explicit=explicit) 90 IF (explicit) THEN 91 CALL get_qs_env(qs_env, do_kpoints=do_kpoints) 92 IF (do_kpoints) THEN 93 CALL do_calculate_band_structure(qs_env) 94 ELSE 95 NULLIFY (qs_env_kp) 96 CALL create_kp_from_gamma(qs_env, qs_env_kp) 97 CALL do_calculate_band_structure(qs_env_kp) 98 CALL qs_env_release(qs_env_kp) 99 END IF 100 END IF 101 102 END SUBROUTINE calculate_band_structure 103 104! ************************************************************************************************** 105!> \brief band structure calculation 106!> \param qs_env ... 107!> \author JGH 108! ************************************************************************************************** 109 SUBROUTINE do_calculate_band_structure(qs_env) 110 TYPE(qs_environment_type), POINTER :: qs_env 111 112 CHARACTER(LEN=*), PARAMETER :: routineN = 'do_calculate_band_structure', & 113 routineP = moduleN//':'//routineN 114 115 CHARACTER(LEN=default_string_length) :: filename, ustr 116 CHARACTER(LEN=default_string_length), & 117 DIMENSION(:), POINTER :: spname, strptr 118 INTEGER :: i_rep, ik, ikk, ikpgr, ip, ispin, iw, & 119 n_ptr, n_rep, nadd, nkp, nmo, npline, & 120 npoints, nspins, unit_nr 121 INTEGER, DIMENSION(2) :: kp_range 122 LOGICAL :: explicit, io_default, my_kpgrp 123 REAL(KIND=dp) :: t1, t2 124 REAL(KIND=dp), DIMENSION(3) :: kpptr 125 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, eigval 126 REAL(kind=dp), DIMENSION(:, :), POINTER :: kpgeneral, kspecial, xkp 127 TYPE(cell_type), POINTER :: cell 128 TYPE(cp_para_env_type), POINTER :: para_env 129 TYPE(dft_control_type), POINTER :: dft_control 130 TYPE(kpoint_env_type), POINTER :: kp 131 TYPE(kpoint_type), POINTER :: kpoint 132 TYPE(section_vals_type), POINTER :: bs_input, kpset 133 134 bs_input => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%BAND_STRUCTURE") 135 CALL section_vals_get(bs_input, explicit=explicit) 136 IF (explicit) THEN 137 CALL section_vals_val_get(bs_input, "FILE_NAME", c_val=filename) 138 CALL section_vals_val_get(bs_input, "ADDED_MOS", i_val=nadd) 139 unit_nr = cp_logger_get_default_io_unit() 140 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 141 CALL get_qs_env(qs_env, cell=cell) 142 kpset => section_vals_get_subs_vals(bs_input, "KPOINT_SET") 143 CALL section_vals_get(kpset, n_repetition=n_rep) 144 IF (unit_nr > 0) THEN 145 WRITE (unit_nr, FMT="(/,T2,A)") "KPOINTS| Band Structure Calculation" 146 WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of K-Point Sets", n_rep 147 IF (nadd /= 0) THEN 148 WRITE (unit_nr, FMT="(T2,A,T71,I10)") "KPOINTS| Number of added Mos/Bands", nadd 149 END IF 150 END IF 151 IF (filename == "") THEN 152 ! use standard output file 153 iw = unit_nr 154 io_default = .TRUE. 155 ELSE 156 io_default = .FALSE. 157 IF (para_env%ionode) THEN 158 CALL open_file(filename, unit_number=iw, file_status="UNKNOWN", file_action="WRITE", & 159 file_position="APPEND") 160 ELSE 161 iw = -1 162 END IF 163 END IF 164 DO i_rep = 1, n_rep 165 t1 = m_walltime() 166 CALL section_vals_val_get(kpset, "NPOINTS", i_rep_section=i_rep, i_val=npline) 167 CALL section_vals_val_get(kpset, "UNITS", i_rep_section=i_rep, c_val=ustr) 168 CALL uppercase(ustr) 169 CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_section=i_rep, n_rep_val=n_ptr) 170 ALLOCATE (kspecial(3, n_ptr)) 171 ALLOCATE (spname(n_ptr)) 172 DO ip = 1, n_ptr 173 CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_section=i_rep, i_rep_val=ip, c_vals=strptr) 174 IF (SIZE(strptr(:), 1) == 4) THEN 175 spname(ip) = strptr(1) 176 READ (strptr(2), "(F12.8)") kpptr(1) 177 READ (strptr(3), "(F12.8)") kpptr(2) 178 READ (strptr(4), "(F12.8)") kpptr(3) 179 ELSE IF (SIZE(strptr(:), 1) == 3) THEN 180 spname(ip) = "not specified" 181 READ (strptr(1), "(F12.8)") kpptr(1) 182 READ (strptr(2), "(F12.8)") kpptr(2) 183 READ (strptr(3), "(F12.8)") kpptr(3) 184 ELSE 185 CPABORT("Input SPECIAL_POINT invalid") 186 END IF 187 SELECT CASE (ustr) 188 CASE ("B_VECTOR") 189 kspecial(1:3, ip) = kpptr(1:3) 190 CASE ("CART_ANGSTROM") 191 kspecial(1:3, ip) = (kpptr(1)*cell%hmat(1, 1:3) + & 192 kpptr(2)*cell%hmat(2, 1:3) + & 193 kpptr(3)*cell%hmat(3, 1:3))/twopi*angstrom 194 CASE ("CART_BOHR") 195 kspecial(1:3, ip) = (kpptr(1)*cell%hmat(1, 1:3) + & 196 kpptr(2)*cell%hmat(2, 1:3) + & 197 kpptr(3)*cell%hmat(3, 1:3))/twopi 198 CASE DEFAULT 199 CPABORT("Unknown Unit for kpoint definition") 200 END SELECT 201 END DO 202 npoints = (n_ptr - 1)*npline + 1 203 ! 204 CPASSERT(npoints >= 1) 205 IF (unit_nr > 0) THEN 206 WRITE (unit_nr, "(T2,A,I4,T71,I10)") "KPOINTS| Number of K-Points in Set ", i_rep, npoints 207 WRITE (unit_nr, "(T2,A)") "KPOINTS| In units of b-vector [2pi/Bohr]" 208 DO ip = 1, n_ptr 209 WRITE (unit_nr, "(T2,A,I4,T31,A15,3F10.4)") "KPOINTS| Special K-Point ", & 210 ip, ADJUSTL(TRIM(spname(ip))), kspecial(1:3, ip) 211 END DO 212 END IF 213 IF (iw > 0 .AND. (iw /= unit_nr)) THEN 214 WRITE (iw, "(A,I8,T30,A,I8)") " SET:", i_rep, " TOTAL POINTS:", npoints 215 DO ip = 1, n_ptr 216 WRITE (iw, "(A,I4,T30,A15,3F12.6)") " POINT", & 217 ip, ADJUSTL(TRIM(spname(ip))), kspecial(1:3, ip) 218 END DO 219 END IF 220 ! initialize environment and calculate MOS 221 ALLOCATE (kpgeneral(3, npoints)) 222 kpgeneral(1:3, 1) = kspecial(1:3, 1) 223 ikk = 1 224 DO ik = 2, n_ptr 225 DO ip = 1, npline 226 ikk = ikk + 1 227 kpgeneral(1:3, ikk) = kspecial(1:3, ik - 1) + & 228 REAL(ip, KIND=dp)/REAL(npline, KIND=dp)* & 229 (kspecial(1:3, ik) - kspecial(1:3, ik - 1)) 230 END DO 231 END DO 232 NULLIFY (kpoint) 233 CALL calculate_kp_orbitals(qs_env, kpoint, "GENERAL", nadd, kpgeneral=kpgeneral) 234 DEALLOCATE (kpgeneral) 235 ! 236 CALL get_qs_env(qs_env, dft_control=dft_control) 237 nspins = dft_control%nspins 238 kp => kpoint%kp_env(1)%kpoint_env 239 CALL get_mo_set(kp%mos(1, 1)%mo_set, nmo=nmo) 240 ALLOCATE (eigval(nmo)) 241 CALL get_kpoint_info(kpoint, nkp=nkp, kp_range=kp_range, xkp=xkp) 242 DO ik = 1, nkp 243 my_kpgrp = (ik >= kp_range(1) .AND. ik <= kp_range(2)) 244 DO ispin = 1, nspins 245 IF (my_kpgrp) THEN 246 ikpgr = ik - kp_range(1) + 1 247 kp => kpoint%kp_env(ikpgr)%kpoint_env 248 CALL get_mo_set(kp%mos(1, ispin)%mo_set, eigenvalues=eigenvalues) 249 eigval(1:nmo) = eigenvalues(1:nmo) 250 ELSE 251 eigval(1:nmo) = 0.0_dp 252 END IF 253 CALL mp_sum(eigval, kpoint%para_env_inter_kp%group) 254 ! output 255 IF (iw > 0) THEN 256 eigval(1:nmo) = eigval(1:nmo)*evolt 257 WRITE (iw, "(T8,A,I5,T20,A,I2,T34,A,3F12.8)") "Nr.", ik, "Spin", ispin, "K-Point", xkp(1:3, ik) 258 WRITE (iw, "(T8,I10)") nmo 259 WRITE (iw, "(T8,4F16.8)") eigval(1:nmo) 260 END IF 261 END DO 262 END DO 263 ! 264 DEALLOCATE (kspecial, spname) 265 DEALLOCATE (eigval) 266 CALL kpoint_release(kpoint) 267 t2 = m_walltime() 268 IF (unit_nr > 0) THEN 269 WRITE (unit_nr, FMT="(T2,A,T67,F14.3)") "KPOINTS| Time for K-Point Line ", t2 - t1 270 END IF 271 END DO 272 ! close output file 273 IF (.NOT. io_default) THEN 274 IF (para_env%ionode) THEN 275 CALL close_file(iw) 276 END IF 277 END IF 278 END IF 279 280 END SUBROUTINE do_calculate_band_structure 281 282! ************************************************************************************************** 283!> \brief diagonalize KS matrices at a set of kpoints 284!> \param qs_env ... 285!> \param kpoint ... 286!> \param scheme ... 287!> \param nadd ... 288!> \param mp_grid ... 289!> \param kpgeneral ... 290!> \param group_size_ext ... 291!> \author JGH 292! ************************************************************************************************** 293 SUBROUTINE calculate_kp_orbitals(qs_env, kpoint, scheme, nadd, mp_grid, kpgeneral, group_size_ext) 294 TYPE(qs_environment_type), POINTER :: qs_env 295 TYPE(kpoint_type), POINTER :: kpoint 296 CHARACTER(LEN=*), INTENT(IN) :: scheme 297 INTEGER, INTENT(IN) :: nadd 298 INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: mp_grid 299 REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), & 300 OPTIONAL :: kpgeneral 301 INTEGER, INTENT(IN), OPTIONAL :: group_size_ext 302 303 CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_kp_orbitals', & 304 routineP = moduleN//':'//routineN 305 306 INTEGER :: i, ic, ik, ikk, ispin, ix, iy, iz, & 307 npoints 308 REAL(KIND=dp), DIMENSION(3) :: kpt_latt 309 REAL(KIND=dp), DIMENSION(3, 3) :: h_inv 310 TYPE(cell_type), POINTER :: cell 311 TYPE(cp_blacs_env_type), POINTER :: blacs_env 312 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools 313 TYPE(cp_fm_type), POINTER :: mo_coeff 314 TYPE(cp_para_env_type), POINTER :: para_env 315 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_s 316 TYPE(dft_control_type), POINTER :: dft_control 317 TYPE(mo_set_p_type), DIMENSION(:, :), POINTER :: mos 318 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 319 POINTER :: sab_nl 320 TYPE(qs_scf_env_type), POINTER :: scf_env 321 TYPE(scf_control_type), POINTER :: scf_control 322 323 CPASSERT(.NOT. ASSOCIATED(kpoint)) 324 325 CALL kpoint_create(kpoint) 326 327 kpoint%kp_scheme = scheme 328 kpoint%symmetry = .FALSE. 329 kpoint%verbose = .FALSE. 330 kpoint%full_grid = .FALSE. 331 kpoint%use_real_wfn = .FALSE. 332 kpoint%eps_geo = 1.e-6_dp 333 IF (PRESENT(group_size_ext)) THEN 334 kpoint%parallel_group_size = group_size_ext 335 ELSE 336 kpoint%parallel_group_size = -1 337 END IF 338 SELECT CASE (scheme) 339 CASE ("GAMMA") 340 kpoint%nkp = 1 341 ALLOCATE (kpoint%xkp(3, 1), kpoint%wkp(1)) 342 kpoint%xkp(1:3, 1) = 0.0_dp 343 kpoint%wkp(1) = 1.0_dp 344 kpoint%symmetry = .TRUE. 345 ALLOCATE (kpoint%kp_sym(1)) 346 NULLIFY (kpoint%kp_sym(1)%kpoint_sym) 347 CALL kpoint_sym_create(kpoint%kp_sym(1)%kpoint_sym) 348 CASE ("MONKHORST-PACK") 349 CPASSERT(PRESENT(mp_grid)) 350 npoints = mp_grid(1)*mp_grid(2)*mp_grid(3) 351 kpoint%nkp_grid(1:3) = mp_grid(1:3) 352 kpoint%full_grid = .TRUE. 353 kpoint%nkp = npoints 354 ALLOCATE (kpoint%xkp(3, npoints), kpoint%wkp(npoints)) 355 kpoint%wkp(:) = 1._dp/REAL(npoints, KIND=dp) 356 CALL get_qs_env(qs_env, cell=cell) 357 CALL get_cell(cell, h_inv=h_inv) 358 i = 0 359 DO ix = 1, mp_grid(1) 360 DO iy = 1, mp_grid(2) 361 DO iz = 1, mp_grid(3) 362 i = i + 1 363 kpt_latt(1) = REAL(2*ix - mp_grid(1) - 1, KIND=dp)/(2._dp*REAL(mp_grid(1), KIND=dp)) 364 kpt_latt(2) = REAL(2*iy - mp_grid(2) - 1, KIND=dp)/(2._dp*REAL(mp_grid(2), KIND=dp)) 365 kpt_latt(3) = REAL(2*iz - mp_grid(3) - 1, KIND=dp)/(2._dp*REAL(mp_grid(3), KIND=dp)) 366 kpoint%xkp(1:3, i) = MATMUL(TRANSPOSE(h_inv), kpt_latt(:)) 367 END DO 368 END DO 369 END DO 370 ! default: no symmetry settings 371 ALLOCATE (kpoint%kp_sym(kpoint%nkp)) 372 DO i = 1, kpoint%nkp 373 NULLIFY (kpoint%kp_sym(i)%kpoint_sym) 374 CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym) 375 END DO 376 CASE ("MACDONALD") 377 CPABORT("MACDONALD not implemented") 378 CASE ("GENERAL") 379 CPASSERT(PRESENT(kpgeneral)) 380 npoints = SIZE(kpgeneral, 2) 381 kpoint%nkp = npoints 382 ALLOCATE (kpoint%xkp(3, npoints), kpoint%wkp(npoints)) 383 kpoint%wkp(:) = 1._dp/REAL(npoints, KIND=dp) 384 kpoint%xkp(1:3, 1:npoints) = kpgeneral(1:3, 1:npoints) 385 ! default: no symmetry settings 386 ALLOCATE (kpoint%kp_sym(kpoint%nkp)) 387 DO i = 1, kpoint%nkp 388 NULLIFY (kpoint%kp_sym(i)%kpoint_sym) 389 CALL kpoint_sym_create(kpoint%kp_sym(i)%kpoint_sym) 390 END DO 391 CASE DEFAULT 392 CPABORT("Unknown kpoint scheme requested") 393 END SELECT 394 395 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) 396 kpoint%para_env => para_env 397 CALL cp_para_env_retain(para_env) 398 kpoint%blacs_env_all => blacs_env 399 CALL cp_blacs_env_retain(blacs_env) 400 CALL kpoint_env_initialize(kpoint) 401 402 CALL kpoint_initialize_mos(kpoint, qs_env%mos, nadd) 403 DO ik = 1, SIZE(kpoint%kp_env) 404 CALL mpools_get(kpoint%mpools, ao_mo_fm_pools=ao_mo_fm_pools) 405 mos => kpoint%kp_env(ik)%kpoint_env%mos 406 ikk = kpoint%kp_range(1) + ik - 1 407 CPASSERT(ASSOCIATED(mos)) 408 DO ispin = 1, SIZE(mos, 2) 409 DO ic = 1, SIZE(mos, 1) 410 CALL get_mo_set(mos(ic, ispin)%mo_set, mo_coeff=mo_coeff) 411 IF (.NOT. ASSOCIATED(mo_coeff)) THEN 412 CALL init_mo_set(mos(ic, ispin)%mo_set, & 413 fm_pool=ao_mo_fm_pools(ispin)%pool, name="kpoints") 414 END IF 415 END DO 416 END DO 417 END DO 418 CALL get_qs_env(qs_env, sab_kp=sab_nl, dft_control=dft_control) 419 CALL kpoint_init_cell_index(kpoint, sab_nl, para_env, dft_control) 420 ! 421 CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, matrix_s_kp=matrix_s, & 422 scf_env=scf_env, scf_control=scf_control) 423 CALL do_general_diag_kp(matrix_ks, matrix_s, kpoint, scf_env, scf_control, .FALSE.) 424 425 END SUBROUTINE calculate_kp_orbitals 426 427! ************************************************************************************************** 428 429END MODULE qs_band_structure 430