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 CP2K+SMEAGOL interface.
10 : !> \author Sergey Chulkov
11 : !> \author Christian Ahart
12 : !> \author Clotilde Cucinotta
13 : ! **************************************************************************************************
14 : MODULE smeagol_interface
15 : USE bibliography, ONLY: Ahart2024, &
16 : Bailey2006, &
17 : cite_reference
18 : USE cell_types, ONLY: cell_type, &
19 : real_to_scaled, &
20 : scaled_to_real
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
23 : USE cp_files, ONLY: close_file, &
24 : open_file
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_get_iter_level_by_name, &
29 : cp_get_iter_nr
30 : USE cp_dbcsr_api, ONLY: dbcsr_copy, &
31 : dbcsr_p_type, &
32 : dbcsr_set
33 : USE input_constants, ONLY: smeagol_bulklead_left, &
34 : smeagol_bulklead_leftright, &
35 : smeagol_bulklead_right, &
36 : smeagol_runtype_bulktransport, &
37 : smeagol_runtype_emtransport
38 : USE kahan_sum, ONLY: accurate_dot_product, &
39 : accurate_sum
40 : USE kinds, ONLY: dp, &
41 : int_8
42 : USE kpoint_types, ONLY: get_kpoint_info, &
43 : kpoint_type
44 : USE message_passing, ONLY: mp_para_env_type
45 : #if defined(__SMEAGOL)
46 : USE mmpi_negf, ONLY: create_communicators_negf, &
47 : destroy_communicators_negf
48 : USE mnegf_interface, ONLY: negf_interface
49 : USE negfmod, ONLY: smeagolglobal_em_nas => em_nas, &
50 : smeagolglobal_em_nau => em_nau, &
51 : smeagolglobal_em_nso => em_nso, &
52 : smeagolglobal_em_nuo => em_nuo, &
53 : smeagolglobal_negfon => negfon
54 : #endif
55 : USE physcon, ONLY: bohr, &
56 : evolt
57 : USE pw_grid_types, ONLY: pw_grid_type
58 : USE pw_types, ONLY: pw_r3d_rs_type
59 : USE qs_matrix_w, ONLY: compute_matrix_w
60 : USE qs_energy_types, ONLY: qs_energy_type
61 : USE qs_environment_types, ONLY: get_qs_env, &
62 : qs_environment_type
63 : USE qs_ks_types, ONLY: qs_ks_env_type, &
64 : set_ks_env
65 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
66 : USE qs_rho_types, ONLY: qs_rho_get, &
67 : qs_rho_type
68 : USE qs_subsys_types, ONLY: qs_subsys_get, &
69 : qs_subsys_type
70 : USE scf_control_types, ONLY: scf_control_type
71 : USE smeagol_control_types, ONLY: smeagol_control_type
72 : USE smeagol_emtoptions, ONLY: ReadOptionsNEGF_DFT, &
73 : emtrans_deallocate_global_arrays, &
74 : emtrans_options, &
75 : reademtr
76 : USE smeagol_matrix_utils, ONLY: convert_dbcsr_to_distributed_siesta, &
77 : convert_distributed_siesta_to_dbcsr, &
78 : siesta_distrib_csc_struct_type, &
79 : siesta_struct_create, &
80 : siesta_struct_release
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 : PRIVATE
85 :
86 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'smeagol_interface'
87 :
88 : PUBLIC :: run_smeagol_bulktrans, run_smeagol_emtrans
89 : PUBLIC :: smeagol_shift_v_hartree
90 :
91 : CONTAINS
92 :
93 : ! **************************************************************************************************
94 : !> \brief Save overlap, Kohn-Sham, electron density, and energy-density matrices of semi-infinite
95 : !> electrodes in SIESTA format.
96 : !> \param qs_env QuickStep environment
97 : ! **************************************************************************************************
98 17637 : SUBROUTINE run_smeagol_bulktrans(qs_env)
99 : TYPE(qs_environment_type), POINTER :: qs_env
100 :
101 : CHARACTER(LEN=*), PARAMETER :: routineN = 'run_smeagol_bulktrans'
102 :
103 : CHARACTER(len=32) :: hst_fmt
104 : INTEGER :: funit, handle, img, ispin, lead_label, &
105 : log_unit, nimages, nspin
106 : INTEGER(kind=int_8) :: ielem
107 : INTEGER, DIMENSION(2) :: max_ij_cell_image
108 17637 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
109 : LOGICAL :: do_kpoints, has_unit_metric, not_regtest
110 : REAL(kind=dp) :: H_to_Ry
111 17637 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: matrix_siesta_1d, nelectrons
112 17637 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: matrix_siesta_2d
113 17637 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: matrix_kp_generic
114 17637 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, matrix_w_kp, &
115 17637 : rho_ao_kp
116 : TYPE(dft_control_type), POINTER :: dft_control
117 : TYPE(kpoint_type), POINTER :: kpoints
118 : TYPE(mp_para_env_type), POINTER :: para_env
119 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
120 17637 : POINTER :: sab_nl
121 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
122 : TYPE(qs_energy_type), POINTER :: energy
123 : TYPE(qs_ks_env_type), POINTER :: ks_env
124 : TYPE(qs_rho_type), POINTER :: rho
125 : TYPE(qs_subsys_type), POINTER :: subsys
126 : TYPE(scf_control_type), POINTER :: scf_control
127 17637 : TYPE(siesta_distrib_csc_struct_type) :: siesta_struct
128 : TYPE(smeagol_control_type), POINTER :: smeagol_control
129 :
130 17637 : CALL get_qs_env(qs_env, dft_control=dft_control)
131 17637 : smeagol_control => dft_control%smeagol_control
132 17637 : IF (.NOT. (smeagol_control%smeagol_enabled .AND. smeagol_control%run_type == smeagol_runtype_bulktransport)) RETURN
133 :
134 4 : CALL timeset(routineN, handle)
135 :
136 4 : log_unit = cp_logger_get_default_io_unit()
137 4 : H_to_Ry = smeagol_control%to_smeagol_energy_units
138 4 : not_regtest = .NOT. smeagol_control%do_regtest
139 :
140 4 : lead_label = smeagol_control%lead_label
141 4 : nspin = dft_control%nspins
142 :
143 4 : NULLIFY (v_hartree_rspace)
144 : CALL get_qs_env(qs_env, energy=energy, para_env=para_env, subsys=subsys, &
145 : scf_control=scf_control, &
146 : do_kpoints=do_kpoints, kpoints=kpoints, &
147 : matrix_ks_kp=matrix_ks_kp, matrix_s_kp=matrix_s_kp, &
148 4 : rho=rho, sab_orb=sab_nl, v_hartree_rspace=v_hartree_rspace)
149 4 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
150 :
151 4 : IF (not_regtest) THEN
152 : ! save average electrostatic potential of the electrode along transport direction
153 2 : CALL write_average_hartree_potential(v_hartree_rspace, smeagol_control%project_name)
154 : END IF
155 :
156 4 : IF (log_unit > 0) THEN
157 2 : WRITE (log_unit, '(/,T2,A,T61,ES20.10E2)') "SMEAGOL| E_FERMI [a.u.] = ", energy%efermi
158 : END IF
159 :
160 4 : IF (do_kpoints) THEN
161 4 : nimages = dft_control%nimages
162 4 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
163 : ELSE
164 0 : nimages = 1
165 0 : ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
166 0 : cell_to_index(0, 0, 0) = 1
167 : ! We do need at least two cell images along transport direction.
168 0 : CPABORT("Please enable k-points")
169 : END IF
170 :
171 : ! largest index of cell images along i and j cell vectors
172 : ! e.g., (2,0) in case of 5 cell images (0,0), (1,0), (-1,0), (2,0), and (-2,0)
173 : ! -1 means to use all available cell images along a particular cell vector.
174 12 : max_ij_cell_image(:) = -1
175 12 : DO img = 1, 2
176 12 : IF (smeagol_control%n_cell_images(img) > 0) THEN
177 0 : max_ij_cell_image(img) = smeagol_control%n_cell_images(img)/2
178 : END IF
179 : END DO
180 :
181 288 : ALLOCATE (matrix_kp_generic(nimages))
182 :
183 : ! compute energy-density (W) matrix. We may need it later to calculate NEGF forces
184 4 : CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric)
185 4 : IF (.NOT. has_unit_metric) THEN
186 4 : CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
187 4 : IF (.NOT. ASSOCIATED(matrix_w_kp)) THEN
188 4 : NULLIFY (matrix_w_kp)
189 4 : CALL get_qs_env(qs_env, ks_env=ks_env)
190 4 : CALL dbcsr_allocate_matrix_set(matrix_w_kp, nspin, nimages)
191 8 : DO ispin = 1, nspin
192 284 : DO img = 1, nimages
193 276 : ALLOCATE (matrix_w_kp(ispin, img)%matrix)
194 276 : CALL dbcsr_copy(matrix_w_kp(ispin, img)%matrix, matrix_s_kp(1, 1)%matrix, name="W MATRIX")
195 280 : CALL dbcsr_set(matrix_w_kp(ispin, img)%matrix, 0.0_dp)
196 : END DO
197 : END DO
198 4 : CALL set_ks_env(ks_env, matrix_w_kp=matrix_w_kp)
199 : END IF
200 : END IF
201 4 : CALL compute_matrix_w(qs_env, calc_forces=.TRUE.)
202 :
203 : ! obtain the sparsity pattern of the overlap matrix
204 280 : DO img = 1, nimages
205 280 : matrix_kp_generic(img)%matrix => matrix_s_kp(1, img)%matrix
206 : END DO
207 :
208 : CALL siesta_struct_create(siesta_struct, matrix_kp_generic, subsys, cell_to_index, &
209 4 : sab_nl, para_env, max_ij_cell_image, do_merge=.FALSE., gather_root=para_env%source)
210 :
211 : ! write 'bulklft.DAT' and 'bulkrgt.DAT' files
212 4 : funit = -1
213 4 : IF (not_regtest .AND. para_env%mepos == para_env%source) THEN
214 : ! I/O process
215 1 : IF (lead_label == smeagol_bulklead_left .OR. lead_label == smeagol_bulklead_leftright) THEN
216 1 : CALL open_file("bulklft.DAT", file_status="REPLACE", file_form="FORMATTED", file_action="WRITE", unit_number=funit)
217 :
218 : CALL write_bulk_dat_file(funit, siesta_struct, smeagol_control%project_name, nspin, &
219 : energy%efermi, scf_control%smear%ELECTRONIC_TEMPERATURE, &
220 1 : H_to_Ry, do_kpoints, max_ij_cell_image)
221 :
222 1 : CALL close_file(funit)
223 : END IF
224 :
225 1 : IF (lead_label == smeagol_bulklead_right .OR. lead_label == smeagol_bulklead_leftright) THEN
226 1 : CALL open_file("bulkrgt.DAT", file_status="REPLACE", file_form="FORMATTED", file_action="WRITE", unit_number=funit)
227 :
228 : CALL write_bulk_dat_file(funit, siesta_struct, smeagol_control%project_name, nspin, &
229 : energy%efermi, scf_control%smear%ELECTRONIC_TEMPERATURE, &
230 1 : H_to_Ry, do_kpoints, max_ij_cell_image)
231 :
232 1 : CALL close_file(funit)
233 : END IF
234 : END IF
235 :
236 : ! write project_name.HST file
237 4 : funit = -1
238 4 : IF (para_env%mepos == para_env%source) THEN
239 6 : ALLOCATE (matrix_siesta_1d(siesta_struct%n_nonzero_elements))
240 8 : ALLOCATE (matrix_siesta_2d(siesta_struct%n_nonzero_elements, nspin))
241 :
242 2 : IF (not_regtest) THEN
243 : CALL open_file(TRIM(smeagol_control%project_name)//".HST", &
244 1 : file_status="REPLACE", file_form="FORMATTED", file_action="WRITE", unit_number=funit)
245 : END IF
246 : ELSE
247 2 : ALLOCATE (matrix_siesta_1d(1))
248 6 : ALLOCATE (matrix_siesta_2d(1, nspin))
249 : END IF
250 :
251 : !DO img = 1, nimages
252 : ! matrix_kp_generic(img)%matrix => matrix_s_kp(1, img)%matrix
253 : !END DO
254 4 : CALL convert_dbcsr_to_distributed_siesta(matrix_siesta_1d, matrix_kp_generic, siesta_struct, para_env)
255 :
256 8 : DO ispin = 1, nspin
257 280 : DO img = 1, nimages
258 280 : matrix_kp_generic(img)%matrix => matrix_ks_kp(ispin, img)%matrix
259 : END DO
260 8 : CALL convert_dbcsr_to_distributed_siesta(matrix_siesta_2d(:, ispin), matrix_kp_generic, siesta_struct, para_env)
261 : END DO
262 : ! As SIESTA's default energy unit is Rydberg, scale the KS-matrix
263 2318878 : matrix_siesta_2d(:, :) = H_to_Ry*matrix_siesta_2d(:, :)
264 :
265 4 : IF (funit > 0) THEN ! not_regtest .AND. para_env%mepos == para_env%source
266 1 : WRITE (hst_fmt, '(A,I0,A)') "(", nspin, "ES26.17E3)"
267 1227637 : DO ielem = 1, siesta_struct%n_nonzero_elements
268 1227636 : WRITE (funit, '(ES26.17E3)') matrix_siesta_1d(ielem)
269 1227637 : WRITE (funit, hst_fmt) matrix_siesta_2d(ielem, :)
270 : END DO
271 :
272 1 : CALL close_file(funit)
273 : END IF
274 :
275 : ! write density matrix
276 8 : DO ispin = 1, nspin
277 280 : DO img = 1, nimages
278 280 : matrix_kp_generic(img)%matrix => rho_ao_kp(ispin, img)%matrix
279 : END DO
280 8 : CALL convert_dbcsr_to_distributed_siesta(matrix_siesta_2d(:, ispin), matrix_kp_generic, siesta_struct, para_env)
281 : END DO
282 :
283 4 : IF (para_env%mepos == para_env%source) THEN
284 6 : ALLOCATE (nelectrons(nspin))
285 4 : DO ispin = 1, nspin
286 4 : nelectrons(ispin) = accurate_dot_product(matrix_siesta_2d(:, ispin), matrix_siesta_1d)
287 : END DO
288 :
289 2 : CPASSERT(log_unit > 0)
290 2 : IF (nspin > 1) THEN
291 0 : WRITE (log_unit, '(T2,A,T61,F20.10)') "SMEAGOL| Number of alpha-spin electrons: ", nelectrons(1)
292 0 : WRITE (log_unit, '(T2,A,T61,F20.10)') "SMEAGOL| Number of beta-spin electrons: ", nelectrons(2)
293 : ELSE
294 2 : WRITE (log_unit, '(T2,A,T61,F20.10)') "SMEAGOL| Number of electrons: ", nelectrons(1)
295 : END IF
296 2 : DEALLOCATE (nelectrons)
297 :
298 2 : IF (not_regtest) THEN
299 : CALL open_file(TRIM(smeagol_control%project_name)//".DM", &
300 1 : file_status="REPLACE", file_form="UNFORMATTED", file_action="WRITE", unit_number=funit)
301 :
302 1 : CALL write_bulk_dm_file(funit, siesta_struct, nspin, matrix_siesta_2d)
303 :
304 1 : CALL close_file(funit)
305 : END IF
306 : END IF
307 :
308 4 : CALL get_qs_env(qs_env, matrix_w_kp=matrix_w_kp)
309 4 : IF (ASSOCIATED(matrix_w_kp)) THEN
310 : ! write energy density matrix
311 8 : DO ispin = 1, nspin
312 280 : DO img = 1, nimages
313 280 : matrix_kp_generic(img)%matrix => matrix_w_kp(ispin, img)%matrix
314 : END DO
315 : CALL convert_dbcsr_to_distributed_siesta(matrix_siesta_2d(:, ispin), matrix_kp_generic, &
316 8 : siesta_struct, para_env)
317 : END DO
318 :
319 4 : IF (not_regtest .AND. para_env%mepos == para_env%source) THEN
320 : CALL open_file(TRIM(smeagol_control%project_name)//".EDM", &
321 1 : file_status="REPLACE", file_form="UNFORMATTED", file_action="WRITE", unit_number=funit)
322 :
323 1 : CALL write_bulk_dm_file(funit, siesta_struct, nspin, matrix_siesta_2d)
324 :
325 1 : CALL close_file(funit)
326 : END IF
327 : END IF
328 :
329 4 : DEALLOCATE (matrix_siesta_2d)
330 4 : DEALLOCATE (matrix_siesta_1d)
331 :
332 4 : DEALLOCATE (matrix_kp_generic)
333 4 : IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
334 :
335 4 : CALL siesta_struct_release(siesta_struct)
336 4 : CALL timestop(handle)
337 35278 : END SUBROUTINE run_smeagol_bulktrans
338 :
339 : ! **************************************************************************************************
340 : !> \brief Write a sparse matrix structure file in SIESTA format. Should be called on I/O MPI process only.
341 : !> \param funit file to write
342 : !> \param siesta_struct sparse matrix structure in SIESTA format
343 : !> \param system_label SMEAGOL project name (first components of file names, i.e. system_label.HST)
344 : !> \param nspin number of spin components
345 : !> \param EFermi Fermi level
346 : !> \param temperature electronic temperature
347 : !> \param H_to_Ry Hartree to Rydberg scale factor
348 : !> \param do_kpoints whether to perform k-point calculation. Should always be enabled as
349 : !> SMEAGOL expects at least 3 cell replicas along the transport direction
350 : !> \param max_ij_cell_image maximum cell-replica indices along i- and j- lattice vectors
351 : !> (perpendicular the transport direction)
352 : ! **************************************************************************************************
353 2 : SUBROUTINE write_bulk_dat_file(funit, siesta_struct, system_label, nspin, EFermi, temperature, &
354 : H_to_Ry, do_kpoints, max_ij_cell_image)
355 : INTEGER, INTENT(in) :: funit
356 : TYPE(siesta_distrib_csc_struct_type), INTENT(in) :: siesta_struct
357 : CHARACTER(len=*), INTENT(in) :: system_label
358 : INTEGER, INTENT(in) :: nspin
359 : REAL(kind=dp), INTENT(in) :: EFermi, temperature, H_to_Ry
360 : LOGICAL, INTENT(in) :: do_kpoints
361 : INTEGER, DIMENSION(2), INTENT(in) :: max_ij_cell_image
362 :
363 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bulk_dat_file'
364 :
365 : INTEGER :: handle, icol, irow, nao_supercell, &
366 : nao_unitcell
367 : INTEGER, DIMENSION(2) :: ncells_siesta
368 :
369 2 : CALL timeset(routineN, handle)
370 :
371 : ! ++ header
372 2 : nao_unitcell = siesta_struct%nrows
373 2 : nao_supercell = siesta_struct%ncols
374 6 : ncells_siesta(1:2) = 2*max_ij_cell_image(1:2) + 1
375 :
376 : ! SMEAGOL expects Temperature and Fermi energy in Rydberg energy units, not in Hartree energy units.
377 : ! It is why these values are doubled.
378 :
379 : WRITE (funit, '(1X,A20,3I12,2ES26.17E3,3I12)') &
380 2 : system_label, nao_unitcell, nspin, siesta_struct%n_nonzero_elements, &
381 4 : H_to_Ry*EFermi, H_to_Ry*temperature, ncells_siesta(1:2), nao_supercell
382 :
383 : ! ++ number of non-zero matrix elements on each row and
384 : ! the index of the first non-zero matrix element on this row -1 in 1D data array.
385 650 : DO irow = 1, nao_unitcell
386 650 : WRITE (funit, '(2I12)') siesta_struct%n_nonzero_cols(irow), siesta_struct%row_offset(irow)
387 : END DO
388 :
389 48602 : DO icol = 1, nao_supercell
390 48602 : WRITE (funit, '(I12)') siesta_struct%indxuo(icol)
391 : END DO
392 :
393 : ! ++ column indices of non-zero matrix elements
394 650 : DO irow = 1, nao_unitcell
395 2455922 : DO icol = 1, siesta_struct%n_nonzero_cols(irow)
396 2455272 : WRITE (funit, '(I12)') siesta_struct%col_index(siesta_struct%row_offset(irow) + icol)
397 :
398 2455920 : IF (do_kpoints) THEN
399 2455272 : WRITE (funit, '(F21.16,5X,F21.16,5X,F21.16)') siesta_struct%xij(:, siesta_struct%row_offset(irow) + icol)
400 : END IF
401 : END DO
402 : END DO
403 :
404 2 : CALL timestop(handle)
405 2 : END SUBROUTINE write_bulk_dat_file
406 :
407 : ! **************************************************************************************************
408 : !> \brief Write an (energy-) density matrix. Should be called on I/O MPI process only.
409 : !> \param funit file to write
410 : !> \param siesta_struct sparse matrix structure in SIESTA format
411 : !> \param nspin number of spin components
412 : !> \param matrix_siesta non-zero matrix elements (1:siesta_struct%n_nonzero_elements, 1:nspin)
413 : ! **************************************************************************************************
414 2 : SUBROUTINE write_bulk_dm_file(funit, siesta_struct, nspin, matrix_siesta)
415 : INTEGER, INTENT(in) :: funit
416 : TYPE(siesta_distrib_csc_struct_type), INTENT(in) :: siesta_struct
417 : INTEGER, INTENT(in) :: nspin
418 : REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: matrix_siesta
419 :
420 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_bulk_dm_file'
421 :
422 : INTEGER :: handle, irow, ispin
423 :
424 2 : CALL timeset(routineN, handle)
425 :
426 : ! ++ number of compressed rows, number of spin components
427 2 : WRITE (funit) siesta_struct%nrows, nspin
428 :
429 : ! ++ number of non-zero matrix elements on each compressed row.
430 : ! The sparsity pattern for alpha- and beta-spin density matrices are identical
431 2 : WRITE (funit) siesta_struct%n_nonzero_cols
432 :
433 : ! ++ column indices of non-zero matrix elements
434 650 : DO irow = 1, siesta_struct%nrows
435 : WRITE (funit) siesta_struct%col_index(siesta_struct%row_offset(irow) + 1: &
436 650 : siesta_struct%row_offset(irow) + siesta_struct%n_nonzero_cols(irow))
437 : END DO
438 :
439 : ! ++ non-zero matrix blocks
440 4 : DO ispin = 1, nspin
441 652 : DO irow = 1, siesta_struct%nrows
442 : WRITE (funit) matrix_siesta(siesta_struct%row_offset(irow) + 1: &
443 650 : siesta_struct%row_offset(irow) + siesta_struct%n_nonzero_cols(irow), ispin)
444 : END DO
445 : END DO
446 :
447 2 : CALL timestop(handle)
448 2 : END SUBROUTINE write_bulk_dm_file
449 :
450 : ! **************************************************************************************************
451 : !> \brief Write the average value of Hartree potential along transport direction.
452 : !> SMEAGOL assumes that the transport direction coincides with z-axis.
453 : !> \param v_hartree_rspace Hartree potential on a real-space 3-D grid
454 : !> \param project_name SMEAGOL project name
455 : !> \note This routine assumes that the lattice vector 'k' coincides with z-axis
456 : ! **************************************************************************************************
457 2 : SUBROUTINE write_average_hartree_potential(v_hartree_rspace, project_name)
458 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
459 : CHARACTER(len=*), INTENT(in) :: project_name
460 :
461 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_average_hartree_potential'
462 :
463 : INTEGER :: funit, handle, iz, lz, uz
464 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: v_hartree_z_average
465 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
466 2 : POINTER :: cr3d
467 : TYPE(pw_grid_type), POINTER :: pw_grid
468 :
469 2 : CALL timeset(routineN, handle)
470 :
471 2 : pw_grid => v_hartree_rspace%pw_grid
472 2 : cr3d => v_hartree_rspace%array
473 :
474 6 : ALLOCATE (v_hartree_z_average(pw_grid%bounds(1, 3):pw_grid%bounds(2, 3)))
475 178 : v_hartree_z_average(:) = 0.0_dp
476 :
477 2 : lz = pw_grid%bounds_local(1, 3)
478 2 : uz = pw_grid%bounds_local(2, 3)
479 :
480 : ! save average electrostatic potential
481 178 : DO iz = lz, uz
482 178 : v_hartree_z_average(iz) = accurate_sum(cr3d(:, :, iz))
483 : END DO
484 2 : CALL pw_grid%para%group%sum(v_hartree_z_average)
485 : v_hartree_z_average(:) = v_hartree_z_average(:)/ &
486 178 : (REAL(pw_grid%npts(1), kind=dp)*REAL(pw_grid%npts(2), kind=dp))
487 :
488 2 : funit = -1
489 2 : IF (pw_grid%para%group%mepos == pw_grid%para%group%source) THEN
490 : CALL open_file(TRIM(ADJUSTL(project_name))//"-VH_AV.dat", &
491 1 : file_status="REPLACE", file_form="FORMATTED", file_action="WRITE", unit_number=funit)
492 1 : WRITE (funit, '(A,T10,A,T25,A)') "#", "z (A)", "V_H average (eV)"
493 89 : DO iz = lz, uz
494 88 : WRITE (funit, '(F20.10,ES20.10E3)') pw_grid%dh(3, 3)*REAL(iz - lz, kind=dp)/bohr, &
495 177 : v_hartree_z_average(iz)/pw_grid%dvol*evolt
496 : END DO
497 1 : CALL close_file(funit)
498 : END IF
499 :
500 2 : DEALLOCATE (v_hartree_z_average)
501 2 : CALL timestop(handle)
502 2 : END SUBROUTINE write_average_hartree_potential
503 :
504 : ! **************************************************************************************************
505 : !> \brief Align Hatree potential of semi-infinite leads to match bulk-transport calculation
506 : !> and apply external electrostatic potential (bias).
507 : !> \param v_hartree_rspace Hartree potential on a real-space 3-D grid [inout]
508 : !> \param cell unit cell
509 : !> \param HartreeLeadsLeft z-coordinate of the left lead
510 : !> \param HartreeLeadsRight z-coordinate of the right lead
511 : !> \param HartreeLeadsBottom average Hartree potential (from bulk-transport calculation)
512 : !> at HartreeLeadsLeft and HartreeLeadsRight
513 : !> \param Vbias the value of external potential to apply
514 : !> \param zleft starting point of external potential drop (initial value 0.5*Vbias)
515 : !> \param zright final point of external potential drop (final value -0.5*Vbias)
516 : !> \param isexplicit_zright whether zright has beed provided explicitly via the input file.
517 : !> If not, use the cell boundary.
518 : !> \param isexplicit_bottom whether the reference Hatree potential for bulk regions has been
519 : !> provided explicitly via the input file. If not, do not align the
520 : !> potential at all (instead of aligning it to 0 which is incorrect).
521 : !> \note This routine assumes that the lattice vector 'k' coincides with z-axis
522 : ! **************************************************************************************************
523 0 : SUBROUTINE smeagol_shift_v_hartree(v_hartree_rspace, cell, &
524 : HartreeLeadsLeft, HartreeLeadsRight, HartreeLeadsBottom, &
525 : Vbias, zleft, zright, isexplicit_zright, isexplicit_bottom)
526 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
527 : TYPE(cell_type), POINTER :: cell
528 : REAL(kind=dp), INTENT(in) :: HartreeLeadsLeft, HartreeLeadsRight, &
529 : HartreeLeadsBottom, Vbias, zleft, &
530 : zright
531 : LOGICAL, INTENT(in) :: isexplicit_zright, isexplicit_bottom
532 :
533 : CHARACTER(LEN=*), PARAMETER :: routineN = 'smeagol_shift_v_hartree'
534 :
535 : INTEGER :: handle, iz, l_right, lz, u_left, uz
536 : REAL(kind=dp) :: v_average_left, v_average_right, &
537 : v_bias_iz, zright_explicit
538 0 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: v_hartree_z_average
539 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
540 0 : POINTER :: cr3d
541 : REAL(kind=dp), DIMENSION(3) :: r, r_pbc
542 : TYPE(pw_grid_type), POINTER :: pw_grid
543 :
544 0 : CALL timeset(routineN, handle)
545 0 : pw_grid => v_hartree_rspace%pw_grid
546 0 : cr3d => v_hartree_rspace%array
547 :
548 0 : zright_explicit = zright
549 0 : IF (.NOT. isexplicit_zright) THEN
550 0 : r_pbc(:) = [0.0_dp, 0.0_dp, 1.0_dp]
551 0 : CALL scaled_to_real(r, r_pbc, cell)
552 0 : zright_explicit = r(3)
553 : END IF
554 :
555 0 : ALLOCATE (v_hartree_z_average(pw_grid%bounds(1, 3):pw_grid%bounds(2, 3)))
556 :
557 0 : lz = pw_grid%bounds_local(1, 3)
558 0 : uz = pw_grid%bounds_local(2, 3)
559 :
560 0 : v_hartree_z_average(:) = 0.0_dp
561 0 : DO iz = lz, uz
562 0 : v_hartree_z_average(iz) = accurate_sum(cr3d(:, :, iz))
563 : END DO
564 :
565 0 : CALL pw_grid%para%group%sum(v_hartree_z_average)
566 : v_hartree_z_average(:) = v_hartree_z_average(:)/ &
567 0 : (REAL(pw_grid%npts(1), kind=dp)*REAL(pw_grid%npts(2), kind=dp))
568 :
569 : ! z-indices of the V_hartree related to the left lead: pw_grid%bounds(1,3) .. u_left
570 0 : r(1:3) = [0.0_dp, 0.0_dp, HartreeLeadsLeft]
571 0 : CALL real_to_scaled(r_pbc, r, cell)
572 0 : u_left = NINT(r_pbc(3)*REAL(pw_grid%npts(3), kind=dp)) + pw_grid%bounds(1, 3)
573 0 : IF (u_left > pw_grid%bounds(2, 3)) u_left = pw_grid%bounds(2, 3)
574 :
575 : ! z-indices of the V_hartree related to the right lead: l_right .. pw_grid%bounds(2, 3)
576 0 : r(1:3) = [0.0_dp, 0.0_dp, HartreeLeadsRight]
577 0 : CALL real_to_scaled(r_pbc, r, cell)
578 0 : l_right = NINT(r_pbc(3)*REAL(pw_grid%npts(3), kind=dp)) + pw_grid%bounds(1, 3)
579 0 : IF (l_right > pw_grid%bounds(2, 3)) l_right = pw_grid%bounds(2, 3)
580 :
581 0 : CPASSERT(u_left <= l_right)
582 :
583 0 : v_average_left = v_hartree_z_average(u_left)
584 0 : v_average_right = v_hartree_z_average(l_right)
585 :
586 : ! align electrostatic potential of leads' regions with ones from bulk transport calculation
587 0 : IF (isexplicit_bottom) THEN
588 0 : v_hartree_z_average(:) = HartreeLeadsBottom*pw_grid%dvol - 0.5_dp*(v_average_left + v_average_right)
589 : ELSE
590 : ! do not align electrostatic potential
591 0 : v_hartree_z_average(:) = 0.0_dp
592 : END IF
593 :
594 : ! external Vbias
595 : ! TO DO: convert zright and zleft to scaled coordinates instead
596 0 : DO iz = lz, uz
597 0 : r_pbc(1) = 0.0_dp
598 0 : r_pbc(2) = 0.0_dp
599 0 : r_pbc(3) = REAL(iz - pw_grid%bounds(1, 3), kind=dp)/REAL(pw_grid%npts(3), kind=dp)
600 0 : CALL scaled_to_real(r, r_pbc, cell)
601 0 : IF (r(3) < zleft) THEN
602 0 : v_bias_iz = 0.5_dp*Vbias
603 0 : ELSE IF (r(3) > zright_explicit) THEN
604 0 : v_bias_iz = -0.5_dp*Vbias
605 : ELSE
606 0 : v_bias_iz = Vbias*(0.5_dp - (r(3) - zleft)/(zright_explicit - zleft))
607 : END IF
608 0 : v_hartree_z_average(iz) = v_hartree_z_average(iz) + v_bias_iz*pw_grid%dvol
609 : END DO
610 :
611 0 : DO iz = lz, uz
612 0 : cr3d(:, :, iz) = cr3d(:, :, iz) + v_hartree_z_average(iz)
613 : END DO
614 :
615 0 : DEALLOCATE (v_hartree_z_average)
616 0 : CALL timestop(handle)
617 0 : END SUBROUTINE smeagol_shift_v_hartree
618 :
619 : ! **************************************************************************************************
620 : !> \brief Run NEGF/SMEAGOL transport calculation.
621 : !> \param qs_env QuickStep environment
622 : !> \param last converged SCF iterations; compute NEGF properties [in]
623 : !> \param iter index of the current iteration [in]
624 : !> \param rho_ao_kp refined electron density; to be mixed with electron density from the previous
625 : !> SCF iteration [out]
626 : ! **************************************************************************************************
627 17637 : SUBROUTINE run_smeagol_emtrans(qs_env, last, iter, rho_ao_kp)
628 : TYPE(qs_environment_type), POINTER :: qs_env
629 : LOGICAL, INTENT(in) :: last
630 : INTEGER, INTENT(in) :: iter
631 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
632 : POINTER :: rho_ao_kp
633 :
634 : CHARACTER(LEN=*), PARAMETER :: routineN = 'run_smeagol_emtrans'
635 :
636 : INTEGER :: handle, img, ispin, md_step, natoms, &
637 : nimages, nspin
638 : INTEGER, DIMENSION(2) :: max_ij_cell_image, n_cell_images
639 17637 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
640 : LOGICAL :: do_kpoints, local_ldos, local_TrCoeff, &
641 : negfon_saved, structure_changed
642 : REAL(kind=dp) :: H_to_Ry
643 17637 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: matrix_s_csc_merged
644 17637 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: Dnew_csc_merged, Enew_csc_merged, &
645 17637 : matrix_ks_csc_merged
646 : TYPE(cell_type), POINTER :: ucell
647 : TYPE(cp_logger_type), POINTER :: logger
648 17637 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: matrix_kp_generic
649 17637 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, matrix_w_kp
650 : TYPE(dft_control_type), POINTER :: dft_control
651 : TYPE(kpoint_type), POINTER :: kpoints
652 : TYPE(mp_para_env_type), POINTER :: para_env
653 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
654 17637 : POINTER :: sab_nl
655 : TYPE(qs_subsys_type), POINTER :: subsys
656 : TYPE(scf_control_type), POINTER :: scf_control
657 17637 : TYPE(siesta_distrib_csc_struct_type) :: siesta_struct_merged
658 : TYPE(smeagol_control_type), POINTER :: smeagol_control
659 :
660 35274 : logger => cp_get_default_logger()
661 :
662 17637 : CALL get_qs_env(qs_env, dft_control=dft_control)
663 :
664 17637 : nspin = dft_control%nspins
665 17637 : smeagol_control => dft_control%smeagol_control
666 17637 : H_to_Ry = smeagol_control%to_smeagol_energy_units
667 :
668 17637 : IF (.NOT. (smeagol_control%smeagol_enabled .AND. smeagol_control%run_type == smeagol_runtype_emtransport)) RETURN
669 0 : CALL timeset(routineN, handle)
670 :
671 0 : NULLIFY (kpoints, matrix_s_kp, matrix_ks_kp, matrix_w_kp, para_env, sab_nl, scf_control, subsys)
672 : #if defined(__SMEAGOL)
673 0 : CALL cite_reference(Ahart2024)
674 0 : CALL cite_reference(Bailey2006)
675 :
676 0 : CPASSERT(ASSOCIATED(smeagol_control%aux))
677 : CALL get_qs_env(qs_env, para_env=para_env, scf_control=scf_control, &
678 : do_kpoints=do_kpoints, kpoints=kpoints, &
679 : matrix_ks_kp=matrix_ks_kp, matrix_s_kp=matrix_s_kp, matrix_w_kp=matrix_w_kp, &
680 0 : sab_orb=sab_nl, subsys=subsys)
681 0 : CALL qs_subsys_get(subsys, cell=ucell, natom=natoms)
682 :
683 0 : IF (do_kpoints) THEN
684 0 : nimages = dft_control%nimages
685 0 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
686 : ELSE
687 0 : nimages = 1
688 0 : ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
689 0 : cell_to_index(0, 0, 0) = 1
690 : END IF
691 :
692 0 : max_ij_cell_image(:) = -1
693 0 : DO img = 1, 2
694 0 : IF (smeagol_control%n_cell_images(img) > 0) THEN
695 0 : max_ij_cell_image(img) = smeagol_control%n_cell_images(img)/2
696 : END IF
697 : END DO
698 :
699 : ! obtain the sparsity pattern of the Kohn-Sham matrix
700 0 : ALLOCATE (matrix_kp_generic(nimages))
701 0 : DO img = 1, nimages
702 0 : matrix_kp_generic(img)%matrix => matrix_s_kp(1, img)%matrix
703 : END DO
704 :
705 : CALL siesta_struct_create(siesta_struct_merged, matrix_kp_generic, subsys, &
706 0 : cell_to_index, sab_nl, para_env, max_ij_cell_image, do_merge=.TRUE., gather_root=-1)
707 :
708 : ! Number of unit cells along (x and y ?) directions
709 0 : n_cell_images(1:2) = 2*max_ij_cell_image(1:2) + 1
710 :
711 0 : ALLOCATE (matrix_s_csc_merged(siesta_struct_merged%n_nonzero_elements))
712 0 : ALLOCATE (matrix_ks_csc_merged(siesta_struct_merged%n_nonzero_elements, nspin))
713 :
714 0 : CALL convert_dbcsr_to_distributed_siesta(matrix_s_csc_merged, matrix_kp_generic, siesta_struct_merged, para_env)
715 :
716 0 : DO ispin = 1, nspin
717 0 : DO img = 1, nimages
718 0 : matrix_kp_generic(img)%matrix => matrix_ks_kp(ispin, img)%matrix
719 : END DO
720 :
721 : CALL convert_dbcsr_to_distributed_siesta(matrix_ks_csc_merged(:, ispin), matrix_kp_generic, &
722 0 : siesta_struct_merged, para_env)
723 : END DO
724 :
725 0 : IF (smeagol_control%aux%md_iter_level > 0) THEN
726 0 : CALL cp_get_iter_nr(logger%iter_info, smeagol_control%aux%md_iter_level, iter_nr=md_step)
727 0 : ELSE IF (smeagol_control%aux%md_iter_level == 0) THEN
728 : ! single-point energy calculation : there is only one MD step in this case
729 0 : md_step = smeagol_control%aux%md_first_step
730 : ELSE
731 : ! first invocation of the subroutine : initialise md_iter_level and md_first_step variables
732 0 : smeagol_control%aux%md_iter_level = cp_get_iter_level_by_name(logger%iter_info, "MD")
733 :
734 0 : IF (smeagol_control%aux%md_iter_level <= 0) &
735 0 : smeagol_control%aux%md_iter_level = cp_get_iter_level_by_name(logger%iter_info, "GEO_OPT")
736 :
737 0 : IF (smeagol_control%aux%md_iter_level <= 0) &
738 0 : smeagol_control%aux%md_iter_level = 0
739 :
740 : ! index of the first GEO_OPT / MD step
741 0 : IF (smeagol_control%aux%md_iter_level > 0) THEN
742 0 : CALL cp_get_iter_nr(logger%iter_info, smeagol_control%aux%md_iter_level, iter_nr=smeagol_control%aux%md_first_step)
743 : ELSE
744 : ! it has already been initialised in read_smeagol_control()
745 0 : smeagol_control%aux%md_first_step = 0
746 : END IF
747 :
748 0 : md_step = smeagol_control%aux%md_first_step
749 : END IF
750 :
751 0 : CALL reademtr(smeagol_control, natoms, gamma_negf=.FALSE.)
752 0 : CALL ReadOptionsNEGF_DFT(smeagol_control, ucell, torqueflag=.FALSE., torquelin=.FALSE.)
753 :
754 : CALL emtrans_options(smeagol_control, &
755 : matrix_s=matrix_s_kp(1, 1)%matrix, para_env=para_env, iter=iter, &
756 : istep=md_step, inicoor=smeagol_control%aux%md_first_step, iv=0, &
757 0 : delta=smeagol_control%aux%delta, nk=kpoints%nkp)
758 :
759 : CALL create_communicators_negf(parent_comm=para_env%get_handle(), &
760 : nprocs_inverse=smeagol_control%aux%nprocs_inverse, &
761 0 : NParallelK=smeagol_control%aux%NParallelK)
762 :
763 0 : ALLOCATE (Dnew_csc_merged(siesta_struct_merged%n_nonzero_elements, nspin))
764 0 : ALLOCATE (Enew_csc_merged(siesta_struct_merged%n_nonzero_elements, nspin))
765 :
766 : ! As SMEAGOL's default energy unit is Rydberg, scale the KS-matrix
767 0 : matrix_ks_csc_merged(:, :) = H_to_Ry*matrix_ks_csc_merged(:, :)
768 :
769 : ! SMEAGOL computes current if EM.OrderN = .false., otherwise the printed current is equal to 0.
770 : ! The following computes current regardless the value of EM.OrderN parameter
771 0 : IF (last) THEN
772 0 : negfon_saved = smeagolglobal_negfon
773 0 : smeagolglobal_negfon = .FALSE.
774 : END IF
775 :
776 0 : IF (last) THEN
777 0 : local_TrCoeff = smeagol_control%aux%TrCoeff
778 0 : local_ldos = .TRUE. ! TO DO: find out initial value of ldos
779 : ELSE
780 0 : local_TrCoeff = .FALSE.
781 0 : local_ldos = .FALSE.
782 : END IF
783 :
784 : ! number of atoms in the unit cell
785 0 : smeagolglobal_em_nau = natoms ! number of atoms in the unit cell
786 :
787 : ! number of atomic orbitals in the unit cell.
788 : ! This global variable defines the size of temporary arrays for (P)DOS calculation.
789 : ! This should be the total number of the atomic orbitals in the unit cell,
790 : ! instead of the number of atomic orbitals local to the current MPI process.
791 0 : smeagolglobal_em_nuo = siesta_struct_merged%ncols/(n_cell_images(1)*n_cell_images(2))
792 :
793 : ! number of atoms in the super cell
794 0 : smeagolglobal_em_nas = SIZE(siesta_struct_merged%xa)
795 :
796 : ! number of atomic orbitals in the super cell
797 0 : smeagolglobal_em_nso = siesta_struct_merged%ncols
798 :
799 : ! The following global variables are only used in writephibin() which is not called by this interface
800 : ! smeagolglobal_em_isa => isa
801 : ! smeagolglobal_em_iaorb => iaorb
802 : ! smeagolglobal_em_iphorb => iphorb
803 :
804 0 : structure_changed = (iter == 1)
805 :
806 : CALL Negf_Interface( &
807 : ! distributed Kohn-Sham, overlap, density, and energy-density matrices in SIESTA format
808 : H=matrix_ks_csc_merged, &
809 : S=matrix_s_csc_merged, &
810 : DM=Dnew_csc_merged, &
811 : Omega=Enew_csc_merged, &
812 : ! interatomic distances for each non-zero matrix element
813 : xij=siesta_struct_merged%xij, &
814 : ! number of atomic orbitals in a supercell
815 : no_s=siesta_struct_merged%ncols, &
816 : ! number of atomic orbitals in the unit cell
817 : no_u=siesta_struct_merged%ncols/(n_cell_images(1)*n_cell_images(2)), &
818 : ! number of AOs local to the given MPI process
819 : no_u_node=siesta_struct_merged%nrows, &
820 : ! atomic coordinates for each atom in the supercell
821 : xa=siesta_struct_merged%xa, &
822 : ! unused dummy variable na_u
823 : na_u=natoms, &
824 : ! number of atoms in the supercell
825 : na_s=SIZE(siesta_struct_merged%xa), &
826 : ! number of spin-components
827 : NspinRealInputMatrix=nspin, &
828 : ! number of non-zero matrix element
829 : maxnh=INT(siesta_struct_merged%n_nonzero_elements), &
830 : ! number of non-zero matrix elements on each locally stored row
831 : numh=siesta_struct_merged%n_nonzero_cols, &
832 : ! offset (index-1) of first non-zero matrix elements on each locally stored row
833 : listhptr=siesta_struct_merged%row_offset, &
834 : ! column number (AO index) for each non-zero matrix element
835 : listh=siesta_struct_merged%col_index, &
836 : ! number of k-points
837 : nkpts=kpoints%nkp, &
838 : ! k-point coordinates
839 : kpoint=kpoints%xkp, &
840 : ! k-point weight
841 : weight_k=kpoints%wkp, &
842 : ! index of equivalent atomic orbital within the primary unit cell for each AO in the supercell
843 : indxuo=siesta_struct_merged%indxuo, &
844 : ! list of atomic indices on which each AO (in the supercell) is centred
845 : iaorb=siesta_struct_merged%iaorb, &
846 : ! GEO_OPT / MD step index
847 : MD_Step=md_step, &
848 : ! GEO_OPT / MD at which SMEAGOL should allocate its internal data structures
849 : inicoor=smeagol_control%aux%md_first_step, &
850 : ! ivv (step over bias, first IV step should always be 0 (hardcoded in SMEAGOL)
851 : ! we iterate over GEO_OPT / MD steps instead of bias steps in order not to overwrite
852 : ! TRC (Transmission coefficients) files
853 : IV_step=md_step - smeagol_control%aux%md_first_step, &
854 : ! applied electrostatic bias
855 : Vb=smeagol_control%aux%VBias*H_to_Ry, &
856 : ! index of the currect SCF iteration
857 : SCF_step=iter, &
858 : ! compute density matrix (.FALSE.) or properties (.TRUE.)
859 : Last_SCF_step=last, &
860 : ! recompute self-energy matrices
861 : StructureChanged=structure_changed, &
862 : ! electronic temperature in Ry
863 : temp=smeagol_control%aux%temperature*H_to_Ry, &
864 : ! number of unit cell replicas along i-, and j- cell verctors
865 : nsc=n_cell_images, &
866 : ! name of SMEAGOL project (partial file name for files created by SMEAGOL)
867 : slabel=smeagol_control%project_name, &
868 : ! number of integration points along real axis, circular path and a line in imaginary space parallel to the real axis
869 : NEnergR=smeagol_control%aux%NEnergR, &
870 : NEnergIC=smeagol_control%aux%NEnergIC, &
871 : NEnergIL=smeagol_control%aux%NEnergIL, &
872 : ! number of poles
873 : NPoles=smeagol_control%aux%NPoles, &
874 : ! small imaginary shift that makes matrix inversion computationally stable
875 : Delta=smeagol_control%aux%deltamin*H_to_Ry, &
876 : ! integration lower bound
877 : EnergLB=smeagol_control%aux%EnergLB*H_to_Ry, &
878 : ! [unused dummy arguments] initial (VInitial) and final (VFinal) voltage
879 : VInitial=smeagol_control%aux%VBias, &
880 : VFinal=smeagol_control%aux%VBias, &
881 : !
882 : SpinCL=smeagol_control%aux%SpinCL, &
883 : ! number of slices for OrderN matrix inversion
884 : NSlices=smeagol_control%aux%NSlices, &
885 : ! whether to compute transmission coefficients (TrCoeff), IETS spectrum (CalcIETS), and local Dnsity-of-States (ldos)
886 : TrCoeff=local_TrCoeff, &
887 : CalcIETS=smeagol_control%aux%CalcIETS, &
888 : ldos=local_ldos, &
889 : ! Transmission coefficients are only computed for certain runtypes (encoded with magic numbers).
890 : ! In case of idyn=0, transmission coefficients are always computed.
891 : idyn=0, &
892 : ! do not compute transmission coefficient for initial GEO_OPT / MD iterations
893 : tmdskip=smeagol_control%aux%tmdskip, &
894 : ! computes transmission coefficients once for each 'tmdsampling' MD iterations
895 0 : tmdsampling=smeagol_control%aux%tmdsampling)
896 :
897 : ! *** Bound-state correction method ***
898 :
899 : ! smeagol_control%bs_add : BS.Add (.FALSE.) ; add bound states
900 : ! smeagol_control%bs_method : BS.Method (0) ; 0 - use effective Hamiltonian; 1 - add small imaginary part to the selfenergies
901 : ! smeagol_control%bssc : BS.SetOccupation (1)
902 : IF (smeagol_control%aux%bs_add .AND. smeagol_control%aux%bs_method == 1 .AND. smeagol_control%aux%bssc == 0 .AND. &
903 0 : kpoints%nkp > 1 .AND. (.NOT. last)) THEN
904 :
905 : CALL Negf_Interface(H=matrix_ks_csc_merged, &
906 : S=matrix_s_csc_merged, &
907 : DM=Dnew_csc_merged, &
908 : Omega=Enew_csc_merged, &
909 : xij=siesta_struct_merged%xij, &
910 : no_s=siesta_struct_merged%ncols, &
911 : no_u=siesta_struct_merged%ncols/(n_cell_images(1)*n_cell_images(2)), &
912 : no_u_node=siesta_struct_merged%nrows, &
913 : xa=siesta_struct_merged%xa, &
914 : ! unused dummy variable
915 : na_u=natoms, &
916 : na_s=SIZE(siesta_struct_merged%xa), &
917 : NspinRealInputMatrix=nspin, &
918 : maxnh=INT(siesta_struct_merged%n_nonzero_elements), &
919 : numh=siesta_struct_merged%n_nonzero_cols, &
920 : listhptr=siesta_struct_merged%row_offset, &
921 : listh=siesta_struct_merged%col_index, &
922 : nkpts=kpoints%nkp, &
923 : kpoint=kpoints%xkp, &
924 : weight_k=kpoints%wkp, &
925 : indxuo=siesta_struct_merged%indxuo, &
926 : iaorb=siesta_struct_merged%iaorb, &
927 : MD_Step=md_step, &
928 : inicoor=smeagol_control%aux%md_first_step, &
929 : IV_step=md_step - smeagol_control%aux%md_first_step, &
930 : Vb=smeagol_control%aux%VBias*H_to_Ry, &
931 : ! This line is the only difference from the first Negf_Interface() call
932 : SCF_step=MERGE(2, 1, structure_changed), &
933 : Last_SCF_step=last, &
934 : StructureChanged=structure_changed, &
935 : temp=smeagol_control%aux%temperature*H_to_Ry, &
936 : nsc=n_cell_images, &
937 : slabel=smeagol_control%project_name, &
938 : NEnergR=smeagol_control%aux%NEnergR, &
939 : NEnergIC=smeagol_control%aux%NEnergIC, &
940 : NEnergIL=smeagol_control%aux%NEnergIL, &
941 : NPoles=smeagol_control%aux%NPoles, &
942 : Delta=smeagol_control%aux%deltamin*H_to_Ry, &
943 : EnergLB=smeagol_control%aux%EnergLB*H_to_Ry, &
944 : ! unused dummy variable
945 : VInitial=smeagol_control%aux%VBias, &
946 : ! unused dummy variable
947 : VFinal=smeagol_control%aux%VBias, &
948 : SpinCL=smeagol_control%aux%SpinCL, &
949 : NSlices=smeagol_control%aux%NSlices, &
950 : TrCoeff=local_TrCoeff, &
951 : CalcIETS=smeagol_control%aux%CalcIETS, &
952 : ldos=local_ldos, &
953 : idyn=0, & !
954 : tmdskip=smeagol_control%aux%tmdskip, &
955 0 : tmdsampling=smeagol_control%aux%tmdsampling)
956 : END IF
957 :
958 : ! Restore ovewriten EM.OrderN parameter
959 0 : IF (last) THEN
960 0 : smeagolglobal_negfon = negfon_saved
961 : END IF
962 :
963 0 : IF (PRESENT(rho_ao_kp)) THEN
964 0 : DO ispin = 1, nspin
965 0 : DO img = 1, nimages
966 : ! To be on a safe size, zeroize the new density matrix first. It is not actually needed.
967 0 : CALL dbcsr_set(rho_ao_kp(ispin, img)%matrix, 0.0_dp)
968 0 : matrix_kp_generic(img)%matrix => rho_ao_kp(ispin, img)%matrix
969 : END DO
970 :
971 : CALL convert_distributed_siesta_to_dbcsr(matrix_kp_generic, Dnew_csc_merged(:, ispin), &
972 0 : siesta_struct_merged, para_env)
973 : END DO
974 :
975 : ! current-induced forces
976 0 : IF (smeagol_control%emforces) THEN
977 0 : Enew_csc_merged(:, :) = (1.0_dp/H_to_Ry)*Enew_csc_merged(:, :)
978 :
979 0 : DO ispin = 1, nspin
980 0 : DO img = 1, nimages
981 : ! To be on a safe size, zeroize the new energy-density matrix first. It is not actually needed.
982 0 : CALL dbcsr_set(matrix_w_kp(ispin, img)%matrix, 0.0_dp)
983 0 : matrix_kp_generic(img)%matrix => matrix_w_kp(ispin, img)%matrix
984 : END DO
985 :
986 : CALL convert_distributed_siesta_to_dbcsr(matrix_kp_generic, Enew_csc_merged(:, ispin), &
987 0 : siesta_struct_merged, para_env)
988 : END DO
989 : END IF
990 : END IF
991 :
992 0 : CALL destroy_communicators_negf()
993 0 : CALL emtrans_deallocate_global_arrays()
994 :
995 0 : DEALLOCATE (Dnew_csc_merged, Enew_csc_merged)
996 0 : DEALLOCATE (matrix_s_csc_merged, matrix_ks_csc_merged)
997 :
998 0 : CALL siesta_struct_release(siesta_struct_merged)
999 :
1000 0 : IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
1001 : #else
1002 : CALL cp_abort(__LOCATION__, &
1003 : "CP2K was compiled with no SMEAGOL support.")
1004 : MARK_USED(last)
1005 : MARK_USED(iter)
1006 : MARK_USED(rho_ao_kp)
1007 : ! local variables
1008 : MARK_USED(cell_to_index)
1009 : MARK_USED(do_kpoints)
1010 : MARK_USED(Dnew_csc_merged)
1011 : MARK_USED(Enew_csc_merged)
1012 : MARK_USED(img)
1013 : MARK_USED(ispin)
1014 : MARK_USED(local_ldos)
1015 : MARK_USED(local_TrCoeff)
1016 : MARK_USED(max_ij_cell_image)
1017 : MARK_USED(matrix_kp_generic)
1018 : MARK_USED(matrix_ks_csc_merged)
1019 : MARK_USED(matrix_s_csc_merged)
1020 : MARK_USED(md_step)
1021 : MARK_USED(natoms)
1022 : MARK_USED(negfon_saved)
1023 : MARK_USED(nimages)
1024 : MARK_USED(n_cell_images)
1025 : MARK_USED(siesta_struct_merged)
1026 : MARK_USED(structure_changed)
1027 : MARK_USED(ucell)
1028 : #endif
1029 :
1030 0 : CALL timestop(handle)
1031 35274 : END SUBROUTINE run_smeagol_emtrans
1032 :
1033 : END MODULE smeagol_interface
|
