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