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