(git:8917686)
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) THEN
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 END IF
297 CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
298 CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
299 CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
300 qs_env_contact=qs_env_contact)
301 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
302 END IF
303 END IF
304 END DO
305
306 ! *** obtain relevant Kohn-Sham matrix blocks for each contact with no separate FORCE_ENV ***
307 is_dft_entire = .false.
308 DO icontact = 1, ncontacts
309 IF (negf_control%contacts(icontact)%force_env_index <= 0) THEN
310 IF (negf_control%contacts(icontact)%read_write_HS) THEN
311 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
312 contact_control=negf_control%contacts(icontact), &
313 sub_env=sub_env, qs_env=qs_env, &
314 eps_geometry=negf_control%eps_geometry)
315 CALL negf_env_contact_read_write_hs(icontact, force_env, para_env, negf_env, sub_env, negf_control, negf_section, &
316 log_unit, is_separate=.false., is_dft_entire=is_dft_entire)
317 ELSE
318 IF (log_unit > 0) THEN
319 WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') "NEGF| Construct the Kohn-Sham matrix for the contact", icontact, &
320 " from the entire system bulk DFT calculation"
321 END IF
322 IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
323 is_dft_entire = .true.
324 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
325 contact_control=negf_control%contacts(icontact), &
326 sub_env=sub_env, qs_env=qs_env, &
327 eps_geometry=negf_control%eps_geometry)
328 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
329 END IF
330 END IF
331 END DO
332
333 ! stage 3: obtain an initial KS-matrix for the scattering region
334 IF (log_unit > 0) THEN
335 WRITE (log_unit, '(/,T2,A,T70)') "NEGF| Construct the Kohn-Sham matrix for the scattering region"
336 END IF
337 IF (negf_control%read_write_HS) THEN
338 CALL negf_env_scatt_read_write_hs(force_env, para_env, negf_env, sub_env, negf_control, negf_section, log_unit, &
339 is_dft_entire=is_dft_entire)
340 ELSE
341 IF (.NOT. is_dft_entire) THEN
342 CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
343 is_dft_entire = .true.
344 END IF
345 ! extract device-related matrix blocks
346 CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
347 END IF
348 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
349
350 negf_control%is_dft_entire = is_dft_entire
351
352 ! electron density mixing;
353 ! the input section below should be consistent with the subroutine create_negf_section()
354 NULLIFY (negf_env%mixing_storage)
355 CALL section_vals_val_get(negf_mixing_section, "METHOD", i_val=negf_env%mixing_method)
356
357 CALL get_qs_env(qs_env, dft_control=dft_control)
358 ALLOCATE (negf_env%mixing_storage)
359 CALL mixing_storage_create(negf_env%mixing_storage, negf_mixing_section, &
360 negf_env%mixing_method, dft_control%qs_control%cutoff)
361
362 CALL timestop(handle)
363 END SUBROUTINE negf_env_create
364
365! **************************************************************************************************
366!> \brief Establish mapping between the primary and the contact force environments
367!> \param contact_env NEGF environment for the given contact (modified on exit)
368!> \param contact_control NEGF control
369!> \param atom_map atomic map
370!> \param eps_geometry accuracy in mapping atoms between different force environments
371!> \param subsys_device QuickStep subsystem of the device force environment
372!> \param subsys_contact QuickStep subsystem of the contact force environment
373!> \author Sergey Chulkov
374! **************************************************************************************************
375 SUBROUTINE negf_env_contact_init_maps(contact_env, contact_control, atom_map, &
376 eps_geometry, subsys_device, subsys_contact)
377 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
378 TYPE(negf_control_contact_type), INTENT(in) :: contact_control
379 TYPE(negf_atom_map_type), ALLOCATABLE, &
380 DIMENSION(:), INTENT(inout) :: atom_map
381 REAL(kind=dp), INTENT(in) :: eps_geometry
382 TYPE(qs_subsys_type), POINTER :: subsys_device, subsys_contact
383
384 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_maps'
385
386 INTEGER :: handle, natoms
387
388 CALL timeset(routinen, handle)
389
390 CALL contact_direction_vector(contact_env%origin, &
391 contact_env%direction_vector, &
392 contact_env%origin_bias, &
393 contact_env%direction_vector_bias, &
394 contact_control%atomlist_screening, &
395 contact_control%atomlist_bulk, &
396 subsys_device)
397
398 contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
399
400 IF (contact_env%direction_axis /= 0) THEN
401 natoms = SIZE(contact_control%atomlist_bulk)
402 ALLOCATE (atom_map(natoms))
403
404 ! map atom listed in 'contact_control%atomlist_bulk' to the corresponding atom/cell replica from the contact force_env
405 CALL negf_map_atomic_indices(atom_map=atom_map, &
406 atom_list=contact_control%atomlist_bulk, &
407 subsys_device=subsys_device, &
408 subsys_contact=subsys_contact, &
409 eps_geometry=eps_geometry)
410
411 ! list atoms from 'contact_control%atomlist_bulk' which belong to
412 ! the primary unit cell of the bulk region for the given contact
413 CALL list_atoms_in_bulk_primary_unit_cell(atomlist_cell0=contact_env%atomlist_cell0, &
414 atom_map_cell0=contact_env%atom_map_cell0, &
415 atomlist_bulk=contact_control%atomlist_bulk, &
416 atom_map=atom_map, &
417 origin=contact_env%origin, &
418 direction_vector=contact_env%direction_vector, &
419 direction_axis=contact_env%direction_axis, &
420 subsys_device=subsys_device)
421
422 ! secondary unit cell
423 CALL list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1=contact_env%atomlist_cell1, &
424 atom_map_cell1=contact_env%atom_map_cell1, &
425 atomlist_bulk=contact_control%atomlist_bulk, &
426 atom_map=atom_map, &
427 origin=contact_env%origin, &
428 direction_vector=contact_env%direction_vector, &
429 direction_axis=contact_env%direction_axis, &
430 subsys_device=subsys_device)
431 END IF
432
433 CALL timestop(handle)
434 END SUBROUTINE negf_env_contact_init_maps
435
436! **************************************************************************************************
437!> \brief Reading and writing of the electrode Hamiltonian and overlap matrices from/to a file.
438!> \param icontact ...
439!> \param el_force_env ...
440!> \param para_env ...
441!> \param negf_env ...
442!> \param sub_env ...
443!> \param negf_control ...
444!> \param negf_section ...
445!> \param log_unit ...
446!> \param is_separate ...
447!> \param is_dft_entire ...
448!> \par History
449!> * 12.2025 created [Dmitry Ryndyk]
450! **************************************************************************************************
451 SUBROUTINE negf_env_contact_read_write_hs(icontact, el_force_env, para_env, negf_env, sub_env, negf_control, &
452 negf_section, log_unit, is_separate, is_dft_entire)
453 INTEGER :: icontact
454 TYPE(force_env_type), POINTER :: el_force_env
455 TYPE(mp_para_env_type), POINTER :: para_env
456 TYPE(negf_env_type), INTENT(inout) :: negf_env
457 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
458 TYPE(negf_control_type), POINTER :: negf_control
459 TYPE(section_vals_type), POINTER :: negf_section
460 INTEGER, INTENT(in) :: log_unit
461 LOGICAL, INTENT(in) :: is_separate
462 LOGICAL, INTENT(inout), OPTIONAL :: is_dft_entire
463
464 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_contact_read_write_hs'
465
466 CHARACTER(len=default_path_length) :: filename_h00_1, filename_h00_2, &
467 filename_h01_1, filename_h01_2, &
468 filename_s00, filename_s01
469 INTEGER :: handle, ispin, ncol, nrow, nspins, &
470 print_unit
471 LOGICAL :: exist, exist_all
472 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: target_m
473 TYPE(cp_fm_struct_type), POINTER :: fm_struct
474 TYPE(cp_logger_type), POINTER :: logger
475 TYPE(dft_control_type), POINTER :: dft_control
476 TYPE(qs_environment_type), POINTER :: qs_env_contact
477 TYPE(qs_subsys_type), POINTER :: subsys
478
479 CALL timeset(routinen, handle)
480 logger => cp_get_default_logger()
481
482 CALL force_env_get(el_force_env, qs_env=qs_env_contact)
483 CALL get_qs_env(qs_env_contact, dft_control=dft_control, subsys=subsys)
484 nspins = dft_control%nspins
485
486 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11)') &
487 "NEGF| Construct the Kohn-Sham matrix for the contact", icontact
488
489 ! Check that the files exist.
490 ! ispin=0 is used to show nspins=1
491 exist_all = .true.
492 IF (para_env%is_source()) THEN
493 CALL negf_restart_file_name(filename_s00, exist, negf_section, logger, icontact, s00=.true.)
494 IF (.NOT. exist) THEN
495 CALL cp_warn(__location__, &
496 "User requested to read the overlap matrix from the file named: "// &
497 trim(filename_s00)//". This file does not exist. The file will be created.")
498 exist_all = .false.
499 END IF
500 CALL negf_restart_file_name(filename_s01, exist, negf_section, logger, icontact, s01=.true.)
501 IF (.NOT. exist) THEN
502 CALL cp_warn(__location__, &
503 "User requested to read the overlap matrix from the file named: "// &
504 trim(filename_s01)//". This file does not exist. The file will be created.")
505 exist_all = .false.
506 END IF
507 IF (nspins == 1) THEN
508 CALL negf_restart_file_name(filename_h00_1, exist, negf_section, logger, icontact, ispin=0, h00=.true.)
509 IF (.NOT. exist) THEN
510 CALL cp_warn(__location__, &
511 "User requested to read the Hamiltonian matrix from the file named: "// &
512 trim(filename_h00_1)//". This file does not exist. The file will be created.")
513 exist_all = .false.
514 END IF
515 CALL negf_restart_file_name(filename_h01_1, exist, negf_section, logger, icontact, ispin=0, h01=.true.)
516 IF (.NOT. exist) THEN
517 CALL cp_warn(__location__, &
518 "User requested to read the Hamiltonian matrix from the file named: "// &
519 trim(filename_h01_1)//". This file does not exist. The file will be created.")
520 exist_all = .false.
521 END IF
522 END IF
523 IF (nspins == 2) THEN
524 CALL negf_restart_file_name(filename_h00_1, exist, negf_section, logger, icontact, ispin=1, h00=.true.)
525 IF (.NOT. exist) THEN
526 CALL cp_warn(__location__, &
527 "User requested to read the Hamiltonian matrix from the file named: "// &
528 trim(filename_h00_1)//". This file does not exist. The file will be created.")
529 exist_all = .false.
530 END IF
531 CALL negf_restart_file_name(filename_h01_1, exist, negf_section, logger, icontact, ispin=1, h01=.true.)
532 IF (.NOT. exist) THEN
533 CALL cp_warn(__location__, &
534 "User requested to read tthe Hamiltonian matrix from the file named: "// &
535 trim(filename_h01_1)//". This file does not exist. The file will be created.")
536 exist_all = .false.
537 END IF
538 CALL negf_restart_file_name(filename_h00_2, exist, negf_section, logger, icontact, ispin=2, h00=.true.)
539 IF (.NOT. exist) THEN
540 CALL cp_warn(__location__, &
541 "User requested to read the Hamiltonian matrix from the file named: "// &
542 trim(filename_h00_2)//". This file does not exist. The file will be created.")
543 exist_all = .false.
544 END IF
545 CALL negf_restart_file_name(filename_h01_2, exist, negf_section, logger, icontact, ispin=2, h01=.true.)
546 IF (.NOT. exist) THEN
547 CALL cp_warn(__location__, &
548 "User requested to read the Hamiltonian matrix from the file named: "// &
549 trim(filename_h01_2)//". This file does not exist. The file will be created.")
550 exist_all = .false.
551 END IF
552 END IF
553 END IF
554 CALL para_env%bcast(exist_all)
555
556 IF (exist_all) THEN
557
558 negf_control%contacts(icontact)%is_restart = .true.
559 IF (log_unit > 0) THEN
560 WRITE (log_unit, '(/,T2,A)') "User requested to read the Hamiltonian and overlap matrices from files."
561 WRITE (log_unit, '(T2,A)') "All restart files exist."
562 END IF
563
564 ! ++ create matrices: s_00, s_01, h_00, h_01
565 IF (para_env%is_source()) THEN
566 CALL open_file(file_name=filename_s00, file_status="OLD", &
567 file_form="FORMATTED", file_action="READ", &
568 file_position="REWIND", unit_number=print_unit)
569 READ (print_unit, *) nrow, ncol
570 CALL close_file(print_unit)
571 END IF
572 CALL para_env%bcast(nrow)
573 CALL para_env%bcast(ncol)
574 NULLIFY (fm_struct)
575 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, context=sub_env%blacs_env)
576 ALLOCATE (negf_env%contacts(icontact)%s_00, negf_env%contacts(icontact)%s_01)
577 CALL cp_fm_create(negf_env%contacts(icontact)%s_00, fm_struct)
578 CALL cp_fm_create(negf_env%contacts(icontact)%s_01, fm_struct)
579 ALLOCATE (negf_env%contacts(icontact)%h_00(nspins), negf_env%contacts(icontact)%h_01(nspins))
580 DO ispin = 1, nspins
581 CALL cp_fm_create(negf_env%contacts(icontact)%h_00(ispin), fm_struct)
582 CALL cp_fm_create(negf_env%contacts(icontact)%h_01(ispin), fm_struct)
583 END DO
584 CALL cp_fm_struct_release(fm_struct)
585
586 ALLOCATE (target_m(nrow, ncol))
587 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s00, target_m)
588 CALL para_env%bcast(target_m)
589 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%s_00, target_m)
590 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_00 is read from "//trim(filename_s00)
591 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s01, target_m)
592 CALL para_env%bcast(target_m)
593 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%s_01, target_m)
594 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_01 is read from "//trim(filename_s01)
595 IF (nspins == 1) THEN
596 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_1, target_m)
597 CALL para_env%bcast(target_m)
598 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
599 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//trim(filename_h00_1)
600 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_1, target_m)
601 CALL para_env%bcast(target_m)
602 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
603 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//trim(filename_h01_1)
604 END IF
605 IF (nspins == 2) THEN
606 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_1, target_m)
607 CALL para_env%bcast(target_m)
608 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
609 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//trim(filename_h00_1)//" for spin 1"
610 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_1, target_m)
611 CALL para_env%bcast(target_m)
612 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
613 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//trim(filename_h01_1)//" for spin 1"
614 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_2, target_m)
615 CALL para_env%bcast(target_m)
616 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(2), target_m)
617 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//trim(filename_h00_2)//" for spin 2"
618 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_2, target_m)
619 CALL para_env%bcast(target_m)
620 CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(2), target_m)
621 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//trim(filename_h01_2)//" for spin 2"
622 END IF
623 DEALLOCATE (target_m)
624
625 ELSE
626
627 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') &
628 "Some restart files do not exist. ALL restart files will be recalculated!"
629
630 IF (is_separate) THEN
631 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') &
632 "Construct the Kohn-Sham matrix from from the separate bulk DFT calculation"
633 CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
634 CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
635 qs_env_contact=qs_env_contact)
636 ELSE
637 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') &
638 "Construct the Kohn-Sham matrix from the entire system bulk DFT calculation"
639 negf_control%contacts(icontact)%read_write_HS = .false.
640 IF (.NOT. is_dft_entire) CALL qs_energies(qs_env_contact, consistent_energies=.false., calc_forces=.false.)
641 CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
642 contact_control=negf_control%contacts(icontact), &
643 sub_env=sub_env, qs_env=qs_env_contact, &
644 eps_geometry=negf_control%eps_geometry)
645 negf_control%contacts(icontact)%read_write_HS = .true.
646 is_dft_entire = .true.
647 END IF
648
649 CALL cp_fm_get_info(negf_env%contacts(icontact)%s_00, nrow_global=nrow)
650 ALLOCATE (target_m(nrow, nrow))
651 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%s_00, target_m)
652 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s00, target_m)
653 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "S_00 is saved to "//trim(filename_s00)
654 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%s_01, target_m)
655 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s01, target_m)
656 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_01 is saved to "//trim(filename_s01)
657 IF (nspins == 1) THEN
658 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
659 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_1, target_m)
660 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//trim(filename_h00_1)
661 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
662 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_1, target_m)
663 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//trim(filename_h01_1)
664 END IF
665 IF (nspins == 2) THEN
666 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
667 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_1, target_m)
668 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//trim(filename_h00_1)//" for spin 1"
669 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
670 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_1, target_m)
671 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//trim(filename_h01_1)//" for spin 1"
672 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(2), target_m)
673 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_2, target_m)
674 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//trim(filename_h00_2)//" for spin 2"
675 CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(2), target_m)
676 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_2, target_m)
677 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//trim(filename_h01_2)//" for spin 2"
678 END IF
679 DEALLOCATE (target_m)
680
681 negf_control%write_common_restart_file = .true.
682
683 END IF
684
685 IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
686
687 CALL timestop(handle)
688 END SUBROUTINE negf_env_contact_read_write_hs
689
690! **************************************************************************************************
691!> \brief Extract relevant matrix blocks for the given contact.
692!> \param contact_env NEGF environment for the contact (modified on exit)
693!> \param sub_env NEGF parallel (sub)group environment
694!> \param qs_env_contact QuickStep environment for the contact force environment
695!> \par History
696!> * 10.2017 created [Sergey Chulkov]
697!> * 10.2025 The subroutine is essentially modified. New functionality of negf_copy_contact_matrix.
698!> [Dmitry Ryndyk]
699! **************************************************************************************************
700 SUBROUTINE negf_env_contact_init_matrices(contact_env, sub_env, qs_env_contact)
701 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
702 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
703 TYPE(qs_environment_type), POINTER :: qs_env_contact
704
705 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_matrices'
706
707 INTEGER :: handle, iatom, ispin, nao, natoms, &
708 nimages, nspins
709 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_list0, atom_list1
710 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
711 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
712 LOGICAL :: do_kpoints
713 TYPE(cp_fm_struct_type), POINTER :: fm_struct
714 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matkp
715 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
716 TYPE(dft_control_type), POINTER :: dft_control
717 TYPE(kpoint_type), POINTER :: kpoints
718 TYPE(mp_para_env_type), POINTER :: para_env
719 TYPE(qs_rho_type), POINTER :: rho_struct
720 TYPE(qs_subsys_type), POINTER :: subsys
721
722 CALL timeset(routinen, handle)
723
724 CALL get_qs_env(qs_env_contact, &
725 dft_control=dft_control, &
726 do_kpoints=do_kpoints, &
727 kpoints=kpoints, &
728 matrix_ks_kp=matrix_ks_kp, &
729 matrix_s_kp=matrix_s_kp, &
730 para_env=para_env, &
731 rho=rho_struct, &
732 subsys=subsys)
733 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
734
735 CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
736
737 natoms = SIZE(contact_env%atomlist_cell0)
738 ALLOCATE (atom_list0(natoms))
739 DO iatom = 1, natoms
740 atom_list0(iatom) = contact_env%atom_map_cell0(iatom)%iatom
741
742 ! with no k-points there is one-to-one correspondence between the primary unit cell
743 ! of the contact force_env and the first contact unit cell of the device force_env
744 IF (sum(abs(contact_env%atom_map_cell0(iatom)%cell(:))) > 0) THEN
745 cpabort("NEGF K-points are not currently supported")
746 END IF
747 END DO
748
749 cpassert(SIZE(contact_env%atomlist_cell1) == natoms)
750 ALLOCATE (atom_list1(natoms))
751 DO iatom = 1, natoms
752 atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
753 END DO
754
755 nspins = dft_control%nspins
756 nimages = dft_control%nimages
757
758 IF (do_kpoints) THEN
759 CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
760 ELSE
761 ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
762 cell_to_index(0, 0, 0) = 1
763 END IF
764
765 ALLOCATE (index_to_cell(3, nimages))
766 CALL invert_cell_to_index(cell_to_index, nimages, index_to_cell)
767 IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
768
769 NULLIFY (fm_struct)
770 nao = number_of_atomic_orbitals(subsys, atom_list0)
771 CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env)
772
773 ! ++ create matrices: s_00, s_01
774 ALLOCATE (contact_env%s_00, contact_env%s_01)
775 CALL cp_fm_create(contact_env%s_00, fm_struct)
776 CALL cp_fm_create(contact_env%s_01, fm_struct)
777
778 ! ++ create matrices: h_00, h_01, rho_00, rho_01
779 ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
780 ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
781 DO ispin = 1, nspins
782 CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
783 CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
784 CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
785 CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
786 END DO
787
788 CALL cp_fm_struct_release(fm_struct)
789
790 ! extract matrices: s_00, s_01
791 matkp => matrix_s_kp(1, :)
792 CALL negf_copy_contact_matrix(fm_cell0=contact_env%s_00, &
793 fm_cell1=contact_env%s_01, &
794 direction_axis=contact_env%direction_axis, &
795 matrix_kp=matkp, &
796 atom_list0=atom_list0, atom_list1=atom_list1, &
797 subsys=subsys, mpi_comm_global=para_env, &
798 kpoints=kpoints)
799
800 ! extract matrices: h_00, h_01, rho_00, rho_01
801 DO ispin = 1, nspins
802 matkp => matrix_ks_kp(ispin, :)
803 CALL negf_copy_contact_matrix(fm_cell0=contact_env%h_00(ispin), &
804 fm_cell1=contact_env%h_01(ispin), &
805 direction_axis=contact_env%direction_axis, &
806 matrix_kp=matkp, &
807 atom_list0=atom_list0, atom_list1=atom_list1, &
808 subsys=subsys, mpi_comm_global=para_env, &
809 kpoints=kpoints)
810
811 matkp => rho_ao_kp(ispin, :)
812 CALL negf_copy_contact_matrix(fm_cell0=contact_env%rho_00(ispin), &
813 fm_cell1=contact_env%rho_01(ispin), &
814 direction_axis=contact_env%direction_axis, &
815 matrix_kp=matkp, &
816 atom_list0=atom_list0, atom_list1=atom_list1, &
817 subsys=subsys, mpi_comm_global=para_env, &
818 kpoints=kpoints)
819 END DO
820
821 DEALLOCATE (atom_list0, atom_list1)
822
823 CALL timestop(handle)
824 END SUBROUTINE negf_env_contact_init_matrices
825
826! **************************************************************************************************
827!> \brief Extract relevant matrix blocks for the given contact using the device's force environment.
828!> \param contact_env NEGF environment for the contact (modified on exit)
829!> \param contact_control NEGF control for the contact
830!> \param sub_env NEGF parallel (sub)group environment
831!> \param qs_env QuickStep environment for the device force environment
832!> \param eps_geometry accuracy in Cartesian coordinates
833!> \author Sergey Chulkov
834! **************************************************************************************************
835 SUBROUTINE negf_env_contact_init_matrices_gamma(contact_env, contact_control, sub_env, qs_env, eps_geometry)
836 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
837 TYPE(negf_control_contact_type), INTENT(in) :: contact_control
838 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
839 TYPE(qs_environment_type), POINTER :: qs_env
840 REAL(kind=dp), INTENT(in) :: eps_geometry
841
842 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_contact_init_matrices_gamma'
843
844 INTEGER :: handle, iatom, icell, ispin, nao_c, &
845 nspins
846 LOGICAL :: do_kpoints
847 REAL(kind=dp), DIMENSION(2) :: r2_origin_cell
848 REAL(kind=dp), DIMENSION(3) :: direction_vector, origin
849 TYPE(cp_fm_struct_type), POINTER :: fm_struct
850 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
851 TYPE(dft_control_type), POINTER :: dft_control
852 TYPE(mp_para_env_type), POINTER :: para_env
853 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
854 TYPE(qs_rho_type), POINTER :: rho_struct
855 TYPE(qs_subsys_type), POINTER :: subsys
856
857 CALL timeset(routinen, handle)
858
859 CALL get_qs_env(qs_env, &
860 dft_control=dft_control, &
861 do_kpoints=do_kpoints, &
862 matrix_ks_kp=matrix_ks_kp, &
863 matrix_s_kp=matrix_s_kp, &
864 para_env=para_env, &
865 rho=rho_struct, &
866 subsys=subsys)
867 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
868
869 IF (do_kpoints) THEN
870 CALL cp_abort(__location__, &
871 "K-points in device region have not been implemented yet.")
872 END IF
873
874 nspins = dft_control%nspins
875
876 nao_c = number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(1)%vector)
877 IF (number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(2)%vector) /= nao_c) THEN
878 CALL cp_abort(__location__, &
879 "Primary and secondary bulk contact cells should be identical "// &
880 "in terms of the number of atoms of each kind, and their basis sets. "// &
881 "No single atom, however, can be shared between these two cells.")
882 END IF
883
884 contact_env%homo_energy = 0.0_dp
885
886 CALL contact_direction_vector(contact_env%origin, &
887 contact_env%direction_vector, &
888 contact_env%origin_bias, &
889 contact_env%direction_vector_bias, &
890 contact_control%atomlist_screening, &
891 contact_control%atomlist_bulk, &
892 subsys)
893
894 contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys, eps_geometry)
895
896 ! choose the primary and secondary contact unit cells
897 CALL qs_subsys_get(subsys, particle_set=particle_set)
898
899 origin = particle_set(contact_control%atomlist_screening(1))%r
900 DO iatom = 2, SIZE(contact_control%atomlist_screening)
901 origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
902 END DO
903 origin = origin/real(SIZE(contact_control%atomlist_screening), kind=dp)
904
905 DO icell = 1, 2
906 direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
907 DO iatom = 2, SIZE(contact_control%atomlist_cell(icell)%vector)
908 direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
909 END DO
910 direction_vector = direction_vector/real(SIZE(contact_control%atomlist_cell(icell)%vector), kind=dp)
911 direction_vector = direction_vector - origin
912 r2_origin_cell(icell) = dot_product(direction_vector, direction_vector)
913 END DO
914
915 IF (abs(r2_origin_cell(1) - r2_origin_cell(2)) < (eps_geometry*eps_geometry)) THEN
916 ! primary and secondary bulk unit cells should not overlap;
917 ! currently we check that they are different by at least one atom that is, indeed, not sufficient.
918 CALL cp_abort(__location__, &
919 "Primary and secondary bulk contact cells should not overlap ")
920 ELSE IF (r2_origin_cell(1) < r2_origin_cell(2)) THEN
921 IF (.NOT. ALLOCATED(contact_env%atomlist_cell0)) THEN
922 ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(1)%vector)))
923 END IF
924 contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
925 IF (.NOT. ALLOCATED(contact_env%atomlist_cell1)) THEN
926 ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(2)%vector)))
927 END IF
928 contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
929 ELSE
930 IF (.NOT. ALLOCATED(contact_env%atomlist_cell0)) THEN
931 ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(2)%vector)))
932 END IF
933 contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
934 IF (.NOT. ALLOCATED(contact_env%atomlist_cell1)) THEN
935 ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(1)%vector)))
936 END IF
937 contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
938 END IF
939 IF (.NOT. contact_control%read_write_HS) THEN
940 NULLIFY (fm_struct)
941 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
942 ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
943 ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
944 DO ispin = 1, nspins
945 CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
946 CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
947 CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
948 CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
949 END DO
950 ALLOCATE (contact_env%s_00, contact_env%s_01)
951 CALL cp_fm_create(contact_env%s_00, fm_struct)
952 CALL cp_fm_create(contact_env%s_01, fm_struct)
953 CALL cp_fm_struct_release(fm_struct)
954
955 DO ispin = 1, nspins
956 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
957 fm=contact_env%h_00(ispin), &
958 atomlist_row=contact_env%atomlist_cell0, &
959 atomlist_col=contact_env%atomlist_cell0, &
960 subsys=subsys, mpi_comm_global=para_env, &
961 do_upper_diag=.true., do_lower=.true.)
962 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
963 fm=contact_env%h_01(ispin), &
964 atomlist_row=contact_env%atomlist_cell0, &
965 atomlist_col=contact_env%atomlist_cell1, &
966 subsys=subsys, mpi_comm_global=para_env, &
967 do_upper_diag=.true., do_lower=.true.)
968
969 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
970 fm=contact_env%rho_00(ispin), &
971 atomlist_row=contact_env%atomlist_cell0, &
972 atomlist_col=contact_env%atomlist_cell0, &
973 subsys=subsys, mpi_comm_global=para_env, &
974 do_upper_diag=.true., do_lower=.true.)
975 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
976 fm=contact_env%rho_01(ispin), &
977 atomlist_row=contact_env%atomlist_cell0, &
978 atomlist_col=contact_env%atomlist_cell1, &
979 subsys=subsys, mpi_comm_global=para_env, &
980 do_upper_diag=.true., do_lower=.true.)
981 END DO
982
983 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
984 fm=contact_env%s_00, &
985 atomlist_row=contact_env%atomlist_cell0, &
986 atomlist_col=contact_env%atomlist_cell0, &
987 subsys=subsys, mpi_comm_global=para_env, &
988 do_upper_diag=.true., do_lower=.true.)
989 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
990 fm=contact_env%s_01, &
991 atomlist_row=contact_env%atomlist_cell0, &
992 atomlist_col=contact_env%atomlist_cell1, &
993 subsys=subsys, mpi_comm_global=para_env, &
994 do_upper_diag=.true., do_lower=.true.)
995 END IF
996 CALL timestop(handle)
997 END SUBROUTINE negf_env_contact_init_matrices_gamma
998
999! **************************************************************************************************
1000!> \brief Reading and writing of the electrode Hamiltonian and overlap matrices from/to a file.
1001!> \param force_env ...
1002!> \param para_env ...
1003!> \param negf_env ...
1004!> \param sub_env ...
1005!> \param negf_control ...
1006!> \param negf_section ...
1007!> \param log_unit ...
1008!> \param is_dft_entire ...
1009!> \par History
1010!> * 01.2026 created [Dmitry Ryndyk]
1011! **************************************************************************************************
1012 SUBROUTINE negf_env_scatt_read_write_hs(force_env, para_env, negf_env, sub_env, negf_control, negf_section, &
1013 log_unit, is_dft_entire)
1014 TYPE(force_env_type), POINTER :: force_env
1015 TYPE(mp_para_env_type), POINTER :: para_env
1016 TYPE(negf_env_type), INTENT(inout) :: negf_env
1017 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1018 TYPE(negf_control_type), POINTER :: negf_control
1019 TYPE(section_vals_type), POINTER :: negf_section
1020 INTEGER, INTENT(in) :: log_unit
1021 LOGICAL, INTENT(inout), OPTIONAL :: is_dft_entire
1022
1023 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_scatt_read_write_hs'
1024
1025 CHARACTER(len=default_path_length) :: filename_h_1, filename_h_2, filename_s
1026 CHARACTER(len=default_path_length), ALLOCATABLE, &
1027 DIMENSION(:) :: filename_hc_1, filename_hc_2, filename_sc
1028 INTEGER :: handle, icontact, ispin, ncol_s, &
1029 ncol_sc, ncontacts, nrow_s, nrow_sc, &
1030 nspins, print_unit
1031 LOGICAL :: exist, exist_all
1032 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: target_m
1033 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1034 TYPE(cp_logger_type), POINTER :: logger
1035 TYPE(dft_control_type), POINTER :: dft_control
1036 TYPE(qs_environment_type), POINTER :: qs_env
1037 TYPE(qs_subsys_type), POINTER :: subsys
1038
1039 CALL timeset(routinen, handle)
1040 logger => cp_get_default_logger()
1041
1042 CALL force_env_get(force_env, qs_env=qs_env)
1043 CALL get_qs_env(qs_env, dft_control=dft_control, subsys=subsys)
1044 ncontacts = SIZE(negf_control%contacts)
1045 nspins = dft_control%nspins
1046 ALLOCATE (filename_sc(ncontacts), filename_hc_1(ncontacts), filename_hc_2(ncontacts))
1047
1048 ! Check that the files exist.
1049 ! ispin=0 is used to show nspins=1
1050 exist_all = .true.
1051 IF (para_env%is_source()) THEN
1052 CALL negf_restart_file_name(filename_s, exist, negf_section, logger, s=.true.)
1053 IF (.NOT. exist) THEN
1054 CALL cp_warn(__location__, &
1055 "User requested to read the overlap matrix from the file named: "// &
1056 trim(filename_s)//". This file does not exist. The file will be created.")
1057 exist_all = .false.
1058 END IF
1059 IF (nspins == 1) THEN
1060 CALL negf_restart_file_name(filename_h_1, exist, negf_section, logger, ispin=0, h=.true.)
1061 IF (.NOT. exist) THEN
1062 CALL cp_warn(__location__, &
1063 "User requested to read the Hamiltonian matrix from the file named: "// &
1064 trim(filename_h_1)//". This file does not exist. The file will be created.")
1065 exist_all = .false.
1066 END IF
1067 END IF
1068 IF (nspins == 2) THEN
1069 CALL negf_restart_file_name(filename_h_1, exist, negf_section, logger, ispin=1, 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_1)//". This file does not exist. The file will be created.")
1074 exist_all = .false.
1075 END IF
1076 CALL negf_restart_file_name(filename_h_2, exist, negf_section, logger, ispin=2, h=.true.)
1077 IF (.NOT. exist) THEN
1078 CALL cp_warn(__location__, &
1079 "User requested to read the Hamiltonian matrix from the file named: "// &
1080 trim(filename_h_2)//". This file does not exist. The file will be created.")
1081 exist_all = .false.
1082 END IF
1083 END IF
1084 DO icontact = 1, ncontacts
1085 CALL negf_restart_file_name(filename_sc(icontact), exist, negf_section, logger, icontact=icontact, sc=.true.)
1086 IF (.NOT. exist) THEN
1087 CALL cp_warn(__location__, &
1088 "User requested to read the overlap matrix from the file named: "// &
1089 trim(filename_sc(icontact))//". This file does not exist. The file will be created.")
1090 exist_all = .false.
1091 END IF
1092 IF (nspins == 1) THEN
1093 CALL negf_restart_file_name(filename_hc_1(icontact), exist, negf_section, logger, icontact=icontact, &
1094 ispin=0, hc=.true.)
1095 IF (.NOT. exist) THEN
1096 CALL cp_warn(__location__, &
1097 "User requested to read the Hamiltonian matrix from the file named: "// &
1098 trim(filename_hc_1(icontact))//". This file does not exist. The file will be created.")
1099 exist_all = .false.
1100 END IF
1101 END IF
1102 IF (nspins == 2) THEN
1103 CALL negf_restart_file_name(filename_hc_1(icontact), exist, negf_section, logger, icontact=icontact, &
1104 ispin=1, hc=.true.)
1105 IF (.NOT. exist) THEN
1106 CALL cp_warn(__location__, &
1107 "User requested to read the Hamiltonian matrix from the file named: "// &
1108 trim(filename_hc_1(icontact))//". This file does not exist. The file will be created.")
1109 exist_all = .false.
1110 END IF
1111 CALL negf_restart_file_name(filename_hc_2(icontact), exist, negf_section, logger, icontact=icontact, &
1112 ispin=2, hc=.true.)
1113 IF (.NOT. exist) THEN
1114 CALL cp_warn(__location__, &
1115 "User requested to read the Hamiltonian matrix from the file named: "// &
1116 trim(filename_hc_2(icontact))//". This file does not exist. The file will be created.")
1117 exist_all = .false.
1118 END IF
1119 END IF
1120 END DO
1121 END IF
1122 CALL para_env%bcast(exist_all)
1123
1124 IF (exist_all) THEN
1125
1126 negf_control%is_restart = .true.
1127
1128 IF (log_unit > 0) THEN
1129 WRITE (log_unit, '(/,T2,A)') "User requested to read the Hamiltonian and overlap matrices from files."
1130 WRITE (log_unit, '(T2,A)') "All restart files exist."
1131 END IF
1132
1133 ! ++ create matrices: s_s, s_sc, h_s, h_sc
1134 IF (para_env%is_source()) THEN
1135 CALL open_file(file_name=filename_s, file_status="OLD", &
1136 file_form="FORMATTED", file_action="READ", &
1137 file_position="REWIND", unit_number=print_unit)
1138 READ (print_unit, *) nrow_s, ncol_s
1139 CALL close_file(print_unit)
1140 END IF
1141 CALL para_env%bcast(nrow_s)
1142 CALL para_env%bcast(ncol_s)
1143 NULLIFY (fm_struct)
1144 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_s, ncol_global=ncol_s, context=sub_env%blacs_env)
1145 ALLOCATE (negf_env%s_s)
1146 CALL cp_fm_create(negf_env%s_s, fm_struct)
1147 ALLOCATE (negf_env%h_s(nspins))
1148 DO ispin = 1, nspins
1149 CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
1150 END DO
1151 CALL cp_fm_struct_release(fm_struct)
1152 ALLOCATE (negf_env%s_sc(ncontacts))
1153 ALLOCATE (negf_env%h_sc(nspins, ncontacts))
1154 DO icontact = 1, ncontacts
1155 IF (para_env%is_source()) THEN
1156 CALL open_file(file_name=filename_sc(icontact), file_status="OLD", &
1157 file_form="FORMATTED", file_action="READ", &
1158 file_position="REWIND", unit_number=print_unit)
1159 READ (print_unit, *) nrow_sc, ncol_sc
1160 CALL close_file(print_unit)
1161 END IF
1162 CALL para_env%bcast(nrow_sc)
1163 CALL para_env%bcast(ncol_sc)
1164 NULLIFY (fm_struct)
1165 CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_sc, ncol_global=ncol_sc, context=sub_env%blacs_env)
1166 CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
1167 DO ispin = 1, nspins
1168 CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
1169 END DO
1170 CALL cp_fm_struct_release(fm_struct)
1171 END DO
1172
1173 ALLOCATE (target_m(nrow_s, ncol_s))
1174 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s, target_m)
1175 CALL para_env%bcast(target_m)
1176 CALL cp_fm_set_submatrix(negf_env%s_s, target_m)
1177 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_s is read from "//trim(filename_s)
1178 IF (nspins == 1) THEN
1179 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_1, target_m)
1180 CALL para_env%bcast(target_m)
1181 CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m)
1182 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//trim(filename_h_1)
1183 END IF
1184 IF (nspins == 2) THEN
1185 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_1, target_m)
1186 CALL para_env%bcast(target_m)
1187 CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m)
1188 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//trim(filename_h_1)//" for spin 1"
1189 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_2, target_m)
1190 CALL para_env%bcast(target_m)
1191 CALL cp_fm_set_submatrix(negf_env%h_s(2), target_m)
1192 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//trim(filename_h_2)//" for spin 2"
1193 END IF
1194 DEALLOCATE (target_m)
1195
1196 DO icontact = 1, ncontacts
1197 ALLOCATE (target_m(nrow_s, ncol_sc))
1198 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_sc(icontact), target_m)
1199 CALL para_env%bcast(target_m)
1200 CALL cp_fm_set_submatrix(negf_env%s_sc(icontact), target_m)
1201 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_sc is read from "//trim(filename_sc(icontact))
1202 IF (nspins == 1) THEN
1203 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_1(icontact), target_m)
1204 CALL para_env%bcast(target_m)
1205 CALL cp_fm_set_submatrix(negf_env%h_sc(1, icontact), target_m)
1206 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//trim(filename_hc_1(icontact))
1207 END IF
1208 IF (nspins == 2) THEN
1209 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_1(icontact), target_m)
1210 CALL para_env%bcast(target_m)
1211 CALL cp_fm_set_submatrix(negf_env%h_sc(1, icontact), target_m)
1212 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//trim(filename_hc_1(icontact))//" for spin 1"
1213 IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_2(icontact), target_m)
1214 CALL para_env%bcast(target_m)
1215 CALL cp_fm_set_submatrix(negf_env%h_sc(2, icontact), target_m)
1216 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//trim(filename_hc_2(icontact))//" for spin 2"
1217 END IF
1218 DEALLOCATE (target_m)
1219
1220 END DO
1221
1222 ELSE
1223
1224 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') &
1225 "Some restart files do not exist. ALL restart files will be recalculated!"
1226
1227 IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.false., calc_forces=.false.)
1228 ! extract device-related matrix blocks
1229 CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
1230 is_dft_entire = .true.
1231
1232 CALL cp_fm_get_info(negf_env%s_s, nrow_global=nrow_s, ncol_global=ncol_s)
1233 ALLOCATE (target_m(nrow_s, ncol_s))
1234 CALL cp_fm_get_submatrix(negf_env%s_s, target_m)
1235 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s, target_m)
1236 IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "S_s is saved to "//trim(filename_s)
1237 IF (nspins == 1) THEN
1238 CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m)
1239 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_1, target_m)
1240 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//trim(filename_h_1)
1241 END IF
1242 IF (nspins == 2) THEN
1243 CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m)
1244 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_1, target_m)
1245 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//trim(filename_h_1)//" for spin 1"
1246 CALL cp_fm_get_submatrix(negf_env%h_s(2), target_m)
1247 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_2, target_m)
1248 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//trim(filename_h_2)//" for spin 2"
1249 END IF
1250 DEALLOCATE (target_m)
1251
1252 DO icontact = 1, ncontacts
1253 CALL cp_fm_get_info(negf_env%contacts(icontact)%s_00, nrow_global=nrow_sc, ncol_global=ncol_sc)
1254 ALLOCATE (target_m(nrow_s, ncol_sc))
1255 CALL cp_fm_get_submatrix(negf_env%s_sc(icontact), target_m)
1256 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_sc(icontact), target_m)
1257 IF (log_unit > 0) WRITE (log_unit, '(T2,A,I3)') &
1258 "S_sc is saved to "//trim(filename_sc(icontact))//" for contact", icontact
1259 IF (nspins == 1) THEN
1260 CALL cp_fm_get_submatrix(negf_env%h_sc(1, icontact), target_m)
1261 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_1(icontact), target_m)
1262 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//trim(filename_hc_1(icontact))
1263 END IF
1264 IF (nspins == 2) THEN
1265 CALL cp_fm_get_submatrix(negf_env%h_sc(1, icontact), target_m)
1266 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_1(icontact), target_m)
1267 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//trim(filename_hc_1(icontact))//" for spin 1"
1268 CALL cp_fm_get_submatrix(negf_env%h_sc(2, icontact), target_m)
1269 IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_2(icontact), target_m)
1270 IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//trim(filename_hc_2(icontact))//" for spin 2"
1271 END IF
1272 DEALLOCATE (target_m)
1273 END DO
1274
1275 negf_control%write_common_restart_file = .true.
1276
1277 END IF
1278
1279 DEALLOCATE (filename_sc, filename_hc_1, filename_hc_2)
1280 CALL timestop(handle)
1281 END SUBROUTINE negf_env_scatt_read_write_hs
1282
1283! **************************************************************************************************
1284!> \brief Extract relevant matrix blocks for the scattering region as well as
1285!> all the scattering -- contact interface regions.
1286!> \param negf_env NEGF environment (modified on exit)
1287!> \param negf_control NEGF control
1288!> \param sub_env NEGF parallel (sub)group environment
1289!> \param qs_env Primary QuickStep environment
1290!> \author Sergey Chulkov
1291! **************************************************************************************************
1292 SUBROUTINE negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
1293 TYPE(negf_env_type), INTENT(inout) :: negf_env
1294 TYPE(negf_control_type), POINTER :: negf_control
1295 TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
1296 TYPE(qs_environment_type), POINTER :: qs_env
1297
1298 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_env_device_init_matrices'
1299
1300 INTEGER :: handle, icontact, ispin, nao_c, nao_s, &
1301 ncontacts, nspins
1302 LOGICAL :: do_kpoints
1303 TYPE(cp_fm_struct_type), POINTER :: fm_struct
1304 TYPE(dbcsr_p_type) :: hmat
1305 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
1306 TYPE(dft_control_type), POINTER :: dft_control
1307 TYPE(mp_para_env_type), POINTER :: para_env
1308 TYPE(pw_env_type), POINTER :: pw_env
1309 TYPE(pw_pool_type), POINTER :: pw_pool
1310 TYPE(pw_r3d_rs_type) :: v_hartree
1311 TYPE(qs_subsys_type), POINTER :: subsys
1312
1313 CALL timeset(routinen, handle)
1314
1315 IF (ALLOCATED(negf_control%atomlist_S_screening)) THEN
1316 CALL get_qs_env(qs_env, &
1317 dft_control=dft_control, &
1318 do_kpoints=do_kpoints, &
1319 matrix_ks_kp=matrix_ks_kp, &
1320 matrix_s_kp=matrix_s_kp, &
1321 para_env=para_env, &
1322 pw_env=pw_env, &
1323 subsys=subsys)
1324 IF (dft_control%qs_control%xtb) CALL rebuild_pw_env(qs_env)
1325 CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
1326
1327 IF (do_kpoints) THEN
1328 CALL cp_abort(__location__, &
1329 "K-points in device region have not been implemented yet.")
1330 END IF
1331
1332 ncontacts = SIZE(negf_control%contacts)
1333 nspins = dft_control%nspins
1334
1335 NULLIFY (fm_struct)
1336 nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening)
1337
1338 ! ++ create matrices: h_s, s_s
1339 NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
1340 ALLOCATE (negf_env%h_s(nspins))
1341
1342 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
1343 ALLOCATE (negf_env%s_s)
1344 CALL cp_fm_create(negf_env%s_s, fm_struct)
1345 DO ispin = 1, nspins
1346 CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
1347 END DO
1348 ALLOCATE (negf_env%v_hartree_s)
1349 CALL cp_fm_create(negf_env%v_hartree_s, fm_struct)
1350 CALL cp_fm_struct_release(fm_struct)
1351
1352 ! ++ create matrices: h_sc, s_sc
1353 ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
1354 DO icontact = 1, ncontacts
1355 nao_c = number_of_atomic_orbitals(subsys, negf_env%contacts(icontact)%atomlist_cell0)
1356 CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
1357
1358 CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
1359
1360 DO ispin = 1, nspins
1361 CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
1362 END DO
1363
1364 CALL cp_fm_struct_release(fm_struct)
1365 END DO
1366
1367 ! extract matrices: h_s, s_s
1368 DO ispin = 1, nspins
1369 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
1370 fm=negf_env%h_s(ispin), &
1371 atomlist_row=negf_control%atomlist_S_screening, &
1372 atomlist_col=negf_control%atomlist_S_screening, &
1373 subsys=subsys, mpi_comm_global=para_env, &
1374 do_upper_diag=.true., do_lower=.true.)
1375 END DO
1376
1377 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
1378 fm=negf_env%s_s, &
1379 atomlist_row=negf_control%atomlist_S_screening, &
1380 atomlist_col=negf_control%atomlist_S_screening, &
1381 subsys=subsys, mpi_comm_global=para_env, &
1382 do_upper_diag=.true., do_lower=.true.)
1383
1384 ! v_hartree_s
1385 NULLIFY (hmat%matrix)
1386 CALL dbcsr_init_p(hmat%matrix)
1387 CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
1388 CALL dbcsr_set(hmat%matrix, 0.0_dp)
1389
1390 CALL pw_pool%create_pw(v_hartree)
1391 CALL negf_env_init_v_hartree(v_hartree, negf_env%contacts, negf_control%contacts)
1392
1393 CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, &
1394 calculate_forces=.false., compute_tau=.false., gapw=.false.)
1395
1396 CALL pw_pool%give_back_pw(v_hartree)
1397
1398 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=hmat%matrix, &
1399 fm=negf_env%v_hartree_s, &
1400 atomlist_row=negf_control%atomlist_S_screening, &
1401 atomlist_col=negf_control%atomlist_S_screening, &
1402 subsys=subsys, mpi_comm_global=para_env, &
1403 do_upper_diag=.true., do_lower=.true.)
1404
1405 CALL dbcsr_deallocate_matrix(hmat%matrix)
1406
1407 ! extract matrices: h_sc, s_sc
1408 DO icontact = 1, ncontacts
1409 DO ispin = 1, nspins
1410 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
1411 fm=negf_env%h_sc(ispin, icontact), &
1412 atomlist_row=negf_control%atomlist_S_screening, &
1413 atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
1414 subsys=subsys, mpi_comm_global=para_env, &
1415 do_upper_diag=.true., do_lower=.true.)
1416 END DO
1417
1418 CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
1419 fm=negf_env%s_sc(icontact), &
1420 atomlist_row=negf_control%atomlist_S_screening, &
1421 atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
1422 subsys=subsys, mpi_comm_global=para_env, &
1423 do_upper_diag=.true., do_lower=.true.)
1424 END DO
1425 END IF
1426
1427 CALL timestop(handle)
1428 END SUBROUTINE negf_env_device_init_matrices
1429
1430! **************************************************************************************************
1431!> \brief Contribution to the Hartree potential related to the external bias voltage.
1432!> \param v_hartree Hartree potential (modified on exit)
1433!> \param contact_env NEGF environment for every contact
1434!> \param contact_control NEGF control for every contact
1435!> \author Sergey Chulkov
1436! **************************************************************************************************
1437 SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control)
1438 TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree
1439 TYPE(negf_env_contact_type), DIMENSION(:), &
1440 INTENT(in) :: contact_env
1441 TYPE(negf_control_contact_type), DIMENSION(:), &
1442 INTENT(in) :: contact_control
1443
1444 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_init_v_hartree'
1445 REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*epsilon(0.0_dp)
1446
1447 INTEGER :: dx, dy, dz, handle, icontact, ix, iy, &
1448 iz, lx, ly, lz, ncontacts, ux, uy, uz
1449 REAL(kind=dp) :: dvol, pot, proj, v1, v2
1450 REAL(kind=dp), DIMENSION(3) :: dirvector_bias, point_coord, &
1451 point_indices, vector
1452
1453 CALL timeset(routinen, handle)
1454
1455 ncontacts = SIZE(contact_env)
1456 cpassert(SIZE(contact_control) == ncontacts)
1457 cpassert(ncontacts == 2)
1458
1459 dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
1460 v1 = contact_control(1)%v_external
1461 v2 = contact_control(2)%v_external
1462
1463 lx = v_hartree%pw_grid%bounds_local(1, 1)
1464 ux = v_hartree%pw_grid%bounds_local(2, 1)
1465 ly = v_hartree%pw_grid%bounds_local(1, 2)
1466 uy = v_hartree%pw_grid%bounds_local(2, 2)
1467 lz = v_hartree%pw_grid%bounds_local(1, 3)
1468 uz = v_hartree%pw_grid%bounds_local(2, 3)
1469
1470 dx = v_hartree%pw_grid%npts(1)/2
1471 dy = v_hartree%pw_grid%npts(2)/2
1472 dz = v_hartree%pw_grid%npts(3)/2
1473
1474 dvol = v_hartree%pw_grid%dvol
1475
1476 DO iz = lz, uz
1477 point_indices(3) = real(iz + dz, kind=dp)
1478 DO iy = ly, uy
1479 point_indices(2) = real(iy + dy, kind=dp)
1480
1481 DO ix = lx, ux
1482 point_indices(1) = real(ix + dx, kind=dp)
1483 point_coord(:) = matmul(v_hartree%pw_grid%dh, point_indices)
1484
1485 vector = point_coord - contact_env(1)%origin_bias
1486 proj = projection_on_direction_vector(vector, dirvector_bias)
1487 IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
1488 ! scattering region
1489 ! proj == 0 we are at the first contact boundary
1490 ! proj == 1 we are at the second contact boundary
1491 IF (proj < 0.0_dp) THEN
1492 proj = 0.0_dp
1493 ELSE IF (proj > 1.0_dp) THEN
1494 proj = 1.0_dp
1495 END IF
1496 pot = v1 + (v2 - v1)*proj
1497 ELSE
1498 pot = 0.0_dp
1499 DO icontact = 1, ncontacts
1500 vector = point_coord - contact_env(icontact)%origin_bias
1501 proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias)
1502
1503 IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
1504 pot = contact_control(icontact)%v_external
1505 EXIT
1506 END IF
1507 END DO
1508 END IF
1509
1510 v_hartree%array(ix, iy, iz) = pot*dvol
1511 END DO
1512 END DO
1513 END DO
1514
1515 CALL timestop(handle)
1516 END SUBROUTINE negf_env_init_v_hartree
1517
1518! **************************************************************************************************
1519!> \brief Detect the axis towards secondary unit cell.
1520!> \param direction_vector direction vector
1521!> \param subsys_contact QuickStep subsystem of the contact force environment
1522!> \param eps_geometry accuracy in mapping atoms between different force environments
1523!> \return direction axis: 0 (undefined), 1 (x), 2(y), 3 (z)
1524!> \par History
1525!> * 08.2017 created [Sergey Chulkov]
1526! **************************************************************************************************
1527 FUNCTION contact_direction_axis(direction_vector, subsys_contact, eps_geometry) RESULT(direction_axis)
1528 REAL(kind=dp), DIMENSION(3), INTENT(in) :: direction_vector
1529 TYPE(qs_subsys_type), POINTER :: subsys_contact
1530 REAL(kind=dp), INTENT(in) :: eps_geometry
1531 INTEGER :: direction_axis
1532
1533 INTEGER :: i, naxes
1534 REAL(kind=dp), DIMENSION(3) :: scaled
1535 TYPE(cell_type), POINTER :: cell
1536
1537 CALL qs_subsys_get(subsys_contact, cell=cell)
1538 CALL real_to_scaled(scaled, direction_vector, cell)
1539
1540 naxes = 0
1541 direction_axis = 0 ! initialize to make GCC<=6 happy
1542
1543 DO i = 1, 3
1544 IF (abs(scaled(i)) > eps_geometry) THEN
1545 IF (scaled(i) > 0.0_dp) THEN
1546 direction_axis = i
1547 ELSE
1548 direction_axis = -i
1549 END IF
1550 naxes = naxes + 1
1551 END IF
1552 END DO
1553
1554 ! direction_vector is not parallel to one of the unit cell's axis
1555 IF (naxes /= 1) direction_axis = 0
1556 END FUNCTION contact_direction_axis
1557
1558! **************************************************************************************************
1559!> \brief Estimate energy of the highest spin-alpha occupied molecular orbital.
1560!> \param homo_energy HOMO energy (initialised on exit)
1561!> \param qs_env QuickStep environment
1562!> \par History
1563!> * 01.2017 created [Sergey Chulkov]
1564! **************************************************************************************************
1565 SUBROUTINE negf_homo_energy_estimate(homo_energy, qs_env)
1566 REAL(kind=dp), INTENT(out) :: homo_energy
1567 TYPE(qs_environment_type), POINTER :: qs_env
1568
1569 CHARACTER(LEN=*), PARAMETER :: routinen = 'negf_homo_energy_estimate'
1570 INTEGER, PARAMETER :: gamma_point = 1
1571
1572 INTEGER :: handle, homo, ikpgr, ikpoint, imo, &
1573 ispin, kplocal, nmo, nspins
1574 INTEGER, DIMENSION(2) :: kp_range
1575 LOGICAL :: do_kpoints
1576 REAL(kind=dp) :: my_homo_energy
1577 REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1578 TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
1579 TYPE(kpoint_type), POINTER :: kpoints
1580 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1581 TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
1582 TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
1583
1584 CALL timeset(routinen, handle)
1585 my_homo_energy = 0.0_dp
1586
1587 CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
1588
1589 IF (do_kpoints) THEN
1590 CALL get_kpoint_info(kpoints, kp_env=kp_env, kp_range=kp_range, para_env_kp=para_env_kp)
1591
1592 ! looking for a processor that holds the gamma point
1593 IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point) THEN
1594 kplocal = kp_range(2) - kp_range(1) + 1
1595
1596 DO ikpgr = 1, kplocal
1597 CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
1598
1599 IF (ikpoint == gamma_point) THEN
1600 ! mos_kp(component, spin), where component = 1 (real), or 2 (imaginary)
1601 CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues) ! mu=fermi_level
1602
1603 my_homo_energy = eigenvalues(homo)
1604 EXIT
1605 END IF
1606 END DO
1607 END IF
1608
1609 CALL para_env%sum(my_homo_energy)
1610 ELSE
1611 ! Hamiltonian of the bulk contact region has been computed without k-points.
1612 ! Try to obtain the HOMO energy assuming there is no OT. We probably should abort here
1613 ! as we do need a second replica of the bulk contact unit cell along transport
1614 ! direction anyway which is not available without k-points.
1615
1616 CALL cp_abort(__location__, &
1617 "It is necessary to use k-points along the transport direction "// &
1618 "for all contact FORCE_EVAL-s")
1619 ! It is necessary to use k-points along the transport direction within all contact FORCE_EVAL-s
1620
1621 nspins = SIZE(mos)
1622
1623 spin_loop: DO ispin = 1, nspins
1624 CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
1625
1626 DO imo = nmo, 1, -1
1627 IF (eigenvalues(imo) /= 0.0_dp) EXIT spin_loop
1628 END DO
1629 END DO spin_loop
1630
1631 IF (imo == 0) THEN
1632 cpabort("Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
1633 END IF
1634
1635 my_homo_energy = eigenvalues(homo)
1636 END IF
1637
1638 homo_energy = my_homo_energy
1639 CALL timestop(handle)
1640 END SUBROUTINE negf_homo_energy_estimate
1641
1642! **************************************************************************************************
1643!> \brief List atoms from the contact's primary unit cell.
1644!> \param atomlist_cell0 list of atoms belonging to the contact's primary unit cell
1645!> (allocate and initialised on exit)
1646!> \param atom_map_cell0 atomic map of atoms from 'atomlist_cell0' (allocate and initialised on exit)
1647!> \param atomlist_bulk list of atoms belonging to the bulk contact region
1648!> \param atom_map atomic map of atoms from 'atomlist_bulk'
1649!> \param origin origin of the contact
1650!> \param direction_vector direction vector of the contact
1651!> \param direction_axis axis towards secondary unit cell
1652!> \param subsys_device QuickStep subsystem of the device force environment
1653!> \par History
1654!> * 08.2017 created [Sergey Chulkov]
1655! **************************************************************************************************
1656 SUBROUTINE list_atoms_in_bulk_primary_unit_cell(atomlist_cell0, atom_map_cell0, atomlist_bulk, atom_map, &
1657 origin, direction_vector, direction_axis, subsys_device)
1658 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell0
1659 TYPE(negf_atom_map_type), ALLOCATABLE, &
1660 DIMENSION(:), INTENT(inout) :: atom_map_cell0
1661 INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1662 TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1663 REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1664 INTEGER, INTENT(in) :: direction_axis
1665 TYPE(qs_subsys_type), POINTER :: subsys_device
1666
1667 CHARACTER(LEN=*), PARAMETER :: routinen = 'list_atoms_in_bulk_primary_unit_cell'
1668
1669 INTEGER :: atom_min, dir_axis_min, &
1670 direction_axis_abs, handle, iatom, &
1671 natoms_bulk, natoms_cell0
1672 REAL(kind=dp) :: proj, proj_min
1673 REAL(kind=dp), DIMENSION(3) :: vector
1674 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1675
1676 CALL timeset(routinen, handle)
1677 CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1678
1679 natoms_bulk = SIZE(atomlist_bulk)
1680 cpassert(SIZE(atom_map, 1) == natoms_bulk)
1681 direction_axis_abs = abs(direction_axis)
1682
1683 ! looking for the nearest atom from the scattering region
1684 proj_min = 1.0_dp
1685 atom_min = 1
1686 DO iatom = 1, natoms_bulk
1687 vector = particle_set(atomlist_bulk(iatom))%r - origin
1688 proj = projection_on_direction_vector(vector, direction_vector)
1689
1690 IF (proj < proj_min) THEN
1691 proj_min = proj
1692 atom_min = iatom
1693 END IF
1694 END DO
1695
1696 dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1697
1698 natoms_cell0 = 0
1699 DO iatom = 1, natoms_bulk
1700 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
1701 natoms_cell0 = natoms_cell0 + 1
1702 END IF
1703 END DO
1704
1705 ALLOCATE (atomlist_cell0(natoms_cell0))
1706 ALLOCATE (atom_map_cell0(natoms_cell0))
1707
1708 natoms_cell0 = 0
1709 DO iatom = 1, natoms_bulk
1710 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
1711 natoms_cell0 = natoms_cell0 + 1
1712 atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
1713 atom_map_cell0(natoms_cell0) = atom_map(iatom)
1714 END IF
1715 END DO
1716
1717 CALL timestop(handle)
1718 END SUBROUTINE list_atoms_in_bulk_primary_unit_cell
1719
1720! **************************************************************************************************
1721!> \brief List atoms from the contact's secondary unit cell.
1722!> \param atomlist_cell1 list of atoms belonging to the contact's secondary unit cell
1723!> (allocate and initialised on exit)
1724!> \param atom_map_cell1 atomic map of atoms from 'atomlist_cell1'
1725!> (allocate and initialised on exit)
1726!> \param atomlist_bulk list of atoms belonging to the bulk contact region
1727!> \param atom_map atomic map of atoms from 'atomlist_bulk'
1728!> \param origin origin of the contact
1729!> \param direction_vector direction vector of the contact
1730!> \param direction_axis axis towards the secondary unit cell
1731!> \param subsys_device QuickStep subsystem of the device force environment
1732!> \par History
1733!> * 11.2017 created [Sergey Chulkov]
1734!> \note Cloned from list_atoms_in_bulk_primary_unit_cell. Will be removed once we can managed to
1735!> maintain consistency between real-space matrices from different force_eval sections.
1736! **************************************************************************************************
1737 SUBROUTINE list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1, atom_map_cell1, atomlist_bulk, atom_map, &
1738 origin, direction_vector, direction_axis, subsys_device)
1739 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell1
1740 TYPE(negf_atom_map_type), ALLOCATABLE, &
1741 DIMENSION(:), INTENT(inout) :: atom_map_cell1
1742 INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1743 TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1744 REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1745 INTEGER, INTENT(in) :: direction_axis
1746 TYPE(qs_subsys_type), POINTER :: subsys_device
1747
1748 CHARACTER(LEN=*), PARAMETER :: routinen = 'list_atoms_in_bulk_secondary_unit_cell'
1749
1750 INTEGER :: atom_min, dir_axis_min, &
1751 direction_axis_abs, handle, iatom, &
1752 natoms_bulk, natoms_cell1, offset
1753 REAL(kind=dp) :: proj, proj_min
1754 REAL(kind=dp), DIMENSION(3) :: vector
1755 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1756
1757 CALL timeset(routinen, handle)
1758 CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1759
1760 natoms_bulk = SIZE(atomlist_bulk)
1761 cpassert(SIZE(atom_map, 1) == natoms_bulk)
1762 direction_axis_abs = abs(direction_axis)
1763 offset = sign(1, direction_axis)
1764
1765 ! looking for the nearest atom from the scattering region
1766 proj_min = 1.0_dp
1767 atom_min = 1
1768 DO iatom = 1, natoms_bulk
1769 vector = particle_set(atomlist_bulk(iatom))%r - origin
1770 proj = projection_on_direction_vector(vector, direction_vector)
1771
1772 IF (proj < proj_min) THEN
1773 proj_min = proj
1774 atom_min = iatom
1775 END IF
1776 END DO
1777
1778 dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1779
1780 natoms_cell1 = 0
1781 DO iatom = 1, natoms_bulk
1782 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
1783 natoms_cell1 = natoms_cell1 + 1
1784 END IF
1785 END DO
1786
1787 ALLOCATE (atomlist_cell1(natoms_cell1))
1788 ALLOCATE (atom_map_cell1(natoms_cell1))
1789
1790 natoms_cell1 = 0
1791 DO iatom = 1, natoms_bulk
1792 IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
1793 natoms_cell1 = natoms_cell1 + 1
1794 atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
1795 atom_map_cell1(natoms_cell1) = atom_map(iatom)
1796 atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
1797 END IF
1798 END DO
1799
1800 CALL timestop(handle)
1801 END SUBROUTINE list_atoms_in_bulk_secondary_unit_cell
1802
1803! **************************************************************************************************
1804!> \brief Release a NEGF environment variable.
1805!> \param negf_env NEGF environment to release
1806!> \par History
1807!> * 01.2017 created [Sergey Chulkov]
1808! **************************************************************************************************
1809 SUBROUTINE negf_env_release(negf_env)
1810 TYPE(negf_env_type), INTENT(inout) :: negf_env
1811
1812 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_release'
1813
1814 INTEGER :: handle, icontact
1815
1816 CALL timeset(routinen, handle)
1817
1818 IF (ALLOCATED(negf_env%contacts)) THEN
1819 DO icontact = SIZE(negf_env%contacts), 1, -1
1820 CALL negf_env_contact_release(negf_env%contacts(icontact))
1821 END DO
1822
1823 DEALLOCATE (negf_env%contacts)
1824 END IF
1825
1826 ! h_s
1827 CALL cp_fm_release(negf_env%h_s)
1828
1829 ! h_sc
1830 CALL cp_fm_release(negf_env%h_sc)
1831
1832 ! s_s
1833 IF (ASSOCIATED(negf_env%s_s)) THEN
1834 CALL cp_fm_release(negf_env%s_s)
1835 DEALLOCATE (negf_env%s_s)
1836 NULLIFY (negf_env%s_s)
1837 END IF
1838
1839 ! s_sc
1840 CALL cp_fm_release(negf_env%s_sc)
1841
1842 ! v_hartree_s
1843 IF (ASSOCIATED(negf_env%v_hartree_s)) THEN
1844 CALL cp_fm_release(negf_env%v_hartree_s)
1845 DEALLOCATE (negf_env%v_hartree_s)
1846 NULLIFY (negf_env%v_hartree_s)
1847 END IF
1848
1849 ! density mixing
1850 IF (ASSOCIATED(negf_env%mixing_storage)) THEN
1851 CALL mixing_storage_release(negf_env%mixing_storage)
1852 DEALLOCATE (negf_env%mixing_storage)
1853 END IF
1854
1855 CALL timestop(handle)
1856 END SUBROUTINE negf_env_release
1857
1858! **************************************************************************************************
1859!> \brief Release a NEGF contact environment variable.
1860!> \param contact_env NEGF contact environment to release
1861! **************************************************************************************************
1862 SUBROUTINE negf_env_contact_release(contact_env)
1863 TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
1864
1865 CHARACTER(len=*), PARAMETER :: routinen = 'negf_env_contact_release'
1866
1867 INTEGER :: handle
1868
1869 CALL timeset(routinen, handle)
1870
1871 ! h_00
1872 CALL cp_fm_release(contact_env%h_00)
1873
1874 ! h_01
1875 CALL cp_fm_release(contact_env%h_01)
1876
1877 ! rho_00
1878 CALL cp_fm_release(contact_env%rho_00)
1879
1880 ! rho_01
1881 CALL cp_fm_release(contact_env%rho_01)
1882
1883 ! s_00
1884 IF (ASSOCIATED(contact_env%s_00)) THEN
1885 CALL cp_fm_release(contact_env%s_00)
1886 DEALLOCATE (contact_env%s_00)
1887 NULLIFY (contact_env%s_00)
1888 END IF
1889
1890 ! s_01
1891 IF (ASSOCIATED(contact_env%s_01)) THEN
1892 CALL cp_fm_release(contact_env%s_01)
1893 DEALLOCATE (contact_env%s_01)
1894 NULLIFY (contact_env%s_01)
1895 END IF
1896
1897 IF (ALLOCATED(contact_env%atomlist_cell0)) DEALLOCATE (contact_env%atomlist_cell0)
1898 IF (ALLOCATED(contact_env%atomlist_cell1)) DEALLOCATE (contact_env%atomlist_cell1)
1899 IF (ALLOCATED(contact_env%atom_map_cell0)) DEALLOCATE (contact_env%atom_map_cell0)
1900
1901 CALL timestop(handle)
1902 END SUBROUTINE negf_env_contact_release
1903
1904END MODULE negf_env_types
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
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:595
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:371
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:1204
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:355
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:1038
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:750
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
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:925
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:344
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:735
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:1052
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:1105
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, gamma_centered)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:529
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:1810
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:306
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:336
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:646
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:246
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:124
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, matrix_vhxc, 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:534
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:401
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:1043
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:137
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:757
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:114
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:97
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:175
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
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:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
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