(git:495eafe)
Loading...
Searching...
No Matches
negf_env_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Environment for NEGF based quantum transport calculations
10! **************************************************************************************************
11MODULE negf_env_types
12 USE cell_types, ONLY: cell_type,&
13 real_to_scaled
14 USE cp_blacs_env, ONLY: cp_blacs_env_type
15 USE cp_control_types, ONLY: dft_control_type
16 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
17 dbcsr_deallocate_matrix,&
18 dbcsr_init_p,&
19 dbcsr_p_type,&
20 dbcsr_set
21 USE cp_files, ONLY: close_file,&
22 open_file
23 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 cp_fm_struct_release,&
25 cp_fm_struct_type
26 USE cp_fm_types, ONLY: cp_fm_create,&
27 cp_fm_get_info,&
28 cp_fm_get_submatrix,&
29 cp_fm_release,&
30 cp_fm_set_submatrix,&
31 cp_fm_type
32 USE cp_log_handling, ONLY: cp_get_default_logger,&
33 cp_logger_type
34 USE force_env_types, ONLY: force_env_get,&
35 force_env_p_type,&
36 force_env_type,&
37 use_qs_force
38 USE input_section_types, ONLY: section_vals_get_subs_vals,&
39 section_vals_type,&
40 section_vals_val_get
41 USE kinds, ONLY: default_path_length,&
42 default_string_length,&
43 dp
44 USE kpoint_types, ONLY: get_kpoint_env,&
45 get_kpoint_info,&
46 kpoint_env_p_type,&
47 kpoint_type
48 USE message_passing, ONLY: mp_para_env_type
49 USE negf_atom_map, ONLY: negf_atom_map_type,&
50 negf_map_atomic_indices
51 USE negf_control_types, ONLY: negf_control_contact_type,&
52 negf_control_type
53 USE negf_io, ONLY: negf_print_matrix_to_file,&
54 negf_read_matrix_from_file,&
55 negf_restart_file_name
56 USE negf_matrix_utils, ONLY: invert_cell_to_index,&
57 negf_copy_contact_matrix,&
58 negf_copy_sym_dbcsr_to_fm_submat,&
59 number_of_atomic_orbitals
60 USE negf_subgroup_types, ONLY: negf_subgroup_env_type
61 USE negf_vectors, ONLY: contact_direction_vector,&
62 projection_on_direction_vector
63 USE particle_types, ONLY: particle_type
64 USE pw_env_types, ONLY: pw_env_get,&
65 pw_env_type
66 USE pw_pool_types, ONLY: pw_pool_type
67 USE pw_types, ONLY: pw_r3d_rs_type
68 USE qs_density_mixing_types, ONLY: mixing_storage_create,&
69 mixing_storage_release,&
70 mixing_storage_type
71 USE qs_energy, ONLY: qs_energies
72 USE qs_energy_init, ONLY: qs_energies_init
73 USE qs_environment_types, ONLY: get_qs_env,&
74 qs_environment_type
75 USE qs_integrate_potential, ONLY: integrate_v_rspace
76 USE qs_mo_types, ONLY: get_mo_set,&
77 mo_set_type
78 USE qs_rho_types, ONLY: qs_rho_get,&
79 qs_rho_type
80 USE qs_scf_post_tb, ONLY: rebuild_pw_env
81 USE qs_subsys_types, ONLY: qs_subsys_get,&
82 qs_subsys_type
83#include "./base/base_uses.f90"
84
85 IMPLICIT NONE
86 PRIVATE
87
88 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_env_types'
89 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .true.
90
91 PUBLIC :: negf_env_type, negf_env_contact_type
92 PUBLIC :: negf_env_create, negf_env_release
93
94! **************************************************************************************************
95!> \brief Contact-specific NEGF environment.
96!> \author Sergey Chulkov
97! **************************************************************************************************
98 TYPE negf_env_contact_type
99 REAL(kind=dp), DIMENSION(3) :: direction_vector = -1.0_dp, origin = -1.0_dp
100 REAL(kind=dp), DIMENSION(3) :: direction_vector_bias = -1.0_dp, origin_bias = -1.0_dp
101 !> an axis towards the secondary contact unit cell which coincides with the transport direction
102 !> 0 (undefined), 1 (+x), 2 (+y), 3 (+z), -1 (-x), -2 (-y), -3 (-z)
103 INTEGER :: direction_axis = -1
104 !> atoms belonging to a primary contact unit cell
105 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_cell0
106 !> atoms belonging to a secondary contact unit cell (will be removed one day ...)
107 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_cell1
108 !> list of equivalent atoms in an appropriate contact force environment
109 TYPE(negf_atom_map_type), ALLOCATABLE, &
110 DIMENSION(:) :: atom_map_cell0, atom_map_cell1
111 !> Fermi energy
112 REAL(kind=dp) :: fermi_energy = 0.0_dp
113 !> energy of the HOMO
114 REAL(kind=dp) :: homo_energy = -1.0_dp
115 !> number of electrons Sp(rho_00,s_00)
116 REAL(kind=dp) :: nelectrons_qs_cell0 = 0.0_dp
117 !> number of electrons Sp(rho_01,s_01)
118 REAL(kind=dp) :: nelectrons_qs_cell1 = 0.0_dp
119 !> diagonal (h_00) and off-diagonal (h_01) blocks of the contact Kohn-Sham matrix ([number_of_spins]).
120 !> The matrix h_01 is of the shape [nao_cell0 x nao_cell1]
121 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h_00, h_01
122 !> diagonal and off-diagonal blocks of the density matrix
123 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: rho_00, rho_01
124 !> diagonal and off-diagonal blocks of the overlap matrix
125 TYPE(cp_fm_type), POINTER :: s_00 => null(), s_01 => null()
126 END TYPE negf_env_contact_type
127
128! **************************************************************************************************
129!> \brief NEGF environment.
130!> \author Sergey Chulkov
131! **************************************************************************************************
132 TYPE negf_env_type
133 !> contact-specific NEGF environments
134 TYPE(negf_env_contact_type), ALLOCATABLE, &
135 DIMENSION(:) :: contacts
136 !> Kohn-Sham matrix of the scattering region
137 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h_s
138 !> Kohn-Sham matrix of the scattering region -- contact interface ([nspins, ncontacts])
139 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: h_sc
140 !> overlap matrix of the scattering region
141 TYPE(cp_fm_type), POINTER :: s_s => null()
142 !> an external Hartree potential in atomic basis set representation
143 TYPE(cp_fm_type), POINTER :: v_hartree_s => null()
144 !> overlap matrix of the scattering region -- contact interface for every contact ([ncontacts])
145 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_sc
146 !> structure needed for density mixing
147 TYPE(mixing_storage_type), POINTER :: mixing_storage => null()
148 !> density mixing method
149 INTEGER :: mixing_method = -1
150 !> number of electrons Sp(rho_s,s_s)
151 REAL(kind=dp) :: nelectrons_ref = 0.0_dp
152 !> number of electrons Sp(rho_s,s_s)
153 REAL(kind=dp) :: nelectrons = 0.0_dp
154 END TYPE negf_env_type
155
156! **************************************************************************************************
157!> \brief Allocatable list of the type 'negf_atom_map_type'.
158!> \author Sergey Chulkov
159! **************************************************************************************************
160 TYPE negf_atom_map_contact_type
161 TYPE(negf_atom_map_type), ALLOCATABLE, DIMENSION(:) :: atom_map
162 END TYPE negf_atom_map_contact_type
163
164CONTAINS
165
166! **************************************************************************************************
167!> \brief Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
168!> \param negf_env NEGF environment to create
169!> \param sub_env NEGF parallel (sub)group environment
170!> \param negf_control NEGF control
171!> \param force_env the primary force environment
172!> \param negf_mixing_section pointer to a mixing section within the NEGF input section
173!> \param log_unit output unit number
174!> \par History
175!> * 01.2017 created [Sergey Chulkov]
176! **************************************************************************************************
177 SUBROUTINE negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
178 TYPE(negf_env_type), INTENT(inout) :: negf_env
179 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
180 TYPE(negf_control_type), POINTER :: negf_control
181 TYPE(force_env_type), POINTER :: force_env
182 TYPE(section_vals_type), POINTER :: negf_mixing_section
183 INTEGER, INTENT(in) :: log_unit
184
185 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_create'
186
187 CHARACTER(len=default_string_length) :: contact_str, force_env_str, &
188 n_force_env_str
189 INTEGER :: handle, icontact, in_use, n_force_env, &
190 ncontacts
191 LOGICAL :: do_kpoints, is_dft_entire
192 TYPE(cp_blacs_env_type), POINTER :: blacs_env
193 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
194 TYPE(dft_control_type), POINTER :: dft_control
195 TYPE(force_env_p_type), DIMENSION(:), POINTER :: sub_force_env
196 TYPE(mp_para_env_type), POINTER :: para_env
197 TYPE(negf_atom_map_contact_type), ALLOCATABLE, &
198 DIMENSION(:) :: map_contact
199 TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
200 TYPE(qs_environment_type), POINTER :: qs_env, qs_env_contact
201 TYPE(qs_subsys_type), POINTER :: subsys, subsys_contact
202 TYPE(section_vals_type), POINTER :: negf_section, root_section
203
204 CALL timeset(routinen, handle)
205
206 ! ensure we have Quickstep enabled for all force_env
207 NULLIFY (sub_force_env)
208 CALL force_env_get(force_env, in_use=in_use, qs_env=qs_env, root_section=root_section, &
209 sub_force_env=sub_force_env)
210
211 IF (ASSOCIATED(sub_force_env)) THEN
212 n_force_env = SIZE(sub_force_env)
213 ELSE
214 n_force_env = 0
215 END IF
216
217 IF (in_use == use_qs_force) THEN
218 DO icontact = 1, n_force_env
219 CALL force_env_get(sub_force_env(icontact)%force_env, in_use=in_use)
220 IF (in_use /= use_qs_force) EXIT
221 END DO
222 END IF
223
224 IF (in_use /= use_qs_force) THEN
225 cpabort("Quickstep is required for NEGF run.")
226 END IF
227
228 ! check that all mentioned FORCE_EVAL sections are actually present
229 ncontacts = SIZE(negf_control%contacts)
230
231 DO icontact = 1, ncontacts
232 IF (negf_control%contacts(icontact)%force_env_index > n_force_env) THEN
233 WRITE (contact_str, '(I11)') icontact
234 WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
235 WRITE (n_force_env_str, '(I11)') n_force_env
236
237 CALL cp_abort(__location__, &
238 "Contact number "//trim(adjustl(contact_str))//" is linked with the FORCE_EVAL section number "// &
239 trim(adjustl(force_env_str))//", however only "//trim(adjustl(n_force_env_str))// &
240 " FORCE_EVAL sections have been found. Note that FORCE_EVAL sections are enumerated from 0"// &
241 " and that the primary (0-th) section must contain all the atoms.")
242 END IF
243 END DO
244
245 ! create basic matrices and neighbour lists for the primary force_env,
246 ! so we know how matrix elements are actually distributed across CPUs.
247 CALL qs_energies_init(qs_env, calc_forces=.false.)
248 CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, &
249 matrix_s_kp=matrix_s_kp, matrix_ks_kp=matrix_ks_kp, &
250 para_env=para_env, subsys=subsys, v_hartree_rspace=v_hartree_rspace)
251
252 negf_section => section_vals_get_subs_vals(root_section, "NEGF")
253
254 IF (do_kpoints) THEN
255 cpabort("k-points are currently not supported for device FORCE_EVAL")
256 END IF
257
258 ! stage 1: map the atoms between the device force_env and all contact force_env-s
259 ALLOCATE (negf_env%contacts(ncontacts))
260 ALLOCATE (map_contact(ncontacts))
261
262 DO icontact = 1, ncontacts
263 IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
264 CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
265 CALL get_qs_env(qs_env_contact, subsys=subsys_contact)
266
267 CALL negf_env_contact_init_maps(contact_env=negf_env%contacts(icontact), &
268 contact_control=negf_control%contacts(icontact), &
269 atom_map=map_contact(icontact)%atom_map, &
270 eps_geometry=negf_control%eps_geometry, &
271 subsys_device=subsys, &
272 subsys_contact=subsys_contact)
273
274 IF (negf_env%contacts(icontact)%direction_axis == 0) THEN
275 WRITE (contact_str, '(I11)') icontact
276 WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
277 CALL cp_abort(__location__, &
278 "One lattice vector of the contact unit cell (FORCE_EVAL section "// &
279 trim(adjustl(force_env_str))//") must be parallel to the direction of the contact "// &
280 trim(adjustl(contact_str))//".")
281 END IF
282 END IF
283 END DO
284
285 ! stage 2: obtain relevant Kohn-Sham matrix blocks for each contact (separate bulk DFT calculation)
286 DO icontact = 1, ncontacts
287 IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
288 IF (negf_control%contacts(icontact)%read_write_HS) THEN
289 CALL negf_env_contact_read_write_hs &
290 (icontact, sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, &
291 para_env, negf_env, sub_env, negf_control, negf_section, log_unit, is_separate=.true.)
292 ELSE
293 IF (log_unit > 0) &
294 WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') "NEGF| Construct the Kohn-Sham matrix for the contact", icontact, &
295 " from the separate bulk DFT calculation"
296 CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
297 CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
298 CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
299 qs_env_contact=qs_env_contact)
300 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
301 END IF
302 END IF
303 END DO
304
305 ! *** obtain relevant Kohn-Sham matrix blocks for each contact with no separate FORCE_ENV ***
306 is_dft_entire = .false.
307 DO icontact = 1, ncontacts
308 IF (negf_control%contacts(icontact)%force_env_index <= 0) THEN
309 IF (negf_control%contacts(icontact)%read_write_HS) THEN
310 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
311 contact_control=negf_control%contacts(icontact), &
312 sub_env=sub_env, qs_env=qs_env, &
313 eps_geometry=negf_control%eps_geometry)
314 CALL negf_env_contact_read_write_hs(icontact, force_env, para_env, negf_env, sub_env, negf_control, negf_section, &
315 log_unit, is_separate=.false., is_dft_entire=is_dft_entire)
316 ELSE
317 IF (log_unit > 0) &
318 WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') "NEGF| Construct the Kohn-Sham matrix for the contact", icontact, &
319 " from the entire system bulk DFT calculation"
320 IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
321 is_dft_entire = .true.
322 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
323 contact_control=negf_control%contacts(icontact), &
324 sub_env=sub_env, qs_env=qs_env, &
325 eps_geometry=negf_control%eps_geometry)
326 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
327 END IF
328 END IF
329 END DO
330
331 ! stage 3: obtain an initial KS-matrix for the scattering region
332 IF (log_unit > 0) &
333 WRITE (log_unit, '(/,T2,A,T70)') "NEGF| Construct the Kohn-Sham matrix for the scattering region"
334 IF (negf_control%read_write_HS) THEN
335 CALL negf_env_scatt_read_write_hs(force_env, para_env, negf_env, sub_env, negf_control, negf_section, log_unit, &
336 is_dft_entire=is_dft_entire)
337 ELSE
338 IF (.NOT. is_dft_entire) THEN
339 CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
340 is_dft_entire = .true.
341 END IF
342 ! extract device-related matrix blocks
343 CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
344 END IF
345 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
346
347 negf_control%is_dft_entire = is_dft_entire
348
349 ! electron density mixing;
350 ! the input section below should be consistent with the subroutine create_negf_section()
351 NULLIFY (negf_env%mixing_storage)
352 CALL section_vals_val_get(negf_mixing_section, "METHOD", i_val=negf_env%mixing_method)
353
354 CALL get_qs_env(qs_env, dft_control=dft_control)
355 ALLOCATE (negf_env%mixing_storage)
356 CALL mixing_storage_create(negf_env%mixing_storage, negf_mixing_section, &
357 negf_env%mixing_method, dft_control%qs_control%cutoff)
358
359 CALL timestop(handle)
360 END SUBROUTINE negf_env_create
361
362! **************************************************************************************************
363!> \brief Establish mapping between the primary and the contact force environments
364!> \param contact_env NEGF environment for the given contact (modified on exit)
365!> \param contact_control NEGF control
366!> \param atom_map atomic map
367!> \param eps_geometry accuracy in mapping atoms between different force environments
368!> \param subsys_device QuickStep subsystem of the device force environment
369!> \param subsys_contact QuickStep subsystem of the contact force environment
370!> \author Sergey Chulkov
371! **************************************************************************************************
372 SUBROUTINE negf_env_contact_init_maps(contact_env, contact_control, atom_map, &
373 eps_geometry, subsys_device, subsys_contact)
374 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
375 TYPE(negf_control_contact_type), INTENT(in) :: contact_control
376 TYPE(negf_atom_map_type), ALLOCATABLE, &
377 DIMENSION(:), INTENT(inout) :: atom_map
378 REAL(kind=dp), INTENT(in) :: eps_geometry
379 TYPE(qs_subsys_type), POINTER :: subsys_device, subsys_contact
380
381 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_maps'
382
383 INTEGER :: handle, natoms
384
385 CALL timeset(routinen, handle)
386
387 CALL contact_direction_vector(contact_env%origin, &
388 contact_env%direction_vector, &
389 contact_env%origin_bias, &
390 contact_env%direction_vector_bias, &
391 contact_control%atomlist_screening, &
392 contact_control%atomlist_bulk, &
393 subsys_device)
394
395 contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
396
397 IF (contact_env%direction_axis /= 0) THEN
398 natoms = SIZE(contact_control%atomlist_bulk)
399 ALLOCATE (atom_map(natoms))
400
401 ! map atom listed in 'contact_control%atomlist_bulk' to the corresponding atom/cell replica from the contact force_env
402 CALL negf_map_atomic_indices(atom_map=atom_map, &
403 atom_list=contact_control%atomlist_bulk, &
404 subsys_device=subsys_device, &
405 subsys_contact=subsys_contact, &
406 eps_geometry=eps_geometry)
407
408 ! list atoms from 'contact_control%atomlist_bulk' which belong to
409 ! the primary unit cell of the bulk region for the given contact
410 CALL list_atoms_in_bulk_primary_unit_cell(atomlist_cell0=contact_env%atomlist_cell0, &
411 atom_map_cell0=contact_env%atom_map_cell0, &
412 atomlist_bulk=contact_control%atomlist_bulk, &
413 atom_map=atom_map, &
414 origin=contact_env%origin, &
415 direction_vector=contact_env%direction_vector, &
416 direction_axis=contact_env%direction_axis, &
417 subsys_device=subsys_device)
418
419 ! secondary unit cell
420 CALL list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1=contact_env%atomlist_cell1, &
421 atom_map_cell1=contact_env%atom_map_cell1, &
422 atomlist_bulk=contact_control%atomlist_bulk, &
423 atom_map=atom_map, &
424 origin=contact_env%origin, &
425 direction_vector=contact_env%direction_vector, &
426 direction_axis=contact_env%direction_axis, &
427 subsys_device=subsys_device)
428 END IF
429
430 CALL timestop(handle)
431 END SUBROUTINE negf_env_contact_init_maps
432
433! **************************************************************************************************
434!> \brief Reading and writing of the electrode Hamiltonian and overlap matrices from/to a file.
435!> \param icontact ...
436!> \param el_force_env ...
437!> \param para_env ...
438!> \param negf_env ...
439!> \param sub_env ...
440!> \param negf_control ...
441!> \param negf_section ...
442!> \param log_unit ...
443!> \param is_separate ...
444!> \param is_dft_entire ...
445!> \par History
446!> * 12.2025 created [Dmitry Ryndyk]
447! **************************************************************************************************
448 SUBROUTINE negf_env_contact_read_write_hs(icontact, el_force_env, para_env, negf_env, sub_env, negf_control, &
449 negf_section, log_unit, is_separate, is_dft_entire)
450 INTEGER :: icontact
451 TYPE(force_env_type), POINTER :: el_force_env
452 TYPE(mp_para_env_type), POINTER :: para_env
453 TYPE(negf_env_type), INTENT(inout) :: negf_env
454 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
455 TYPE(negf_control_type), POINTER :: negf_control
456 TYPE(section_vals_type), POINTER :: negf_section
457 INTEGER, INTENT(in) :: log_unit
458 LOGICAL, INTENT(in) :: is_separate
459 LOGICAL, INTENT(inout), OPTIONAL :: is_dft_entire
460
461 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_contact_read_write_hs'
462
463 CHARACTER(len=default_path_length) :: filename_h00_1, filename_h00_2, &
464 filename_h01_1, filename_h01_2, &
465 filename_s00, filename_s01
466 INTEGER :: handle, ispin, ncol, nrow, nspins, &
467 print_unit
468 LOGICAL :: exist, exist_all
469 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: target_m
470 TYPE(cp_fm_struct_type), POINTER :: fm_struct
471 TYPE(cp_logger_type), POINTER :: logger
472 TYPE(dft_control_type), POINTER :: dft_control
473 TYPE(qs_environment_type), POINTER :: qs_env_contact
474 TYPE(qs_subsys_type), POINTER :: subsys
475
476 CALL timeset(routinen, handle)
477 logger => cp_get_default_logger()
478
479 CALL force_env_get(el_force_env, qs_env=qs_env_contact)
480 CALL get_qs_env(qs_env_contact, dft_control=dft_control, subsys=subsys)
481 nspins = dft_control%nspins
482
483 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11)') &
484 "NEGF| Construct the Kohn-Sham matrix for the contact", icontact
485
486 ! Check that the files exist.
487 ! ispin=0 is used to show nspins=1
488 exist_all = .true.
489 IF (para_env%is_source()) THEN
490 CALL negf_restart_file_name(filename_s00, exist, negf_section, logger, icontact, s00=.true.)
491 IF (.NOT. exist) THEN
492 CALL cp_warn(__location__, &
493 "User requested to read the overlap matrix from the file named: "// &
494 trim(filename_s00)//". This file does not exist. The file will be created.")
495 exist_all = .false.
496 END IF
497 CALL negf_restart_file_name(filename_s01, exist, negf_section, logger, icontact, s01=.true.)
498 IF (.NOT. exist) THEN
499 CALL cp_warn(__location__, &
500 "User requested to read the overlap matrix from the file named: "// &
501 trim(filename_s01)//". This file does not exist. The file will be created.")
502 exist_all = .false.
503 END IF
504 IF (nspins == 1) THEN
505 CALL negf_restart_file_name(filename_h00_1, exist, negf_section, logger, icontact, ispin=0, h00=.true.)
506 IF (.NOT. exist) THEN
507 CALL cp_warn(__location__, &
508 "User requested to read the Hamiltonian matrix from the file named: "// &
509 trim(filename_h00_1)//". This file does not exist. The file will be created.")
510 exist_all = .false.
511 END IF
512 CALL negf_restart_file_name(filename_h01_1, exist, negf_section, logger, icontact, ispin=0, h01=.true.)
513 IF (.NOT. exist) THEN
514 CALL cp_warn(__location__, &
515 "User requested to read the Hamiltonian matrix from the file named: "// &
516 trim(filename_h01_1)//". This file does not exist. The file will be created.")
517 exist_all = .false.
518 END IF
519 END IF
520 IF (nspins == 2) THEN
521 CALL negf_restart_file_name(filename_h00_1, exist, negf_section, logger, icontact, ispin=1, h00=.true.)
522 IF (.NOT. exist) THEN
523 CALL cp_warn(__location__, &
524 "User requested to read the Hamiltonian matrix from the file named: "// &
525 trim(filename_h00_1)//". This file does not exist. The file will be created.")
526 exist_all = .false.
527 END IF
528 CALL negf_restart_file_name(filename_h01_1, exist, negf_section, logger, icontact, ispin=1, h01=.true.)
529 IF (.NOT. exist) THEN
530 CALL cp_warn(__location__, &
531 "User requested to read tthe Hamiltonian matrix from the file named: "// &
532 trim(filename_h01_1)//". This file does not exist. The file will be created.")
533 exist_all = .false.
534 END IF
535 CALL negf_restart_file_name(filename_h00_2, exist, negf_section, logger, icontact, ispin=2, h00=.true.)
536 IF (.NOT. exist) THEN
537 CALL cp_warn(__location__, &
538 "User requested to read the Hamiltonian matrix from the file named: "// &
539 trim(filename_h00_2)//". This file does not exist. The file will be created.")
540 exist_all = .false.
541 END IF
542 CALL negf_restart_file_name(filename_h01_2, exist, negf_section, logger, icontact, ispin=2, h01=.true.)
543 IF (.NOT. exist) THEN
544 CALL cp_warn(__location__, &
545 "User requested to read the Hamiltonian matrix from the file named: "// &
546 trim(filename_h01_2)//". This file does not exist. The file will be created.")
547 exist_all = .false.
548 END IF
549 END IF
550 END IF
551 CALL para_env%bcast(exist_all)
552
553 IF (exist_all) THEN
554
555 negf_control%contacts(icontact)%is_restart = .true.
556 IF (log_unit > 0) THEN
557 WRITE (log_unit, '(/,T2,A)') "User requested to read the Hamiltonian and overlap matrices from files."
558 WRITE (log_unit, '(T2,A)') "All restart files exist."
559 END IF
560
561 ! ++ create matrices: s_00, s_01, h_00, h_01
562 IF (para_env%is_source()) THEN
563 CALL open_file(file_name=filename_s00, file_status="OLD", &
564 file_form="FORMATTED", file_action="READ", &
565 file_position="REWIND", unit_number=print_unit)
566 READ (print_unit, *) nrow, ncol
567 CALL close_file(print_unit)
568 END IF
569 CALL para_env%bcast(nrow)
570 CALL para_env%bcast(ncol)
571 NULLIFY (fm_struct)
572 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, context=sub_env%blacs_env)
573 ALLOCATE (negf_env%contacts(icontact)%s_00, negf_env%contacts(icontact)%s_01)
574 CALL cp_fm_create(negf_env%contacts(icontact)%s_00, fm_struct)
575 CALL cp_fm_create(negf_env%contacts(icontact)%s_01, fm_struct)
576 ALLOCATE (negf_env%contacts(icontact)%h_00(nspins), negf_env%contacts(icontact)%h_01(nspins))
577 DO ispin = 1, nspins
578 CALL cp_fm_create(negf_env%contacts(icontact)%h_00(ispin), fm_struct)
579 CALL cp_fm_create(negf_env%contacts(icontact)%h_01(ispin), fm_struct)
580 END DO
581 CALL cp_fm_struct_release(fm_struct)
582
583 ALLOCATE (target_m(nrow, ncol))
584 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s00, target_m)
585 CALL para_env%bcast(target_m)
586 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%s_00, target_m)
587 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_00 is read from "//trim(filename_s00)
588 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s01, target_m)
589 CALL para_env%bcast(target_m)
590 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%s_01, target_m)
591 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_01 is read from "//trim(filename_s01)
592 IF (nspins == 1) THEN
593 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_1, target_m)
594 CALL para_env%bcast(target_m)
595 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
596 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//trim(filename_h00_1)
597 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_1, target_m)
598 CALL para_env%bcast(target_m)
599 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
600 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//trim(filename_h01_1)
601 END IF
602 IF (nspins == 2) THEN
603 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_1, target_m)
604 CALL para_env%bcast(target_m)
605 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
606 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//trim(filename_h00_1)//" for spin 1"
607 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_1, target_m)
608 CALL para_env%bcast(target_m)
609 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
610 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//trim(filename_h01_1)//" for spin 1"
611 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_2, target_m)
612 CALL para_env%bcast(target_m)
613 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(2), target_m)
614 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//trim(filename_h00_2)//" for spin 2"
615 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_2, target_m)
616 CALL para_env%bcast(target_m)
617 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(2), target_m)
618 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//trim(filename_h01_2)//" for spin 2"
619 END IF
620 DEALLOCATE (target_m)
621
622 ELSE
623
624 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') &
625 "Some restart files do not exist. ALL restart files will be recalculated!"
626
627 IF (is_separate) THEN
628 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') &
629 "Construct the Kohn-Sham matrix from from the separate bulk DFT calculation"
630 CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
631 CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
632 qs_env_contact=qs_env_contact)
633 ELSE
634 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') &
635 "Construct the Kohn-Sham matrix from the entire system bulk DFT calculation"
636 negf_control%contacts(icontact)%read_write_HS = .false.
637 IF (.NOT. is_dft_entire) CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
638 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
639 contact_control=negf_control%contacts(icontact), &
640 sub_env=sub_env, qs_env=qs_env_contact, &
641 eps_geometry=negf_control%eps_geometry)
642 negf_control%contacts(icontact)%read_write_HS = .true.
643 is_dft_entire = .true.
644 END IF
645
646 CALL cp_fm_get_info(negf_env%contacts(icontact)%s_00, nrow_global=nrow)
647 ALLOCATE (target_m(nrow, nrow))
648 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%s_00, target_m)
649 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s00, target_m)
650 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "S_00 is saved to "//trim(filename_s00)
651 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%s_01, target_m)
652 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s01, target_m)
653 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_01 is saved to "//trim(filename_s01)
654 IF (nspins == 1) THEN
655 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
656 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_1, target_m)
657 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//trim(filename_h00_1)
658 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
659 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_1, target_m)
660 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//trim(filename_h01_1)
661 END IF
662 IF (nspins == 2) THEN
663 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
664 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_1, target_m)
665 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//trim(filename_h00_1)//" for spin 1"
666 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
667 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_1, target_m)
668 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//trim(filename_h01_1)//" for spin 1"
669 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(2), target_m)
670 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_2, target_m)
671 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//trim(filename_h00_2)//" for spin 2"
672 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(2), target_m)
673 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_2, target_m)
674 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//trim(filename_h01_2)//" for spin 2"
675 END IF
676 DEALLOCATE (target_m)
677
678 negf_control%write_common_restart_file = .true.
679
680 END IF
681
682 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
683
684 CALL timestop(handle)
685 END SUBROUTINE negf_env_contact_read_write_hs
686
687! **************************************************************************************************
688!> \brief Extract relevant matrix blocks for the given contact.
689!> \param contact_env NEGF environment for the contact (modified on exit)
690!> \param sub_env NEGF parallel (sub)group environment
691!> \param qs_env_contact QuickStep environment for the contact force environment
692!> \par History
693!> * 10.2017 created [Sergey Chulkov]
694!> * 10.2025 The subroutine is essentially modified. New functionality of negf_copy_contact_matrix.
695!> [Dmitry Ryndyk]
696! **************************************************************************************************
697 SUBROUTINE negf_env_contact_init_matrices(contact_env, sub_env, qs_env_contact)
698 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
699 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
700 TYPE(qs_environment_type), POINTER :: qs_env_contact
701
702 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_matrices'
703
704 INTEGER :: handle, iatom, ispin, nao, natoms, &
705 nimages, nspins
706 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_list0, atom_list1
707 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
708 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
709 LOGICAL :: do_kpoints
710 TYPE(cp_fm_struct_type), POINTER :: fm_struct
711 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matkp
712 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
713 TYPE(dft_control_type), POINTER :: dft_control
714 TYPE(kpoint_type), POINTER :: kpoints
715 TYPE(mp_para_env_type), POINTER :: para_env
716 TYPE(qs_rho_type), POINTER :: rho_struct
717 TYPE(qs_subsys_type), POINTER :: subsys
718
719 CALL timeset(routinen, handle)
720
721 CALL get_qs_env(qs_env_contact, &
722 dft_control=dft_control, &
723 do_kpoints=do_kpoints, &
724 kpoints=kpoints, &
725 matrix_ks_kp=matrix_ks_kp, &
726 matrix_s_kp=matrix_s_kp, &
727 para_env=para_env, &
728 rho=rho_struct, &
729 subsys=subsys)
730 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
731
732 CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
733
734 natoms = SIZE(contact_env%atomlist_cell0)
735 ALLOCATE (atom_list0(natoms))
736 DO iatom = 1, natoms
737 atom_list0(iatom) = contact_env%atom_map_cell0(iatom)%iatom
738
739 ! with no k-points there is one-to-one correspondence between the primary unit cell
740 ! of the contact force_env and the first contact unit cell of the device force_env
741 IF (sum(abs(contact_env%atom_map_cell0(iatom)%cell(:))) > 0) THEN
742 cpabort("NEGF K-points are not currently supported")
743 END IF
744 END DO
745
746 cpassert(SIZE(contact_env%atomlist_cell1) == natoms)
747 ALLOCATE (atom_list1(natoms))
748 DO iatom = 1, natoms
749 atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
750 END DO
751
752 nspins = dft_control%nspins
753 nimages = dft_control%nimages
754
755 IF (do_kpoints) THEN
756 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
757 ELSE
758 ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
759 cell_to_index(0, 0, 0) = 1
760 END IF
761
762 ALLOCATE (index_to_cell(3, nimages))
763 CALL invert_cell_to_index(cell_to_index, nimages, index_to_cell)
764 IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
765
766 NULLIFY (fm_struct)
767 nao = number_of_atomic_orbitals(subsys, atom_list0)
768 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env)
769
770 ! ++ create matrices: s_00, s_01
771 ALLOCATE (contact_env%s_00, contact_env%s_01)
772 CALL cp_fm_create(contact_env%s_00, fm_struct)
773 CALL cp_fm_create(contact_env%s_01, fm_struct)
774
775 ! ++ create matrices: h_00, h_01, rho_00, rho_01
776 ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
777 ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
778 DO ispin = 1, nspins
779 CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
780 CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
781 CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
782 CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
783 END DO
784
785 CALL cp_fm_struct_release(fm_struct)
786
787 ! extract matrices: s_00, s_01
788 matkp => matrix_s_kp(1, :)
789 CALL negf_copy_contact_matrix(fm_cell0=contact_env%s_00, &
790 fm_cell1=contact_env%s_01, &
791 direction_axis=contact_env%direction_axis, &
792 matrix_kp=matkp, &
793 atom_list0=atom_list0, atom_list1=atom_list1, &
794 subsys=subsys, mpi_comm_global=para_env, &
795 kpoints=kpoints)
796
797 ! extract matrices: h_00, h_01, rho_00, rho_01
798 DO ispin = 1, nspins
799 matkp => matrix_ks_kp(ispin, :)
800 CALL negf_copy_contact_matrix(fm_cell0=contact_env%h_00(ispin), &
801 fm_cell1=contact_env%h_01(ispin), &
802 direction_axis=contact_env%direction_axis, &
803 matrix_kp=matkp, &
804 atom_list0=atom_list0, atom_list1=atom_list1, &
805 subsys=subsys, mpi_comm_global=para_env, &
806 kpoints=kpoints)
807
808 matkp => rho_ao_kp(ispin, :)
809 CALL negf_copy_contact_matrix(fm_cell0=contact_env%rho_00(ispin), &
810 fm_cell1=contact_env%rho_01(ispin), &
811 direction_axis=contact_env%direction_axis, &
812 matrix_kp=matkp, &
813 atom_list0=atom_list0, atom_list1=atom_list1, &
814 subsys=subsys, mpi_comm_global=para_env, &
815 kpoints=kpoints)
816 END DO
817
818 DEALLOCATE (atom_list0, atom_list1)
819
820 CALL timestop(handle)
821 END SUBROUTINE negf_env_contact_init_matrices
822
823! **************************************************************************************************
824!> \brief Extract relevant matrix blocks for the given contact using the device's force environment.
825!> \param contact_env NEGF environment for the contact (modified on exit)
826!> \param contact_control NEGF control for the contact
827!> \param sub_env NEGF parallel (sub)group environment
828!> \param qs_env QuickStep environment for the device force environment
829!> \param eps_geometry accuracy in Cartesian coordinates
830!> \author Sergey Chulkov
831! **************************************************************************************************
832 SUBROUTINE negf_env_contact_init_matrices_gamma(contact_env, contact_control, sub_env, qs_env, eps_geometry)
833 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
834 TYPE(negf_control_contact_type), INTENT(in) :: contact_control
835 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
836 TYPE(qs_environment_type), POINTER :: qs_env
837 REAL(kind=dp), INTENT(in) :: eps_geometry
838
839 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_matrices_gamma'
840
841 INTEGER :: handle, iatom, icell, ispin, nao_c, &
842 nspins
843 LOGICAL :: do_kpoints
844 REAL(kind=dp), DIMENSION(2) :: r2_origin_cell
845 REAL(kind=dp), DIMENSION(3) :: direction_vector, origin
846 TYPE(cp_fm_struct_type), POINTER :: fm_struct
847 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
848 TYPE(dft_control_type), POINTER :: dft_control
849 TYPE(mp_para_env_type), POINTER :: para_env
850 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
851 TYPE(qs_rho_type), POINTER :: rho_struct
852 TYPE(qs_subsys_type), POINTER :: subsys
853
854 CALL timeset(routinen, handle)
855
856 CALL get_qs_env(qs_env, &
857 dft_control=dft_control, &
858 do_kpoints=do_kpoints, &
859 matrix_ks_kp=matrix_ks_kp, &
860 matrix_s_kp=matrix_s_kp, &
861 para_env=para_env, &
862 rho=rho_struct, &
863 subsys=subsys)
864 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
865
866 IF (do_kpoints) THEN
867 CALL cp_abort(__location__, &
868 "K-points in device region have not been implemented yet.")
869 END IF
870
871 nspins = dft_control%nspins
872
873 nao_c = number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(1)%vector)
874 IF (number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(2)%vector) /= nao_c) THEN
875 CALL cp_abort(__location__, &
876 "Primary and secondary bulk contact cells should be identical "// &
877 "in terms of the number of atoms of each kind, and their basis sets. "// &
878 "No single atom, however, can be shared between these two cells.")
879 END IF
880
881 contact_env%homo_energy = 0.0_dp
882
883 CALL contact_direction_vector(contact_env%origin, &
884 contact_env%direction_vector, &
885 contact_env%origin_bias, &
886 contact_env%direction_vector_bias, &
887 contact_control%atomlist_screening, &
888 contact_control%atomlist_bulk, &
889 subsys)
890
891 contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys, eps_geometry)
892
893 ! choose the primary and secondary contact unit cells
894 CALL qs_subsys_get(subsys, particle_set=particle_set)
895
896 origin = particle_set(contact_control%atomlist_screening(1))%r
897 DO iatom = 2, SIZE(contact_control%atomlist_screening)
898 origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
899 END DO
900 origin = origin/real(SIZE(contact_control%atomlist_screening), kind=dp)
901
902 DO icell = 1, 2
903 direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
904 DO iatom = 2, SIZE(contact_control%atomlist_cell(icell)%vector)
905 direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
906 END DO
907 direction_vector = direction_vector/real(SIZE(contact_control%atomlist_cell(icell)%vector), kind=dp)
908 direction_vector = direction_vector - origin
909 r2_origin_cell(icell) = dot_product(direction_vector, direction_vector)
910 END DO
911
912 IF (abs(r2_origin_cell(1) - r2_origin_cell(2)) < (eps_geometry*eps_geometry)) THEN
913 ! primary and secondary bulk unit cells should not overlap;
914 ! currently we check that they are different by at least one atom that is, indeed, not sufficient.
915 CALL cp_abort(__location__, &
916 "Primary and secondary bulk contact cells should not overlap ")
917 ELSE IF (r2_origin_cell(1) < r2_origin_cell(2)) THEN
918 IF (.NOT. ALLOCATED(contact_env%atomlist_cell0)) &
919 ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(1)%vector)))
920 contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
921 IF (.NOT. ALLOCATED(contact_env%atomlist_cell1)) &
922 ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(2)%vector)))
923 contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
924 ELSE
925 IF (.NOT. ALLOCATED(contact_env%atomlist_cell0)) &
926 ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(2)%vector)))
927 contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
928 IF (.NOT. ALLOCATED(contact_env%atomlist_cell1)) &
929 ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(1)%vector)))
930 contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
931 END IF
932 IF (.NOT. contact_control%read_write_HS) THEN
933 NULLIFY (fm_struct)
934 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
935 ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
936 ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
937 DO ispin = 1, nspins
938 CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
939 CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
940 CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
941 CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
942 END DO
943 ALLOCATE (contact_env%s_00, contact_env%s_01)
944 CALL cp_fm_create(contact_env%s_00, fm_struct)
945 CALL cp_fm_create(contact_env%s_01, fm_struct)
946 CALL cp_fm_struct_release(fm_struct)
947
948 DO ispin = 1, nspins
949 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
950 fm=contact_env%h_00(ispin), &
951 atomlist_row=contact_env%atomlist_cell0, &
952 atomlist_col=contact_env%atomlist_cell0, &
953 subsys=subsys, mpi_comm_global=para_env, &
954 do_upper_diag=.true., do_lower=.true.)
955 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
956 fm=contact_env%h_01(ispin), &
957 atomlist_row=contact_env%atomlist_cell0, &
958 atomlist_col=contact_env%atomlist_cell1, &
959 subsys=subsys, mpi_comm_global=para_env, &
960 do_upper_diag=.true., do_lower=.true.)
961
962 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
963 fm=contact_env%rho_00(ispin), &
964 atomlist_row=contact_env%atomlist_cell0, &
965 atomlist_col=contact_env%atomlist_cell0, &
966 subsys=subsys, mpi_comm_global=para_env, &
967 do_upper_diag=.true., do_lower=.true.)
968 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
969 fm=contact_env%rho_01(ispin), &
970 atomlist_row=contact_env%atomlist_cell0, &
971 atomlist_col=contact_env%atomlist_cell1, &
972 subsys=subsys, mpi_comm_global=para_env, &
973 do_upper_diag=.true., do_lower=.true.)
974 END DO
975
976 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
977 fm=contact_env%s_00, &
978 atomlist_row=contact_env%atomlist_cell0, &
979 atomlist_col=contact_env%atomlist_cell0, &
980 subsys=subsys, mpi_comm_global=para_env, &
981 do_upper_diag=.true., do_lower=.true.)
982 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
983 fm=contact_env%s_01, &
984 atomlist_row=contact_env%atomlist_cell0, &
985 atomlist_col=contact_env%atomlist_cell1, &
986 subsys=subsys, mpi_comm_global=para_env, &
987 do_upper_diag=.true., do_lower=.true.)
988 END IF
989 CALL timestop(handle)
990 END SUBROUTINE negf_env_contact_init_matrices_gamma
991
992! **************************************************************************************************
993!> \brief Reading and writing of the electrode Hamiltonian and overlap matrices from/to a file.
994!> \param force_env ...
995!> \param para_env ...
996!> \param negf_env ...
997!> \param sub_env ...
998!> \param negf_control ...
999!> \param negf_section ...
1000!> \param log_unit ...
1001!> \param is_dft_entire ...
1002!> \par History
1003!> * 01.2026 created [Dmitry Ryndyk]
1004! **************************************************************************************************
1005 SUBROUTINE negf_env_scatt_read_write_hs(force_env, para_env, negf_env, sub_env, negf_control, negf_section, &
1006 log_unit, is_dft_entire)
1007 TYPE(force_env_type), POINTER :: force_env
1008 TYPE(mp_para_env_type), POINTER :: para_env
1009 TYPE(negf_env_type), INTENT(inout) :: negf_env
1010 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1011 TYPE(negf_control_type), POINTER :: negf_control
1012 TYPE(section_vals_type), POINTER :: negf_section
1013 INTEGER, INTENT(in) :: log_unit
1014 LOGICAL, INTENT(inout), OPTIONAL :: is_dft_entire
1015
1016 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_scatt_read_write_hs'
1017
1018 CHARACTER(len=default_path_length) :: filename_h_1, filename_h_2, filename_s
1019 CHARACTER(len=default_path_length), ALLOCATABLE, &
1020 DIMENSION(:) :: filename_hc_1, filename_hc_2, filename_sc
1021 INTEGER :: handle, icontact, ispin, ncol_s, &
1022 ncol_sc, ncontacts, nrow_s, nrow_sc, &
1023 nspins, print_unit
1024 LOGICAL :: exist, exist_all
1025 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: target_m
1026 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1027 TYPE(cp_logger_type), POINTER :: logger
1028 TYPE(dft_control_type), POINTER :: dft_control
1029 TYPE(qs_environment_type), POINTER :: qs_env
1030 TYPE(qs_subsys_type), POINTER :: subsys
1031
1032 CALL timeset(routinen, handle)
1033 logger => cp_get_default_logger()
1034
1035 CALL force_env_get(force_env, qs_env=qs_env)
1036 CALL get_qs_env(qs_env, dft_control=dft_control, subsys=subsys)
1037 ncontacts = SIZE(negf_control%contacts)
1038 nspins = dft_control%nspins
1039 ALLOCATE (filename_sc(ncontacts), filename_hc_1(ncontacts), filename_hc_2(ncontacts))
1040
1041 ! Check that the files exist.
1042 ! ispin=0 is used to show nspins=1
1043 exist_all = .true.
1044 IF (para_env%is_source()) THEN
1045 CALL negf_restart_file_name(filename_s, exist, negf_section, logger, s=.true.)
1046 IF (.NOT. exist) THEN
1047 CALL cp_warn(__location__, &
1048 "User requested to read the overlap matrix from the file named: "// &
1049 trim(filename_s)//". This file does not exist. The file will be created.")
1050 exist_all = .false.
1051 END IF
1052 IF (nspins == 1) THEN
1053 CALL negf_restart_file_name(filename_h_1, exist, negf_section, logger, ispin=0, h=.true.)
1054 IF (.NOT. exist) THEN
1055 CALL cp_warn(__location__, &
1056 "User requested to read the Hamiltonian matrix from the file named: "// &
1057 trim(filename_h_1)//". This file does not exist. The file will be created.")
1058 exist_all = .false.
1059 END IF
1060 END IF
1061 IF (nspins == 2) THEN
1062 CALL negf_restart_file_name(filename_h_1, exist, negf_section, logger, ispin=1, h=.true.)
1063 IF (.NOT. exist) THEN
1064 CALL cp_warn(__location__, &
1065 "User requested to read the Hamiltonian matrix from the file named: "// &
1066 trim(filename_h_1)//". This file does not exist. The file will be created.")
1067 exist_all = .false.
1068 END IF
1069 CALL negf_restart_file_name(filename_h_2, exist, negf_section, logger, ispin=2, h=.true.)
1070 IF (.NOT. exist) THEN
1071 CALL cp_warn(__location__, &
1072 "User requested to read the Hamiltonian matrix from the file named: "// &
1073 trim(filename_h_2)//". This file does not exist. The file will be created.")
1074 exist_all = .false.
1075 END IF
1076 END IF
1077 DO icontact = 1, ncontacts
1078 CALL negf_restart_file_name(filename_sc(icontact), exist, negf_section, logger, icontact=icontact, sc=.true.)
1079 IF (.NOT. exist) THEN
1080 CALL cp_warn(__location__, &
1081 "User requested to read the overlap matrix from the file named: "// &
1082 trim(filename_sc(icontact))//". This file does not exist. The file will be created.")
1083 exist_all = .false.
1084 END IF
1085 IF (nspins == 1) THEN
1086 CALL negf_restart_file_name(filename_hc_1(icontact), exist, negf_section, logger, icontact=icontact, &
1087 ispin=0, hc=.true.)
1088 IF (.NOT. exist) THEN
1089 CALL cp_warn(__location__, &
1090 "User requested to read the Hamiltonian matrix from the file named: "// &
1091 trim(filename_hc_1(icontact))//". This file does not exist. The file will be created.")
1092 exist_all = .false.
1093 END IF
1094 END IF
1095 IF (nspins == 2) THEN
1096 CALL negf_restart_file_name(filename_hc_1(icontact), exist, negf_section, logger, icontact=icontact, &
1097 ispin=1, hc=.true.)
1098 IF (.NOT. exist) THEN
1099 CALL cp_warn(__location__, &
1100 "User requested to read the Hamiltonian matrix from the file named: "// &
1101 trim(filename_hc_1(icontact))//". This file does not exist. The file will be created.")
1102 exist_all = .false.
1103 END IF
1104 CALL negf_restart_file_name(filename_hc_2(icontact), exist, negf_section, logger, icontact=icontact, &
1105 ispin=2, hc=.true.)
1106 IF (.NOT. exist) THEN
1107 CALL cp_warn(__location__, &
1108 "User requested to read the Hamiltonian matrix from the file named: "// &
1109 trim(filename_hc_2(icontact))//". This file does not exist. The file will be created.")
1110 exist_all = .false.
1111 END IF
1112 END IF
1113 END DO
1114 END IF
1115 CALL para_env%bcast(exist_all)
1116
1117 IF (exist_all) THEN
1118
1119 negf_control%is_restart = .true.
1120
1121 IF (log_unit > 0) THEN
1122 WRITE (log_unit, '(/,T2,A)') "User requested to read the Hamiltonian and overlap matrices from files."
1123 WRITE (log_unit, '(T2,A)') "All restart files exist."
1124 END IF
1125
1126 ! ++ create matrices: s_s, s_sc, h_s, h_sc
1127 IF (para_env%is_source()) THEN
1128 CALL open_file(file_name=filename_s, file_status="OLD", &
1129 file_form="FORMATTED", file_action="READ", &
1130 file_position="REWIND", unit_number=print_unit)
1131 READ (print_unit, *) nrow_s, ncol_s
1132 CALL close_file(print_unit)
1133 END IF
1134 CALL para_env%bcast(nrow_s)
1135 CALL para_env%bcast(ncol_s)
1136 NULLIFY (fm_struct)
1137 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_s, ncol_global=ncol_s, context=sub_env%blacs_env)
1138 ALLOCATE (negf_env%s_s)
1139 CALL cp_fm_create(negf_env%s_s, fm_struct)
1140 ALLOCATE (negf_env%h_s(nspins))
1141 DO ispin = 1, nspins
1142 CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
1143 END DO
1144 CALL cp_fm_struct_release(fm_struct)
1145 ALLOCATE (negf_env%s_sc(ncontacts))
1146 ALLOCATE (negf_env%h_sc(nspins, ncontacts))
1147 DO icontact = 1, ncontacts
1148 IF (para_env%is_source()) THEN
1149 CALL open_file(file_name=filename_sc(icontact), file_status="OLD", &
1150 file_form="FORMATTED", file_action="READ", &
1151 file_position="REWIND", unit_number=print_unit)
1152 READ (print_unit, *) nrow_sc, ncol_sc
1153 CALL close_file(print_unit)
1154 END IF
1155 CALL para_env%bcast(nrow_sc)
1156 CALL para_env%bcast(ncol_sc)
1157 NULLIFY (fm_struct)
1158 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_sc, ncol_global=ncol_sc, context=sub_env%blacs_env)
1159 CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
1160 DO ispin = 1, nspins
1161 CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
1162 END DO
1163 CALL cp_fm_struct_release(fm_struct)
1164 END DO
1165
1166 ALLOCATE (target_m(nrow_s, ncol_s))
1167 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s, target_m)
1168 CALL para_env%bcast(target_m)
1169 CALL cp_fm_set_submatrix(negf_env%s_s, target_m)
1170 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_s is read from "//trim(filename_s)
1171 IF (nspins == 1) THEN
1172 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_1, target_m)
1173 CALL para_env%bcast(target_m)
1174 CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m)
1175 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//trim(filename_h_1)
1176 END IF
1177 IF (nspins == 2) THEN
1178 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_1, target_m)
1179 CALL para_env%bcast(target_m)
1180 CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m)
1181 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//trim(filename_h_1)//" for spin 1"
1182 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_2, target_m)
1183 CALL para_env%bcast(target_m)
1184 CALL cp_fm_set_submatrix(negf_env%h_s(2), target_m)
1185 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//trim(filename_h_2)//" for spin 2"
1186 END IF
1187 DEALLOCATE (target_m)
1188
1189 DO icontact = 1, ncontacts
1190 ALLOCATE (target_m(nrow_s, ncol_sc))
1191 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_sc(icontact), target_m)
1192 CALL para_env%bcast(target_m)
1193 CALL cp_fm_set_submatrix(negf_env%s_sc(icontact), target_m)
1194 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_sc is read from "//trim(filename_sc(icontact))
1195 IF (nspins == 1) THEN
1196 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_1(icontact), target_m)
1197 CALL para_env%bcast(target_m)
1198 CALL cp_fm_set_submatrix(negf_env%h_sc(1, icontact), target_m)
1199 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//trim(filename_hc_1(icontact))
1200 END IF
1201 IF (nspins == 2) THEN
1202 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_1(icontact), target_m)
1203 CALL para_env%bcast(target_m)
1204 CALL cp_fm_set_submatrix(negf_env%h_sc(1, icontact), target_m)
1205 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//trim(filename_hc_1(icontact))//" for spin 1"
1206 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_2(icontact), target_m)
1207 CALL para_env%bcast(target_m)
1208 CALL cp_fm_set_submatrix(negf_env%h_sc(2, icontact), target_m)
1209 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//trim(filename_hc_2(icontact))//" for spin 2"
1210 END IF
1211 DEALLOCATE (target_m)
1212
1213 END DO
1214
1215 ELSE
1216
1217 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') &
1218 "Some restart files do not exist. ALL restart files will be recalculated!"
1219
1220 IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
1221 ! extract device-related matrix blocks
1222 CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
1223 is_dft_entire = .true.
1224
1225 CALL cp_fm_get_info(negf_env%s_s, nrow_global=nrow_s, ncol_global=ncol_s)
1226 ALLOCATE (target_m(nrow_s, ncol_s))
1227 CALL cp_fm_get_submatrix(negf_env%s_s, target_m)
1228 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s, target_m)
1229 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "S_s is saved to "//trim(filename_s)
1230 IF (nspins == 1) THEN
1231 CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m)
1232 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_1, target_m)
1233 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//trim(filename_h_1)
1234 END IF
1235 IF (nspins == 2) THEN
1236 CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m)
1237 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_1, target_m)
1238 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//trim(filename_h_1)//" for spin 1"
1239 CALL cp_fm_get_submatrix(negf_env%h_s(2), target_m)
1240 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_2, target_m)
1241 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//trim(filename_h_2)//" for spin 2"
1242 END IF
1243 DEALLOCATE (target_m)
1244
1245 DO icontact = 1, ncontacts
1246 CALL cp_fm_get_info(negf_env%contacts(icontact)%s_00, nrow_global=nrow_sc, ncol_global=ncol_sc)
1247 ALLOCATE (target_m(nrow_s, ncol_sc))
1248 CALL cp_fm_get_submatrix(negf_env%s_sc(icontact), target_m)
1249 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_sc(icontact), target_m)
1250 IF (log_unit > 0) WRITE (log_unit, '(T2,A,I3)') &
1251 "S_sc is saved to "//trim(filename_sc(icontact))//" for contact", icontact
1252 IF (nspins == 1) THEN
1253 CALL cp_fm_get_submatrix(negf_env%h_sc(1, icontact), target_m)
1254 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_1(icontact), target_m)
1255 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//trim(filename_hc_1(icontact))
1256 END IF
1257 IF (nspins == 2) THEN
1258 CALL cp_fm_get_submatrix(negf_env%h_sc(1, icontact), target_m)
1259 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_1(icontact), target_m)
1260 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//trim(filename_hc_1(icontact))//" for spin 1"
1261 CALL cp_fm_get_submatrix(negf_env%h_sc(2, icontact), target_m)
1262 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_2(icontact), target_m)
1263 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//trim(filename_hc_2(icontact))//" for spin 2"
1264 END IF
1265 DEALLOCATE (target_m)
1266 END DO
1267
1268 negf_control%write_common_restart_file = .true.
1269
1270 END IF
1271
1272 DEALLOCATE (filename_sc, filename_hc_1, filename_hc_2)
1273 CALL timestop(handle)
1274 END SUBROUTINE negf_env_scatt_read_write_hs
1275
1276! **************************************************************************************************
1277!> \brief Extract relevant matrix blocks for the scattering region as well as
1278!> all the scattering -- contact interface regions.
1279!> \param negf_env NEGF environment (modified on exit)
1280!> \param negf_control NEGF control
1281!> \param sub_env NEGF parallel (sub)group environment
1282!> \param qs_env Primary QuickStep environment
1283!> \author Sergey Chulkov
1284! **************************************************************************************************
1285 SUBROUTINE negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
1286 TYPE(negf_env_type), INTENT(inout) :: negf_env
1287 TYPE(negf_control_type), POINTER :: negf_control
1288 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1289 TYPE(qs_environment_type), POINTER :: qs_env
1290
1291 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_device_init_matrices'
1292
1293 INTEGER :: handle, icontact, ispin, nao_c, nao_s, &
1294 ncontacts, nspins
1295 LOGICAL :: do_kpoints
1296 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1297 TYPE(dbcsr_p_type) :: hmat
1298 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
1299 TYPE(dft_control_type), POINTER :: dft_control
1300 TYPE(mp_para_env_type), POINTER :: para_env
1301 TYPE(pw_env_type), POINTER :: pw_env
1302 TYPE(pw_pool_type), POINTER :: pw_pool
1303 TYPE(pw_r3d_rs_type) :: v_hartree
1304 TYPE(qs_subsys_type), POINTER :: subsys
1305
1306 CALL timeset(routinen, handle)
1307
1308 IF (ALLOCATED(negf_control%atomlist_S_screening)) THEN
1309 CALL get_qs_env(qs_env, &
1310 dft_control=dft_control, &
1311 do_kpoints=do_kpoints, &
1312 matrix_ks_kp=matrix_ks_kp, &
1313 matrix_s_kp=matrix_s_kp, &
1314 para_env=para_env, &
1315 pw_env=pw_env, &
1316 subsys=subsys)
1317 IF (dft_control%qs_control%xtb) CALL rebuild_pw_env(qs_env)
1318 CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
1319
1320 IF (do_kpoints) THEN
1321 CALL cp_abort(__location__, &
1322 "K-points in device region have not been implemented yet.")
1323 END IF
1324
1325 ncontacts = SIZE(negf_control%contacts)
1326 nspins = dft_control%nspins
1327
1328 NULLIFY (fm_struct)
1329 nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening)
1330
1331 ! ++ create matrices: h_s, s_s
1332 NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
1333 ALLOCATE (negf_env%h_s(nspins))
1334
1335 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
1336 ALLOCATE (negf_env%s_s)
1337 CALL cp_fm_create(negf_env%s_s, fm_struct)
1338 DO ispin = 1, nspins
1339 CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
1340 END DO
1341 ALLOCATE (negf_env%v_hartree_s)
1342 CALL cp_fm_create(negf_env%v_hartree_s, fm_struct)
1343 CALL cp_fm_struct_release(fm_struct)
1344
1345 ! ++ create matrices: h_sc, s_sc
1346 ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
1347 DO icontact = 1, ncontacts
1348 nao_c = number_of_atomic_orbitals(subsys, negf_env%contacts(icontact)%atomlist_cell0)
1349 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
1350
1351 CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
1352
1353 DO ispin = 1, nspins
1354 CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
1355 END DO
1356
1357 CALL cp_fm_struct_release(fm_struct)
1358 END DO
1359
1360 ! extract matrices: h_s, s_s
1361 DO ispin = 1, nspins
1362 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
1363 fm=negf_env%h_s(ispin), &
1364 atomlist_row=negf_control%atomlist_S_screening, &
1365 atomlist_col=negf_control%atomlist_S_screening, &
1366 subsys=subsys, mpi_comm_global=para_env, &
1367 do_upper_diag=.true., do_lower=.true.)
1368 END DO
1369
1370 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
1371 fm=negf_env%s_s, &
1372 atomlist_row=negf_control%atomlist_S_screening, &
1373 atomlist_col=negf_control%atomlist_S_screening, &
1374 subsys=subsys, mpi_comm_global=para_env, &
1375 do_upper_diag=.true., do_lower=.true.)
1376
1377 ! v_hartree_s
1378 NULLIFY (hmat%matrix)
1379 CALL dbcsr_init_p(hmat%matrix)
1380 CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
1381 CALL dbcsr_set(hmat%matrix, 0.0_dp)
1382
1383 CALL pw_pool%create_pw(v_hartree)
1384 CALL negf_env_init_v_hartree(v_hartree, negf_env%contacts, negf_control%contacts)
1385
1386 CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, &
1387 calculate_forces=.false., compute_tau=.false., gapw=.false.)
1388
1389 CALL pw_pool%give_back_pw(v_hartree)
1390
1391 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=hmat%matrix, &
1392 fm=negf_env%v_hartree_s, &
1393 atomlist_row=negf_control%atomlist_S_screening, &
1394 atomlist_col=negf_control%atomlist_S_screening, &
1395 subsys=subsys, mpi_comm_global=para_env, &
1396 do_upper_diag=.true., do_lower=.true.)
1397
1398 CALL dbcsr_deallocate_matrix(hmat%matrix)
1399
1400 ! extract matrices: h_sc, s_sc
1401 DO icontact = 1, ncontacts
1402 DO ispin = 1, nspins
1403 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
1404 fm=negf_env%h_sc(ispin, icontact), &
1405 atomlist_row=negf_control%atomlist_S_screening, &
1406 atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
1407 subsys=subsys, mpi_comm_global=para_env, &
1408 do_upper_diag=.true., do_lower=.true.)
1409 END DO
1410
1411 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
1412 fm=negf_env%s_sc(icontact), &
1413 atomlist_row=negf_control%atomlist_S_screening, &
1414 atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
1415 subsys=subsys, mpi_comm_global=para_env, &
1416 do_upper_diag=.true., do_lower=.true.)
1417 END DO
1418 END IF
1419
1420 CALL timestop(handle)
1421 END SUBROUTINE negf_env_device_init_matrices
1422
1423! **************************************************************************************************
1424!> \brief Contribution to the Hartree potential related to the external bias voltage.
1425!> \param v_hartree Hartree potential (modified on exit)
1426!> \param contact_env NEGF environment for every contact
1427!> \param contact_control NEGF control for every contact
1428!> \author Sergey Chulkov
1429! **************************************************************************************************
1430 SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control)
1431 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree
1432 TYPE(negf_env_contact_type), DIMENSION(:), &
1433 INTENT(in) :: contact_env
1434 TYPE(negf_control_contact_type), DIMENSION(:), &
1435 INTENT(in) :: contact_control
1436
1437 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_init_v_hartree'
1438 REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*epsilon(0.0_dp)
1439
1440 INTEGER :: dx, dy, dz, handle, icontact, ix, iy, &
1441 iz, lx, ly, lz, ncontacts, ux, uy, uz
1442 REAL(kind=dp) :: dvol, pot, proj, v1, v2
1443 REAL(kind=dp), DIMENSION(3) :: dirvector_bias, point_coord, &
1444 point_indices, vector
1445
1446 CALL timeset(routinen, handle)
1447
1448 ncontacts = SIZE(contact_env)
1449 cpassert(SIZE(contact_control) == ncontacts)
1450 cpassert(ncontacts == 2)
1451
1452 dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
1453 v1 = contact_control(1)%v_external
1454 v2 = contact_control(2)%v_external
1455
1456 lx = v_hartree%pw_grid%bounds_local(1, 1)
1457 ux = v_hartree%pw_grid%bounds_local(2, 1)
1458 ly = v_hartree%pw_grid%bounds_local(1, 2)
1459 uy = v_hartree%pw_grid%bounds_local(2, 2)
1460 lz = v_hartree%pw_grid%bounds_local(1, 3)
1461 uz = v_hartree%pw_grid%bounds_local(2, 3)
1462
1463 dx = v_hartree%pw_grid%npts(1)/2
1464 dy = v_hartree%pw_grid%npts(2)/2
1465 dz = v_hartree%pw_grid%npts(3)/2
1466
1467 dvol = v_hartree%pw_grid%dvol
1468
1469 DO iz = lz, uz
1470 point_indices(3) = real(iz + dz, kind=dp)
1471 DO iy = ly, uy
1472 point_indices(2) = real(iy + dy, kind=dp)
1473
1474 DO ix = lx, ux
1475 point_indices(1) = real(ix + dx, kind=dp)
1476 point_coord(:) = matmul(v_hartree%pw_grid%dh, point_indices)
1477
1478 vector = point_coord - contact_env(1)%origin_bias
1479 proj = projection_on_direction_vector(vector, dirvector_bias)
1480 IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
1481 ! scattering region
1482 ! proj == 0 we are at the first contact boundary
1483 ! proj == 1 we are at the second contact boundary
1484 IF (proj < 0.0_dp) THEN
1485 proj = 0.0_dp
1486 ELSE IF (proj > 1.0_dp) THEN
1487 proj = 1.0_dp
1488 END IF
1489 pot = v1 + (v2 - v1)*proj
1490 ELSE
1491 pot = 0.0_dp
1492 DO icontact = 1, ncontacts
1493 vector = point_coord - contact_env(icontact)%origin_bias
1494 proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias)
1495
1496 IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
1497 pot = contact_control(icontact)%v_external
1498 EXIT
1499 END IF
1500 END DO
1501 END IF
1502
1503 v_hartree%array(ix, iy, iz) = pot*dvol
1504 END DO
1505 END DO
1506 END DO
1507
1508 CALL timestop(handle)
1509 END SUBROUTINE negf_env_init_v_hartree
1510
1511! **************************************************************************************************
1512!> \brief Detect the axis towards secondary unit cell.
1513!> \param direction_vector direction vector
1514!> \param subsys_contact QuickStep subsystem of the contact force environment
1515!> \param eps_geometry accuracy in mapping atoms between different force environments
1516!> \return direction axis: 0 (undefined), 1 (x), 2(y), 3 (z)
1517!> \par History
1518!> * 08.2017 created [Sergey Chulkov]
1519! **************************************************************************************************
1520 FUNCTION contact_direction_axis(direction_vector, subsys_contact, eps_geometry) RESULT(direction_axis)
1521 REAL(kind=dp), DIMENSION(3), INTENT(in) :: direction_vector
1522 TYPE(qs_subsys_type), POINTER :: subsys_contact
1523 REAL(kind=dp), INTENT(in) :: eps_geometry
1524 INTEGER :: direction_axis
1525
1526 INTEGER :: i, naxes
1527 REAL(kind=dp), DIMENSION(3) :: scaled
1528 TYPE(cell_type), POINTER :: cell
1529
1530 CALL qs_subsys_get(subsys_contact, cell=cell)
1531 CALL real_to_scaled(scaled, direction_vector, cell)
1532
1533 naxes = 0
1534 direction_axis = 0 ! initialize to make GCC<=6 happy
1535
1536 DO i = 1, 3
1537 IF (abs(scaled(i)) > eps_geometry) THEN
1538 IF (scaled(i) > 0.0_dp) THEN
1539 direction_axis = i
1540 ELSE
1541 direction_axis = -i
1542 END IF
1543 naxes = naxes + 1
1544 END IF
1545 END DO
1546
1547 ! direction_vector is not parallel to one of the unit cell's axis
1548 IF (naxes /= 1) direction_axis = 0
1549 END FUNCTION contact_direction_axis
1550
1551! **************************************************************************************************
1552!> \brief Estimate energy of the highest spin-alpha occupied molecular orbital.
1553!> \param homo_energy HOMO energy (initialised on exit)
1554!> \param qs_env QuickStep environment
1555!> \par History
1556!> * 01.2017 created [Sergey Chulkov]
1557! **************************************************************************************************
1558 SUBROUTINE negf_homo_energy_estimate(homo_energy, qs_env)
1559 REAL(kind=dp), INTENT(out) :: homo_energy
1560 TYPE(qs_environment_type), POINTER :: qs_env
1561
1562 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_homo_energy_estimate'
1563 INTEGER, PARAMETER :: gamma_point = 1
1564
1565 INTEGER :: handle, homo, ikpgr, ikpoint, imo, &
1566 ispin, kplocal, nmo, nspins
1567 INTEGER, DIMENSION(2) :: kp_range
1568 LOGICAL :: do_kpoints
1569 REAL(kind=dp) :: my_homo_energy
1570 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1571 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
1572 TYPE(kpoint_type), POINTER :: kpoints
1573 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1574 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
1575 TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
1576
1577 CALL timeset(routinen, handle)
1578 my_homo_energy = 0.0_dp
1579
1580 CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
1581
1582 IF (do_kpoints) THEN
1583 CALL get_kpoint_info(kpoints, kp_env=kp_env, kp_range=kp_range, para_env_kp=para_env_kp)
1584
1585 ! looking for a processor that holds the gamma point
1586 IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point) THEN
1587 kplocal = kp_range(2) - kp_range(1) + 1
1588
1589 DO ikpgr = 1, kplocal
1590 CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
1591
1592 IF (ikpoint == gamma_point) THEN
1593 ! mos_kp(component, spin), where component = 1 (real), or 2 (imaginary)
1594 CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues) ! mu=fermi_level
1595
1596 my_homo_energy = eigenvalues(homo)
1597 EXIT
1598 END IF
1599 END DO
1600 END IF
1601
1602 CALL para_env%sum(my_homo_energy)
1603 ELSE
1604 ! Hamiltonian of the bulk contact region has been computed without k-points.
1605 ! Try to obtain the HOMO energy assuming there is no OT. We probably should abort here
1606 ! as we do need a second replica of the bulk contact unit cell along transport
1607 ! direction anyway which is not available without k-points.
1608
1609 CALL cp_abort(__location__, &
1610 "It is necessary to use k-points along the transport direction "// &
1611 "for all contact FORCE_EVAL-s")
1612 ! It is necessary to use k-points along the transport direction within all contact FORCE_EVAL-s
1613
1614 nspins = SIZE(mos)
1615
1616 spin_loop: DO ispin = 1, nspins
1617 CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
1618
1619 DO imo = nmo, 1, -1
1620 IF (eigenvalues(imo) /= 0.0_dp) EXIT spin_loop
1621 END DO
1622 END DO spin_loop
1623
1624 IF (imo == 0) THEN
1625 cpabort("Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
1626 END IF
1627
1628 my_homo_energy = eigenvalues(homo)
1629 END IF
1630
1631 homo_energy = my_homo_energy
1632 CALL timestop(handle)
1633 END SUBROUTINE negf_homo_energy_estimate
1634
1635! **************************************************************************************************
1636!> \brief List atoms from the contact's primary unit cell.
1637!> \param atomlist_cell0 list of atoms belonging to the contact's primary unit cell
1638!> (allocate and initialised on exit)
1639!> \param atom_map_cell0 atomic map of atoms from 'atomlist_cell0' (allocate and initialised on exit)
1640!> \param atomlist_bulk list of atoms belonging to the bulk contact region
1641!> \param atom_map atomic map of atoms from 'atomlist_bulk'
1642!> \param origin origin of the contact
1643!> \param direction_vector direction vector of the contact
1644!> \param direction_axis axis towards secondary unit cell
1645!> \param subsys_device QuickStep subsystem of the device force environment
1646!> \par History
1647!> * 08.2017 created [Sergey Chulkov]
1648! **************************************************************************************************
1649 SUBROUTINE list_atoms_in_bulk_primary_unit_cell(atomlist_cell0, atom_map_cell0, atomlist_bulk, atom_map, &
1650 origin, direction_vector, direction_axis, subsys_device)
1651 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell0
1652 TYPE(negf_atom_map_type), ALLOCATABLE, &
1653 DIMENSION(:), INTENT(inout) :: atom_map_cell0
1654 INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1655 TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1656 REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1657 INTEGER, INTENT(in) :: direction_axis
1658 TYPE(qs_subsys_type), POINTER :: subsys_device
1659
1660 CHARACTER(LEN=*), PARAMETER :: routinen = 'list_atoms_in_bulk_primary_unit_cell'
1661
1662 INTEGER :: atom_min, dir_axis_min, &
1663 direction_axis_abs, handle, iatom, &
1664 natoms_bulk, natoms_cell0
1665 REAL(kind=dp) :: proj, proj_min
1666 REAL(kind=dp), DIMENSION(3) :: vector
1667 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1668
1669 CALL timeset(routinen, handle)
1670 CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1671
1672 natoms_bulk = SIZE(atomlist_bulk)
1673 cpassert(SIZE(atom_map, 1) == natoms_bulk)
1674 direction_axis_abs = abs(direction_axis)
1675
1676 ! looking for the nearest atom from the scattering region
1677 proj_min = 1.0_dp
1678 atom_min = 1
1679 DO iatom = 1, natoms_bulk
1680 vector = particle_set(atomlist_bulk(iatom))%r - origin
1681 proj = projection_on_direction_vector(vector, direction_vector)
1682
1683 IF (proj < proj_min) THEN
1684 proj_min = proj
1685 atom_min = iatom
1686 END IF
1687 END DO
1688
1689 dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1690
1691 natoms_cell0 = 0
1692 DO iatom = 1, natoms_bulk
1693 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) &
1694 natoms_cell0 = natoms_cell0 + 1
1695 END DO
1696
1697 ALLOCATE (atomlist_cell0(natoms_cell0))
1698 ALLOCATE (atom_map_cell0(natoms_cell0))
1699
1700 natoms_cell0 = 0
1701 DO iatom = 1, natoms_bulk
1702 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
1703 natoms_cell0 = natoms_cell0 + 1
1704 atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
1705 atom_map_cell0(natoms_cell0) = atom_map(iatom)
1706 END IF
1707 END DO
1708
1709 CALL timestop(handle)
1710 END SUBROUTINE list_atoms_in_bulk_primary_unit_cell
1711
1712! **************************************************************************************************
1713!> \brief List atoms from the contact's secondary unit cell.
1714!> \param atomlist_cell1 list of atoms belonging to the contact's secondary unit cell
1715!> (allocate and initialised on exit)
1716!> \param atom_map_cell1 atomic map of atoms from 'atomlist_cell1'
1717!> (allocate and initialised on exit)
1718!> \param atomlist_bulk list of atoms belonging to the bulk contact region
1719!> \param atom_map atomic map of atoms from 'atomlist_bulk'
1720!> \param origin origin of the contact
1721!> \param direction_vector direction vector of the contact
1722!> \param direction_axis axis towards the secondary unit cell
1723!> \param subsys_device QuickStep subsystem of the device force environment
1724!> \par History
1725!> * 11.2017 created [Sergey Chulkov]
1726!> \note Cloned from list_atoms_in_bulk_primary_unit_cell. Will be removed once we can managed to
1727!> maintain consistency between real-space matrices from different force_eval sections.
1728! **************************************************************************************************
1729 SUBROUTINE list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1, atom_map_cell1, atomlist_bulk, atom_map, &
1730 origin, direction_vector, direction_axis, subsys_device)
1731 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell1
1732 TYPE(negf_atom_map_type), ALLOCATABLE, &
1733 DIMENSION(:), INTENT(inout) :: atom_map_cell1
1734 INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1735 TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1736 REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1737 INTEGER, INTENT(in) :: direction_axis
1738 TYPE(qs_subsys_type), POINTER :: subsys_device
1739
1740 CHARACTER(LEN=*), PARAMETER :: routinen = 'list_atoms_in_bulk_secondary_unit_cell'
1741
1742 INTEGER :: atom_min, dir_axis_min, &
1743 direction_axis_abs, handle, iatom, &
1744 natoms_bulk, natoms_cell1, offset
1745 REAL(kind=dp) :: proj, proj_min
1746 REAL(kind=dp), DIMENSION(3) :: vector
1747 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1748
1749 CALL timeset(routinen, handle)
1750 CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1751
1752 natoms_bulk = SIZE(atomlist_bulk)
1753 cpassert(SIZE(atom_map, 1) == natoms_bulk)
1754 direction_axis_abs = abs(direction_axis)
1755 offset = sign(1, direction_axis)
1756
1757 ! looking for the nearest atom from the scattering region
1758 proj_min = 1.0_dp
1759 atom_min = 1
1760 DO iatom = 1, natoms_bulk
1761 vector = particle_set(atomlist_bulk(iatom))%r - origin
1762 proj = projection_on_direction_vector(vector, direction_vector)
1763
1764 IF (proj < proj_min) THEN
1765 proj_min = proj
1766 atom_min = iatom
1767 END IF
1768 END DO
1769
1770 dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1771
1772 natoms_cell1 = 0
1773 DO iatom = 1, natoms_bulk
1774 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) &
1775 natoms_cell1 = natoms_cell1 + 1
1776 END DO
1777
1778 ALLOCATE (atomlist_cell1(natoms_cell1))
1779 ALLOCATE (atom_map_cell1(natoms_cell1))
1780
1781 natoms_cell1 = 0
1782 DO iatom = 1, natoms_bulk
1783 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
1784 natoms_cell1 = natoms_cell1 + 1
1785 atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
1786 atom_map_cell1(natoms_cell1) = atom_map(iatom)
1787 atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
1788 END IF
1789 END DO
1790
1791 CALL timestop(handle)
1792 END SUBROUTINE list_atoms_in_bulk_secondary_unit_cell
1793
1794! **************************************************************************************************
1795!> \brief Release a NEGF environment variable.
1796!> \param negf_env NEGF environment to release
1797!> \par History
1798!> * 01.2017 created [Sergey Chulkov]
1799! **************************************************************************************************
1800 SUBROUTINE negf_env_release(negf_env)
1801 TYPE(negf_env_type), INTENT(inout) :: negf_env
1802
1803 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_release'
1804
1805 INTEGER :: handle, icontact
1806
1807 CALL timeset(routinen, handle)
1808
1809 IF (ALLOCATED(negf_env%contacts)) THEN
1810 DO icontact = SIZE(negf_env%contacts), 1, -1
1811 CALL negf_env_contact_release(negf_env%contacts(icontact))
1812 END DO
1813
1814 DEALLOCATE (negf_env%contacts)
1815 END IF
1816
1817 ! h_s
1818 CALL cp_fm_release(negf_env%h_s)
1819
1820 ! h_sc
1821 CALL cp_fm_release(negf_env%h_sc)
1822
1823 ! s_s
1824 IF (ASSOCIATED(negf_env%s_s)) THEN
1825 CALL cp_fm_release(negf_env%s_s)
1826 DEALLOCATE (negf_env%s_s)
1827 NULLIFY (negf_env%s_s)
1828 END IF
1829
1830 ! s_sc
1831 CALL cp_fm_release(negf_env%s_sc)
1832
1833 ! v_hartree_s
1834 IF (ASSOCIATED(negf_env%v_hartree_s)) THEN
1835 CALL cp_fm_release(negf_env%v_hartree_s)
1836 DEALLOCATE (negf_env%v_hartree_s)
1837 NULLIFY (negf_env%v_hartree_s)
1838 END IF
1839
1840 ! density mixing
1841 IF (ASSOCIATED(negf_env%mixing_storage)) THEN
1842 CALL mixing_storage_release(negf_env%mixing_storage)
1843 DEALLOCATE (negf_env%mixing_storage)
1844 END IF
1845
1846 CALL timestop(handle)
1847 END SUBROUTINE negf_env_release
1848
1849! **************************************************************************************************
1850!> \brief Release a NEGF contact environment variable.
1851!> \param contact_env NEGF contact environment to release
1852! **************************************************************************************************
1853 SUBROUTINE negf_env_contact_release(contact_env)
1854 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
1855
1856 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_contact_release'
1857
1858 INTEGER :: handle
1859
1860 CALL timeset(routinen, handle)
1861
1862 ! h_00
1863 CALL cp_fm_release(contact_env%h_00)
1864
1865 ! h_01
1866 CALL cp_fm_release(contact_env%h_01)
1867
1868 ! rho_00
1869 CALL cp_fm_release(contact_env%rho_00)
1870
1871 ! rho_01
1872 CALL cp_fm_release(contact_env%rho_01)
1873
1874 ! s_00
1875 IF (ASSOCIATED(contact_env%s_00)) THEN
1876 CALL cp_fm_release(contact_env%s_00)
1877 DEALLOCATE (contact_env%s_00)
1878 NULLIFY (contact_env%s_00)
1879 END IF
1880
1881 ! s_01
1882 IF (ASSOCIATED(contact_env%s_01)) THEN
1883 CALL cp_fm_release(contact_env%s_01)
1884 DEALLOCATE (contact_env%s_01)
1885 NULLIFY (contact_env%s_01)
1886 END IF
1887
1888 IF (ALLOCATED(contact_env%atomlist_cell0)) DEALLOCATE (contact_env%atomlist_cell0)
1889 IF (ALLOCATED(contact_env%atomlist_cell1)) DEALLOCATE (contact_env%atomlist_cell1)
1890 IF (ALLOCATED(contact_env%atom_map_cell0)) DEALLOCATE (contact_env%atom_map_cell0)
1891
1892 CALL timestop(handle)
1893 END SUBROUTINE negf_env_contact_release
1894
1895END MODULE negf_env_types
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:491
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:311
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:122
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:789
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:972
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:342
force_env_types::use_qs_force
integer, parameter, public use_qs_force
Definition force_env_types.F:83
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_env
subroutine, public get_kpoint_env(kpoint_env, nkpoint, wkp, xkp, is_local, mos)
Get information from a single kpoint environment.
Definition kpoint_types.F:1010
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:501
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
negf_atom_map
Map atoms between various force environments.
Definition negf_atom_map.F:13
negf_atom_map::negf_map_atomic_indices
subroutine, public negf_map_atomic_indices(atom_map, atom_list, subsys_device, subsys_contact, eps_geometry)
Map atoms in the cell 'subsys_device' listed in 'atom_list' with the corresponding atoms in the cell ...
Definition negf_atom_map.F:80
negf_control_types
Input control types for NEGF based quantum transport calculations.
Definition negf_control_types.F:12
negf_env_types
Environment for NEGF based quantum transport calculations.
Definition negf_env_types.F:11
negf_env_types::negf_env_release
subroutine, public negf_env_release(negf_env)
Release a NEGF environment variable.
Definition negf_env_types.F:1801
negf_env_types::negf_env_create
subroutine, public negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
Definition negf_env_types.F:178
negf_io
Routines for reading and writing NEGF restart files.
Definition negf_io.F:12
negf_io::negf_restart_file_name
subroutine, public negf_restart_file_name(filename, exist, negf_section, logger, icontact, ispin, h00, h01, s00, s01, h, s, hc, sc, h_scf)
Checks if the restart file exists and returns the filename.
Definition negf_io.F:61
negf_io::negf_print_matrix_to_file
subroutine, public negf_print_matrix_to_file(filename, matrix)
Prints full matrix to a file.
Definition negf_io.F:279
negf_io::negf_read_matrix_from_file
subroutine, public negf_read_matrix_from_file(filename, matrix)
Reads full matrix from a file.
Definition negf_io.F:309
negf_matrix_utils
Helper routines to manipulate with matrices.
Definition negf_matrix_utils.F:11
negf_matrix_utils::number_of_atomic_orbitals
integer function, public number_of_atomic_orbitals(subsys, atom_list)
Compute the number of atomic orbitals of the given set of atoms.
Definition negf_matrix_utils.F:66
negf_matrix_utils::invert_cell_to_index
subroutine, public invert_cell_to_index(cell_to_index, nimages, index_to_cell)
Invert cell_to_index mapping between unit cells and DBCSR matrix images.
Definition negf_matrix_utils.F:636
negf_matrix_utils::negf_copy_contact_matrix
subroutine, public negf_copy_contact_matrix(fm_cell0, fm_cell1, direction_axis, matrix_kp, atom_list0, atom_list1, subsys, mpi_comm_global, kpoints)
Driver routine to extract diagonal and off-diagonal blocks from a symmetric DBCSR matrix.
Definition negf_matrix_utils.F:305
negf_matrix_utils::negf_copy_sym_dbcsr_to_fm_submat
subroutine, public negf_copy_sym_dbcsr_to_fm_submat(matrix, fm, atomlist_row, atomlist_col, subsys, mpi_comm_global, do_upper_diag, do_lower)
Extract part of the DBCSR matrix based on selected atoms and copy it into a dense matrix.
Definition negf_matrix_utils.F:192
negf_subgroup_types
Environment for NEGF based quantum transport calculations.
Definition negf_subgroup_types.F:13
negf_vectors
Routines to deal with vectors in 3-D real space.
Definition negf_vectors.F:12
negf_vectors::projection_on_direction_vector
pure real(kind=dp) function, public projection_on_direction_vector(vector, vector0)
project the 'vector' onto the direction 'vector0'. Both vectors should have the same origin.
Definition negf_vectors.F:134
negf_vectors::contact_direction_vector
subroutine, public contact_direction_vector(origin, direction_vector, origin_bias, direction_vector_bias, atomlist_screening, atomlist_bulk, subsys)
compute direction vector of the given contact
Definition negf_vectors.F:44
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_density_mixing_types
module that contains the definitions of the scf types
Definition qs_density_mixing_types.F:14
qs_density_mixing_types::mixing_storage_release
subroutine, public mixing_storage_release(mixing_store)
releases a mixing_storage
Definition qs_density_mixing_types.F:230
qs_density_mixing_types::mixing_storage_create
subroutine, public mixing_storage_create(mixing_store, mixing_section, mixing_method, ecut)
creates a mixing_storage
Definition qs_density_mixing_types.F:119
qs_energy_init
Utility subroutine for qs energy calculation.
Definition qs_energy_init.F:14
qs_energy_init::qs_energies_init
subroutine, public qs_energies_init(qs_env, calc_forces)
Refactoring of qs_energies_scf. Driver routine for the initial setup and calculations for a qs energy...
Definition qs_energy_init.F:92
qs_energy
Perform a QUICKSTEP wavefunction optimization (single point)
Definition qs_energy.F:14
qs_energy::qs_energies
subroutine, public qs_energies(qs_env, consistent_energies, calc_forces)
Driver routine for QUICKSTEP single point wavefunction optimization.
Definition qs_energy.F:70
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:530
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:399
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_post_tb
Does all kind of post scf calculations for DFTB.
Definition qs_scf_post_tb.F:15
qs_scf_post_tb::rebuild_pw_env
subroutine, public rebuild_pw_env(qs_env)
...
Definition qs_scf_post_tb.F:1037
qs_subsys_types
types that represent a quickstep subsys
Definition qs_subsys_types.F:12
qs_subsys_types::qs_subsys_get
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Definition qs_subsys_types.F:134
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:654
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
force_env_types::force_env_p_type
allows for the creation of an array of force_env
Definition force_env_types.F:174
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_env_p_type
Definition kpoint_types.F:96
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:166
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
negf_atom_map::negf_atom_map_type
Structure that maps the given atom in the sourse FORCE_EVAL section with another atom from the target...
Definition negf_atom_map.F:38
negf_control_types::negf_control_contact_type
Input parameters related to a single contact.
Definition negf_control_types.F:47
negf_control_types::negf_control_type
Input parameters related to the NEGF run.
Definition negf_control_types.F:79
negf_env_types::negf_env_contact_type
Contact-specific NEGF environment.
Definition negf_env_types.F:98
negf_env_types::negf_env_type
NEGF environment.
Definition negf_env_types.F:132
negf_subgroup_types::negf_subgroup_env_type
Parallel (sub)group environment.
Definition negf_subgroup_types.F:36
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_density_mixing_types::mixing_storage_type
Definition qs_density_mixing_types.F:58
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_subsys_types::qs_subsys_type
Definition qs_subsys_types.F:56