Line data Source code
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 : ! **************************************************************************************************
11 : MODULE 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 = 0.0_dp
116 : !> number of electrons Sp(rho_01,s_01)
117 : REAL(kind=dp) :: nelectrons_qs_cell1 = 0.0_dp
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 = 0.0_dp
151 : !> number of electrons Sp(rho_s,s_s)
152 : REAL(kind=dp) :: nelectrons = 0.0_dp
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 :
163 : CONTAINS
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 4 : 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 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
193 : TYPE(dft_control_type), POINTER :: dft_control
194 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
204 :
205 : ! ensure we have Quickstep enabled for all force_env
206 4 : NULLIFY (sub_force_env)
207 : CALL force_env_get(force_env, in_use=in_use, qs_env=qs_env, root_section=root_section, &
208 4 : sub_force_env=sub_force_env)
209 :
210 4 : IF (ASSOCIATED(sub_force_env)) THEN
211 2 : n_force_env = SIZE(sub_force_env)
212 : ELSE
213 2 : n_force_env = 0
214 : END IF
215 :
216 4 : IF (in_use == use_qs_force) THEN
217 8 : DO icontact = 1, n_force_env
218 4 : CALL force_env_get(sub_force_env(icontact)%force_env, in_use=in_use)
219 8 : IF (in_use /= use_qs_force) EXIT
220 : END DO
221 : END IF
222 :
223 4 : IF (in_use /= use_qs_force) THEN
224 0 : CPABORT("Quickstep is required for NEGF run.")
225 : END IF
226 :
227 : ! check that all mentioned FORCE_EVAL sections are actually present
228 4 : ncontacts = SIZE(negf_control%contacts)
229 :
230 12 : DO icontact = 1, ncontacts
231 12 : IF (negf_control%contacts(icontact)%force_env_index > n_force_env) THEN
232 0 : WRITE (contact_str, '(I11)') icontact
233 0 : WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
234 0 : 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 0 : " 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 4 : 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 4 : para_env=para_env, subsys=subsys, v_hartree_rspace=v_hartree_rspace)
250 :
251 4 : negf_section => section_vals_get_subs_vals(root_section, "NEGF")
252 :
253 4 : IF (do_kpoints) THEN
254 0 : 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 80 : ALLOCATE (negf_env%contacts(ncontacts))
259 20 : ALLOCATE (map_contact(ncontacts))
260 :
261 12 : DO icontact = 1, ncontacts
262 12 : IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
263 4 : CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
264 4 : 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 4 : subsys_contact=subsys_contact)
272 :
273 4 : IF (negf_env%contacts(icontact)%direction_axis == 0) THEN
274 0 : WRITE (contact_str, '(I11)') icontact
275 0 : 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 0 : 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 12 : DO icontact = 1, ncontacts
286 12 : IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
287 4 : 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 0 : para_env, negf_env, sub_env, negf_control, negf_section, log_unit, is_separate=.TRUE.)
291 : ELSE
292 4 : IF (log_unit > 0) &
293 2 : WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') "NEGF| Construct the Kohn-Sham matrix for the contact", icontact, &
294 4 : " from the separate bulk DFT calculation"
295 4 : CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
296 4 : 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 4 : qs_env_contact=qs_env_contact)
299 4 : 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 4 : is_dft_entire = .FALSE.
306 12 : DO icontact = 1, ncontacts
307 12 : IF (negf_control%contacts(icontact)%force_env_index <= 0) THEN
308 4 : 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 0 : 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 0 : log_unit, is_separate=.FALSE., is_dft_entire=is_dft_entire)
315 : ELSE
316 4 : IF (log_unit > 0) &
317 2 : WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') "NEGF| Construct the Kohn-Sham matrix for the contact", icontact, &
318 4 : " from the entire system bulk DFT calculation"
319 4 : IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.)
320 4 : 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 4 : eps_geometry=negf_control%eps_geometry)
325 4 : 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 4 : IF (log_unit > 0) &
332 2 : WRITE (log_unit, '(/,T2,A,T70)') "NEGF| Construct the Kohn-Sham matrix for the scattering region"
333 4 : 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 0 : is_dft_entire=is_dft_entire)
336 : ELSE
337 4 : IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.)
338 : ! extract device-related matrix blocks
339 4 : CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
340 : END IF
341 4 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
342 :
343 4 : 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 4 : NULLIFY (negf_env%mixing_storage)
348 4 : CALL section_vals_val_get(negf_mixing_section, "METHOD", i_val=negf_env%mixing_method)
349 :
350 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
351 16 : ALLOCATE (negf_env%mixing_storage)
352 : CALL mixing_storage_create(negf_env%mixing_storage, negf_mixing_section, &
353 4 : negf_env%mixing_method, dft_control%qs_control%cutoff)
354 :
355 4 : CALL timestop(handle)
356 16 : 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 4 : 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 4 : 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 4 : subsys_device)
390 :
391 4 : contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
392 :
393 4 : IF (contact_env%direction_axis /= 0) THEN
394 4 : natoms = SIZE(contact_control%atomlist_bulk)
395 56 : 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 4 : 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 4 : 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 4 : subsys_device=subsys_device)
424 : END IF
425 :
426 4 : CALL timestop(handle)
427 4 : 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 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
473 0 : logger => cp_get_default_logger()
474 :
475 0 : CALL force_env_get(el_force_env, qs_env=qs_env_contact)
476 0 : CALL get_qs_env(qs_env_contact, dft_control=dft_control, subsys=subsys)
477 0 : nspins = dft_control%nspins
478 :
479 0 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11)') &
480 0 : "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 0 : exist_all = .TRUE.
485 0 : IF (para_env%is_source()) THEN
486 0 : CALL negf_restart_file_name(filename_s00, exist, negf_section, logger, icontact, s00=.TRUE.)
487 0 : IF (.NOT. exist) THEN
488 : CALL cp_warn(__LOCATION__, &
489 : "User requested to read the overlap matrix from the file named: "// &
490 0 : TRIM(filename_s00)//". This file does not exist. The file will be created.")
491 0 : exist_all = .FALSE.
492 : END IF
493 0 : CALL negf_restart_file_name(filename_s01, exist, negf_section, logger, icontact, s01=.TRUE.)
494 0 : IF (.NOT. exist) THEN
495 : CALL cp_warn(__LOCATION__, &
496 : "User requested to read the overlap matrix from the file named: "// &
497 0 : TRIM(filename_s01)//". This file does not exist. The file will be created.")
498 0 : exist_all = .FALSE.
499 : END IF
500 0 : IF (nspins == 1) THEN
501 0 : CALL negf_restart_file_name(filename_h00_1, exist, negf_section, logger, icontact, ispin=0, h00=.TRUE.)
502 0 : IF (.NOT. exist) THEN
503 : CALL cp_warn(__LOCATION__, &
504 : "User requested to read the Hamiltonian matrix from the file named: "// &
505 0 : TRIM(filename_h00_1)//". This file does not exist. The file will be created.")
506 0 : exist_all = .FALSE.
507 : END IF
508 0 : CALL negf_restart_file_name(filename_h01_1, exist, negf_section, logger, icontact, ispin=0, h01=.TRUE.)
509 0 : IF (.NOT. exist) THEN
510 : CALL cp_warn(__LOCATION__, &
511 : "User requested to read the Hamiltonian matrix from the file named: "// &
512 0 : TRIM(filename_h01_1)//". This file does not exist. The file will be created.")
513 0 : exist_all = .FALSE.
514 : END IF
515 : END IF
516 0 : IF (nspins == 2) THEN
517 0 : CALL negf_restart_file_name(filename_h00_1, exist, negf_section, logger, icontact, ispin=1, h00=.TRUE.)
518 0 : IF (.NOT. exist) THEN
519 : CALL cp_warn(__LOCATION__, &
520 : "User requested to read the Hamiltonian matrix from the file named: "// &
521 0 : TRIM(filename_h00_1)//". This file does not exist. The file will be created.")
522 0 : exist_all = .FALSE.
523 : END IF
524 0 : CALL negf_restart_file_name(filename_h01_1, exist, negf_section, logger, icontact, ispin=1, h01=.TRUE.)
525 0 : IF (.NOT. exist) THEN
526 : CALL cp_warn(__LOCATION__, &
527 : "User requested to read tthe Hamiltonian matrix from the file named: "// &
528 0 : TRIM(filename_h01_1)//". This file does not exist. The file will be created.")
529 0 : exist_all = .FALSE.
530 : END IF
531 0 : CALL negf_restart_file_name(filename_h00_2, exist, negf_section, logger, icontact, ispin=2, h00=.TRUE.)
532 0 : IF (.NOT. exist) THEN
533 : CALL cp_warn(__LOCATION__, &
534 : "User requested to read the Hamiltonian matrix from the file named: "// &
535 0 : TRIM(filename_h00_2)//". This file does not exist. The file will be created.")
536 0 : exist_all = .FALSE.
537 : END IF
538 0 : CALL negf_restart_file_name(filename_h01_2, exist, negf_section, logger, icontact, ispin=2, h01=.TRUE.)
539 0 : IF (.NOT. exist) THEN
540 : CALL cp_warn(__LOCATION__, &
541 : "User requested to read the Hamiltonian matrix from the file named: "// &
542 0 : TRIM(filename_h01_2)//". This file does not exist. The file will be created.")
543 0 : exist_all = .FALSE.
544 : END IF
545 : END IF
546 : END IF
547 0 : CALL para_env%bcast(exist_all)
548 :
549 0 : IF (exist_all) THEN
550 :
551 0 : negf_control%contacts(icontact)%is_restart = .TRUE.
552 :
553 0 : 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 0 : 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 0 : file_position="REWIND", unit_number=print_unit)
560 0 : READ (print_unit, *) nrow, ncol
561 0 : CALL close_file(print_unit)
562 : END IF
563 0 : CALL para_env%bcast(nrow)
564 0 : CALL para_env%bcast(ncol)
565 0 : NULLIFY (fm_struct)
566 0 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow, ncol_global=ncol, context=sub_env%blacs_env)
567 0 : ALLOCATE (negf_env%contacts(icontact)%s_00, negf_env%contacts(icontact)%s_01)
568 0 : CALL cp_fm_create(negf_env%contacts(icontact)%s_00, fm_struct)
569 0 : CALL cp_fm_create(negf_env%contacts(icontact)%s_01, fm_struct)
570 0 : ALLOCATE (negf_env%contacts(icontact)%h_00(nspins), negf_env%contacts(icontact)%h_01(nspins))
571 0 : DO ispin = 1, nspins
572 0 : CALL cp_fm_create(negf_env%contacts(icontact)%h_00(ispin), fm_struct)
573 0 : CALL cp_fm_create(negf_env%contacts(icontact)%h_01(ispin), fm_struct)
574 : END DO
575 0 : CALL cp_fm_struct_release(fm_struct)
576 :
577 0 : ALLOCATE (target_m(nrow, ncol))
578 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s00, target_m)
579 0 : CALL para_env%bcast(target_m)
580 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%s_00, target_m)
581 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_00 is read from "//TRIM(filename_s00)
582 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s01, target_m)
583 0 : CALL para_env%bcast(target_m)
584 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%s_01, target_m)
585 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_01 is read from "//TRIM(filename_s01)
586 0 : IF (nspins == 1) THEN
587 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_1, target_m)
588 0 : CALL para_env%bcast(target_m)
589 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
590 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//TRIM(filename_h00_1)
591 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_1, target_m)
592 0 : CALL para_env%bcast(target_m)
593 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
594 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//TRIM(filename_H01_1)
595 : END IF
596 0 : IF (nspins == 2) THEN
597 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_1, target_m)
598 0 : CALL para_env%bcast(target_m)
599 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
600 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//TRIM(filename_h00_1)//" for spin 1"
601 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_1, target_m)
602 0 : CALL para_env%bcast(target_m)
603 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
604 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is read from "//TRIM(filename_H01_1)//" for spin 1"
605 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h00_2, target_m)
606 0 : CALL para_env%bcast(target_m)
607 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_00(2), target_m)
608 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is read from "//TRIM(filename_h00_2)//" for spin 2"
609 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h01_2, target_m)
610 0 : CALL para_env%bcast(target_m)
611 0 : CALL cp_fm_set_submatrix(negf_env%contacts(icontact)%h_01(2), target_m)
612 0 : 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 0 : DEALLOCATE (target_m)
615 :
616 : ELSE
617 :
618 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') &
619 0 : "Some restart files do not exist. ALL restart files will be recalculated!"
620 :
621 0 : IF (is_separate) THEN
622 0 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') &
623 0 : "Construct the Kohn-Sham matrix from from the separate bulk DFT calculation"
624 0 : 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 0 : qs_env_contact=qs_env_contact)
627 : ELSE
628 0 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,A,T70,I11,/,A)') &
629 0 : "Construct the Kohn-Sham matrix from the entire system bulk DFT calculation"
630 0 : negf_control%contacts(icontact)%read_write_HS = .FALSE.
631 0 : 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 0 : eps_geometry=negf_control%eps_geometry)
636 0 : negf_control%contacts(icontact)%read_write_HS = .TRUE.
637 0 : is_dft_entire = .TRUE.
638 : END IF
639 :
640 0 : CALL cp_fm_get_info(negf_env%contacts(icontact)%s_00, nrow_global=nrow)
641 0 : ALLOCATE (target_m(nrow, nrow))
642 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%s_00, target_m)
643 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s00, target_m)
644 0 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "S_00 is saved to "//TRIM(filename_s00)
645 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%s_01, target_m)
646 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s01, target_m)
647 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_01 is saved to "//TRIM(filename_s01)
648 0 : IF (nspins == 1) THEN
649 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
650 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_1, target_m)
651 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//TRIM(filename_h00_1)
652 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
653 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_1, target_m)
654 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//TRIM(filename_h01_1)
655 : END IF
656 0 : IF (nspins == 2) THEN
657 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(1), target_m)
658 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_1, target_m)
659 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//TRIM(filename_h00_1)//" for spin 1"
660 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(1), target_m)
661 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_1, target_m)
662 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_01 is saved to "//TRIM(filename_h01_1)//" for spin 1"
663 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_00(2), target_m)
664 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h00_2, target_m)
665 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_00 is saved to "//TRIM(filename_h00_2)//" for spin 2"
666 0 : CALL cp_fm_get_submatrix(negf_env%contacts(icontact)%h_01(2), target_m)
667 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h01_2, target_m)
668 0 : 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 0 : DEALLOCATE (target_m)
671 :
672 0 : negf_control%write_common_restart_file = .TRUE.
673 :
674 : END IF
675 :
676 0 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
677 :
678 0 : CALL timestop(handle)
679 0 : 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 4 : 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 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_list0, atom_list1
701 4 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
702 4 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
703 : LOGICAL :: do_kpoints
704 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
705 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matkp
706 4 : 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 4 : 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 4 : subsys=subsys)
724 4 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
725 :
726 4 : CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
727 :
728 4 : natoms = SIZE(contact_env%atomlist_cell0)
729 12 : ALLOCATE (atom_list0(natoms))
730 20 : DO iatom = 1, natoms
731 16 : 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 68 : IF (SUM(ABS(contact_env%atom_map_cell0(iatom)%cell(:))) > 0) THEN
736 0 : CPABORT("NEGF K-points are not currently supported")
737 : END IF
738 : END DO
739 :
740 4 : CPASSERT(SIZE(contact_env%atomlist_cell1) == natoms)
741 8 : ALLOCATE (atom_list1(natoms))
742 20 : DO iatom = 1, natoms
743 20 : atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
744 : END DO
745 :
746 4 : nspins = dft_control%nspins
747 4 : nimages = dft_control%nimages
748 :
749 4 : IF (do_kpoints) THEN
750 4 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
751 : ELSE
752 0 : ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
753 0 : cell_to_index(0, 0, 0) = 1
754 : END IF
755 :
756 12 : ALLOCATE (index_to_cell(3, nimages))
757 4 : CALL invert_cell_to_index(cell_to_index, nimages, index_to_cell)
758 4 : IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
759 :
760 4 : NULLIFY (fm_struct)
761 4 : nao = number_of_atomic_orbitals(subsys, atom_list0)
762 4 : 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 4 : ALLOCATE (contact_env%s_00, contact_env%s_01)
766 4 : CALL cp_fm_create(contact_env%s_00, fm_struct)
767 4 : CALL cp_fm_create(contact_env%s_01, fm_struct)
768 :
769 : ! ++ create matrices: h_00, h_01, rho_00, rho_01
770 24 : ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
771 20 : ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
772 8 : DO ispin = 1, nspins
773 4 : CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
774 4 : CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
775 4 : CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
776 8 : CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
777 : END DO
778 :
779 4 : CALL cp_fm_struct_release(fm_struct)
780 :
781 : ! extract matrices: s_00, s_01
782 4 : 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 4 : kpoints=kpoints)
790 :
791 : ! extract matrices: h_00, h_01, rho_00, rho_01
792 8 : DO ispin = 1, nspins
793 4 : 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 4 : kpoints=kpoints)
801 :
802 4 : 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 8 : kpoints=kpoints)
810 : END DO
811 :
812 4 : DEALLOCATE (atom_list0, atom_list1)
813 :
814 4 : CALL timestop(handle)
815 8 : 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 4 : 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 4 : 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 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
845 : TYPE(qs_rho_type), POINTER :: rho_struct
846 : TYPE(qs_subsys_type), POINTER :: subsys
847 :
848 4 : 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 4 : subsys=subsys)
858 4 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
859 :
860 4 : IF (do_kpoints) THEN
861 : CALL cp_abort(__LOCATION__, &
862 0 : "K-points in device region have not been implemented yet.")
863 : END IF
864 :
865 4 : nspins = dft_control%nspins
866 :
867 4 : nao_c = number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(1)%vector)
868 4 : 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 0 : "No single atom, however, can be shared between these two cells.")
873 : END IF
874 :
875 4 : 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 4 : subsys)
884 :
885 4 : 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 4 : CALL qs_subsys_get(subsys, particle_set=particle_set)
889 :
890 16 : origin = particle_set(contact_control%atomlist_screening(1))%r
891 16 : DO iatom = 2, SIZE(contact_control%atomlist_screening)
892 52 : origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
893 : END DO
894 16 : origin = origin/REAL(SIZE(contact_control%atomlist_screening), kind=dp)
895 :
896 12 : DO icell = 1, 2
897 32 : direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
898 32 : DO iatom = 2, SIZE(contact_control%atomlist_cell(icell)%vector)
899 104 : direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
900 : END DO
901 32 : direction_vector = direction_vector/REAL(SIZE(contact_control%atomlist_cell(icell)%vector), kind=dp)
902 32 : direction_vector = direction_vector - origin
903 36 : r2_origin_cell(icell) = DOT_PRODUCT(direction_vector, direction_vector)
904 : END DO
905 :
906 4 : 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 0 : "Primary and secondary bulk contact cells should not overlap ")
911 4 : ELSE IF (r2_origin_cell(1) < r2_origin_cell(2)) THEN
912 2 : IF (.NOT. ALLOCATED(contact_env%atomlist_cell0)) &
913 6 : ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(1)%vector)))
914 10 : contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
915 2 : IF (.NOT. ALLOCATED(contact_env%atomlist_cell1)) &
916 6 : ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(2)%vector)))
917 10 : contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
918 : ELSE
919 2 : IF (.NOT. ALLOCATED(contact_env%atomlist_cell0)) &
920 6 : ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(2)%vector)))
921 10 : contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
922 2 : IF (.NOT. ALLOCATED(contact_env%atomlist_cell1)) &
923 6 : ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(1)%vector)))
924 10 : contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
925 : END IF
926 4 : IF (.NOT. contact_control%read_write_HS) THEN
927 4 : NULLIFY (fm_struct)
928 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
929 28 : ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
930 28 : ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
931 8 : DO ispin = 1, nspins
932 4 : CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
933 4 : CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
934 4 : CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
935 8 : CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
936 : END DO
937 4 : ALLOCATE (contact_env%s_00, contact_env%s_01)
938 4 : CALL cp_fm_create(contact_env%s_00, fm_struct)
939 4 : CALL cp_fm_create(contact_env%s_01, fm_struct)
940 4 : CALL cp_fm_struct_release(fm_struct)
941 :
942 8 : 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 4 : 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 4 : 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 4 : 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 8 : 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 4 : 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 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
982 : END IF
983 4 : CALL timestop(handle)
984 4 : 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 0 : 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 0 : 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 0 : 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 0 : CALL timeset(routineN, handle)
1027 0 : logger => cp_get_default_logger()
1028 :
1029 0 : CALL force_env_get(force_env, qs_env=qs_env)
1030 0 : CALL get_qs_env(qs_env, dft_control=dft_control, subsys=subsys)
1031 0 : ncontacts = SIZE(negf_control%contacts)
1032 0 : nspins = dft_control%nspins
1033 0 : 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 0 : exist_all = .TRUE.
1038 0 : IF (para_env%is_source()) THEN
1039 0 : CALL negf_restart_file_name(filename_s, exist, negf_section, logger, s=.TRUE.)
1040 0 : IF (.NOT. exist) THEN
1041 : CALL cp_warn(__LOCATION__, &
1042 : "User requested to read the overlap matrix from the file named: "// &
1043 0 : TRIM(filename_s)//". This file does not exist. The file will be created.")
1044 0 : exist_all = .FALSE.
1045 : END IF
1046 0 : IF (nspins == 1) THEN
1047 0 : CALL negf_restart_file_name(filename_h_1, exist, negf_section, logger, ispin=0, h=.TRUE.)
1048 0 : IF (.NOT. exist) THEN
1049 : CALL cp_warn(__LOCATION__, &
1050 : "User requested to read the Hamiltonian matrix from the file named: "// &
1051 0 : TRIM(filename_h_1)//". This file does not exist. The file will be created.")
1052 0 : exist_all = .FALSE.
1053 : END IF
1054 : END IF
1055 0 : IF (nspins == 2) THEN
1056 0 : CALL negf_restart_file_name(filename_h_1, exist, negf_section, logger, ispin=1, h=.TRUE.)
1057 0 : IF (.NOT. exist) THEN
1058 : CALL cp_warn(__LOCATION__, &
1059 : "User requested to read the Hamiltonian matrix from the file named: "// &
1060 0 : TRIM(filename_h_1)//". This file does not exist. The file will be created.")
1061 0 : exist_all = .FALSE.
1062 : END IF
1063 0 : CALL negf_restart_file_name(filename_h_2, exist, negf_section, logger, ispin=2, h=.TRUE.)
1064 0 : IF (.NOT. exist) THEN
1065 : CALL cp_warn(__LOCATION__, &
1066 : "User requested to read the Hamiltonian matrix from the file named: "// &
1067 0 : TRIM(filename_h_2)//". This file does not exist. The file will be created.")
1068 0 : exist_all = .FALSE.
1069 : END IF
1070 : END IF
1071 0 : DO icontact = 1, ncontacts
1072 0 : CALL negf_restart_file_name(filename_sc(icontact), exist, negf_section, logger, icontact=icontact, sc=.TRUE.)
1073 0 : IF (.NOT. exist) THEN
1074 : CALL cp_warn(__LOCATION__, &
1075 : "User requested to read the overlap matrix from the file named: "// &
1076 0 : TRIM(filename_sc(icontact))//". This file does not exist. The file will be created.")
1077 0 : exist_all = .FALSE.
1078 : END IF
1079 0 : IF (nspins == 1) THEN
1080 : CALL negf_restart_file_name(filename_hc_1(icontact), exist, negf_section, logger, icontact=icontact, &
1081 0 : ispin=0, hc=.TRUE.)
1082 0 : IF (.NOT. exist) THEN
1083 : CALL cp_warn(__LOCATION__, &
1084 : "User requested to read the Hamiltonian matrix from the file named: "// &
1085 0 : TRIM(filename_hc_1(icontact))//". This file does not exist. The file will be created.")
1086 0 : exist_all = .FALSE.
1087 : END IF
1088 : END IF
1089 0 : IF (nspins == 2) THEN
1090 : CALL negf_restart_file_name(filename_hc_1(icontact), exist, negf_section, logger, icontact=icontact, &
1091 0 : ispin=1, hc=.TRUE.)
1092 0 : IF (.NOT. exist) THEN
1093 : CALL cp_warn(__LOCATION__, &
1094 : "User requested to read the Hamiltonian matrix from the file named: "// &
1095 0 : TRIM(filename_hc_1(icontact))//". This file does not exist. The file will be created.")
1096 0 : exist_all = .FALSE.
1097 : END IF
1098 : CALL negf_restart_file_name(filename_hc_2(icontact), exist, negf_section, logger, icontact=icontact, &
1099 0 : ispin=2, hc=.TRUE.)
1100 0 : IF (.NOT. exist) THEN
1101 : CALL cp_warn(__LOCATION__, &
1102 : "User requested to read the Hamiltonian matrix from the file named: "// &
1103 0 : TRIM(filename_hc_2(icontact))//". This file does not exist. The file will be created.")
1104 0 : exist_all = .FALSE.
1105 : END IF
1106 : END IF
1107 : END DO
1108 : END IF
1109 0 : CALL para_env%bcast(exist_all)
1110 :
1111 0 : IF (exist_all) THEN
1112 :
1113 0 : negf_control%is_restart = .TRUE.
1114 :
1115 0 : 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 0 : 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 0 : file_position="REWIND", unit_number=print_unit)
1122 0 : READ (print_unit, *) nrow_s, ncol_s
1123 0 : CALL close_file(print_unit)
1124 : END IF
1125 0 : CALL para_env%bcast(nrow_s)
1126 0 : CALL para_env%bcast(ncol_s)
1127 0 : NULLIFY (fm_struct)
1128 0 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_s, ncol_global=ncol_s, context=sub_env%blacs_env)
1129 0 : ALLOCATE (negf_env%s_s)
1130 0 : CALL cp_fm_create(negf_env%s_s, fm_struct)
1131 0 : ALLOCATE (negf_env%h_s(nspins))
1132 0 : DO ispin = 1, nspins
1133 0 : CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
1134 : END DO
1135 0 : CALL cp_fm_struct_release(fm_struct)
1136 0 : ALLOCATE (negf_env%s_sc(ncontacts))
1137 0 : ALLOCATE (negf_env%h_sc(nspins, ncontacts))
1138 0 : DO icontact = 1, ncontacts
1139 0 : 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 0 : file_position="REWIND", unit_number=print_unit)
1143 0 : READ (print_unit, *) nrow_sc, ncol_sc
1144 0 : CALL close_file(print_unit)
1145 : END IF
1146 0 : CALL para_env%bcast(nrow_sc)
1147 0 : CALL para_env%bcast(ncol_sc)
1148 0 : NULLIFY (fm_struct)
1149 0 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_sc, ncol_global=ncol_sc, context=sub_env%blacs_env)
1150 0 : CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
1151 0 : DO ispin = 1, nspins
1152 0 : CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
1153 : END DO
1154 0 : CALL cp_fm_struct_release(fm_struct)
1155 : END DO
1156 :
1157 0 : ALLOCATE (target_m(nrow_s, ncol_s))
1158 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_s, target_m)
1159 0 : CALL para_env%bcast(target_m)
1160 0 : CALL cp_fm_set_submatrix(negf_env%s_s, target_m)
1161 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_s is read from "//TRIM(filename_s)
1162 0 : IF (nspins == 1) THEN
1163 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_1, target_m)
1164 0 : CALL para_env%bcast(target_m)
1165 0 : CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m)
1166 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//TRIM(filename_h_1)
1167 : END IF
1168 0 : IF (nspins == 2) THEN
1169 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_1, target_m)
1170 0 : CALL para_env%bcast(target_m)
1171 0 : CALL cp_fm_set_submatrix(negf_env%h_s(1), target_m)
1172 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is read from "//TRIM(filename_h_1)//" for spin 1"
1173 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_h_2, target_m)
1174 0 : CALL para_env%bcast(target_m)
1175 0 : CALL cp_fm_set_submatrix(negf_env%h_s(2), target_m)
1176 0 : 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 0 : DEALLOCATE (target_m)
1179 :
1180 0 : DO icontact = 1, ncontacts
1181 0 : ALLOCATE (target_m(nrow_s, ncol_sc))
1182 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_sc(icontact), target_m)
1183 0 : CALL para_env%bcast(target_m)
1184 0 : CALL cp_fm_set_submatrix(negf_env%s_sc(icontact), target_m)
1185 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "S_sc is read from "//TRIM(filename_sc(icontact))
1186 0 : IF (nspins == 1) THEN
1187 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_1(icontact), target_m)
1188 0 : CALL para_env%bcast(target_m)
1189 0 : CALL cp_fm_set_submatrix(negf_env%h_sc(1, icontact), target_m)
1190 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//TRIM(filename_hc_1(icontact))
1191 : END IF
1192 0 : IF (nspins == 2) THEN
1193 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_1(icontact), target_m)
1194 0 : CALL para_env%bcast(target_m)
1195 0 : CALL cp_fm_set_submatrix(negf_env%h_sc(1, icontact), target_m)
1196 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is read from "//TRIM(filename_hc_1(icontact))//" for spin 1"
1197 0 : IF (para_env%is_source()) CALL negf_read_matrix_from_file(filename_hc_2(icontact), target_m)
1198 0 : CALL para_env%bcast(target_m)
1199 0 : CALL cp_fm_set_submatrix(negf_env%h_sc(2, icontact), target_m)
1200 0 : 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 0 : DEALLOCATE (target_m)
1203 :
1204 : END DO
1205 :
1206 : ELSE
1207 :
1208 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') &
1209 0 : "Some restart files do not exist. ALL restart files will be recalculated!"
1210 :
1211 0 : IF (.NOT. is_dft_entire) CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.)
1212 : ! extract device-related matrix blocks
1213 0 : CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
1214 0 : is_dft_entire = .TRUE.
1215 :
1216 0 : CALL cp_fm_get_info(negf_env%s_s, nrow_global=nrow_s, ncol_global=ncol_s)
1217 0 : ALLOCATE (target_m(nrow_s, ncol_s))
1218 0 : CALL cp_fm_get_submatrix(negf_env%s_s, target_m)
1219 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_s, target_m)
1220 0 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,A)') "S_s is saved to "//TRIM(filename_s)
1221 0 : IF (nspins == 1) THEN
1222 0 : CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m)
1223 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_1, target_m)
1224 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//TRIM(filename_h_1)
1225 : END IF
1226 0 : IF (nspins == 2) THEN
1227 0 : CALL cp_fm_get_submatrix(negf_env%h_s(1), target_m)
1228 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_1, target_m)
1229 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_s is saved to "//TRIM(filename_h_1)//" for spin 1"
1230 0 : CALL cp_fm_get_submatrix(negf_env%h_s(2), target_m)
1231 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_h_2, target_m)
1232 0 : 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 0 : DEALLOCATE (target_m)
1235 :
1236 0 : DO icontact = 1, ncontacts
1237 0 : CALL cp_fm_get_info(negf_env%contacts(icontact)%s_00, nrow_global=nrow_sc, ncol_global=ncol_sc)
1238 0 : ALLOCATE (target_m(nrow_s, ncol_sc))
1239 0 : CALL cp_fm_get_submatrix(negf_env%s_sc(icontact), target_m)
1240 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_sc(icontact), target_m)
1241 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A,I3)') &
1242 0 : "S_sc is saved to "//TRIM(filename_sc(icontact))//" for contact", icontact
1243 0 : IF (nspins == 1) THEN
1244 0 : CALL cp_fm_get_submatrix(negf_env%h_sc(1, icontact), target_m)
1245 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_1(icontact), target_m)
1246 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//TRIM(filename_hc_1(icontact))
1247 : END IF
1248 0 : IF (nspins == 2) THEN
1249 0 : CALL cp_fm_get_submatrix(negf_env%h_sc(1, icontact), target_m)
1250 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_1(icontact), target_m)
1251 0 : IF (log_unit > 0) WRITE (log_unit, '(T2,A)') "H_sc is saved to "//TRIM(filename_hc_1(icontact))//" for spin 1"
1252 0 : CALL cp_fm_get_submatrix(negf_env%h_sc(2, icontact), target_m)
1253 0 : IF (para_env%is_source()) CALL negf_print_matrix_to_file(filename_hc_2(icontact), target_m)
1254 0 : 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 0 : DEALLOCATE (target_m)
1257 : END DO
1258 :
1259 0 : negf_control%write_common_restart_file = .TRUE.
1260 :
1261 : END IF
1262 :
1263 0 : DEALLOCATE (filename_sc, filename_hc_1, filename_hc_2)
1264 0 : CALL timestop(handle)
1265 0 : 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 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
1298 :
1299 4 : 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 4 : subsys=subsys)
1308 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
1309 :
1310 4 : IF (do_kpoints) THEN
1311 : CALL cp_abort(__LOCATION__, &
1312 0 : "K-points in device region have not been implemented yet.")
1313 : END IF
1314 :
1315 4 : ncontacts = SIZE(negf_control%contacts)
1316 4 : nspins = dft_control%nspins
1317 :
1318 4 : NULLIFY (fm_struct)
1319 4 : nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening)
1320 :
1321 : ! ++ create matrices: h_s, s_s
1322 4 : NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
1323 16 : ALLOCATE (negf_env%h_s(nspins))
1324 :
1325 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
1326 4 : ALLOCATE (negf_env%s_s)
1327 4 : CALL cp_fm_create(negf_env%s_s, fm_struct)
1328 8 : DO ispin = 1, nspins
1329 8 : CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
1330 : END DO
1331 4 : ALLOCATE (negf_env%v_hartree_s)
1332 4 : CALL cp_fm_create(negf_env%v_hartree_s, fm_struct)
1333 4 : CALL cp_fm_struct_release(fm_struct)
1334 :
1335 : ! ++ create matrices: h_sc, s_sc
1336 48 : ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
1337 12 : DO icontact = 1, ncontacts
1338 8 : nao_c = number_of_atomic_orbitals(subsys, negf_env%contacts(icontact)%atomlist_cell0)
1339 8 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
1340 :
1341 8 : CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
1342 :
1343 16 : DO ispin = 1, nspins
1344 16 : CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
1345 : END DO
1346 :
1347 12 : CALL cp_fm_struct_release(fm_struct)
1348 : END DO
1349 :
1350 : ! extract matrices: h_s, s_s
1351 8 : 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 8 : 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 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
1366 :
1367 : ! v_hartree_s
1368 4 : NULLIFY (hmat%matrix)
1369 4 : CALL dbcsr_init_p(hmat%matrix)
1370 4 : CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
1371 4 : CALL dbcsr_set(hmat%matrix, 0.0_dp)
1372 :
1373 4 : CALL pw_pool%create_pw(v_hartree)
1374 4 : 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 4 : calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE.)
1378 :
1379 4 : 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 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
1387 :
1388 4 : CALL dbcsr_deallocate_matrix(hmat%matrix)
1389 :
1390 : ! extract matrices: h_sc, s_sc
1391 12 : DO icontact = 1, ncontacts
1392 16 : 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 16 : 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 12 : do_upper_diag=.TRUE., do_lower=.TRUE.)
1407 : END DO
1408 : END IF
1409 :
1410 4 : CALL timestop(handle)
1411 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
1437 :
1438 4 : ncontacts = SIZE(contact_env)
1439 4 : CPASSERT(SIZE(contact_control) == ncontacts)
1440 4 : CPASSERT(ncontacts == 2)
1441 :
1442 16 : dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
1443 4 : v1 = contact_control(1)%v_external
1444 4 : v2 = contact_control(2)%v_external
1445 :
1446 4 : lx = v_hartree%pw_grid%bounds_local(1, 1)
1447 4 : ux = v_hartree%pw_grid%bounds_local(2, 1)
1448 4 : ly = v_hartree%pw_grid%bounds_local(1, 2)
1449 4 : uy = v_hartree%pw_grid%bounds_local(2, 2)
1450 4 : lz = v_hartree%pw_grid%bounds_local(1, 3)
1451 4 : uz = v_hartree%pw_grid%bounds_local(2, 3)
1452 :
1453 4 : dx = v_hartree%pw_grid%npts(1)/2
1454 4 : dy = v_hartree%pw_grid%npts(2)/2
1455 4 : dz = v_hartree%pw_grid%npts(3)/2
1456 :
1457 4 : dvol = v_hartree%pw_grid%dvol
1458 :
1459 1028 : DO iz = lz, uz
1460 1024 : point_indices(3) = REAL(iz + dz, kind=dp)
1461 41988 : DO iy = ly, uy
1462 40960 : point_indices(2) = REAL(iy + dy, kind=dp)
1463 :
1464 861184 : DO ix = lx, ux
1465 819200 : point_indices(1) = REAL(ix + dx, kind=dp)
1466 10649600 : point_coord(:) = MATMUL(v_hartree%pw_grid%dh, point_indices)
1467 :
1468 3276800 : vector = point_coord - contact_env(1)%origin_bias
1469 819200 : proj = projection_on_direction_vector(vector, dirvector_bias)
1470 819200 : 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 320000 : IF (proj < 0.0_dp) THEN
1475 : proj = 0.0_dp
1476 320000 : ELSE IF (proj > 1.0_dp) THEN
1477 0 : proj = 1.0_dp
1478 : END IF
1479 320000 : pot = v1 + (v2 - v1)*proj
1480 : ELSE
1481 790400 : pot = 0.0_dp
1482 790400 : DO icontact = 1, ncontacts
1483 3046400 : vector = point_coord - contact_env(icontact)%origin_bias
1484 761600 : proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias)
1485 :
1486 790400 : IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
1487 470400 : pot = contact_control(icontact)%v_external
1488 470400 : EXIT
1489 : END IF
1490 : END DO
1491 : END IF
1492 :
1493 860160 : v_hartree%array(ix, iy, iz) = pot*dvol
1494 : END DO
1495 : END DO
1496 : END DO
1497 :
1498 4 : CALL timestop(handle)
1499 4 : 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 8 : 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 8 : CALL qs_subsys_get(subsys_contact, cell=cell)
1521 8 : CALL real_to_scaled(scaled, direction_vector, cell)
1522 :
1523 8 : naxes = 0
1524 8 : direction_axis = 0 ! initialize to make GCC<=6 happy
1525 :
1526 32 : DO i = 1, 3
1527 32 : IF (ABS(scaled(i)) > eps_geometry) THEN
1528 8 : IF (scaled(i) > 0.0_dp) THEN
1529 : direction_axis = i
1530 : ELSE
1531 4 : direction_axis = -i
1532 : END IF
1533 8 : 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 8 : IF (naxes /= 1) direction_axis = 0
1539 8 : 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 4 : 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 4 : REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
1561 4 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
1562 : TYPE(kpoint_type), POINTER :: kpoints
1563 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1564 4 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
1565 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
1566 :
1567 4 : CALL timeset(routineN, handle)
1568 4 : my_homo_energy = 0.0_dp
1569 :
1570 4 : CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
1571 :
1572 4 : IF (do_kpoints) THEN
1573 4 : 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 4 : IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point) THEN
1577 2 : kplocal = kp_range(2) - kp_range(1) + 1
1578 :
1579 2 : DO ikpgr = 1, kplocal
1580 2 : CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
1581 :
1582 2 : IF (ikpoint == gamma_point) THEN
1583 : ! mos_kp(component, spin), where component = 1 (real), or 2 (imaginary)
1584 2 : CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues) ! mu=fermi_level
1585 :
1586 2 : my_homo_energy = eigenvalues(homo)
1587 2 : EXIT
1588 : END IF
1589 : END DO
1590 : END IF
1591 :
1592 4 : 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 0 : "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 0 : nspins = SIZE(mos)
1605 :
1606 0 : spin_loop: DO ispin = 1, nspins
1607 0 : CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
1608 :
1609 0 : DO imo = nmo, 1, -1
1610 0 : IF (eigenvalues(imo) /= 0.0_dp) EXIT spin_loop
1611 : END DO
1612 : END DO spin_loop
1613 :
1614 0 : IF (imo == 0) THEN
1615 0 : CPABORT("Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
1616 : END IF
1617 :
1618 0 : my_homo_energy = eigenvalues(homo)
1619 : END IF
1620 :
1621 4 : homo_energy = my_homo_energy
1622 4 : CALL timestop(handle)
1623 4 : 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 4 : 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 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1658 :
1659 4 : CALL timeset(routineN, handle)
1660 4 : CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1661 :
1662 4 : natoms_bulk = SIZE(atomlist_bulk)
1663 4 : CPASSERT(SIZE(atom_map, 1) == natoms_bulk)
1664 4 : direction_axis_abs = ABS(direction_axis)
1665 :
1666 : ! looking for the nearest atom from the scattering region
1667 4 : proj_min = 1.0_dp
1668 4 : atom_min = 1
1669 36 : DO iatom = 1, natoms_bulk
1670 128 : vector = particle_set(atomlist_bulk(iatom))%r - origin
1671 32 : proj = projection_on_direction_vector(vector, direction_vector)
1672 :
1673 36 : IF (proj < proj_min) THEN
1674 16 : proj_min = proj
1675 16 : atom_min = iatom
1676 : END IF
1677 : END DO
1678 :
1679 4 : dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1680 :
1681 4 : natoms_cell0 = 0
1682 36 : DO iatom = 1, natoms_bulk
1683 32 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) &
1684 20 : natoms_cell0 = natoms_cell0 + 1
1685 : END DO
1686 :
1687 12 : ALLOCATE (atomlist_cell0(natoms_cell0))
1688 40 : ALLOCATE (atom_map_cell0(natoms_cell0))
1689 :
1690 4 : natoms_cell0 = 0
1691 36 : DO iatom = 1, natoms_bulk
1692 36 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
1693 16 : natoms_cell0 = natoms_cell0 + 1
1694 16 : atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
1695 16 : atom_map_cell0(natoms_cell0) = atom_map(iatom)
1696 : END IF
1697 : END DO
1698 :
1699 4 : CALL timestop(handle)
1700 4 : 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 4 : 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 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1738 :
1739 4 : CALL timeset(routineN, handle)
1740 4 : CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1741 :
1742 4 : natoms_bulk = SIZE(atomlist_bulk)
1743 4 : CPASSERT(SIZE(atom_map, 1) == natoms_bulk)
1744 4 : direction_axis_abs = ABS(direction_axis)
1745 4 : offset = SIGN(1, direction_axis)
1746 :
1747 : ! looking for the nearest atom from the scattering region
1748 4 : proj_min = 1.0_dp
1749 4 : atom_min = 1
1750 36 : DO iatom = 1, natoms_bulk
1751 128 : vector = particle_set(atomlist_bulk(iatom))%r - origin
1752 32 : proj = projection_on_direction_vector(vector, direction_vector)
1753 :
1754 36 : IF (proj < proj_min) THEN
1755 16 : proj_min = proj
1756 16 : atom_min = iatom
1757 : END IF
1758 : END DO
1759 :
1760 4 : dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1761 :
1762 4 : natoms_cell1 = 0
1763 36 : DO iatom = 1, natoms_bulk
1764 32 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) &
1765 20 : natoms_cell1 = natoms_cell1 + 1
1766 : END DO
1767 :
1768 12 : ALLOCATE (atomlist_cell1(natoms_cell1))
1769 40 : ALLOCATE (atom_map_cell1(natoms_cell1))
1770 :
1771 4 : natoms_cell1 = 0
1772 36 : DO iatom = 1, natoms_bulk
1773 36 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
1774 16 : natoms_cell1 = natoms_cell1 + 1
1775 16 : atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
1776 16 : atom_map_cell1(natoms_cell1) = atom_map(iatom)
1777 16 : atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
1778 : END IF
1779 : END DO
1780 :
1781 4 : CALL timestop(handle)
1782 4 : 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 4 : 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 4 : CALL timeset(routineN, handle)
1798 :
1799 4 : IF (ALLOCATED(negf_env%contacts)) THEN
1800 12 : DO icontact = SIZE(negf_env%contacts), 1, -1
1801 12 : CALL negf_env_contact_release(negf_env%contacts(icontact))
1802 : END DO
1803 :
1804 12 : DEALLOCATE (negf_env%contacts)
1805 : END IF
1806 :
1807 : ! h_s
1808 4 : CALL cp_fm_release(negf_env%h_s)
1809 :
1810 : ! h_sc
1811 4 : CALL cp_fm_release(negf_env%h_sc)
1812 :
1813 : ! s_s
1814 4 : IF (ASSOCIATED(negf_env%s_s)) THEN
1815 4 : CALL cp_fm_release(negf_env%s_s)
1816 4 : DEALLOCATE (negf_env%s_s)
1817 : NULLIFY (negf_env%s_s)
1818 : END IF
1819 :
1820 : ! s_sc
1821 4 : CALL cp_fm_release(negf_env%s_sc)
1822 :
1823 : ! v_hartree_s
1824 4 : IF (ASSOCIATED(negf_env%v_hartree_s)) THEN
1825 4 : CALL cp_fm_release(negf_env%v_hartree_s)
1826 4 : DEALLOCATE (negf_env%v_hartree_s)
1827 : NULLIFY (negf_env%v_hartree_s)
1828 : END IF
1829 :
1830 : ! density mixing
1831 4 : IF (ASSOCIATED(negf_env%mixing_storage)) THEN
1832 4 : CALL mixing_storage_release(negf_env%mixing_storage)
1833 4 : DEALLOCATE (negf_env%mixing_storage)
1834 : END IF
1835 :
1836 4 : CALL timestop(handle)
1837 4 : 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 8 : 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 8 : CALL timeset(routineN, handle)
1851 :
1852 : ! h_00
1853 8 : CALL cp_fm_release(contact_env%h_00)
1854 :
1855 : ! h_01
1856 8 : CALL cp_fm_release(contact_env%h_01)
1857 :
1858 : ! rho_00
1859 8 : CALL cp_fm_release(contact_env%rho_00)
1860 :
1861 : ! rho_01
1862 8 : CALL cp_fm_release(contact_env%rho_01)
1863 :
1864 : ! s_00
1865 8 : IF (ASSOCIATED(contact_env%s_00)) THEN
1866 8 : CALL cp_fm_release(contact_env%s_00)
1867 8 : DEALLOCATE (contact_env%s_00)
1868 : NULLIFY (contact_env%s_00)
1869 : END IF
1870 :
1871 : ! s_01
1872 8 : IF (ASSOCIATED(contact_env%s_01)) THEN
1873 8 : CALL cp_fm_release(contact_env%s_01)
1874 8 : DEALLOCATE (contact_env%s_01)
1875 : NULLIFY (contact_env%s_01)
1876 : END IF
1877 :
1878 8 : IF (ALLOCATED(contact_env%atomlist_cell0)) DEALLOCATE (contact_env%atomlist_cell0)
1879 8 : IF (ALLOCATED(contact_env%atomlist_cell1)) DEALLOCATE (contact_env%atomlist_cell1)
1880 8 : IF (ALLOCATED(contact_env%atom_map_cell0)) DEALLOCATE (contact_env%atom_map_cell0)
1881 :
1882 8 : CALL timestop(handle)
1883 8 : END SUBROUTINE negf_env_contact_release
1884 :
1885 0 : END MODULE negf_env_types
|
