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 and writing of projected density of states 8!> The DOS is computed per angular momentum and per kind 9!> \par History 10!> - 11!> \author Marcella (29.02.2008,MK) 12! ************************************************************************************************** 13MODULE qs_pdos 14 USE atomic_kind_types, ONLY: atomic_kind_type,& 15 get_atomic_kind,& 16 get_atomic_kind_set 17 USE basis_set_types, ONLY: get_gto_basis_set,& 18 gto_basis_set_type 19 USE cell_types, ONLY: cell_type,& 20 pbc 21 USE cp_blacs_env, ONLY: cp_blacs_env_type 22 USE cp_control_types, ONLY: dft_control_type 23 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm 24 USE cp_fm_diag, ONLY: cp_fm_power 25 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 26 cp_fm_struct_release,& 27 cp_fm_struct_type 28 USE cp_fm_types, ONLY: cp_fm_create,& 29 cp_fm_get_info,& 30 cp_fm_get_submatrix,& 31 cp_fm_init_random,& 32 cp_fm_release,& 33 cp_fm_type 34 USE cp_gemm_interface, ONLY: cp_gemm 35 USE cp_log_handling, ONLY: cp_get_default_logger,& 36 cp_logger_get_default_io_unit,& 37 cp_logger_type,& 38 cp_to_string 39 USE cp_output_handling, ONLY: cp_p_file,& 40 cp_print_key_finished_output,& 41 cp_print_key_should_output,& 42 cp_print_key_unit_nr 43 USE cp_para_types, ONLY: cp_para_env_type 44 USE dbcsr_api, ONLY: dbcsr_p_type 45 USE input_section_types, ONLY: section_vals_get,& 46 section_vals_get_subs_vals,& 47 section_vals_type,& 48 section_vals_val_get 49 USE kinds, ONLY: default_string_length,& 50 dp 51 USE memory_utilities, ONLY: reallocate 52 USE message_passing, ONLY: mp_sum 53 USE orbital_pointers, ONLY: nso,& 54 nsoset 55 USE orbital_symbols, ONLY: l_sym,& 56 sgf_symbol 57 USE particle_types, ONLY: particle_type 58 USE preconditioner_types, ONLY: preconditioner_type 59 USE pw_env_types, ONLY: pw_env_get,& 60 pw_env_type 61 USE pw_pool_types, ONLY: pw_pool_create_pw,& 62 pw_pool_give_back_pw,& 63 pw_pool_p_type,& 64 pw_pool_type 65 USE pw_types, ONLY: COMPLEXDATA1D,& 66 REALDATA3D,& 67 REALSPACE,& 68 RECIPROCALSPACE,& 69 pw_p_type 70 USE qs_collocate_density, ONLY: calculate_wavefunction 71 USE qs_environment_types, ONLY: get_qs_env,& 72 qs_environment_type 73 USE qs_kind_types, ONLY: get_qs_kind,& 74 get_qs_kind_set,& 75 qs_kind_type 76 USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues 77 USE qs_mo_types, ONLY: get_mo_set,& 78 mo_set_type 79 USE qs_ot_eigensolver, ONLY: ot_eigensolver 80 USE scf_control_types, ONLY: scf_control_type 81#include "./base/base_uses.f90" 82 83 IMPLICIT NONE 84 85 PRIVATE 86 87 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_pdos' 88 89! ************************************************************************************************** 90 ! *** Public subroutines *** 91 92 PUBLIC :: calculate_projected_dos 93 94 TYPE ldos_type 95 INTEGER :: maxl, nlist 96 LOGICAL :: separate_components 97 INTEGER, DIMENSION(:), POINTER :: list_index 98 REAL(KIND=dp), DIMENSION(:, :), & 99 POINTER :: pdos_array 100 END TYPE ldos_type 101 102 TYPE r_ldos_type 103 INTEGER :: nlist, npoints 104 INTEGER, DIMENSION(:, :), POINTER :: index_grid_local 105 INTEGER, DIMENSION(:), POINTER :: list_index 106 REAL(KIND=dp), DIMENSION(:), POINTER :: x_range, y_range, z_range 107 REAL(KIND=dp), DIMENSION(:), POINTER :: eval_range 108 REAL(KIND=dp), DIMENSION(:), & 109 POINTER :: pdos_array 110 END TYPE r_ldos_type 111 112 TYPE ldos_p_type 113 TYPE(ldos_type), POINTER :: ldos 114 END TYPE ldos_p_type 115 116 TYPE r_ldos_p_type 117 TYPE(r_ldos_type), POINTER :: ldos 118 END TYPE r_ldos_p_type 119CONTAINS 120 121! ************************************************************************************************** 122!> \brief Compute and write projected density of states 123!> \param mo_set ... 124!> \param atomic_kind_set ... 125!> \param qs_kind_set ... 126!> \param particle_set ... 127!> \param qs_env ... 128!> \param dft_section ... 129!> \param ispin ... 130!> \param xas_mittle ... 131!> \param external_matrix_shalf ... 132!> \date 26.02.2008 133!> \par History: 134!> - Added optional external matrix_shalf to avoid recomputing it (A. Bussy, 09.2019) 135!> \par Variables 136!> - 137!> - 138!> \author MI 139!> \version 1.0 140! ************************************************************************************************** 141 SUBROUTINE calculate_projected_dos(mo_set, atomic_kind_set, qs_kind_set, particle_set, qs_env, & 142 dft_section, ispin, xas_mittle, external_matrix_shalf) 143 144 TYPE(mo_set_type), POINTER :: mo_set 145 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 146 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 147 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 148 TYPE(qs_environment_type), POINTER :: qs_env 149 TYPE(section_vals_type), POINTER :: dft_section 150 INTEGER, INTENT(IN), OPTIONAL :: ispin 151 CHARACTER(LEN=default_string_length), INTENT(IN), & 152 OPTIONAL :: xas_mittle 153 TYPE(cp_fm_type), OPTIONAL, POINTER :: external_matrix_shalf 154 155 CHARACTER(len=*), PARAMETER :: routineN = 'calculate_projected_dos', & 156 routineP = moduleN//':'//routineN 157 158 CHARACTER(LEN=16) :: fmtstr2 159 CHARACTER(LEN=27) :: fmtstr1 160 CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:, :, :) :: tmp_str 161 CHARACTER(LEN=default_string_length) :: kind_name, my_act, my_mittle, my_pos, & 162 spin(2) 163 CHARACTER(LEN=default_string_length), & 164 ALLOCATABLE, DIMENSION(:) :: ldos_index, r_ldos_index 165 INTEGER :: handle, homo, i, iatom, ikind, il, ildos, im, imo, in_x, in_y, in_z, ir, irow, & 166 iset, isgf, ishell, iso, iterstep, iw, j, jx, jy, jz, k, lcomponent, lshell, maxl, & 167 maxlgto, my_spin, n_dependent, n_r_ldos, n_rep, nao, natom, ncol_global, nkind, nldos, & 168 nlumo, nmo, np_tot, npoints, nrow_global, nset, nsgf, nvirt, out_each, output_unit 169 INTEGER, ALLOCATABLE, DIMENSION(:) :: firstrow 170 INTEGER, DIMENSION(:), POINTER :: list, nshell 171 INTEGER, DIMENSION(:, :), POINTER :: bo, l 172 LOGICAL :: append, calc_matsh, do_ldos, do_r_ldos, & 173 do_virt, ionode, separate_components, & 174 should_output 175 LOGICAL, DIMENSION(:, :), POINTER :: read_r 176 REAL(KIND=dp) :: dh(3, 3), dvol, e_fermi, r(3), r_vec(3), & 177 ratom(3) 178 REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, evals_virt, & 179 occupation_numbers 180 REAL(KIND=dp), DIMENSION(:, :), POINTER :: vecbuffer 181 REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pdos_array 182 TYPE(cell_type), POINTER :: cell 183 TYPE(cp_blacs_env_type), POINTER :: context 184 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 185 TYPE(cp_fm_type), POINTER :: matrix_shalf, matrix_shalfc, & 186 matrix_work, mo_coeff, mo_virt 187 TYPE(cp_logger_type), POINTER :: logger 188 TYPE(cp_para_env_type), POINTER :: para_env 189 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: s_matrix 190 TYPE(dft_control_type), POINTER :: dft_control 191 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 192 TYPE(ldos_p_type), DIMENSION(:), POINTER :: ldos_p 193 TYPE(pw_env_type), POINTER :: pw_env 194 TYPE(pw_p_type) :: wf_g, wf_r 195 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 196 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 197 TYPE(r_ldos_p_type), DIMENSION(:), POINTER :: r_ldos_p 198 TYPE(section_vals_type), POINTER :: ldos_section 199 200 NULLIFY (logger) 201 logger => cp_get_default_logger() 202 ionode = logger%para_env%ionode 203 should_output = BTEST(cp_print_key_should_output(logger%iter_info, dft_section, & 204 "PRINT%PDOS"), cp_p_file) 205 output_unit = cp_logger_get_default_io_unit(logger) 206 207 spin(1) = "ALPHA" 208 spin(2) = "BETA" 209 IF ((.NOT. should_output)) RETURN 210 211 NULLIFY (context, s_matrix, orb_basis_set, para_env, pdos_array) 212 NULLIFY (eigenvalues, fm_struct_tmp, matrix_work, matrix_shalf, mo_coeff, vecbuffer) 213 NULLIFY (ldos_section, list, cell, pw_env, auxbas_pw_pool, mo_virt, evals_virt) 214 NULLIFY (occupation_numbers, ldos_p, r_ldos_p, dft_control, occupation_numbers) 215 216 CALL timeset(routineN, handle) 217 iterstep = logger%iter_info%iteration(logger%iter_info%n_rlevel) 218 219 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') & 220 " Calculate PDOS at iteration step ", iterstep 221 CALL get_qs_env(qs_env=qs_env, & 222 matrix_s=s_matrix) 223 224 CALL get_atomic_kind_set(atomic_kind_set, natom=natom) 225 CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf, maxlgto=maxlgto) 226 nkind = SIZE(atomic_kind_set) 227 228 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, homo=homo, nao=nao, nmo=nmo, & 229 mu=e_fermi) 230 CALL cp_fm_get_info(mo_coeff, & 231 context=context, para_env=para_env, & 232 nrow_global=nrow_global, & 233 ncol_global=ncol_global) 234 235 CALL section_vals_val_get(dft_section, "PRINT%PDOS%OUT_EACH_MO", i_val=out_each) 236 IF (out_each == -1) out_each = nao + 1 237 CALL section_vals_val_get(dft_section, "PRINT%PDOS%NLUMO", i_val=nlumo) 238 IF (nlumo == -1) nlumo = nao - homo 239 do_virt = (nlumo > (nmo - homo)) 240 nvirt = nlumo - (nmo - homo) 241 ! Generate virtual orbitals 242 IF (do_virt) THEN 243 IF (PRESENT(ispin)) THEN 244 my_spin = ispin 245 ELSE 246 my_spin = 1 247 END IF 248 249 CALL generate_virtual_mo(qs_env, mo_set, evals_virt, mo_virt, nvirt, ispin=my_spin) 250 ELSE 251 NULLIFY (evals_virt, mo_virt) 252 nvirt = 0 253 END IF 254 255 calc_matsh = .TRUE. 256 IF (PRESENT(external_matrix_shalf)) calc_matsh = .FALSE. 257 258 ! Create S^1/2 : from sparse to full matrix, if no external available 259 IF (calc_matsh) THEN 260 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & 261 nrow_global=nrow_global, ncol_global=nrow_global) 262 CALL cp_fm_create(matrix_shalf, fm_struct_tmp, name="matrix_shalf") 263 CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_work") 264 CALL cp_fm_struct_release(fm_struct_tmp) 265 CALL copy_dbcsr_to_fm(s_matrix(1)%matrix, matrix_shalf) 266 CALL cp_fm_power(matrix_shalf, matrix_work, 0.5_dp, EPSILON(0.0_dp), n_dependent) 267 CALL cp_fm_release(matrix_work) 268 ELSE 269 matrix_shalf => external_matrix_shalf 270 END IF 271 272 ! Multiply S^(1/2) time the mOS coefficients to get orthonormalized MOS 273 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & 274 nrow_global=nrow_global, ncol_global=ncol_global) 275 CALL cp_fm_create(matrix_shalfc, fm_struct_tmp, name="matrix_shalfc") 276 CALL cp_gemm("N", "N", nrow_global, ncol_global, nrow_global, & 277 1.0_dp, matrix_shalf, mo_coeff, 0.0_dp, matrix_shalfc) 278 CALL cp_fm_struct_release(fm_struct_tmp) 279 280 IF (do_virt) THEN 281 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T14,I10,T27,A))') & 282 " Compute ", nvirt, " additional unoccupied KS orbitals" 283 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & 284 nrow_global=nrow_global, ncol_global=nvirt) 285 CALL cp_fm_create(matrix_work, fm_struct_tmp, name="matrix_shalfc") 286 CALL cp_gemm("N", "N", nrow_global, nvirt, nrow_global, & 287 1.0_dp, matrix_shalf, mo_virt, 0.0_dp, matrix_work) 288 CALL cp_fm_struct_release(fm_struct_tmp) 289 END IF 290 291 IF (calc_matsh) CALL cp_fm_release(matrix_shalf) 292 ! Array to store the PDOS per kind and angular momentum 293 do_ldos = .FALSE. 294 ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%LDOS") 295 296 CALL section_vals_get(ldos_section, n_repetition=nldos) 297 IF (nldos > 0) THEN 298 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') & 299 " Prepare the list of atoms for LDOS. Number of lists: ", nldos 300 do_ldos = .TRUE. 301 ALLOCATE (ldos_p(nldos)) 302 ALLOCATE (ldos_index(nldos)) 303 DO ildos = 1, nldos 304 WRITE (ldos_index(ildos), '(I0)') ildos 305 ALLOCATE (ldos_p(ildos)%ldos) 306 NULLIFY (ldos_p(ildos)%ldos%pdos_array) 307 NULLIFY (ldos_p(ildos)%ldos%list_index) 308 309 CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep) 310 IF (n_rep > 0) THEN 311 ldos_p(ildos)%ldos%nlist = 0 312 DO ir = 1, n_rep 313 NULLIFY (list) 314 CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, & 315 i_vals=list) 316 IF (ASSOCIATED(list)) THEN 317 CALL reallocate(ldos_p(ildos)%ldos%list_index, 1, ldos_p(ildos)%ldos%nlist + SIZE(list)) 318 DO i = 1, SIZE(list) 319 ldos_p(ildos)%ldos%list_index(i + ldos_p(ildos)%ldos%nlist) = list(i) 320 END DO 321 ldos_p(ildos)%ldos%nlist = ldos_p(ildos)%ldos%nlist + SIZE(list) 322 END IF 323 END DO 324 ELSE 325 ! stop, LDOS without list of atoms is not implemented 326 END IF 327 328 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') & 329 " List ", ildos, " contains ", ldos_p(ildos)%ldos%nlist, " atoms" 330 CALL section_vals_val_get(ldos_section, "COMPONENTS", i_rep_section=ildos, & 331 l_val=ldos_p(ildos)%ldos%separate_components) 332 IF (ldos_p(ildos)%ldos%separate_components) THEN 333 ALLOCATE (ldos_p(ildos)%ldos%pdos_array(nsoset(maxlgto), nmo + nvirt)) 334 ELSE 335 ALLOCATE (ldos_p(ildos)%ldos%pdos_array(0:maxlgto, nmo + nvirt)) 336 END IF 337 ldos_p(ildos)%ldos%pdos_array = 0.0_dp 338 ldos_p(ildos)%ldos%maxl = -1 339 340 END DO 341 END IF 342 343 do_r_ldos = .FALSE. 344 ldos_section => section_vals_get_subs_vals(dft_section, "PRINT%PDOS%R_LDOS") 345 CALL section_vals_get(ldos_section, n_repetition=n_r_ldos) 346 IF (n_r_ldos > 0) THEN 347 do_r_ldos = .TRUE. 348 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T3,A,T61,I10))') & 349 " Prepare the list of points for R_LDOS. Number of lists: ", n_r_ldos 350 ALLOCATE (r_ldos_p(n_r_ldos)) 351 ALLOCATE (r_ldos_index(n_r_ldos)) 352 CALL get_qs_env(qs_env=qs_env, & 353 cell=cell, & 354 dft_control=dft_control, & 355 pw_env=pw_env) 356 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 357 pw_pools=pw_pools) 358 359 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 360 use_data=REALDATA3D, & 361 in_space=REALSPACE) 362 CALL pw_pool_create_pw(auxbas_pw_pool, wf_g%pw, & 363 use_data=COMPLEXDATA1D, & 364 in_space=RECIPROCALSPACE) 365 ALLOCATE (read_r(4, n_r_ldos)) 366 DO ildos = 1, n_r_ldos 367 WRITE (r_ldos_index(ildos), '(I0)') ildos 368 ALLOCATE (r_ldos_p(ildos)%ldos) 369 NULLIFY (r_ldos_p(ildos)%ldos%pdos_array) 370 NULLIFY (r_ldos_p(ildos)%ldos%list_index) 371 372 CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, n_rep_val=n_rep) 373 IF (n_rep > 0) THEN 374 r_ldos_p(ildos)%ldos%nlist = 0 375 DO ir = 1, n_rep 376 NULLIFY (list) 377 CALL section_vals_val_get(ldos_section, "LIST", i_rep_section=ildos, i_rep_val=ir, & 378 i_vals=list) 379 IF (ASSOCIATED(list)) THEN 380 CALL reallocate(r_ldos_p(ildos)%ldos%list_index, 1, r_ldos_p(ildos)%ldos%nlist + SIZE(list)) 381 DO i = 1, SIZE(list) 382 r_ldos_p(ildos)%ldos%list_index(i + r_ldos_p(ildos)%ldos%nlist) = list(i) 383 END DO 384 r_ldos_p(ildos)%ldos%nlist = r_ldos_p(ildos)%ldos%nlist + SIZE(list) 385 END IF 386 END DO 387 ELSE 388 ! stop, LDOS without list of atoms is not implemented 389 END IF 390 391 ALLOCATE (r_ldos_p(ildos)%ldos%pdos_array(nmo + nvirt)) 392 r_ldos_p(ildos)%ldos%pdos_array = 0.0_dp 393 read_r(1:3, ildos) = .FALSE. 394 CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, explicit=read_r(1, ildos)) 395 IF (read_r(1, ildos)) THEN 396 CALL section_vals_val_get(ldos_section, "XRANGE", i_rep_section=ildos, r_vals= & 397 r_ldos_p(ildos)%ldos%x_range) 398 ELSE 399 ALLOCATE (r_ldos_p(ildos)%ldos%x_range(2)) 400 r_ldos_p(ildos)%ldos%x_range = 0.0_dp 401 END IF 402 CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, explicit=read_r(2, ildos)) 403 IF (read_r(2, ildos)) THEN 404 CALL section_vals_val_get(ldos_section, "YRANGE", i_rep_section=ildos, r_vals= & 405 r_ldos_p(ildos)%ldos%y_range) 406 ELSE 407 ALLOCATE (r_ldos_p(ildos)%ldos%y_range(2)) 408 r_ldos_p(ildos)%ldos%y_range = 0.0_dp 409 END IF 410 CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, explicit=read_r(3, ildos)) 411 IF (read_r(3, ildos)) THEN 412 CALL section_vals_val_get(ldos_section, "ZRANGE", i_rep_section=ildos, r_vals= & 413 r_ldos_p(ildos)%ldos%z_range) 414 ELSE 415 ALLOCATE (r_ldos_p(ildos)%ldos%z_range(2)) 416 r_ldos_p(ildos)%ldos%z_range = 0.0_dp 417 END IF 418 419 CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, explicit=read_r(4, ildos)) 420 IF (read_r(4, ildos)) THEN 421 CALL section_vals_val_get(ldos_section, "ERANGE", i_rep_section=ildos, & 422 r_vals=r_ldos_p(ildos)%ldos%eval_range) 423 ELSE 424 ALLOCATE (r_ldos_p(ildos)%ldos%eval_range(2)) 425 r_ldos_p(ildos)%ldos%eval_range(1) = -HUGE(0.0_dp) 426 r_ldos_p(ildos)%ldos%eval_range(2) = +HUGE(0.0_dp) 427 END IF 428 429 bo => wf_r%pw%pw_grid%bounds_local 430 dh = wf_r%pw%pw_grid%dh 431 dvol = wf_r%pw%pw_grid%dvol 432 np_tot = wf_r%pw%pw_grid%npts(1)*wf_r%pw%pw_grid%npts(2)*wf_r%pw%pw_grid%npts(3) 433 ALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local(3, np_tot)) 434 435 r_ldos_p(ildos)%ldos%npoints = 0 436 DO jz = bo(1, 3), bo(2, 3) 437 DO jy = bo(1, 2), bo(2, 2) 438 DO jx = bo(1, 1), bo(2, 1) 439 !compute the position of the grid point 440 i = jx - wf_r%pw%pw_grid%bounds(1, 1) 441 j = jy - wf_r%pw%pw_grid%bounds(1, 2) 442 k = jz - wf_r%pw%pw_grid%bounds(1, 3) 443 r(3) = k*dh(3, 3) + j*dh(3, 2) + i*dh(3, 1) 444 r(2) = k*dh(2, 3) + j*dh(2, 2) + i*dh(2, 1) 445 r(1) = k*dh(1, 3) + j*dh(1, 2) + i*dh(1, 1) 446 447 DO il = 1, r_ldos_p(ildos)%ldos%nlist 448 iatom = r_ldos_p(ildos)%ldos%list_index(il) 449 ratom = particle_set(iatom)%r 450 r_vec = pbc(ratom, r, cell) 451 IF (cell%orthorhombic) THEN 452 IF (cell%perd(1) == 0) r_vec(1) = MODULO(r_vec(1), cell%hmat(1, 1)) 453 IF (cell%perd(2) == 0) r_vec(2) = MODULO(r_vec(2), cell%hmat(2, 2)) 454 IF (cell%perd(3) == 0) r_vec(3) = MODULO(r_vec(3), cell%hmat(3, 3)) 455 END IF 456 457 in_x = 0 458 in_y = 0 459 in_z = 0 460 IF (r_ldos_p(ildos)%ldos%x_range(1) /= 0.0_dp) THEN 461 IF (r_vec(1) > r_ldos_p(ildos)%ldos%x_range(1) .AND. & 462 r_vec(1) < r_ldos_p(ildos)%ldos%x_range(2)) THEN 463 in_x = 1 464 END IF 465 ELSE 466 in_x = 1 467 END IF 468 IF (r_ldos_p(ildos)%ldos%y_range(1) /= 0.0_dp) THEN 469 IF (r_vec(2) > r_ldos_p(ildos)%ldos%y_range(1) .AND. & 470 r_vec(2) < r_ldos_p(ildos)%ldos%y_range(2)) THEN 471 in_y = 1 472 END IF 473 ELSE 474 in_y = 1 475 END IF 476 IF (r_ldos_p(ildos)%ldos%z_range(1) /= 0.0_dp) THEN 477 IF (r_vec(3) > r_ldos_p(ildos)%ldos%z_range(1) .AND. & 478 r_vec(3) < r_ldos_p(ildos)%ldos%z_range(2)) THEN 479 in_z = 1 480 END IF 481 ELSE 482 in_z = 1 483 END IF 484 IF (in_x*in_y*in_z > 0) THEN 485 r_ldos_p(ildos)%ldos%npoints = r_ldos_p(ildos)%ldos%npoints + 1 486 r_ldos_p(ildos)%ldos%index_grid_local(1, r_ldos_p(ildos)%ldos%npoints) = jx 487 r_ldos_p(ildos)%ldos%index_grid_local(2, r_ldos_p(ildos)%ldos%npoints) = jy 488 r_ldos_p(ildos)%ldos%index_grid_local(3, r_ldos_p(ildos)%ldos%npoints) = jz 489 EXIT 490 END IF 491 END DO 492 END DO 493 END DO 494 END DO 495 CALL reallocate(r_ldos_p(ildos)%ldos%index_grid_local, 1, 3, 1, r_ldos_p(ildos)%ldos%npoints) 496 npoints = r_ldos_p(ildos)%ldos%npoints 497 CALL mp_sum(npoints, para_env%group) 498 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='((T10,A,T18,I6,T25,A,T36,I10,A))') & 499 " List ", ildos, " contains ", npoints, " grid points" 500 END DO 501 END IF 502 503 CALL section_vals_val_get(dft_section, "PRINT%PDOS%COMPONENTS", l_val=separate_components) 504 IF (separate_components) THEN 505 ALLOCATE (pdos_array(nsoset(maxlgto), nkind, nmo + nvirt)) 506 ELSE 507 ALLOCATE (pdos_array(0:maxlgto, nkind, nmo + nvirt)) 508 END IF 509 IF (do_virt) THEN 510 ALLOCATE (eigenvalues(nmo + nvirt)) 511 eigenvalues(1:nmo) = mo_set%eigenvalues(1:nmo) 512 eigenvalues(nmo + 1:nmo + nvirt) = evals_virt(1:nvirt) 513 ALLOCATE (occupation_numbers(nmo + nvirt)) 514 occupation_numbers(:) = 0.0_dp 515 occupation_numbers(1:nmo) = mo_set%occupation_numbers(1:nmo) 516 ELSE 517 eigenvalues => mo_set%eigenvalues 518 occupation_numbers => mo_set%occupation_numbers 519 END IF 520 521 pdos_array = 0.0_dp 522 nao = mo_set%nao 523 ALLOCATE (vecbuffer(1, nao)) 524 vecbuffer = 0.0_dp 525 ALLOCATE (firstrow(natom)) 526 firstrow = 0 527 528 !Adjust energy range for r_ldos 529 DO ildos = 1, n_r_ldos 530 IF (eigenvalues(1) > r_ldos_p(ildos)%ldos%eval_range(1)) & 531 r_ldos_p(ildos)%ldos%eval_range(1) = eigenvalues(1) 532 IF (eigenvalues(nmo + nvirt) < r_ldos_p(ildos)%ldos%eval_range(2)) & 533 r_ldos_p(ildos)%ldos%eval_range(2) = eigenvalues(nmo + nvirt) 534 END DO 535 536 IF (output_unit > 0) WRITE (UNIT=output_unit, FMT='(/,(T15,A))') & 537 "---- PDOS: start iteration on the KS states --- " 538 539 DO imo = 1, nmo + nvirt 540 541 IF (output_unit > 0 .AND. MOD(imo, out_each) == 0) WRITE (UNIT=output_unit, FMT='((T20,A,I10))') & 542 " KS state index : ", imo 543 ! Extract the eigenvector from the distributed full matrix 544 IF (imo > nmo) THEN 545 CALL cp_fm_get_submatrix(matrix_work, vecbuffer, 1, imo - nmo, & 546 nao, 1, transpose=.TRUE.) 547 ELSE 548 CALL cp_fm_get_submatrix(matrix_shalfc, vecbuffer, 1, imo, & 549 nao, 1, transpose=.TRUE.) 550 END IF 551 552 ! Calculate the pdos for all the kinds 553 irow = 1 554 DO iatom = 1, natom 555 firstrow(iatom) = irow 556 NULLIFY (orb_basis_set) 557 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) 558 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 559 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 560 nset=nset, & 561 nshell=nshell, & 562 l=l, maxl=maxl) 563 IF (separate_components) THEN 564 isgf = 1 565 DO iset = 1, nset 566 DO ishell = 1, nshell(iset) 567 lshell = l(ishell, iset) 568 DO iso = 1, nso(lshell) 569 lcomponent = nsoset(lshell - 1) + iso 570 pdos_array(lcomponent, ikind, imo) = & 571 pdos_array(lcomponent, ikind, imo) + & 572 vecbuffer(1, irow)*vecbuffer(1, irow) 573 irow = irow + 1 574 END DO ! iso 575 END DO ! ishell 576 END DO ! iset 577 ELSE 578 isgf = 1 579 DO iset = 1, nset 580 DO ishell = 1, nshell(iset) 581 lshell = l(ishell, iset) 582 DO iso = 1, nso(lshell) 583 pdos_array(lshell, ikind, imo) = & 584 pdos_array(lshell, ikind, imo) + & 585 vecbuffer(1, irow)*vecbuffer(1, irow) 586 irow = irow + 1 587 END DO ! iso 588 END DO ! ishell 589 END DO ! iset 590 END IF 591 END DO ! iatom 592 593 ! Calculate the pdos for all the lists 594 DO ildos = 1, nldos 595 DO il = 1, ldos_p(ildos)%ldos%nlist 596 iatom = ldos_p(ildos)%ldos%list_index(il) 597 598 irow = firstrow(iatom) 599 NULLIFY (orb_basis_set) 600 CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind) 601 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 602 603 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, & 604 nset=nset, & 605 nshell=nshell, & 606 l=l, maxl=maxl) 607 ldos_p(ildos)%ldos%maxl = MAX(ldos_p(ildos)%ldos%maxl, maxl) 608 IF (ldos_p(ildos)%ldos%separate_components) THEN 609 isgf = 1 610 DO iset = 1, nset 611 DO ishell = 1, nshell(iset) 612 lshell = l(ishell, iset) 613 DO iso = 1, nso(lshell) 614 lcomponent = nsoset(lshell - 1) + iso 615 ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) = & 616 ldos_p(ildos)%ldos%pdos_array(lcomponent, imo) + & 617 vecbuffer(1, irow)*vecbuffer(1, irow) 618 irow = irow + 1 619 END DO ! iso 620 END DO ! ishell 621 END DO ! iset 622 ELSE 623 isgf = 1 624 DO iset = 1, nset 625 DO ishell = 1, nshell(iset) 626 lshell = l(ishell, iset) 627 DO iso = 1, nso(lshell) 628 ldos_p(ildos)%ldos%pdos_array(lshell, imo) = & 629 ldos_p(ildos)%ldos%pdos_array(lshell, imo) + & 630 vecbuffer(1, irow)*vecbuffer(1, irow) 631 irow = irow + 1 632 END DO ! iso 633 END DO ! ishell 634 END DO ! iset 635 END IF 636 END DO !il 637 END DO !ildos 638 639 ! Calculate the DOS projected in a given volume in real space 640 DO ildos = 1, n_r_ldos 641 IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. & 642 r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN 643 644 IF (imo > nmo) THEN 645 CALL calculate_wavefunction(mo_virt, imo - nmo, & 646 wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, & 647 pw_env) 648 ELSE 649 CALL calculate_wavefunction(mo_coeff, imo, & 650 wf_r, wf_g, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, & 651 pw_env) 652 END IF 653 r_ldos_p(ildos)%ldos%pdos_array(imo) = 0.0_dp 654 DO il = 1, r_ldos_p(ildos)%ldos%npoints 655 j = j + 1 656 jx = r_ldos_p(ildos)%ldos%index_grid_local(1, il) 657 jy = r_ldos_p(ildos)%ldos%index_grid_local(2, il) 658 jz = r_ldos_p(ildos)%ldos%index_grid_local(3, il) 659 r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo) + & 660 wf_r%pw%cr3d(jx, jy, jz)*wf_r%pw%cr3d(jx, jy, jz) 661 END DO 662 r_ldos_p(ildos)%ldos%pdos_array(imo) = r_ldos_p(ildos)%ldos%pdos_array(imo)*dvol 663 END IF 664 END DO 665 END DO ! imo 666 667 CALL cp_fm_release(matrix_shalfc) 668 DEALLOCATE (vecbuffer) 669 670 CALL section_vals_val_get(dft_section, "PRINT%PDOS%APPEND", l_val=append) 671 IF (append .AND. iterstep > 1) THEN 672 my_pos = "APPEND" 673 ELSE 674 my_pos = "REWIND" 675 END IF 676 my_act = "WRITE" 677 DO ikind = 1, nkind 678 679 NULLIFY (orb_basis_set) 680 CALL get_atomic_kind(atomic_kind_set(ikind), name=kind_name) 681 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set) 682 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, maxl=maxl) 683 684 ! basis none has no associated maxl, and no pdos 685 IF (maxl < 0) CYCLE 686 687 IF (PRESENT(ispin)) THEN 688 IF (PRESENT(xas_mittle)) THEN 689 my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind))) 690 ELSE 691 my_mittle = TRIM(spin(ispin))//"_k"//TRIM(ADJUSTL(cp_to_string(ikind))) 692 END IF 693 my_spin = ispin 694 ELSE 695 my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind))) 696 my_spin = 1 697 END IF 698 699 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", & 700 extension=".pdos", file_position=my_pos, file_action=my_act, & 701 file_form="FORMATTED", middle_name=TRIM(my_mittle)) 702 IF (iw > 0) THEN 703 704 fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))" 705 fmtstr2 = "(A42, (10X,A8))" 706 IF (separate_components) THEN 707 WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(maxl) 708 WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(maxl) 709 ELSE 710 WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") maxl + 1 711 WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") maxl + 1 712 END IF 713 714 WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,A)") & 715 "# Projected DOS for atomic kind "//TRIM(kind_name)//" at iteration step i = ", & 716 iterstep, ", E(Fermi) = ", e_fermi, " a.u." 717 IF (separate_components) THEN 718 ALLOCATE (tmp_str(0:0, 0:maxl, -maxl:maxl)) 719 tmp_str = "" 720 DO j = 0, maxl 721 DO i = -j, j 722 tmp_str(0, j, i) = sgf_symbol(0, j, i) 723 END DO 724 END DO 725 726 WRITE (UNIT=iw, FMT=fmtstr2) & 727 "# MO Eigenvalue [a.u.] Occupation", & 728 ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, maxl) 729 DO imo = 1, nmo + nvirt 730 WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), & 731 (pdos_array(lshell, ikind, imo), lshell=1, nsoset(maxl)) 732 END DO 733 DEALLOCATE (tmp_str) 734 ELSE 735 WRITE (UNIT=iw, FMT=fmtstr2) & 736 "# MO Eigenvalue [a.u.] Occupation", & 737 (TRIM(l_sym(il)), il=0, maxl) 738 DO imo = 1, nmo + nvirt 739 WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), & 740 (pdos_array(lshell, ikind, imo), lshell=0, maxl) 741 END DO 742 END IF 743 END IF 744 CALL cp_print_key_finished_output(iw, logger, dft_section, & 745 "PRINT%PDOS") 746 747 END DO ! ikind 748 749 ! write the pdos for the lists, each ona different file, 750 ! the filenames are indexed with the list number 751 DO ildos = 1, nldos 752 ! basis none has no associated maxl, and no pdos 753 IF (ldos_p(ildos)%ldos%maxl > 0) THEN 754 755 IF (PRESENT(ispin)) THEN 756 IF (PRESENT(xas_mittle)) THEN 757 my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos)) 758 ELSE 759 my_mittle = TRIM(spin(ispin))//"_list"//TRIM(ldos_index(ildos)) 760 END IF 761 my_spin = ispin 762 ELSE 763 my_mittle = "list"//TRIM(ldos_index(ildos)) 764 my_spin = 1 765 END IF 766 767 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", & 768 extension=".pdos", file_position=my_pos, file_action=my_act, & 769 file_form="FORMATTED", middle_name=TRIM(my_mittle)) 770 IF (iw > 0) THEN 771 772 fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))" 773 fmtstr2 = "(A42, (10X,A8))" 774 IF (ldos_p(ildos)%ldos%separate_components) THEN 775 WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl) 776 WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") nsoset(ldos_p(ildos)%ldos%maxl) 777 ELSE 778 WRITE (UNIT=fmtstr1(15:16), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 1 779 WRITE (UNIT=fmtstr2(6:7), FMT="(I2)") ldos_p(ildos)%ldos%maxl + 1 780 END IF 781 782 WRITE (UNIT=iw, FMT="(A,I0,A,I0,A,I0,A,F12.6,A)") & 783 "# Projected DOS for list ", ildos, " of ", ldos_p(ildos)%ldos%nlist, & 784 " atoms, at iteration step i = ", iterstep, & 785 ", E(Fermi) = ", e_fermi, " a.u." 786 IF (ldos_p(ildos)%ldos%separate_components) THEN 787 ALLOCATE (tmp_str(0:0, 0:ldos_p(ildos)%ldos%maxl, -ldos_p(ildos)%ldos%maxl:ldos_p(ildos)%ldos%maxl)) 788 tmp_str = "" 789 DO j = 0, ldos_p(ildos)%ldos%maxl 790 DO i = -j, j 791 tmp_str(0, j, i) = sgf_symbol(0, j, i) 792 END DO 793 END DO 794 795 WRITE (UNIT=iw, FMT=fmtstr2) & 796 "# MO Eigenvalue [a.u.] Occupation", & 797 ((TRIM(tmp_str(0, il, im)), im=-il, il), il=0, ldos_p(ildos)%ldos%maxl) 798 DO imo = 1, nmo + nvirt 799 WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), & 800 (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=1, nsoset(ldos_p(ildos)%ldos%maxl)) 801 END DO 802 DEALLOCATE (tmp_str) 803 ELSE 804 WRITE (UNIT=iw, FMT=fmtstr2) & 805 "# MO Eigenvalue [a.u.] Occupation", & 806 (TRIM(l_sym(il)), il=0, ldos_p(ildos)%ldos%maxl) 807 DO imo = 1, nmo + nvirt 808 WRITE (UNIT=iw, FMT=fmtstr1) imo, eigenvalues(imo), occupation_numbers(imo), & 809 (ldos_p(ildos)%ldos%pdos_array(lshell, imo), lshell=0, ldos_p(ildos)%ldos%maxl) 810 END DO 811 END IF 812 END IF 813 CALL cp_print_key_finished_output(iw, logger, dft_section, & 814 "PRINT%PDOS") 815 END IF ! maxl>0 816 END DO ! ildos 817 818 ! write the pdos for the lists, each ona different file, 819 ! the filenames are indexed with the list number 820 DO ildos = 1, n_r_ldos 821 822 npoints = r_ldos_p(ildos)%ldos%npoints 823 CALL mp_sum(npoints, para_env%group) 824 CALL mp_sum(np_tot, para_env%group) 825 CALL mp_sum(r_ldos_p(ildos)%ldos%pdos_array, para_env%group) 826 IF (PRESENT(ispin)) THEN 827 IF (PRESENT(xas_mittle)) THEN 828 my_mittle = TRIM(xas_mittle)//TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos)) 829 ELSE 830 my_mittle = TRIM(spin(ispin))//"_r_list"//TRIM(r_ldos_index(ildos)) 831 END IF 832 my_spin = ispin 833 ELSE 834 my_mittle = "r_list"//TRIM(r_ldos_index(ildos)) 835 my_spin = 1 836 END IF 837 838 iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PDOS", & 839 extension=".pdos", file_position=my_pos, file_action=my_act, & 840 file_form="FORMATTED", middle_name=TRIM(my_mittle)) 841 IF (iw > 0) THEN 842 fmtstr1 = "(I8,2X,2F16.6, (2X,F16.8))" 843 fmtstr2 = "(A42, (10X,A8))" 844 845 WRITE (UNIT=iw, FMT="(A,I0,A,F12.6,F12.6,A,F12.6,A)") & 846 "# Projected DOS in real space, using ", npoints, & 847 " points of the grid, and eval in the range", r_ldos_p(ildos)%ldos%eval_range(1:2), & 848 " Hartree, E(Fermi) = ", e_fermi, " a.u." 849 WRITE (UNIT=iw, FMT="(A)") & 850 "# MO Eigenvalue [a.u.] Occupation LDOS" 851 DO imo = 1, nmo + nvirt 852 IF (r_ldos_p(ildos)%ldos%eval_range(1) <= eigenvalues(imo) .AND. & 853 r_ldos_p(ildos)%ldos%eval_range(2) >= eigenvalues(imo)) THEN 854 WRITE (UNIT=iw, FMT="(I8,2X,2F16.6,E20.10,E20.10)") imo, eigenvalues(imo), occupation_numbers(imo), & 855 r_ldos_p(ildos)%ldos%pdos_array(imo), r_ldos_p(ildos)%ldos%pdos_array(imo)*np_tot 856 END IF 857 END DO 858 859 END IF 860 CALL cp_print_key_finished_output(iw, logger, dft_section, & 861 "PRINT%PDOS") 862 END DO 863 864 ! deallocate local variables 865 DEALLOCATE (pdos_array) 866 DEALLOCATE (firstrow) 867 IF (do_ldos) THEN 868 DO ildos = 1, nldos 869 DEALLOCATE (ldos_p(ildos)%ldos%pdos_array) 870 DEALLOCATE (ldos_p(ildos)%ldos%list_index) 871 DEALLOCATE (ldos_p(ildos)%ldos) 872 END DO 873 DEALLOCATE (ldos_p) 874 DEALLOCATE (ldos_index) 875 END IF 876 IF (do_r_ldos) THEN 877 DO ildos = 1, n_r_ldos 878 DEALLOCATE (r_ldos_p(ildos)%ldos%index_grid_local) 879 DEALLOCATE (r_ldos_p(ildos)%ldos%pdos_array) 880 DEALLOCATE (r_ldos_p(ildos)%ldos%list_index) 881 IF (.NOT. read_r(1, ildos)) THEN 882 DEALLOCATE (r_ldos_p(ildos)%ldos%x_range) 883 END IF 884 IF (.NOT. read_r(2, ildos)) THEN 885 DEALLOCATE (r_ldos_p(ildos)%ldos%y_range) 886 END IF 887 IF (.NOT. read_r(3, ildos)) THEN 888 DEALLOCATE (r_ldos_p(ildos)%ldos%z_range) 889 END IF 890 IF (.NOT. read_r(4, ildos)) THEN 891 DEALLOCATE (r_ldos_p(ildos)%ldos%eval_range) 892 END IF 893 DEALLOCATE (r_ldos_p(ildos)%ldos) 894 END DO 895 DEALLOCATE (read_r) 896 DEALLOCATE (r_ldos_p) 897 DEALLOCATE (r_ldos_index) 898 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 899 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_g%pw) 900 END IF 901 IF (do_virt) THEN 902 DEALLOCATE (evals_virt) 903 CALL cp_fm_release(mo_virt) 904 CALL cp_fm_release(matrix_work) 905 DEALLOCATE (eigenvalues) 906 DEALLOCATE (occupation_numbers) 907 END IF 908 909 CALL timestop(handle) 910 911 END SUBROUTINE calculate_projected_dos 912 913! ************************************************************************************************** 914!> \brief Compute additional virtual states starting from the available MOS 915!> \param qs_env ... 916!> \param mo_set ... 917!> \param evals_virt ... 918!> \param mo_virt ... 919!> \param nvirt ... 920!> \param ispin ... 921!> \date 08.03.2008 922!> \par History: 923!> - 924!> \par Variables 925!> - 926!> - 927!> \author MI 928!> \version 1.0 929! ************************************************************************************************** 930 931 SUBROUTINE generate_virtual_mo(qs_env, mo_set, evals_virt, mo_virt, & 932 nvirt, ispin) 933 934 TYPE(qs_environment_type), POINTER :: qs_env 935 TYPE(mo_set_type), POINTER :: mo_set 936 REAL(KIND=dp), DIMENSION(:), POINTER :: evals_virt 937 TYPE(cp_fm_type), POINTER :: mo_virt 938 INTEGER, INTENT(IN) :: nvirt, ispin 939 940 CHARACTER(len=*), PARAMETER :: routineN = 'generate_virtual_mo', & 941 routineP = moduleN//':'//routineN 942 943 INTEGER :: nmo, nrow_global 944 TYPE(cp_blacs_env_type), POINTER :: context 945 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 946 TYPE(cp_fm_type), POINTER :: mo_coeff 947 TYPE(cp_para_env_type), POINTER :: para_env 948 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, s_matrix 949 TYPE(preconditioner_type), POINTER :: local_preconditioner 950 TYPE(scf_control_type), POINTER :: scf_control 951 952 NULLIFY (evals_virt, mo_virt) 953 ALLOCATE (evals_virt(nvirt)) 954 955 CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, & 956 scf_control=scf_control) 957 CALL get_mo_set(mo_set=mo_set, mo_coeff=mo_coeff, nmo=nmo) 958 CALL cp_fm_get_info(mo_coeff, context=context, para_env=para_env, & 959 nrow_global=nrow_global) 960 961 CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, & 962 nrow_global=nrow_global, ncol_global=nvirt) 963 CALL cp_fm_create(mo_virt, fm_struct_tmp, name="virtual") 964 CALL cp_fm_struct_release(fm_struct_tmp) 965 CALL cp_fm_init_random(mo_virt, nvirt) 966 967 NULLIFY (local_preconditioner) 968 969 CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, & 970 matrix_c_fm=mo_virt, matrix_orthogonal_space_fm=mo_coeff, & 971 eps_gradient=scf_control%eps_lumos, & 972 preconditioner=local_preconditioner, & 973 iter_max=scf_control%max_iter_lumos, & 974 size_ortho_space=nmo) 975 976 CALL calculate_subspace_eigenvalues(mo_virt, ks_matrix(ispin)%matrix, & 977 evals_virt) 978 979 END SUBROUTINE generate_virtual_mo 980 981END MODULE qs_pdos 982 983