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