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