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 STM image as post processing of an electronic 8!> structure calculation, 9!> \par History 10!> Started as a copy from the code in qs_scf_post 11!> \author Joost VandeVondele 7.2008, MI 02.2009 12! ************************************************************************************************** 13MODULE stm_images 14 USE cp_array_utils, ONLY: cp_1d_r_p_type 15 USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t 16 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale 17 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 18 cp_fm_struct_release,& 19 cp_fm_struct_type 20 USE cp_fm_types, ONLY: cp_fm_create,& 21 cp_fm_p_type,& 22 cp_fm_release,& 23 cp_fm_to_fm,& 24 cp_fm_type 25 USE cp_log_handling, ONLY: cp_get_default_logger,& 26 cp_logger_get_default_io_unit,& 27 cp_logger_type 28 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 29 cp_print_key_unit_nr 30 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube 31 USE dbcsr_api, ONLY: dbcsr_copy,& 32 dbcsr_deallocate_matrix,& 33 dbcsr_p_type,& 34 dbcsr_set,& 35 dbcsr_type 36 USE input_section_types, ONLY: section_get_ivals,& 37 section_vals_type,& 38 section_vals_val_get 39 USE kinds, ONLY: default_path_length,& 40 default_string_length,& 41 dp 42 USE particle_list_types, ONLY: particle_list_type 43 USE pw_env_types, ONLY: pw_env_get,& 44 pw_env_type 45 USE pw_pool_types, ONLY: pw_pool_create_pw,& 46 pw_pool_give_back_pw,& 47 pw_pool_p_type,& 48 pw_pool_type 49 USE pw_types, ONLY: COMPLEXDATA1D,& 50 REALDATA3D,& 51 REALSPACE,& 52 RECIPROCALSPACE,& 53 pw_p_type 54 USE qs_collocate_density, ONLY: calculate_rho_elec 55 USE qs_environment_types, ONLY: get_qs_env,& 56 qs_environment_type 57 USE qs_ks_types, ONLY: qs_ks_env_type 58 USE qs_mo_types, ONLY: get_mo_set,& 59 mo_set_p_type 60 USE qs_rho_types, ONLY: qs_rho_get,& 61 qs_rho_type 62#include "./base/base_uses.f90" 63 64 IMPLICIT NONE 65 PRIVATE 66 67 ! Global parameters 68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'stm_images' 69 PUBLIC :: th_stm_image 70 71CONTAINS 72! ************************************************************************************************** 73!> \brief Driver for the calculation of STM image, as post processing of a 74!> ground-state electronic structure calculation. 75!> \param qs_env ... 76!> \param stm_section ... 77!> \param particles ... 78!> \param unoccupied_orbs ... 79!> \param unoccupied_evals ... 80!> \param 81!> \par History 82!> 02.2009 Created [MI] 83!> \author MI 84!> \note 85!> The Tersoff-Hamman 86!> approximation is applied, occupied and a sufficient number of 87!> unoccupied eigenstates are needed (depending on the given Bias potential) 88!> and should be computed in advance. Unoccupied states are calculated 89!> before enetering this module when NLUMO =/ 0 90! ************************************************************************************************** 91 92 SUBROUTINE th_stm_image(qs_env, stm_section, particles, unoccupied_orbs, & 93 unoccupied_evals) 94 95 TYPE(qs_environment_type), POINTER :: qs_env 96 TYPE(section_vals_type), POINTER :: stm_section 97 TYPE(particle_list_type), POINTER :: particles 98 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: unoccupied_orbs 99 TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: unoccupied_evals 100 101 CHARACTER(len=*), PARAMETER :: routineN = 'th_stm_image', routineP = moduleN//':'//routineN 102 103 INTEGER :: handle, irep, ispin, n_rep, ndim, nmo, & 104 nspin, output_unit 105 INTEGER, DIMENSION(:), POINTER :: nadd_unocc, stm_th_torb 106 LOGICAL :: append_cube, use_ref_energy 107 REAL(KIND=dp) :: efermi, ref_energy 108 REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues, mo_occ, stm_biases 109 TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:) :: evals, occupation 110 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: mo_arrays 111 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 112 TYPE(cp_fm_type), POINTER :: mo_coeff 113 TYPE(cp_logger_type), POINTER :: logger 114 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao 115 TYPE(dbcsr_type), POINTER :: stm_density_ao 116 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 117 TYPE(pw_env_type), POINTER :: pw_env 118 TYPE(pw_p_type) :: wf_g, wf_r 119 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools 120 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 121 TYPE(qs_ks_env_type), POINTER :: ks_env 122 TYPE(qs_rho_type), POINTER :: rho 123 124 CALL timeset(routineN, handle) 125 logger => cp_get_default_logger() 126 output_unit = cp_logger_get_default_io_unit(logger) 127 128 NULLIFY (ks_env, mos, rho, rho_ao, pw_env, stm_th_torb, fm_struct_tmp) 129 NULLIFY (auxbas_pw_pool, pw_pools, stm_density_ao, mo_coeff) 130 131 CALL get_qs_env(qs_env, & 132 ks_env=ks_env, & 133 mos=mos, & 134 rho=rho, & 135 pw_env=pw_env) 136 137 CALL qs_rho_get(rho, rho_ao=rho_ao) 138 139 CALL section_vals_val_get(stm_section, "APPEND", l_val=append_cube) 140 CALL section_vals_val_get(stm_section, "BIAS", r_vals=stm_biases) 141 CALL section_vals_val_get(stm_section, "REF_ENERGY", r_val=ref_energy, explicit=use_ref_energy) 142 CALL section_vals_val_get(stm_section, "TH_TORB", n_rep_val=n_rep) 143 IF (n_rep == 0) THEN 144 ALLOCATE (stm_th_torb(1)) 145 stm_th_torb(1) = 0 146 ELSE 147 ALLOCATE (stm_th_torb(n_rep)) 148 DO irep = 1, n_rep 149 CALL section_vals_val_get(stm_section, "TH_TORB", & 150 i_rep_val=irep, i_val=stm_th_torb(irep)) 151 END DO 152 END IF 153 154 ALLOCATE (stm_density_ao) 155 CALL dbcsr_copy(stm_density_ao, rho_ao(1)%matrix, & 156 name="stm_density_ao") 157 158 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, & 159 pw_pools=pw_pools) 160 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 161 use_data=REALDATA3D, & 162 in_space=REALSPACE) 163 CALL pw_pool_create_pw(auxbas_pw_pool, wf_g%pw, & 164 use_data=COMPLEXDATA1D, & 165 in_space=RECIPROCALSPACE) 166 167 nspin = SIZE(mos, 1) 168 ALLOCATE (nadd_unocc(nspin)) 169 nadd_unocc = 0 170 IF (ASSOCIATED(unoccupied_orbs)) THEN 171 DO ispin = 1, nspin 172 nadd_unocc(ispin) = SIZE(unoccupied_evals(ispin)%array) 173 END DO 174 END IF 175 176 ALLOCATE (mo_arrays(nspin)) 177 ALLOCATE (evals(nspin)) 178 ALLOCATE (occupation(nspin)) 179 DO ispin = 1, nspin 180 IF (nadd_unocc(ispin) == 0) THEN 181 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, & 182 eigenvalues=mo_eigenvalues, nmo=nmo, mu=efermi, occupation_numbers=mo_occ) 183 mo_arrays(ispin)%matrix => mo_coeff 184 evals(ispin)%array => mo_eigenvalues 185 occupation(ispin)%array => mo_occ 186 ELSE 187 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff, & 188 eigenvalues=mo_eigenvalues, nmo=nmo, mu=efermi, occupation_numbers=mo_occ) 189 ndim = nmo + nadd_unocc(ispin) 190 ALLOCATE (evals(ispin)%array(ndim)) 191 evals(ispin)%array(1:nmo) = mo_eigenvalues(1:nmo) 192 evals(ispin)%array(1 + nmo:ndim) = unoccupied_evals(ispin)%array(1:nadd_unocc(ispin)) 193 ALLOCATE (occupation(ispin)%array(ndim)) 194 occupation(ispin)%array(1:nmo) = mo_occ(1:nmo) 195 occupation(ispin)%array(1 + nmo:ndim) = 0.0_dp 196 CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=ndim, & 197 template_fmstruct=mo_coeff%matrix_struct) 198 CALL cp_fm_create(mo_arrays(ispin)%matrix, fm_struct_tmp, name="mo_arrays") 199 CALL cp_fm_struct_release(fm_struct_tmp) 200 CALL cp_fm_to_fm(mo_coeff, mo_arrays(ispin)%matrix, nmo, 1, 1) 201 CALL cp_fm_to_fm(unoccupied_orbs(ispin)%matrix, mo_arrays(ispin)%matrix, & 202 nadd_unocc(ispin), 1, nmo + 1) 203 END IF 204 ENDDO 205 IF (use_ref_energy) efermi = ref_energy 206 207 CALL stm_cubes(ks_env, stm_section, stm_density_ao, wf_r, wf_g, mo_arrays, evals, & 208 occupation, efermi, stm_biases, stm_th_torb, particles, & 209 output_unit, append_cube) 210 DO ispin = 1, nspin 211 IF (nadd_unocc(ispin) > 0) THEN 212 DEALLOCATE (evals(ispin)%array) 213 DEALLOCATE (occupation(ispin)%array) 214 CALL cp_fm_release(mo_arrays(ispin)%matrix) 215 END IF 216 END DO 217 DEALLOCATE (mo_arrays) 218 DEALLOCATE (evals) 219 DEALLOCATE (occupation) 220 221 CALL dbcsr_deallocate_matrix(stm_density_ao) 222 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 223 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_g%pw) 224 225 DEALLOCATE (stm_th_torb) 226 DEALLOCATE (nadd_unocc) 227 228 CALL timestop(handle) 229 230 END SUBROUTINE th_stm_image 231 232! ************************************************************************************************** 233!> \brief computes a simple approximation to the tunneling current for STM 234!> \param ks_env ... 235!> \param stm_section ... 236!> \param stm_density_ao ... 237!> \param wf_r ... 238!> \param wf_g ... 239!> \param mo_arrays ... 240!> \param evals ... 241!> \param occupation ... 242!> \param efermi ... 243!> \param stm_biases ... 244!> \param stm_th_torb ... 245!> \param particles ... 246!> \param output_unit ... 247!> \param append_cube ... 248!> \param 249!> \par History 250!> 7.2008 Created [Joost VandeVondele] 251!> 07.2009 modified MI 252!> \author Joost VandeVondele 253!> \note 254!> requires the MOs that are passed to be eigenstates, and energy ordered 255! ************************************************************************************************** 256 SUBROUTINE stm_cubes(ks_env, stm_section, stm_density_ao, wf_r, wf_g, mo_arrays, evals, & 257 occupation, efermi, stm_biases, stm_th_torb, particles, & 258 output_unit, append_cube) 259 260 TYPE(qs_ks_env_type), POINTER :: ks_env 261 TYPE(section_vals_type), POINTER :: stm_section 262 TYPE(dbcsr_type), POINTER :: stm_density_ao 263 TYPE(pw_p_type) :: wf_r, wf_g 264 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(IN) :: mo_arrays 265 TYPE(cp_1d_r_p_type), DIMENSION(:), INTENT(IN) :: evals, occupation 266 REAL(KIND=dp) :: efermi 267 REAL(KIND=dp), DIMENSION(:), POINTER :: stm_biases 268 INTEGER, DIMENSION(:), POINTER :: stm_th_torb 269 TYPE(particle_list_type), POINTER :: particles 270 INTEGER, INTENT(IN) :: output_unit 271 LOGICAL, INTENT(IN) :: append_cube 272 273 CHARACTER(LEN=*), DIMENSION(0:9), PARAMETER :: & 274 torb_string = (/" s", " px", " py", " pz", "dxy", "dyz", "dzx", "dx2", "dy2", "dz2"/) 275 CHARACTER(len=*), PARAMETER :: routineN = 'stm_cubes', routineP = moduleN//':'//routineN 276 277 CHARACTER(LEN=default_path_length) :: filename 278 CHARACTER(LEN=default_string_length) :: my_pos, oname, title 279 INTEGER :: handle, i, ibias, imo, iorb, ispin, & 280 istates, nmo, nspin, nstates(2), & 281 state_start(2), unit_nr 282 LOGICAL :: mpi_io 283 REAL(KIND=dp) :: alpha, total_rho 284 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: occ_tot 285 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp 286 TYPE(cp_fm_type), POINTER :: matrix_v, matrix_vf 287 TYPE(cp_logger_type), POINTER :: logger 288 289 CALL timeset(routineN, handle) 290 291 logger => cp_get_default_logger() 292 NULLIFY (fm_struct_tmp) 293 294 nspin = SIZE(mo_arrays) 295 296 IF (output_unit > 0) WRITE (output_unit, '(T2,A)') "" 297 IF (output_unit > 0) WRITE (output_unit, '(T2,A,F12.6, A)') "STM : Reference energy ", efermi, " a.u. " 298 DO ibias = 1, SIZE(stm_biases) 299 300 IF (output_unit > 0) WRITE (output_unit, '(T2,A)') "" 301 IF (output_unit > 0) WRITE (output_unit, '(T2,A,F16.6)') & 302 "Preparing for STM image at bias [a.u.] ", stm_biases(ibias) 303 304 istates = 0 305 nstates = 0 306 state_start = 0 307 DO ispin = 1, nspin 308 IF (stm_biases(ibias) < 0.0_dp) THEN 309 nmo = SIZE(evals(ispin)%array) 310 DO imo = 1, nmo 311 IF (evals(ispin)%array(imo) > (efermi + stm_biases(ibias)) .AND. & 312 evals(ispin)%array(imo) <= efermi) THEN 313 IF (nstates(ispin) == 0) state_start(ispin) = imo 314 nstates(ispin) = nstates(ispin) + 1 315 END IF 316 END DO 317 IF ((output_unit > 0) .AND. evals(ispin)%array(1) > efermi + stm_biases(ibias)) & 318 WRITE (output_unit, '(T4,A)') "Warning: EFermi+bias below lowest computed occupied MO" 319 ELSE 320 nmo = SIZE(evals(ispin)%array) 321 DO imo = 1, nmo 322 IF (evals(ispin)%array(imo) <= (efermi + stm_biases(ibias)) .AND. & 323 evals(ispin)%array(imo) > efermi) THEN 324 IF (nstates(ispin) == 0) state_start(ispin) = imo 325 nstates(ispin) = nstates(ispin) + 1 326 END IF 327 END DO 328 IF ((output_unit > 0) .AND. evals(ispin)%array(nmo) < efermi + stm_biases(ibias)) & 329 WRITE (output_unit, '(T4,A)') "Warning: E-Fermi+bias above highest computed unoccupied MO" 330 ENDIF 331 istates = istates + nstates(ispin) 332 ENDDO 333 IF ((output_unit > 0)) WRITE (output_unit, '(T4,A,I0,A)') "Using a total of ", istates, " states" 334 IF (istates == 0) CYCLE 335 336 CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=istates, & 337 template_fmstruct=mo_arrays(1)%matrix%matrix_struct) 338 CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v") 339 CALL cp_fm_create(matrix_vf, fm_struct_tmp, name="matrix_vf") 340 CALL cp_fm_struct_release(fm_struct_tmp) 341 342 ALLOCATE (occ_tot(istates)) 343 344 ! we sum both alpha and beta electrons together for this density of states 345 istates = 0 346 alpha = 1.0_dp 347 IF (nspin == 1) alpha = 2.0_dp 348 DO ispin = 1, nspin 349 CALL cp_fm_to_fm(mo_arrays(ispin)%matrix, matrix_v, nstates(ispin), state_start(ispin), istates + 1) 350 CALL cp_fm_to_fm(mo_arrays(ispin)%matrix, matrix_vf, nstates(ispin), state_start(ispin), istates + 1) 351 IF (stm_biases(ibias) < 0.0_dp) THEN 352 occ_tot(istates + 1:istates + nstates(ispin)) = & 353 occupation(ispin)%array(state_start(ispin):state_start(ispin) - 1 + nstates(ispin)) 354 ELSE 355 occ_tot(istates + 1:istates + nstates(ispin)) = & 356 alpha - occupation(ispin)%array(state_start(ispin):state_start(ispin) - 1 + nstates(ispin)) 357 END IF 358 istates = istates + nstates(ispin) 359 ENDDO 360 361 CALL cp_fm_column_scale(matrix_vf, occ_tot(1:istates)) 362 alpha = 1.0_dp 363 364 CALL dbcsr_set(stm_density_ao, 0.0_dp) 365 CALL cp_dbcsr_plus_fm_fm_t(stm_density_ao, matrix_v=matrix_v, matrix_g=matrix_vf, ncol=istates, & 366 alpha=alpha) 367 368 DO i = 1, SIZE(stm_th_torb) 369 iorb = stm_th_torb(i) 370 CALL calculate_rho_elec(matrix_p=stm_density_ao, & 371 rho=wf_r, rho_gspace=wf_g, total_rho=total_rho, & 372 ks_env=ks_env, der_type=iorb) 373 374 oname = torb_string(iorb) 375! fname = "STM_"//TRIM(torb_string(iorb)) 376 WRITE (filename, '(a4,I2.2,a1,I5.5)') "STM_d", iorb, "_", ibias 377 my_pos = "REWIND" 378 IF (append_cube) THEN 379 my_pos = "APPEND" 380 END IF 381 382 mpi_io = .TRUE. 383 unit_nr = cp_print_key_unit_nr(logger, stm_section, extension=".cube", & 384 middle_name=TRIM(filename), file_position=my_pos, file_action="WRITE", & 385 log_filename=.FALSE., mpi_io=mpi_io) 386 WRITE (title, '(A,I0,A,I0,A,F16.8)') "STM cube ", ibias, " wfn deriv. ", iorb, " at bias ", stm_biases(ibias) 387 CALL cp_pw_to_cube(wf_r%pw, unit_nr, title, particles=particles, & 388 stride=section_get_ivals(stm_section, "STRIDE"), zero_tails=.TRUE., & 389 mpi_io=mpi_io) 390 391 CALL cp_print_key_finished_output(unit_nr, logger, stm_section, mpi_io=mpi_io) 392 END DO 393 394 CALL cp_fm_release(matrix_v) 395 CALL cp_fm_release(matrix_vf) 396 DEALLOCATE (occ_tot) 397 398 ENDDO 399 400 CALL timestop(handle) 401 402 END SUBROUTINE stm_cubes 403 404END MODULE stm_images 405