Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_fm_struct, ONLY: cp_fm_struct_create,&
22 : cp_fm_struct_release,&
23 : cp_fm_struct_type
24 : USE cp_fm_types, ONLY: cp_fm_create,&
25 : cp_fm_release,&
26 : cp_fm_type
27 : USE force_env_types, ONLY: force_env_get,&
28 : force_env_p_type,&
29 : force_env_type,&
30 : use_qs_force
31 : USE input_section_types, ONLY: section_vals_type,&
32 : section_vals_val_get
33 : USE kinds, ONLY: default_string_length,&
34 : dp
35 : USE kpoint_types, ONLY: get_kpoint_env,&
36 : get_kpoint_info,&
37 : kpoint_env_p_type,&
38 : kpoint_type
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE negf_atom_map, ONLY: negf_atom_map_type,&
41 : negf_map_atomic_indices
42 : USE negf_control_types, ONLY: negf_control_contact_type,&
43 : negf_control_type
44 : USE negf_matrix_utils, ONLY: invert_cell_to_index,&
45 : negf_copy_contact_matrix,&
46 : negf_copy_sym_dbcsr_to_fm_submat,&
47 : number_of_atomic_orbitals
48 : USE negf_subgroup_types, ONLY: negf_subgroup_env_type
49 : USE negf_vectors, ONLY: contact_direction_vector,&
50 : projection_on_direction_vector
51 : USE particle_types, ONLY: particle_type
52 : USE pw_env_types, ONLY: pw_env_get,&
53 : pw_env_type
54 : USE pw_pool_types, ONLY: pw_pool_type
55 : USE pw_types, ONLY: pw_r3d_rs_type
56 : USE qs_density_mixing_types, ONLY: mixing_storage_create,&
57 : mixing_storage_release,&
58 : mixing_storage_type
59 : USE qs_energy, ONLY: qs_energies
60 : USE qs_energy_init, ONLY: qs_energies_init
61 : USE qs_environment_types, ONLY: get_qs_env,&
62 : qs_environment_type
63 : USE qs_integrate_potential, ONLY: integrate_v_rspace
64 : USE qs_mo_types, ONLY: get_mo_set,&
65 : mo_set_type
66 : USE qs_rho_types, ONLY: qs_rho_get,&
67 : qs_rho_type
68 : USE qs_subsys_types, ONLY: qs_subsys_get,&
69 : qs_subsys_type
70 : #include "./base/base_uses.f90"
71 :
72 : IMPLICIT NONE
73 : PRIVATE
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_env_types'
76 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
77 :
78 : PUBLIC :: negf_env_type, negf_env_contact_type
79 : PUBLIC :: negf_env_create, negf_env_release
80 :
81 : ! **************************************************************************************************
82 : !> \brief Contact-specific NEGF environment.
83 : !> \author Sergey Chulkov
84 : ! **************************************************************************************************
85 : TYPE negf_env_contact_type
86 : REAL(kind=dp), DIMENSION(3) :: direction_vector = -1.0_dp, origin = -1.0_dp
87 : REAL(kind=dp), DIMENSION(3) :: direction_vector_bias = -1.0_dp, origin_bias = -1.0_dp
88 : !> an axis towards the secondary contact unit cell which coincides with the transport direction
89 : !> 0 (undefined), 1 (+x), 2 (+y), 3 (+z), -1 (-x), -2 (-y), -3 (-z)
90 : INTEGER :: direction_axis = -1
91 : !> atoms belonging to a primary contact unit cell
92 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_cell0
93 : !> atoms belonging to a secondary contact unit cell (will be removed one day ...)
94 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_cell1
95 : !> list of equivalent atoms in an appropriate contact force environment
96 : TYPE(negf_atom_map_type), ALLOCATABLE, &
97 : DIMENSION(:) :: atom_map_cell0, atom_map_cell1
98 : !> energy of the HOMO
99 : REAL(kind=dp) :: homo_energy = -1.0_dp
100 : !> diagonal (h_00) and off-diagonal (h_01) blocks of the contact Kohn-Sham matrix ([number_of_spins]).
101 : !> The matrix h_01 is of the shape [nao_cell0 x nao_cell1]
102 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h_00, h_01
103 : !> diagonal and off-diagonal blocks of the density matrix
104 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: rho_00, rho_01
105 : !> diagonal and off-diagonal blocks of the overlap matrix
106 : TYPE(cp_fm_type), POINTER :: s_00 => null(), s_01 => null()
107 : END TYPE negf_env_contact_type
108 :
109 : ! **************************************************************************************************
110 : !> \brief NEGF environment.
111 : !> \author Sergey Chulkov
112 : ! **************************************************************************************************
113 : TYPE negf_env_type
114 : !> contact-specific NEGF environments
115 : TYPE(negf_env_contact_type), ALLOCATABLE, &
116 : DIMENSION(:) :: contacts
117 : !> Kohn-Sham matrix of the scattering region
118 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: h_s
119 : !> Kohn-Sham matrix of the scattering region -- contact interface ([nspins, ncontacts])
120 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: h_sc
121 : !> overlap matrix of the scattering region
122 : TYPE(cp_fm_type), POINTER :: s_s => null()
123 : !> an external Hartree potential in atomic basis set representation
124 : TYPE(cp_fm_type), POINTER :: v_hartree_s => null()
125 : !> overlap matrix of the scattering region -- contact interface for every contact ([ncontacts])
126 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: s_sc
127 : !> structure needed for density mixing
128 : TYPE(mixing_storage_type), POINTER :: mixing_storage => NULL()
129 : !> density mixing method
130 : INTEGER :: mixing_method = -1
131 : END TYPE negf_env_type
132 :
133 : ! **************************************************************************************************
134 : !> \brief Allocatable list of the type 'negf_atom_map_type'.
135 : !> \author Sergey Chulkov
136 : ! **************************************************************************************************
137 : TYPE negf_atom_map_contact_type
138 : TYPE(negf_atom_map_type), ALLOCATABLE, DIMENSION(:) :: atom_map
139 : END TYPE negf_atom_map_contact_type
140 :
141 : CONTAINS
142 :
143 : ! **************************************************************************************************
144 : !> \brief Create a new NEGF environment and compute the relevant Kohn-Sham matrices.
145 : !> \param negf_env NEGF environment to create
146 : !> \param sub_env NEGF parallel (sub)group environment
147 : !> \param negf_control NEGF control
148 : !> \param force_env the primary force environment
149 : !> \param negf_mixing_section pointer to a mixing section within the NEGF input section
150 : !> \param log_unit output unit number
151 : !> \par History
152 : !> * 01.2017 created [Sergey Chulkov]
153 : ! **************************************************************************************************
154 4 : SUBROUTINE negf_env_create(negf_env, sub_env, negf_control, force_env, negf_mixing_section, log_unit)
155 : TYPE(negf_env_type), INTENT(inout) :: negf_env
156 : TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
157 : TYPE(negf_control_type), POINTER :: negf_control
158 : TYPE(force_env_type), POINTER :: force_env
159 : TYPE(section_vals_type), POINTER :: negf_mixing_section
160 : INTEGER, INTENT(in) :: log_unit
161 :
162 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_create'
163 :
164 : CHARACTER(len=default_string_length) :: contact_str, force_env_str, &
165 : n_force_env_str
166 : INTEGER :: handle, icontact, in_use, n_force_env, &
167 : ncontacts
168 : LOGICAL :: do_kpoints
169 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
170 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
171 : TYPE(dft_control_type), POINTER :: dft_control
172 4 : TYPE(force_env_p_type), DIMENSION(:), POINTER :: sub_force_env
173 : TYPE(negf_atom_map_contact_type), ALLOCATABLE, &
174 4 : DIMENSION(:) :: map_contact
175 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_rspace
176 : TYPE(qs_environment_type), POINTER :: qs_env, qs_env_contact
177 : TYPE(qs_subsys_type), POINTER :: subsys, subsys_contact
178 :
179 4 : CALL timeset(routineN, handle)
180 :
181 : ! ensure we have Quickstep enabled for all force_env
182 4 : NULLIFY (sub_force_env)
183 4 : CALL force_env_get(force_env, in_use=in_use, qs_env=qs_env, sub_force_env=sub_force_env)
184 :
185 4 : IF (ASSOCIATED(sub_force_env)) THEN
186 2 : n_force_env = SIZE(sub_force_env)
187 : ELSE
188 2 : n_force_env = 0
189 : END IF
190 :
191 4 : IF (in_use == use_qs_force) THEN
192 8 : DO icontact = 1, n_force_env
193 4 : CALL force_env_get(sub_force_env(icontact)%force_env, in_use=in_use)
194 8 : IF (in_use /= use_qs_force) EXIT
195 : END DO
196 : END IF
197 :
198 4 : IF (in_use /= use_qs_force) THEN
199 0 : CPABORT("Quickstep is required for NEGF run.")
200 : END IF
201 :
202 : ! check that all mentioned FORCE_EVAL sections are actually present
203 4 : ncontacts = SIZE(negf_control%contacts)
204 :
205 12 : DO icontact = 1, ncontacts
206 12 : IF (negf_control%contacts(icontact)%force_env_index > n_force_env) THEN
207 0 : WRITE (contact_str, '(I11)') icontact
208 0 : WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
209 0 : WRITE (n_force_env_str, '(I11)') n_force_env
210 :
211 : CALL cp_abort(__LOCATION__, &
212 : "Contact number "//TRIM(ADJUSTL(contact_str))//" is linked with the FORCE_EVAL section number "// &
213 : TRIM(ADJUSTL(force_env_str))//", however only "//TRIM(ADJUSTL(n_force_env_str))// &
214 : " FORCE_EVAL sections have been found. Note that FORCE_EVAL sections are enumerated from 0"// &
215 0 : " and that the primary (0-th) section must contain all the atoms.")
216 : END IF
217 : END DO
218 :
219 : ! create basic matrices and neighbour lists for the primary force_env,
220 : ! so we know how matrix elements are actually distributed across CPUs.
221 4 : CALL qs_energies_init(qs_env, calc_forces=.FALSE.)
222 : CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, &
223 : matrix_s_kp=matrix_s_kp, matrix_ks_kp=matrix_ks_kp, &
224 4 : subsys=subsys, v_hartree_rspace=v_hartree_rspace)
225 :
226 4 : IF (do_kpoints) THEN
227 0 : CPABORT("k-points are currently not supported for device FORCE_EVAL")
228 : END IF
229 :
230 : ! stage 1: map the atoms between the device force_env and all contact force_env-s
231 80 : ALLOCATE (negf_env%contacts(ncontacts))
232 20 : ALLOCATE (map_contact(ncontacts))
233 :
234 12 : DO icontact = 1, ncontacts
235 12 : IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
236 4 : CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
237 4 : CALL get_qs_env(qs_env_contact, subsys=subsys_contact)
238 :
239 : CALL negf_env_contact_init_maps(contact_env=negf_env%contacts(icontact), &
240 : contact_control=negf_control%contacts(icontact), &
241 : atom_map=map_contact(icontact)%atom_map, &
242 : eps_geometry=negf_control%eps_geometry, &
243 : subsys_device=subsys, &
244 4 : subsys_contact=subsys_contact)
245 :
246 4 : IF (negf_env%contacts(icontact)%direction_axis == 0) THEN
247 0 : WRITE (contact_str, '(I11)') icontact
248 0 : WRITE (force_env_str, '(I11)') negf_control%contacts(icontact)%force_env_index
249 : CALL cp_abort(__LOCATION__, &
250 : "One lattice vector of the contact unit cell (FORCE_EVAL section "// &
251 : TRIM(ADJUSTL(force_env_str))//") must be parallel to the direction of the contact "// &
252 0 : TRIM(ADJUSTL(contact_str))//".")
253 : END IF
254 : END IF
255 : END DO
256 :
257 : ! stage 2: obtain relevant Kohn-Sham matrix blocks for each contact (separate bulk DFT calculation)
258 12 : DO icontact = 1, ncontacts
259 12 : IF (negf_control%contacts(icontact)%force_env_index > 0) THEN
260 4 : IF (log_unit > 0) &
261 2 : WRITE (log_unit, '(/,T2,A,T70,I11)') "NEGF| Construct the Kohn-Sham matrix for the contact", icontact
262 :
263 4 : CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env_contact)
264 :
265 4 : CALL qs_energies(qs_env_contact, consistent_energies=.FALSE., calc_forces=.FALSE.)
266 :
267 : CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, &
268 4 : qs_env_contact=qs_env_contact)
269 4 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
270 : END IF
271 : END DO
272 :
273 : ! obtain an initial KS-matrix for the scattering region
274 4 : IF (log_unit > 0) &
275 2 : WRITE (log_unit, '(/,T2,A,T70)') "NEGF| Construct the Kohn-Sham matrix for the entire system"
276 4 : CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.)
277 4 : IF (log_unit > 0) WRITE (log_unit, '(/,T2,79("-"))')
278 :
279 : ! *** obtain relevant Kohn-Sham matrix blocks for each contact with no separate FORCE_ENV ***
280 12 : DO icontact = 1, ncontacts
281 12 : IF (negf_control%contacts(icontact)%force_env_index <= 0) THEN
282 : CALL negf_env_contact_init_matrices_gamma(contact_env=negf_env%contacts(icontact), &
283 : contact_control=negf_control%contacts(icontact), &
284 : sub_env=sub_env, qs_env=qs_env, &
285 4 : eps_geometry=negf_control%eps_geometry)
286 : END IF
287 : END DO
288 :
289 : ! extract device-related matrix blocks
290 4 : CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
291 :
292 : ! electron density mixing;
293 : ! the input section below should be consistent with the subroutine create_negf_section()
294 4 : NULLIFY (negf_env%mixing_storage)
295 4 : CALL section_vals_val_get(negf_mixing_section, "METHOD", i_val=negf_env%mixing_method)
296 :
297 4 : CALL get_qs_env(qs_env, dft_control=dft_control)
298 16 : ALLOCATE (negf_env%mixing_storage)
299 : CALL mixing_storage_create(negf_env%mixing_storage, negf_mixing_section, &
300 4 : negf_env%mixing_method, dft_control%qs_control%cutoff)
301 :
302 4 : CALL timestop(handle)
303 16 : END SUBROUTINE negf_env_create
304 :
305 : ! **************************************************************************************************
306 : !> \brief Establish mapping between the primary and the contact force environments
307 : !> \param contact_env NEGF environment for the given contact (modified on exit)
308 : !> \param contact_control NEGF control
309 : !> \param atom_map atomic map
310 : !> \param eps_geometry accuracy in mapping atoms between different force environments
311 : !> \param subsys_device QuickStep subsystem of the device force environment
312 : !> \param subsys_contact QuickStep subsystem of the contact force environment
313 : !> \author Sergey Chulkov
314 : ! **************************************************************************************************
315 4 : SUBROUTINE negf_env_contact_init_maps(contact_env, contact_control, atom_map, &
316 : eps_geometry, subsys_device, subsys_contact)
317 : TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
318 : TYPE(negf_control_contact_type), INTENT(in) :: contact_control
319 : TYPE(negf_atom_map_type), ALLOCATABLE, &
320 : DIMENSION(:), INTENT(inout) :: atom_map
321 : REAL(kind=dp), INTENT(in) :: eps_geometry
322 : TYPE(qs_subsys_type), POINTER :: subsys_device, subsys_contact
323 :
324 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_contact_init_maps'
325 :
326 : INTEGER :: handle, natoms
327 :
328 4 : CALL timeset(routineN, handle)
329 :
330 : CALL contact_direction_vector(contact_env%origin, &
331 : contact_env%direction_vector, &
332 : contact_env%origin_bias, &
333 : contact_env%direction_vector_bias, &
334 : contact_control%atomlist_screening, &
335 : contact_control%atomlist_bulk, &
336 4 : subsys_device)
337 :
338 4 : contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys_contact, eps_geometry)
339 :
340 4 : IF (contact_env%direction_axis /= 0) THEN
341 4 : natoms = SIZE(contact_control%atomlist_bulk)
342 56 : ALLOCATE (atom_map(natoms))
343 :
344 : ! map atom listed in 'contact_control%atomlist_bulk' to the corresponding atom/cell replica from the contact force_env
345 : CALL negf_map_atomic_indices(atom_map=atom_map, &
346 : atom_list=contact_control%atomlist_bulk, &
347 : subsys_device=subsys_device, &
348 : subsys_contact=subsys_contact, &
349 4 : eps_geometry=eps_geometry)
350 :
351 : ! list atoms from 'contact_control%atomlist_bulk' which belong to
352 : ! the primary unit cell of the bulk region for the given contact
353 : CALL list_atoms_in_bulk_primary_unit_cell(atomlist_cell0=contact_env%atomlist_cell0, &
354 : atom_map_cell0=contact_env%atom_map_cell0, &
355 : atomlist_bulk=contact_control%atomlist_bulk, &
356 : atom_map=atom_map, &
357 : origin=contact_env%origin, &
358 : direction_vector=contact_env%direction_vector, &
359 : direction_axis=contact_env%direction_axis, &
360 4 : subsys_device=subsys_device)
361 :
362 : ! secondary unit cell
363 : CALL list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1=contact_env%atomlist_cell1, &
364 : atom_map_cell1=contact_env%atom_map_cell1, &
365 : atomlist_bulk=contact_control%atomlist_bulk, &
366 : atom_map=atom_map, &
367 : origin=contact_env%origin, &
368 : direction_vector=contact_env%direction_vector, &
369 : direction_axis=contact_env%direction_axis, &
370 4 : subsys_device=subsys_device)
371 : END IF
372 :
373 4 : CALL timestop(handle)
374 4 : END SUBROUTINE negf_env_contact_init_maps
375 :
376 : ! **************************************************************************************************
377 : !> \brief Extract relevant matrix blocks for the given contact.
378 : !> \param contact_env NEGF environment for the contact (modified on exit)
379 : !> \param sub_env NEGF parallel (sub)group environment
380 : !> \param qs_env_contact QuickStep environment for the contact force environment
381 : !> \par History
382 : !> * 10.2017 created [Sergey Chulkov]
383 : !> * 10.2025 The subroutine is essentially modified. New functionality of negf_copy_contact_matrix.
384 : !> [Dmitry Ryndyk]
385 : ! **************************************************************************************************
386 4 : SUBROUTINE negf_env_contact_init_matrices(contact_env, sub_env, qs_env_contact)
387 : TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
388 : TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
389 : TYPE(qs_environment_type), POINTER :: qs_env_contact
390 :
391 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_contact_init_matrices'
392 :
393 : INTEGER :: handle, iatom, ispin, nao, natoms, &
394 : nimages, nspins
395 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_list0, atom_list1
396 4 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: index_to_cell
397 4 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
398 : LOGICAL :: do_kpoints
399 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
400 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matkp
401 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
402 : TYPE(dft_control_type), POINTER :: dft_control
403 : TYPE(kpoint_type), POINTER :: kpoints
404 : TYPE(mp_para_env_type), POINTER :: para_env
405 : TYPE(qs_rho_type), POINTER :: rho_struct
406 : TYPE(qs_subsys_type), POINTER :: subsys
407 :
408 4 : CALL timeset(routineN, handle)
409 :
410 : CALL get_qs_env(qs_env_contact, &
411 : dft_control=dft_control, &
412 : do_kpoints=do_kpoints, &
413 : kpoints=kpoints, &
414 : matrix_ks_kp=matrix_ks_kp, &
415 : matrix_s_kp=matrix_s_kp, &
416 : para_env=para_env, &
417 : rho=rho_struct, &
418 4 : subsys=subsys)
419 4 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
420 :
421 4 : CALL negf_homo_energy_estimate(contact_env%homo_energy, qs_env_contact)
422 :
423 4 : natoms = SIZE(contact_env%atomlist_cell0)
424 12 : ALLOCATE (atom_list0(natoms))
425 20 : DO iatom = 1, natoms
426 16 : atom_list0(iatom) = contact_env%atom_map_cell0(iatom)%iatom
427 :
428 : ! with no k-points there is one-to-one correspondence between the primary unit cell
429 : ! of the contact force_env and the first contact unit cell of the device force_env
430 68 : IF (SUM(ABS(contact_env%atom_map_cell0(iatom)%cell(:))) > 0) THEN
431 0 : CPABORT("NEGF K-points are not currently supported")
432 : END IF
433 : END DO
434 :
435 4 : CPASSERT(SIZE(contact_env%atomlist_cell1) == natoms)
436 8 : ALLOCATE (atom_list1(natoms))
437 20 : DO iatom = 1, natoms
438 20 : atom_list1(iatom) = contact_env%atom_map_cell1(iatom)%iatom
439 : END DO
440 :
441 4 : nspins = dft_control%nspins
442 4 : nimages = dft_control%nimages
443 :
444 4 : IF (do_kpoints) THEN
445 4 : CALL get_kpoint_info(kpoints, cell_to_index=cell_to_index)
446 : ELSE
447 0 : ALLOCATE (cell_to_index(0:0, 0:0, 0:0))
448 0 : cell_to_index(0, 0, 0) = 1
449 : END IF
450 :
451 12 : ALLOCATE (index_to_cell(3, nimages))
452 4 : CALL invert_cell_to_index(cell_to_index, nimages, index_to_cell)
453 4 : IF (.NOT. do_kpoints) DEALLOCATE (cell_to_index)
454 :
455 4 : NULLIFY (fm_struct)
456 4 : nao = number_of_atomic_orbitals(subsys, atom_list0)
457 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=sub_env%blacs_env)
458 :
459 : ! ++ create matrices: s_00, s_01
460 4 : ALLOCATE (contact_env%s_00, contact_env%s_01)
461 4 : CALL cp_fm_create(contact_env%s_00, fm_struct)
462 4 : CALL cp_fm_create(contact_env%s_01, fm_struct)
463 :
464 : ! ++ create matrices: h_00, h_01, rho_00, rho_01
465 24 : ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
466 20 : ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
467 8 : DO ispin = 1, nspins
468 4 : CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
469 4 : CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
470 4 : CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
471 8 : CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
472 : END DO
473 :
474 4 : CALL cp_fm_struct_release(fm_struct)
475 :
476 : ! extract matrices: s_00, s_01
477 4 : matkp => matrix_s_kp(1, :)
478 : CALL negf_copy_contact_matrix(fm_cell0=contact_env%s_00, &
479 : fm_cell1=contact_env%s_01, &
480 : direction_axis=contact_env%direction_axis, &
481 : matrix_kp=matkp, &
482 : atom_list0=atom_list0, atom_list1=atom_list1, &
483 : subsys=subsys, mpi_comm_global=para_env, &
484 4 : kpoints=kpoints)
485 :
486 : ! extract matrices: h_00, h_01, rho_00, rho_01
487 8 : DO ispin = 1, nspins
488 4 : matkp => matrix_ks_kp(ispin, :)
489 : CALL negf_copy_contact_matrix(fm_cell0=contact_env%h_00(ispin), &
490 : fm_cell1=contact_env%h_01(ispin), &
491 : direction_axis=contact_env%direction_axis, &
492 : matrix_kp=matkp, &
493 : atom_list0=atom_list0, atom_list1=atom_list1, &
494 : subsys=subsys, mpi_comm_global=para_env, &
495 4 : kpoints=kpoints)
496 :
497 4 : matkp => rho_ao_kp(ispin, :)
498 : CALL negf_copy_contact_matrix(fm_cell0=contact_env%rho_00(ispin), &
499 : fm_cell1=contact_env%rho_01(ispin), &
500 : direction_axis=contact_env%direction_axis, &
501 : matrix_kp=matkp, &
502 : atom_list0=atom_list0, atom_list1=atom_list1, &
503 : subsys=subsys, mpi_comm_global=para_env, &
504 8 : kpoints=kpoints)
505 : END DO
506 :
507 4 : DEALLOCATE (atom_list0, atom_list1)
508 :
509 4 : CALL timestop(handle)
510 8 : END SUBROUTINE negf_env_contact_init_matrices
511 :
512 : ! **************************************************************************************************
513 : !> \brief Extract relevant matrix blocks for the given contact using the device's force environment.
514 : !> \param contact_env NEGF environment for the contact (modified on exit)
515 : !> \param contact_control NEGF control for the contact
516 : !> \param sub_env NEGF parallel (sub)group environment
517 : !> \param qs_env QuickStep environment for the device force environment
518 : !> \param eps_geometry accuracy in Cartesian coordinates
519 : !> \author Sergey Chulkov
520 : ! **************************************************************************************************
521 4 : SUBROUTINE negf_env_contact_init_matrices_gamma(contact_env, contact_control, sub_env, qs_env, eps_geometry)
522 : TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
523 : TYPE(negf_control_contact_type), INTENT(in) :: contact_control
524 : TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
525 : TYPE(qs_environment_type), POINTER :: qs_env
526 : REAL(kind=dp), INTENT(in) :: eps_geometry
527 :
528 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_contact_init_matrices_gamma'
529 :
530 : INTEGER :: handle, iatom, icell, ispin, nao_c, &
531 : nspins
532 : LOGICAL :: do_kpoints
533 : REAL(kind=dp), DIMENSION(2) :: r2_origin_cell
534 : REAL(kind=dp), DIMENSION(3) :: direction_vector, origin
535 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
536 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp, rho_ao_kp
537 : TYPE(dft_control_type), POINTER :: dft_control
538 : TYPE(mp_para_env_type), POINTER :: para_env
539 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
540 : TYPE(qs_rho_type), POINTER :: rho_struct
541 : TYPE(qs_subsys_type), POINTER :: subsys
542 :
543 4 : CALL timeset(routineN, handle)
544 :
545 : CALL get_qs_env(qs_env, &
546 : dft_control=dft_control, &
547 : do_kpoints=do_kpoints, &
548 : matrix_ks_kp=matrix_ks_kp, &
549 : matrix_s_kp=matrix_s_kp, &
550 : para_env=para_env, &
551 : rho=rho_struct, &
552 4 : subsys=subsys)
553 4 : CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp)
554 :
555 4 : IF (do_kpoints) THEN
556 : CALL cp_abort(__LOCATION__, &
557 0 : "K-points in device region have not been implemented yet.")
558 : END IF
559 :
560 4 : nspins = dft_control%nspins
561 :
562 4 : nao_c = number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(1)%vector)
563 4 : IF (number_of_atomic_orbitals(subsys, contact_control%atomlist_cell(2)%vector) /= nao_c) THEN
564 : CALL cp_abort(__LOCATION__, &
565 : "Primary and secondary bulk contact cells should be identical "// &
566 : "in terms of the number of atoms of each kind, and their basis sets. "// &
567 0 : "No single atom, however, can be shared between these two cells.")
568 : END IF
569 :
570 4 : contact_env%homo_energy = 0.0_dp
571 :
572 : CALL contact_direction_vector(contact_env%origin, &
573 : contact_env%direction_vector, &
574 : contact_env%origin_bias, &
575 : contact_env%direction_vector_bias, &
576 : contact_control%atomlist_screening, &
577 : contact_control%atomlist_bulk, &
578 4 : subsys)
579 :
580 4 : contact_env%direction_axis = contact_direction_axis(contact_env%direction_vector, subsys, eps_geometry)
581 :
582 : ! choose the primary and secondary contact unit cells
583 4 : CALL qs_subsys_get(subsys, particle_set=particle_set)
584 :
585 16 : origin = particle_set(contact_control%atomlist_screening(1))%r
586 16 : DO iatom = 2, SIZE(contact_control%atomlist_screening)
587 52 : origin = origin + particle_set(contact_control%atomlist_screening(iatom))%r
588 : END DO
589 16 : origin = origin/REAL(SIZE(contact_control%atomlist_screening), kind=dp)
590 :
591 12 : DO icell = 1, 2
592 32 : direction_vector = particle_set(contact_control%atomlist_cell(icell)%vector(1))%r
593 32 : DO iatom = 2, SIZE(contact_control%atomlist_cell(icell)%vector)
594 104 : direction_vector = direction_vector + particle_set(contact_control%atomlist_cell(icell)%vector(iatom))%r
595 : END DO
596 32 : direction_vector = direction_vector/REAL(SIZE(contact_control%atomlist_cell(icell)%vector), kind=dp)
597 32 : direction_vector = direction_vector - origin
598 36 : r2_origin_cell(icell) = DOT_PRODUCT(direction_vector, direction_vector)
599 : END DO
600 :
601 4 : IF (ABS(r2_origin_cell(1) - r2_origin_cell(2)) < (eps_geometry*eps_geometry)) THEN
602 : ! primary and secondary bulk unit cells should not overlap;
603 : ! currently we check that they are different by at least one atom that is, indeed, not sufficient.
604 : CALL cp_abort(__LOCATION__, &
605 0 : "Primary and secondary bulk contact cells should not overlap ")
606 4 : ELSE IF (r2_origin_cell(1) < r2_origin_cell(2)) THEN
607 6 : ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(1)%vector)))
608 10 : contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(1)%vector(:)
609 6 : ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(2)%vector)))
610 10 : contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(2)%vector(:)
611 : ELSE
612 6 : ALLOCATE (contact_env%atomlist_cell0(SIZE(contact_control%atomlist_cell(2)%vector)))
613 10 : contact_env%atomlist_cell0(:) = contact_control%atomlist_cell(2)%vector(:)
614 6 : ALLOCATE (contact_env%atomlist_cell1(SIZE(contact_control%atomlist_cell(1)%vector)))
615 10 : contact_env%atomlist_cell1(:) = contact_control%atomlist_cell(1)%vector(:)
616 : END IF
617 :
618 4 : NULLIFY (fm_struct)
619 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao_c, ncol_global=nao_c, context=sub_env%blacs_env)
620 28 : ALLOCATE (contact_env%h_00(nspins), contact_env%h_01(nspins))
621 28 : ALLOCATE (contact_env%rho_00(nspins), contact_env%rho_01(nspins))
622 8 : DO ispin = 1, nspins
623 4 : CALL cp_fm_create(contact_env%h_00(ispin), fm_struct)
624 4 : CALL cp_fm_create(contact_env%h_01(ispin), fm_struct)
625 4 : CALL cp_fm_create(contact_env%rho_00(ispin), fm_struct)
626 8 : CALL cp_fm_create(contact_env%rho_01(ispin), fm_struct)
627 : END DO
628 4 : ALLOCATE (contact_env%s_00, contact_env%s_01)
629 4 : CALL cp_fm_create(contact_env%s_00, fm_struct)
630 4 : CALL cp_fm_create(contact_env%s_01, fm_struct)
631 4 : CALL cp_fm_struct_release(fm_struct)
632 :
633 8 : DO ispin = 1, nspins
634 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
635 : fm=contact_env%h_00(ispin), &
636 : atomlist_row=contact_env%atomlist_cell0, &
637 : atomlist_col=contact_env%atomlist_cell0, &
638 : subsys=subsys, mpi_comm_global=para_env, &
639 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
640 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
641 : fm=contact_env%h_01(ispin), &
642 : atomlist_row=contact_env%atomlist_cell0, &
643 : atomlist_col=contact_env%atomlist_cell1, &
644 : subsys=subsys, mpi_comm_global=para_env, &
645 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
646 :
647 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
648 : fm=contact_env%rho_00(ispin), &
649 : atomlist_row=contact_env%atomlist_cell0, &
650 : atomlist_col=contact_env%atomlist_cell0, &
651 : subsys=subsys, mpi_comm_global=para_env, &
652 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
653 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_kp(ispin, 1)%matrix, &
654 : fm=contact_env%rho_01(ispin), &
655 : atomlist_row=contact_env%atomlist_cell0, &
656 : atomlist_col=contact_env%atomlist_cell1, &
657 : subsys=subsys, mpi_comm_global=para_env, &
658 8 : do_upper_diag=.TRUE., do_lower=.TRUE.)
659 : END DO
660 :
661 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
662 : fm=contact_env%s_00, &
663 : atomlist_row=contact_env%atomlist_cell0, &
664 : atomlist_col=contact_env%atomlist_cell0, &
665 : subsys=subsys, mpi_comm_global=para_env, &
666 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
667 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
668 : fm=contact_env%s_01, &
669 : atomlist_row=contact_env%atomlist_cell0, &
670 : atomlist_col=contact_env%atomlist_cell1, &
671 : subsys=subsys, mpi_comm_global=para_env, &
672 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
673 :
674 4 : CALL timestop(handle)
675 4 : END SUBROUTINE negf_env_contact_init_matrices_gamma
676 :
677 : ! **************************************************************************************************
678 : !> \brief Extract relevant matrix blocks for the scattering region as well as
679 : !> all the scattering -- contact interface regions.
680 : !> \param negf_env NEGF environment (modified on exit)
681 : !> \param negf_control NEGF control
682 : !> \param sub_env NEGF parallel (sub)group environment
683 : !> \param qs_env Primary QuickStep environment
684 : !> \author Sergey Chulkov
685 : ! **************************************************************************************************
686 4 : SUBROUTINE negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env)
687 : TYPE(negf_env_type), INTENT(inout) :: negf_env
688 : TYPE(negf_control_type), POINTER :: negf_control
689 : TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env
690 : TYPE(qs_environment_type), POINTER :: qs_env
691 :
692 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_env_device_init_matrices'
693 :
694 : INTEGER :: handle, icontact, ispin, nao_c, nao_s, &
695 : ncontacts, nspins
696 : LOGICAL :: do_kpoints
697 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
698 : TYPE(dbcsr_p_type) :: hmat
699 4 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp
700 : TYPE(dft_control_type), POINTER :: dft_control
701 : TYPE(mp_para_env_type), POINTER :: para_env
702 : TYPE(pw_env_type), POINTER :: pw_env
703 : TYPE(pw_pool_type), POINTER :: pw_pool
704 : TYPE(pw_r3d_rs_type) :: v_hartree
705 : TYPE(qs_subsys_type), POINTER :: subsys
706 :
707 4 : CALL timeset(routineN, handle)
708 :
709 4 : IF (ALLOCATED(negf_control%atomlist_S_screening)) THEN
710 : CALL get_qs_env(qs_env, &
711 : dft_control=dft_control, &
712 : do_kpoints=do_kpoints, &
713 : matrix_ks_kp=matrix_ks_kp, &
714 : matrix_s_kp=matrix_s_kp, &
715 : para_env=para_env, &
716 : pw_env=pw_env, &
717 4 : subsys=subsys)
718 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
719 :
720 4 : IF (do_kpoints) THEN
721 : CALL cp_abort(__LOCATION__, &
722 0 : "K-points in device region have not been implemented yet.")
723 : END IF
724 :
725 4 : ncontacts = SIZE(negf_control%contacts)
726 4 : nspins = dft_control%nspins
727 :
728 4 : NULLIFY (fm_struct)
729 4 : nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening)
730 :
731 : ! ++ create matrices: h_s, s_s
732 4 : NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct)
733 16 : ALLOCATE (negf_env%h_s(nspins))
734 :
735 4 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env)
736 4 : ALLOCATE (negf_env%s_s)
737 4 : CALL cp_fm_create(negf_env%s_s, fm_struct)
738 8 : DO ispin = 1, nspins
739 8 : CALL cp_fm_create(negf_env%h_s(ispin), fm_struct)
740 : END DO
741 4 : ALLOCATE (negf_env%v_hartree_s)
742 4 : CALL cp_fm_create(negf_env%v_hartree_s, fm_struct)
743 4 : CALL cp_fm_struct_release(fm_struct)
744 :
745 : ! ++ create matrices: h_sc, s_sc
746 48 : ALLOCATE (negf_env%h_sc(nspins, ncontacts), negf_env%s_sc(ncontacts))
747 12 : DO icontact = 1, ncontacts
748 8 : nao_c = number_of_atomic_orbitals(subsys, negf_env%contacts(icontact)%atomlist_cell0)
749 8 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_c, context=sub_env%blacs_env)
750 :
751 8 : CALL cp_fm_create(negf_env%s_sc(icontact), fm_struct)
752 :
753 16 : DO ispin = 1, nspins
754 16 : CALL cp_fm_create(negf_env%h_sc(ispin, icontact), fm_struct)
755 : END DO
756 :
757 12 : CALL cp_fm_struct_release(fm_struct)
758 : END DO
759 :
760 : ! extract matrices: h_s, s_s
761 8 : DO ispin = 1, nspins
762 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
763 : fm=negf_env%h_s(ispin), &
764 : atomlist_row=negf_control%atomlist_S_screening, &
765 : atomlist_col=negf_control%atomlist_S_screening, &
766 : subsys=subsys, mpi_comm_global=para_env, &
767 8 : do_upper_diag=.TRUE., do_lower=.TRUE.)
768 : END DO
769 :
770 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
771 : fm=negf_env%s_s, &
772 : atomlist_row=negf_control%atomlist_S_screening, &
773 : atomlist_col=negf_control%atomlist_S_screening, &
774 : subsys=subsys, mpi_comm_global=para_env, &
775 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
776 :
777 : ! v_hartree_s
778 4 : NULLIFY (hmat%matrix)
779 4 : CALL dbcsr_init_p(hmat%matrix)
780 4 : CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix)
781 4 : CALL dbcsr_set(hmat%matrix, 0.0_dp)
782 :
783 4 : CALL pw_pool%create_pw(v_hartree)
784 4 : CALL negf_env_init_v_hartree(v_hartree, negf_env%contacts, negf_control%contacts)
785 :
786 : CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, &
787 4 : calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE.)
788 :
789 4 : CALL pw_pool%give_back_pw(v_hartree)
790 :
791 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=hmat%matrix, &
792 : fm=negf_env%v_hartree_s, &
793 : atomlist_row=negf_control%atomlist_S_screening, &
794 : atomlist_col=negf_control%atomlist_S_screening, &
795 : subsys=subsys, mpi_comm_global=para_env, &
796 4 : do_upper_diag=.TRUE., do_lower=.TRUE.)
797 :
798 4 : CALL dbcsr_deallocate_matrix(hmat%matrix)
799 :
800 : ! extract matrices: h_sc, s_sc
801 12 : DO icontact = 1, ncontacts
802 16 : DO ispin = 1, nspins
803 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_kp(ispin, 1)%matrix, &
804 : fm=negf_env%h_sc(ispin, icontact), &
805 : atomlist_row=negf_control%atomlist_S_screening, &
806 : atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
807 : subsys=subsys, mpi_comm_global=para_env, &
808 16 : do_upper_diag=.TRUE., do_lower=.TRUE.)
809 : END DO
810 :
811 : CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_s_kp(1, 1)%matrix, &
812 : fm=negf_env%s_sc(icontact), &
813 : atomlist_row=negf_control%atomlist_S_screening, &
814 : atomlist_col=negf_env%contacts(icontact)%atomlist_cell0, &
815 : subsys=subsys, mpi_comm_global=para_env, &
816 12 : do_upper_diag=.TRUE., do_lower=.TRUE.)
817 : END DO
818 : END IF
819 :
820 4 : CALL timestop(handle)
821 4 : END SUBROUTINE negf_env_device_init_matrices
822 :
823 : ! **************************************************************************************************
824 : !> \brief Contribution to the Hartree potential related to the external bias voltage.
825 : !> \param v_hartree Hartree potential (modified on exit)
826 : !> \param contact_env NEGF environment for every contact
827 : !> \param contact_control NEGF control for every contact
828 : !> \author Sergey Chulkov
829 : ! **************************************************************************************************
830 4 : SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control)
831 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_hartree
832 : TYPE(negf_env_contact_type), DIMENSION(:), &
833 : INTENT(in) :: contact_env
834 : TYPE(negf_control_contact_type), DIMENSION(:), &
835 : INTENT(in) :: contact_control
836 :
837 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_init_v_hartree'
838 : REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp)
839 :
840 : INTEGER :: dx, dy, dz, handle, icontact, ix, iy, &
841 : iz, lx, ly, lz, ncontacts, ux, uy, uz
842 : REAL(kind=dp) :: dvol, pot, proj, v1, v2
843 : REAL(kind=dp), DIMENSION(3) :: dirvector_bias, point_coord, &
844 : point_indices, vector
845 :
846 4 : CALL timeset(routineN, handle)
847 :
848 4 : ncontacts = SIZE(contact_env)
849 4 : CPASSERT(SIZE(contact_control) == ncontacts)
850 4 : CPASSERT(ncontacts == 2)
851 :
852 16 : dirvector_bias = contact_env(2)%origin_bias - contact_env(1)%origin_bias
853 4 : v1 = contact_control(1)%v_external
854 4 : v2 = contact_control(2)%v_external
855 :
856 4 : lx = v_hartree%pw_grid%bounds_local(1, 1)
857 4 : ux = v_hartree%pw_grid%bounds_local(2, 1)
858 4 : ly = v_hartree%pw_grid%bounds_local(1, 2)
859 4 : uy = v_hartree%pw_grid%bounds_local(2, 2)
860 4 : lz = v_hartree%pw_grid%bounds_local(1, 3)
861 4 : uz = v_hartree%pw_grid%bounds_local(2, 3)
862 :
863 4 : dx = v_hartree%pw_grid%npts(1)/2
864 4 : dy = v_hartree%pw_grid%npts(2)/2
865 4 : dz = v_hartree%pw_grid%npts(3)/2
866 :
867 4 : dvol = v_hartree%pw_grid%dvol
868 :
869 1764 : DO iz = lz, uz
870 1760 : point_indices(3) = REAL(iz + dz, kind=dp)
871 112644 : DO iy = ly, uy
872 110880 : point_indices(2) = REAL(iy + dy, kind=dp)
873 :
874 3605360 : DO ix = lx, ux
875 3492720 : point_indices(1) = REAL(ix + dx, kind=dp)
876 45405360 : point_coord(:) = MATMUL(v_hartree%pw_grid%dh, point_indices)
877 :
878 13970880 : vector = point_coord - contact_env(1)%origin_bias
879 3492720 : proj = projection_on_direction_vector(vector, dirvector_bias)
880 3492720 : IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
881 : ! scattering region
882 : ! proj == 0 we are at the first contact boundary
883 : ! proj == 1 we are at the second contact boundary
884 1373274 : IF (proj < 0.0_dp) THEN
885 : proj = 0.0_dp
886 1373274 : ELSE IF (proj > 1.0_dp) THEN
887 0 : proj = 1.0_dp
888 : END IF
889 1373274 : pot = v1 + (v2 - v1)*proj
890 : ELSE
891 3357774 : pot = 0.0_dp
892 3357774 : DO icontact = 1, ncontacts
893 12954816 : vector = point_coord - contact_env(icontact)%origin_bias
894 3238704 : proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias)
895 :
896 3357774 : IF (proj + threshold >= 0.0_dp .AND. proj - threshold <= 1.0_dp) THEN
897 2000376 : pot = contact_control(icontact)%v_external
898 2000376 : EXIT
899 : END IF
900 : END DO
901 : END IF
902 :
903 3603600 : v_hartree%array(ix, iy, iz) = pot*dvol
904 : END DO
905 : END DO
906 : END DO
907 :
908 4 : CALL timestop(handle)
909 4 : END SUBROUTINE negf_env_init_v_hartree
910 :
911 : ! **************************************************************************************************
912 : !> \brief Detect the axis towards secondary unit cell.
913 : !> \param direction_vector direction vector
914 : !> \param subsys_contact QuickStep subsystem of the contact force environment
915 : !> \param eps_geometry accuracy in mapping atoms between different force environments
916 : !> \return direction axis: 0 (undefined), 1 (x), 2(y), 3 (z)
917 : !> \par History
918 : !> * 08.2017 created [Sergey Chulkov]
919 : ! **************************************************************************************************
920 8 : FUNCTION contact_direction_axis(direction_vector, subsys_contact, eps_geometry) RESULT(direction_axis)
921 : REAL(kind=dp), DIMENSION(3), INTENT(in) :: direction_vector
922 : TYPE(qs_subsys_type), POINTER :: subsys_contact
923 : REAL(kind=dp), INTENT(in) :: eps_geometry
924 : INTEGER :: direction_axis
925 :
926 : INTEGER :: i, naxes
927 : REAL(kind=dp), DIMENSION(3) :: scaled
928 : TYPE(cell_type), POINTER :: cell
929 :
930 8 : CALL qs_subsys_get(subsys_contact, cell=cell)
931 8 : CALL real_to_scaled(scaled, direction_vector, cell)
932 :
933 8 : naxes = 0
934 8 : direction_axis = 0 ! initialize to make GCC<=6 happy
935 :
936 32 : DO i = 1, 3
937 32 : IF (ABS(scaled(i)) > eps_geometry) THEN
938 8 : IF (scaled(i) > 0.0_dp) THEN
939 : direction_axis = i
940 : ELSE
941 4 : direction_axis = -i
942 : END IF
943 8 : naxes = naxes + 1
944 : END IF
945 : END DO
946 :
947 : ! direction_vector is not parallel to one of the unit cell's axis
948 8 : IF (naxes /= 1) direction_axis = 0
949 8 : END FUNCTION contact_direction_axis
950 :
951 : ! **************************************************************************************************
952 : !> \brief Estimate energy of the highest spin-alpha occupied molecular orbital.
953 : !> \param homo_energy HOMO energy (initialised on exit)
954 : !> \param qs_env QuickStep environment
955 : !> \par History
956 : !> * 01.2017 created [Sergey Chulkov]
957 : ! **************************************************************************************************
958 4 : SUBROUTINE negf_homo_energy_estimate(homo_energy, qs_env)
959 : REAL(kind=dp), INTENT(out) :: homo_energy
960 : TYPE(qs_environment_type), POINTER :: qs_env
961 :
962 : CHARACTER(LEN=*), PARAMETER :: routineN = 'negf_homo_energy_estimate'
963 : INTEGER, PARAMETER :: gamma_point = 1
964 :
965 : INTEGER :: handle, homo, ikpgr, ikpoint, imo, &
966 : ispin, kplocal, nmo, nspins
967 : INTEGER, DIMENSION(2) :: kp_range
968 : LOGICAL :: do_kpoints
969 : REAL(kind=dp) :: my_homo_energy
970 4 : REAL(kind=dp), DIMENSION(:), POINTER :: eigenvalues
971 4 : TYPE(kpoint_env_p_type), DIMENSION(:), POINTER :: kp_env
972 : TYPE(kpoint_type), POINTER :: kpoints
973 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
974 4 : TYPE(mo_set_type), DIMENSION(:, :), POINTER :: mos_kp
975 : TYPE(mp_para_env_type), POINTER :: para_env, para_env_kp
976 :
977 4 : CALL timeset(routineN, handle)
978 4 : my_homo_energy = 0.0_dp
979 :
980 4 : CALL get_qs_env(qs_env, para_env=para_env, mos=mos, kpoints=kpoints, do_kpoints=do_kpoints)
981 :
982 4 : IF (do_kpoints) THEN
983 4 : CALL get_kpoint_info(kpoints, kp_env=kp_env, kp_range=kp_range, para_env_kp=para_env_kp)
984 :
985 : ! looking for a processor that holds the gamma point
986 4 : IF (para_env_kp%mepos == 0 .AND. kp_range(1) <= gamma_point .AND. kp_range(2) >= gamma_point) THEN
987 2 : kplocal = kp_range(2) - kp_range(1) + 1
988 :
989 2 : DO ikpgr = 1, kplocal
990 2 : CALL get_kpoint_env(kp_env(ikpgr)%kpoint_env, nkpoint=ikpoint, mos=mos_kp)
991 :
992 2 : IF (ikpoint == gamma_point) THEN
993 : ! mos_kp(component, spin), where component = 1 (real), or 2 (imaginary)
994 2 : CALL get_mo_set(mos_kp(1, 1), homo=homo, eigenvalues=eigenvalues) ! mu=fermi_level
995 :
996 2 : my_homo_energy = eigenvalues(homo)
997 2 : EXIT
998 : END IF
999 : END DO
1000 : END IF
1001 :
1002 4 : CALL para_env%sum(my_homo_energy)
1003 : ELSE
1004 : ! Hamiltonian of the bulk contact region has been computed without k-points.
1005 : ! Try to obtain the HOMO energy assuming there is no OT. We probably should abort here
1006 : ! as we do need a second replica of the bulk contact unit cell along transport
1007 : ! direction anyway which is not available without k-points.
1008 :
1009 0 : CPABORT("It is strongly advised to use k-points within all contact FORCE_EVAL-s")
1010 :
1011 0 : nspins = SIZE(mos)
1012 :
1013 0 : spin_loop: DO ispin = 1, nspins
1014 0 : CALL get_mo_set(mos(ispin), homo=homo, nmo=nmo, eigenvalues=eigenvalues)
1015 :
1016 0 : DO imo = nmo, 1, -1
1017 0 : IF (eigenvalues(imo) /= 0.0_dp) EXIT spin_loop
1018 : END DO
1019 : END DO spin_loop
1020 :
1021 0 : IF (imo == 0) THEN
1022 0 : CPABORT("Orbital transformation (OT) for contact FORCE_EVAL-s is not supported")
1023 : END IF
1024 :
1025 0 : my_homo_energy = eigenvalues(homo)
1026 : END IF
1027 :
1028 4 : homo_energy = my_homo_energy
1029 4 : CALL timestop(handle)
1030 4 : END SUBROUTINE negf_homo_energy_estimate
1031 :
1032 : ! **************************************************************************************************
1033 : !> \brief List atoms from the contact's primary unit cell.
1034 : !> \param atomlist_cell0 list of atoms belonging to the contact's primary unit cell
1035 : !> (allocate and initialised on exit)
1036 : !> \param atom_map_cell0 atomic map of atoms from 'atomlist_cell0' (allocate and initialised on exit)
1037 : !> \param atomlist_bulk list of atoms belonging to the bulk contact region
1038 : !> \param atom_map atomic map of atoms from 'atomlist_bulk'
1039 : !> \param origin origin of the contact
1040 : !> \param direction_vector direction vector of the contact
1041 : !> \param direction_axis axis towards secondary unit cell
1042 : !> \param subsys_device QuickStep subsystem of the device force environment
1043 : !> \par History
1044 : !> * 08.2017 created [Sergey Chulkov]
1045 : ! **************************************************************************************************
1046 4 : SUBROUTINE list_atoms_in_bulk_primary_unit_cell(atomlist_cell0, atom_map_cell0, atomlist_bulk, atom_map, &
1047 : origin, direction_vector, direction_axis, subsys_device)
1048 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell0
1049 : TYPE(negf_atom_map_type), ALLOCATABLE, &
1050 : DIMENSION(:), INTENT(inout) :: atom_map_cell0
1051 : INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1052 : TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1053 : REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1054 : INTEGER, INTENT(in) :: direction_axis
1055 : TYPE(qs_subsys_type), POINTER :: subsys_device
1056 :
1057 : CHARACTER(LEN=*), PARAMETER :: routineN = 'list_atoms_in_bulk_primary_unit_cell'
1058 :
1059 : INTEGER :: atom_min, dir_axis_min, &
1060 : direction_axis_abs, handle, iatom, &
1061 : natoms_bulk, natoms_cell0
1062 : REAL(kind=dp) :: proj, proj_min
1063 : REAL(kind=dp), DIMENSION(3) :: vector
1064 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1065 :
1066 4 : CALL timeset(routineN, handle)
1067 4 : CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1068 :
1069 4 : natoms_bulk = SIZE(atomlist_bulk)
1070 4 : CPASSERT(SIZE(atom_map, 1) == natoms_bulk)
1071 4 : direction_axis_abs = ABS(direction_axis)
1072 :
1073 : ! looking for the nearest atom from the scattering region
1074 4 : proj_min = 1.0_dp
1075 4 : atom_min = 1
1076 36 : DO iatom = 1, natoms_bulk
1077 128 : vector = particle_set(atomlist_bulk(iatom))%r - origin
1078 32 : proj = projection_on_direction_vector(vector, direction_vector)
1079 :
1080 36 : IF (proj < proj_min) THEN
1081 16 : proj_min = proj
1082 16 : atom_min = iatom
1083 : END IF
1084 : END DO
1085 :
1086 4 : dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1087 :
1088 4 : natoms_cell0 = 0
1089 36 : DO iatom = 1, natoms_bulk
1090 32 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) &
1091 20 : natoms_cell0 = natoms_cell0 + 1
1092 : END DO
1093 :
1094 12 : ALLOCATE (atomlist_cell0(natoms_cell0))
1095 40 : ALLOCATE (atom_map_cell0(natoms_cell0))
1096 :
1097 4 : natoms_cell0 = 0
1098 36 : DO iatom = 1, natoms_bulk
1099 36 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min) THEN
1100 16 : natoms_cell0 = natoms_cell0 + 1
1101 16 : atomlist_cell0(natoms_cell0) = atomlist_bulk(iatom)
1102 16 : atom_map_cell0(natoms_cell0) = atom_map(iatom)
1103 : END IF
1104 : END DO
1105 :
1106 4 : CALL timestop(handle)
1107 4 : END SUBROUTINE list_atoms_in_bulk_primary_unit_cell
1108 :
1109 : ! **************************************************************************************************
1110 : !> \brief List atoms from the contact's secondary unit cell.
1111 : !> \param atomlist_cell1 list of atoms belonging to the contact's secondary unit cell
1112 : !> (allocate and initialised on exit)
1113 : !> \param atom_map_cell1 atomic map of atoms from 'atomlist_cell1'
1114 : !> (allocate and initialised on exit)
1115 : !> \param atomlist_bulk list of atoms belonging to the bulk contact region
1116 : !> \param atom_map atomic map of atoms from 'atomlist_bulk'
1117 : !> \param origin origin of the contact
1118 : !> \param direction_vector direction vector of the contact
1119 : !> \param direction_axis axis towards the secondary unit cell
1120 : !> \param subsys_device QuickStep subsystem of the device force environment
1121 : !> \par History
1122 : !> * 11.2017 created [Sergey Chulkov]
1123 : !> \note Cloned from list_atoms_in_bulk_primary_unit_cell. Will be removed once we can managed to
1124 : !> maintain consistency between real-space matrices from different force_eval sections.
1125 : ! **************************************************************************************************
1126 4 : SUBROUTINE list_atoms_in_bulk_secondary_unit_cell(atomlist_cell1, atom_map_cell1, atomlist_bulk, atom_map, &
1127 : origin, direction_vector, direction_axis, subsys_device)
1128 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(inout) :: atomlist_cell1
1129 : TYPE(negf_atom_map_type), ALLOCATABLE, &
1130 : DIMENSION(:), INTENT(inout) :: atom_map_cell1
1131 : INTEGER, DIMENSION(:), INTENT(in) :: atomlist_bulk
1132 : TYPE(negf_atom_map_type), DIMENSION(:), INTENT(in) :: atom_map
1133 : REAL(kind=dp), DIMENSION(3), INTENT(in) :: origin, direction_vector
1134 : INTEGER, INTENT(in) :: direction_axis
1135 : TYPE(qs_subsys_type), POINTER :: subsys_device
1136 :
1137 : CHARACTER(LEN=*), PARAMETER :: routineN = 'list_atoms_in_bulk_secondary_unit_cell'
1138 :
1139 : INTEGER :: atom_min, dir_axis_min, &
1140 : direction_axis_abs, handle, iatom, &
1141 : natoms_bulk, natoms_cell1, offset
1142 : REAL(kind=dp) :: proj, proj_min
1143 : REAL(kind=dp), DIMENSION(3) :: vector
1144 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1145 :
1146 4 : CALL timeset(routineN, handle)
1147 4 : CALL qs_subsys_get(subsys_device, particle_set=particle_set)
1148 :
1149 4 : natoms_bulk = SIZE(atomlist_bulk)
1150 4 : CPASSERT(SIZE(atom_map, 1) == natoms_bulk)
1151 4 : direction_axis_abs = ABS(direction_axis)
1152 4 : offset = SIGN(1, direction_axis)
1153 :
1154 : ! looking for the nearest atom from the scattering region
1155 4 : proj_min = 1.0_dp
1156 4 : atom_min = 1
1157 36 : DO iatom = 1, natoms_bulk
1158 128 : vector = particle_set(atomlist_bulk(iatom))%r - origin
1159 32 : proj = projection_on_direction_vector(vector, direction_vector)
1160 :
1161 36 : IF (proj < proj_min) THEN
1162 16 : proj_min = proj
1163 16 : atom_min = iatom
1164 : END IF
1165 : END DO
1166 :
1167 4 : dir_axis_min = atom_map(atom_min)%cell(direction_axis_abs)
1168 :
1169 4 : natoms_cell1 = 0
1170 36 : DO iatom = 1, natoms_bulk
1171 32 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) &
1172 20 : natoms_cell1 = natoms_cell1 + 1
1173 : END DO
1174 :
1175 12 : ALLOCATE (atomlist_cell1(natoms_cell1))
1176 40 : ALLOCATE (atom_map_cell1(natoms_cell1))
1177 :
1178 4 : natoms_cell1 = 0
1179 36 : DO iatom = 1, natoms_bulk
1180 36 : IF (atom_map(iatom)%cell(direction_axis_abs) == dir_axis_min + offset) THEN
1181 16 : natoms_cell1 = natoms_cell1 + 1
1182 16 : atomlist_cell1(natoms_cell1) = atomlist_bulk(iatom)
1183 16 : atom_map_cell1(natoms_cell1) = atom_map(iatom)
1184 16 : atom_map_cell1(natoms_cell1)%cell(direction_axis_abs) = dir_axis_min
1185 : END IF
1186 : END DO
1187 :
1188 4 : CALL timestop(handle)
1189 4 : END SUBROUTINE list_atoms_in_bulk_secondary_unit_cell
1190 :
1191 : ! **************************************************************************************************
1192 : !> \brief Release a NEGF environment variable.
1193 : !> \param negf_env NEGF environment to release
1194 : !> \par History
1195 : !> * 01.2017 created [Sergey Chulkov]
1196 : ! **************************************************************************************************
1197 4 : SUBROUTINE negf_env_release(negf_env)
1198 : TYPE(negf_env_type), INTENT(inout) :: negf_env
1199 :
1200 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_release'
1201 :
1202 : INTEGER :: handle, icontact
1203 :
1204 4 : CALL timeset(routineN, handle)
1205 :
1206 4 : IF (ALLOCATED(negf_env%contacts)) THEN
1207 12 : DO icontact = SIZE(negf_env%contacts), 1, -1
1208 12 : CALL negf_env_contact_release(negf_env%contacts(icontact))
1209 : END DO
1210 :
1211 12 : DEALLOCATE (negf_env%contacts)
1212 : END IF
1213 :
1214 : ! h_s
1215 4 : CALL cp_fm_release(negf_env%h_s)
1216 :
1217 : ! h_sc
1218 4 : CALL cp_fm_release(negf_env%h_sc)
1219 :
1220 : ! s_s
1221 4 : IF (ASSOCIATED(negf_env%s_s)) THEN
1222 4 : CALL cp_fm_release(negf_env%s_s)
1223 4 : DEALLOCATE (negf_env%s_s)
1224 : NULLIFY (negf_env%s_s)
1225 : END IF
1226 :
1227 : ! s_sc
1228 4 : CALL cp_fm_release(negf_env%s_sc)
1229 :
1230 : ! v_hartree_s
1231 4 : IF (ASSOCIATED(negf_env%v_hartree_s)) THEN
1232 4 : CALL cp_fm_release(negf_env%v_hartree_s)
1233 4 : DEALLOCATE (negf_env%v_hartree_s)
1234 : NULLIFY (negf_env%v_hartree_s)
1235 : END IF
1236 :
1237 : ! density mixing
1238 4 : IF (ASSOCIATED(negf_env%mixing_storage)) THEN
1239 4 : CALL mixing_storage_release(negf_env%mixing_storage)
1240 4 : DEALLOCATE (negf_env%mixing_storage)
1241 : END IF
1242 :
1243 4 : CALL timestop(handle)
1244 4 : END SUBROUTINE negf_env_release
1245 :
1246 : ! **************************************************************************************************
1247 : !> \brief Release a NEGF contact environment variable.
1248 : !> \param contact_env NEGF contact environment to release
1249 : ! **************************************************************************************************
1250 8 : SUBROUTINE negf_env_contact_release(contact_env)
1251 : TYPE(negf_env_contact_type), INTENT(inout) :: contact_env
1252 :
1253 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_contact_release'
1254 :
1255 : INTEGER :: handle
1256 :
1257 8 : CALL timeset(routineN, handle)
1258 :
1259 : ! h_00
1260 8 : CALL cp_fm_release(contact_env%h_00)
1261 :
1262 : ! h_01
1263 8 : CALL cp_fm_release(contact_env%h_01)
1264 :
1265 : ! rho_00
1266 8 : CALL cp_fm_release(contact_env%rho_00)
1267 :
1268 : ! rho_01
1269 8 : CALL cp_fm_release(contact_env%rho_01)
1270 :
1271 : ! s_00
1272 8 : IF (ASSOCIATED(contact_env%s_00)) THEN
1273 8 : CALL cp_fm_release(contact_env%s_00)
1274 8 : DEALLOCATE (contact_env%s_00)
1275 : NULLIFY (contact_env%s_00)
1276 : END IF
1277 :
1278 : ! s_01
1279 8 : IF (ASSOCIATED(contact_env%s_01)) THEN
1280 8 : CALL cp_fm_release(contact_env%s_01)
1281 8 : DEALLOCATE (contact_env%s_01)
1282 : NULLIFY (contact_env%s_01)
1283 : END IF
1284 :
1285 8 : IF (ALLOCATED(contact_env%atomlist_cell0)) DEALLOCATE (contact_env%atomlist_cell0)
1286 8 : IF (ALLOCATED(contact_env%atomlist_cell1)) DEALLOCATE (contact_env%atomlist_cell1)
1287 8 : IF (ALLOCATED(contact_env%atom_map_cell0)) DEALLOCATE (contact_env%atom_map_cell0)
1288 :
1289 8 : CALL timestop(handle)
1290 8 : END SUBROUTINE negf_env_contact_release
1291 0 : END MODULE negf_env_types
|
