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