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 Input control types for NEGF based quantum transport calculations
10 : ! **************************************************************************************************
11 :
12 : MODULE negf_control_types
13 : USE cp_subsys_types, ONLY: cp_subsys_get,&
14 : cp_subsys_type
15 : USE input_constants, ONLY: negf_run
16 : USE input_section_types, ONLY: section_vals_get,&
17 : section_vals_get_subs_vals,&
18 : section_vals_type,&
19 : section_vals_val_get
20 : USE kinds, ONLY: default_string_length,&
21 : dp
22 : USE mathconstants, ONLY: pi
23 : USE molecule_kind_types, ONLY: get_molecule_kind,&
24 : molecule_kind_type
25 : USE molecule_types, ONLY: get_molecule,&
26 : molecule_type
27 : USE negf_alloc_types, ONLY: negf_allocatable_ivector
28 : USE particle_types, ONLY: particle_type
29 : USE physcon, ONLY: kelvin
30 : USE string_utilities, ONLY: integer_to_string
31 : USE util, ONLY: sort
32 : #include "./base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_control_types'
38 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
39 :
40 : PUBLIC :: negf_control_type, negf_control_contact_type
41 : PUBLIC :: negf_control_create, negf_control_release, read_negf_control
42 :
43 : ! **************************************************************************************************
44 : !> \brief Input parameters related to a single contact.
45 : !> \author Sergey Chulkov
46 : ! **************************************************************************************************
47 : TYPE negf_control_contact_type
48 : !> atoms belonging to bulk and screening regions
49 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_bulk, atomlist_screening
50 : !> atom belonging to the primary and secondary bulk unit cells
51 : TYPE(negf_allocatable_ivector), ALLOCATABLE, &
52 : DIMENSION(:) :: atomlist_cell
53 : !> index of the sub_force_env which should be used for bulk calculation
54 : INTEGER :: force_env_index = -1
55 : !> contact Fermi level needs to be computed prior NEGF run
56 : LOGICAL :: compute_fermi_level = .FALSE.
57 : !> when computing contact Fermi level, use the energy given in 'fermi_level' (instead of HOMO)
58 : !> (instead of the HOMO energy) as a starting point
59 : LOGICAL :: refine_fermi_level = .FALSE.
60 : !> Fermi level
61 : REAL(kind=dp) :: fermi_level = -1.0_dp
62 : !> temperature [in a.u.]
63 : REAL(kind=dp) :: temperature = -1.0_dp
64 : !> applied electric potential
65 : REAL(kind=dp) :: v_external = 0.0_dp
66 : END TYPE negf_control_contact_type
67 :
68 : ! **************************************************************************************************
69 : !> \brief Input parameters related to the NEGF run.
70 : !> \author Sergey Chulkov
71 : ! **************************************************************************************************
72 : TYPE negf_control_type
73 : !> input options for every contact
74 : TYPE(negf_control_contact_type), ALLOCATABLE, &
75 : DIMENSION(:) :: contacts
76 : !> atoms belonging to the scattering region
77 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_S
78 : !> atoms belonging to the scattering region as well as atoms belonging to
79 : !> screening regions of all the contacts
80 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_S_screening
81 : !> do not keep contact self-energy matrices
82 : LOGICAL :: disable_cache = .FALSE.
83 : !> convergence criteria for adaptive integration methods
84 : REAL(kind=dp) :: conv_density = -1.0_dp
85 : !> convergence criteria for iterative Lopez-Sancho algorithm
86 : REAL(kind=dp) :: conv_green = -1.0_dp
87 : !> convergence criteria for self-consistent iterations
88 : REAL(kind=dp) :: conv_scf = -1.0_dp
89 : !> accuracy in mapping atoms between different force environments
90 : REAL(kind=dp) :: eps_geometry = -1.0_dp
91 : !> applied bias [in a.u.]
92 : REAL(kind=dp) :: v_bias = -1.0_dp
93 : !> integration lower bound [in a.u.]
94 : REAL(kind=dp) :: energy_lbound = -1.0_dp
95 : !> infinitesimal offset along the imaginary axis [in a.u.]
96 : REAL(kind=dp) :: eta = -1.0_dp
97 : !> initial guess to determine the actual Fermi level of bulk contacts [in a.u.]
98 : REAL(kind=dp) :: homo_lumo_gap = -1.0_dp
99 : !> number of residuals (poles of the Fermi function)
100 : INTEGER :: delta_npoles = -1
101 : !> offset along the x-axis away from the poles of the Fermi function [in units of kT]
102 : INTEGER :: gamma_kT = -1
103 : !> integration method
104 : INTEGER :: integr_method = -1
105 : !> minimal number of grid points along the closed contour
106 : INTEGER :: integr_min_points = -1
107 : !> maximal number of grid points along the closed contour
108 : INTEGER :: integr_max_points = -1
109 : !> maximal number of SCF iterations
110 : INTEGER :: max_scf = -1
111 : !> minimal number of MPI processes to be used to compute Green's function per energy point
112 : INTEGER :: nprocs = -1
113 : !> shift in Hartree potential [in a.u.]
114 : REAL(kind=dp) :: v_shift = -1.0_dp
115 : !> initial offset to determine the correct shift in Hartree potential [in a.u.]
116 : REAL(kind=dp) :: v_shift_offset = -1.0_dp
117 : !> maximal number of iteration to determine the shift in Hartree potential
118 : INTEGER :: v_shift_maxiters = -1
119 : END TYPE negf_control_type
120 :
121 : PRIVATE :: read_negf_atomlist
122 :
123 : CONTAINS
124 :
125 : ! **************************************************************************************************
126 : !> \brief allocate control options for Non-equilibrium Green's Function calculation
127 : !> \param negf_control an object to create
128 : !> \par History
129 : !> * 02.2017 created [Sergey Chulkov]
130 : ! **************************************************************************************************
131 8 : SUBROUTINE negf_control_create(negf_control)
132 : TYPE(negf_control_type), POINTER :: negf_control
133 :
134 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_control_create'
135 :
136 : INTEGER :: handle
137 :
138 4 : CPASSERT(.NOT. ASSOCIATED(negf_control))
139 4 : CALL timeset(routineN, handle)
140 :
141 4 : ALLOCATE (negf_control)
142 :
143 4 : CALL timestop(handle)
144 4 : END SUBROUTINE negf_control_create
145 :
146 : ! **************************************************************************************************
147 : !> \brief release memory allocated for NEGF control options
148 : !> \param negf_control an object to release
149 : !> \par History
150 : !> * 02.2017 created [Sergey Chulkov]
151 : ! **************************************************************************************************
152 4 : SUBROUTINE negf_control_release(negf_control)
153 : TYPE(negf_control_type), POINTER :: negf_control
154 :
155 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_control_release'
156 :
157 : INTEGER :: handle, i, j
158 :
159 4 : CALL timeset(routineN, handle)
160 :
161 4 : IF (ASSOCIATED(negf_control)) THEN
162 4 : IF (ALLOCATED(negf_control%atomlist_S)) DEALLOCATE (negf_control%atomlist_S)
163 4 : IF (ALLOCATED(negf_control%atomlist_S_screening)) DEALLOCATE (negf_control%atomlist_S_screening)
164 :
165 4 : IF (ALLOCATED(negf_control%contacts)) THEN
166 12 : DO i = SIZE(negf_control%contacts), 1, -1
167 8 : IF (ALLOCATED(negf_control%contacts(i)%atomlist_bulk)) &
168 8 : DEALLOCATE (negf_control%contacts(i)%atomlist_bulk)
169 :
170 8 : IF (ALLOCATED(negf_control%contacts(i)%atomlist_screening)) &
171 8 : DEALLOCATE (negf_control%contacts(i)%atomlist_screening)
172 :
173 12 : IF (ALLOCATED(negf_control%contacts(i)%atomlist_cell)) THEN
174 12 : DO j = SIZE(negf_control%contacts(i)%atomlist_cell), 1, -1
175 8 : IF (ALLOCATED(negf_control%contacts(i)%atomlist_cell(j)%vector)) &
176 12 : DEALLOCATE (negf_control%contacts(i)%atomlist_cell(j)%vector)
177 : END DO
178 12 : DEALLOCATE (negf_control%contacts(i)%atomlist_cell)
179 : END IF
180 : END DO
181 :
182 12 : DEALLOCATE (negf_control%contacts)
183 : END IF
184 :
185 4 : DEALLOCATE (negf_control)
186 : END IF
187 :
188 4 : CALL timestop(handle)
189 4 : END SUBROUTINE negf_control_release
190 :
191 : ! **************************************************************************************************
192 : !> \brief Read NEGF input parameters.
193 : !> \param negf_control NEGF control parameters
194 : !> \param input root input section
195 : !> \param subsys subsystem environment
196 : ! **************************************************************************************************
197 4 : SUBROUTINE read_negf_control(negf_control, input, subsys)
198 : TYPE(negf_control_type), POINTER :: negf_control
199 : TYPE(section_vals_type), POINTER :: input
200 : TYPE(cp_subsys_type), POINTER :: subsys
201 :
202 : CHARACTER(len=*), PARAMETER :: routineN = 'read_negf_control'
203 :
204 : CHARACTER(len=default_string_length) :: contact_id_str, eta_current_str, eta_max_str, &
205 : npoles_current_str, npoles_min_str, temp_current_str, temp_min_str
206 : INTEGER :: delta_npoles_min, handle, i2_rep, i_rep, &
207 : n2_rep, n_rep, natoms_current, &
208 : natoms_total, run_type
209 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
210 : LOGICAL :: do_negf, is_explicit
211 : REAL(kind=dp) :: eta_max, temp_current, temp_min
212 : TYPE(section_vals_type), POINTER :: cell_section, contact_section, &
213 : negf_section, region_section
214 :
215 4 : CALL timeset(routineN, handle)
216 :
217 4 : CALL section_vals_val_get(input, "GLOBAL%RUN_TYPE", i_val=run_type)
218 4 : do_negf = run_type == negf_run
219 :
220 4 : negf_section => section_vals_get_subs_vals(input, "NEGF")
221 :
222 4 : contact_section => section_vals_get_subs_vals(negf_section, "CONTACT")
223 4 : CALL section_vals_get(contact_section, n_repetition=n_rep, explicit=is_explicit)
224 4 : IF ((.NOT. is_explicit) .AND. do_negf) THEN
225 : CALL cp_abort(__LOCATION__, &
226 0 : "At least one contact is needed for NEGF calculation.")
227 : END IF
228 :
229 20 : ALLOCATE (negf_control%contacts(n_rep))
230 12 : DO i_rep = 1, n_rep
231 8 : region_section => section_vals_get_subs_vals(contact_section, "SCREENING_REGION", i_rep_section=i_rep)
232 8 : CALL section_vals_get(region_section, explicit=is_explicit)
233 :
234 8 : IF ((.NOT. is_explicit) .AND. do_negf) THEN
235 0 : WRITE (contact_id_str, '(I11)') i_rep
236 : CALL cp_abort(__LOCATION__, &
237 0 : "The screening region must be defined for the contact "//TRIM(ADJUSTL(contact_id_str))//".")
238 : END IF
239 :
240 8 : IF (is_explicit) THEN
241 8 : CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_screening, region_section, 1, subsys)
242 : END IF
243 :
244 8 : region_section => section_vals_get_subs_vals(contact_section, "BULK_REGION", i_rep_section=i_rep)
245 :
246 8 : CALL section_vals_get(region_section, explicit=is_explicit)
247 :
248 8 : IF ((.NOT. is_explicit) .AND. do_negf) THEN
249 0 : WRITE (contact_id_str, '(I11)') i_rep
250 : CALL cp_abort(__LOCATION__, &
251 0 : "The bulk region must be defined for the contact "//TRIM(ADJUSTL(contact_id_str))//".")
252 : END IF
253 :
254 8 : IF (is_explicit) THEN
255 8 : CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_bulk, region_section, 1, subsys)
256 : END IF
257 :
258 : CALL section_vals_val_get(contact_section, "FORCE_EVAL_SECTION", &
259 : i_val=negf_control%contacts(i_rep)%force_env_index, &
260 8 : i_rep_section=i_rep)
261 :
262 8 : cell_section => section_vals_get_subs_vals(region_section, "CELL")
263 8 : CALL section_vals_get(cell_section, n_repetition=n2_rep, explicit=is_explicit)
264 :
265 8 : IF (((.NOT. is_explicit) .OR. n2_rep /= 2) .AND. negf_control%contacts(i_rep)%force_env_index <= 0 .AND. do_negf) THEN
266 0 : WRITE (contact_id_str, '(I11)') i_rep
267 : CALL cp_abort(__LOCATION__, &
268 : "You must either provide indices of atoms belonging to two adjacent bulk unit cells "// &
269 : "(BULK_REGION/CELL) for the contact, or the index of the FORCE_EVAL section (FORCE_EVAL_SECTION) "// &
270 : "which will be used to construct Kohn-Sham matrix for the bulk contact "// &
271 0 : TRIM(ADJUSTL(contact_id_str))//".")
272 : END IF
273 :
274 8 : IF (is_explicit .AND. n2_rep > 0) THEN
275 20 : ALLOCATE (negf_control%contacts(i_rep)%atomlist_cell(n2_rep))
276 :
277 12 : DO i2_rep = 1, n2_rep
278 12 : CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_cell(i2_rep)%vector, cell_section, i2_rep, subsys)
279 : END DO
280 : END IF
281 :
282 : CALL section_vals_val_get(contact_section, "REFINE_FERMI_LEVEL", &
283 : l_val=negf_control%contacts(i_rep)%refine_fermi_level, &
284 8 : i_rep_section=i_rep)
285 :
286 : CALL section_vals_val_get(contact_section, "FERMI_LEVEL", &
287 : r_val=negf_control%contacts(i_rep)%fermi_level, &
288 8 : i_rep_section=i_rep, explicit=is_explicit)
289 : negf_control%contacts(i_rep)%compute_fermi_level = (.NOT. is_explicit) .OR. &
290 8 : negf_control%contacts(i_rep)%refine_fermi_level
291 :
292 8 : IF (do_negf .AND. negf_control%contacts(i_rep)%force_env_index <= 0 .AND. &
293 : (.NOT. (is_explicit .OR. negf_control%contacts(i_rep)%refine_fermi_level))) THEN
294 0 : WRITE (contact_id_str, '(I11)') i_rep
295 : CALL cp_warn(__LOCATION__, &
296 : "There is no way to reasonably guess the Fermi level for the bulk contact "// &
297 : TRIM(ADJUSTL(contact_id_str))//" without running a separate bulk DFT calculation first. "// &
298 : "Therefore, 0.0 Hartree will be used as an initial guess. It is strongly advised to enable "// &
299 : "the REFINE_FERMI_LEVEL switch and to provide an initial guess using the FERMI_LEVEL keyword. "// &
300 0 : "Alternatively, a bulk FORCE_EVAL_SECTION can be set up.")
301 : END IF
302 :
303 : CALL section_vals_val_get(contact_section, "TEMPERATURE", &
304 : r_val=negf_control%contacts(i_rep)%temperature, &
305 8 : i_rep_section=i_rep)
306 8 : IF (negf_control%contacts(i_rep)%temperature <= 0.0_dp) THEN
307 0 : CALL cp_abort(__LOCATION__, "Electronic temperature must be > 0")
308 : END IF
309 :
310 : CALL section_vals_val_get(contact_section, "ELECTRIC_POTENTIAL", &
311 : r_val=negf_control%contacts(i_rep)%v_external, &
312 44 : i_rep_section=i_rep)
313 : END DO
314 :
315 4 : region_section => section_vals_get_subs_vals(negf_section, "SCATTERING_REGION")
316 4 : CALL section_vals_get(region_section, explicit=is_explicit)
317 4 : IF (is_explicit) THEN
318 4 : CALL read_negf_atomlist(negf_control%atomlist_S, region_section, 1, subsys)
319 : END IF
320 :
321 4 : CALL section_vals_val_get(negf_section, "DISABLE_CACHE", l_val=negf_control%disable_cache)
322 :
323 4 : CALL section_vals_val_get(negf_section, "EPS_DENSITY", r_val=negf_control%conv_density)
324 4 : CALL section_vals_val_get(negf_section, "EPS_GREEN", r_val=negf_control%conv_green)
325 4 : CALL section_vals_val_get(negf_section, "EPS_SCF", r_val=negf_control%conv_scf)
326 :
327 4 : CALL section_vals_val_get(negf_section, "EPS_GEO", r_val=negf_control%eps_geometry)
328 :
329 4 : CALL section_vals_val_get(negf_section, "ENERGY_LBOUND", r_val=negf_control%energy_lbound)
330 4 : CALL section_vals_val_get(negf_section, "ETA", r_val=negf_control%eta)
331 4 : CALL section_vals_val_get(negf_section, "HOMO_LUMO_GAP", r_val=negf_control%homo_lumo_gap)
332 4 : CALL section_vals_val_get(negf_section, "DELTA_NPOLES", i_val=negf_control%delta_npoles)
333 4 : CALL section_vals_val_get(negf_section, "GAMMA_KT", i_val=negf_control%gamma_kT)
334 :
335 4 : CALL section_vals_val_get(negf_section, "INTEGRATION_METHOD", i_val=negf_control%integr_method)
336 4 : CALL section_vals_val_get(negf_section, "INTEGRATION_MIN_POINTS", i_val=negf_control%integr_min_points)
337 4 : CALL section_vals_val_get(negf_section, "INTEGRATION_MAX_POINTS", i_val=negf_control%integr_max_points)
338 :
339 4 : IF (negf_control%integr_max_points < negf_control%integr_min_points) &
340 0 : negf_control%integr_max_points = negf_control%integr_min_points
341 :
342 4 : CALL section_vals_val_get(negf_section, "MAX_SCF", i_val=negf_control%max_scf)
343 :
344 4 : CALL section_vals_val_get(negf_section, "NPROC_POINT", i_val=negf_control%nprocs)
345 :
346 4 : CALL section_vals_val_get(negf_section, "V_SHIFT", r_val=negf_control%v_shift)
347 4 : CALL section_vals_val_get(negf_section, "V_SHIFT_OFFSET", r_val=negf_control%v_shift_offset)
348 4 : CALL section_vals_val_get(negf_section, "V_SHIFT_MAX_ITERS", i_val=negf_control%v_shift_maxiters)
349 :
350 : ! check consistency
351 4 : IF (negf_control%eta < 0.0_dp) THEN
352 0 : CALL cp_abort(__LOCATION__, "ETA must be >= 0")
353 : END IF
354 :
355 4 : IF (n_rep > 0) THEN
356 16 : delta_npoles_min = NINT(0.5_dp*(negf_control%eta/(pi*MAXVAL(negf_control%contacts(:)%temperature)) + 1.0_dp))
357 : ELSE
358 0 : delta_npoles_min = 1
359 : END IF
360 :
361 4 : IF (negf_control%delta_npoles < delta_npoles_min) THEN
362 0 : IF (n_rep > 0) THEN
363 0 : eta_max = REAL(2*negf_control%delta_npoles - 1, kind=dp)*pi*MAXVAL(negf_control%contacts(:)%temperature)
364 0 : temp_current = MAXVAL(negf_control%contacts(:)%temperature)*kelvin
365 0 : temp_min = negf_control%eta/(pi*REAL(2*negf_control%delta_npoles - 1, kind=dp))*kelvin
366 :
367 0 : WRITE (eta_current_str, '(ES11.4E2)') negf_control%eta
368 0 : WRITE (eta_max_str, '(ES11.4E2)') eta_max
369 0 : WRITE (npoles_current_str, '(I11)') negf_control%delta_npoles
370 0 : WRITE (npoles_min_str, '(I11)') delta_npoles_min
371 0 : WRITE (temp_current_str, '(F11.3)') temp_current
372 0 : WRITE (temp_min_str, '(F11.3)') temp_min
373 :
374 : CALL cp_abort(__LOCATION__, &
375 : "Parameter DELTA_NPOLES must be at least "//TRIM(ADJUSTL(npoles_min_str))// &
376 : " (instead of "//TRIM(ADJUSTL(npoles_current_str))// &
377 : ") for given TEMPERATURE ("//TRIM(ADJUSTL(temp_current_str))// &
378 : " K) and ETA ("//TRIM(ADJUSTL(eta_current_str))// &
379 : "). Alternatively you can increase TEMPERATURE above "//TRIM(ADJUSTL(temp_min_str))// &
380 : " K, or decrease ETA below "//TRIM(ADJUSTL(eta_max_str))// &
381 : ". Please keep in mind that very tight ETA may result in dramatical precision loss"// &
382 0 : " due to inversion of ill-conditioned matrices.")
383 : ELSE
384 : ! no leads have been defined, so calculation will abort anyway
385 0 : negf_control%delta_npoles = delta_npoles_min
386 : END IF
387 : END IF
388 :
389 : ! expand scattering region by adding atoms from contact screening regions
390 4 : n_rep = SIZE(negf_control%contacts)
391 4 : IF (ALLOCATED(negf_control%atomlist_S)) THEN
392 4 : natoms_total = SIZE(negf_control%atomlist_S)
393 : ELSE
394 0 : natoms_total = 0
395 : END IF
396 :
397 12 : DO i_rep = 1, n_rep
398 12 : IF (ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) THEN
399 : IF (ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) &
400 8 : natoms_total = natoms_total + SIZE(negf_control%contacts(i_rep)%atomlist_screening)
401 : END IF
402 : END DO
403 :
404 4 : IF (natoms_total > 0) THEN
405 12 : ALLOCATE (negf_control%atomlist_S_screening(natoms_total))
406 4 : IF (ALLOCATED(negf_control%atomlist_S)) THEN
407 4 : natoms_total = SIZE(negf_control%atomlist_S)
408 20 : negf_control%atomlist_S_screening(1:natoms_total) = negf_control%atomlist_S(1:natoms_total)
409 : ELSE
410 0 : natoms_total = 0
411 : END IF
412 :
413 12 : DO i_rep = 1, n_rep
414 12 : IF (ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) THEN
415 8 : natoms_current = SIZE(negf_control%contacts(i_rep)%atomlist_screening)
416 :
417 : negf_control%atomlist_S_screening(natoms_total + 1:natoms_total + natoms_current) = &
418 40 : negf_control%contacts(i_rep)%atomlist_screening(1:natoms_current)
419 :
420 8 : natoms_total = natoms_total + natoms_current
421 : END IF
422 : END DO
423 :
424 : ! sort and remove duplicated atoms
425 12 : ALLOCATE (inds(natoms_total))
426 4 : CALL sort(negf_control%atomlist_S_screening, natoms_total, inds)
427 4 : DEALLOCATE (inds)
428 :
429 4 : natoms_current = 1
430 48 : DO i_rep = natoms_current + 1, natoms_total
431 48 : IF (negf_control%atomlist_S_screening(i_rep) /= negf_control%atomlist_S_screening(natoms_current)) THEN
432 44 : natoms_current = natoms_current + 1
433 44 : negf_control%atomlist_S_screening(natoms_current) = negf_control%atomlist_S_screening(i_rep)
434 : END IF
435 : END DO
436 :
437 4 : IF (natoms_current < natoms_total) THEN
438 0 : CALL MOVE_ALLOC(negf_control%atomlist_S_screening, inds)
439 :
440 0 : ALLOCATE (negf_control%atomlist_S_screening(natoms_current))
441 0 : negf_control%atomlist_S_screening(1:natoms_current) = inds(1:natoms_current)
442 0 : DEALLOCATE (inds)
443 : END IF
444 : END IF
445 :
446 4 : IF (do_negf .AND. SIZE(negf_control%contacts) > 2) THEN
447 : CALL cp_abort(__LOCATION__, &
448 0 : "General case (> 2 contacts) has not been implemented yet")
449 : END IF
450 :
451 4 : CALL timestop(handle)
452 12 : END SUBROUTINE read_negf_control
453 :
454 : ! **************************************************************************************************
455 : !> \brief Read region-specific list of atoms.
456 : !> \param atomlist list of atoms
457 : !> \param input_section input section which contains 'LIST' and 'MOLNAME' keywords
458 : !> \param i_rep_section repetition index of the input_section
459 : !> \param subsys subsystem environment
460 : ! **************************************************************************************************
461 28 : SUBROUTINE read_negf_atomlist(atomlist, input_section, i_rep_section, subsys)
462 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(out) :: atomlist
463 : TYPE(section_vals_type), POINTER :: input_section
464 : INTEGER, INTENT(in) :: i_rep_section
465 : TYPE(cp_subsys_type), POINTER :: subsys
466 :
467 : CHARACTER(len=*), PARAMETER :: routineN = 'read_negf_atomlist'
468 :
469 : CHARACTER(len=default_string_length) :: index_str, natoms_str
470 : CHARACTER(len=default_string_length), &
471 28 : DIMENSION(:), POINTER :: cptr
472 : INTEGER :: first_atom, handle, iatom, ikind, imol, iname, irep, last_atom, natoms_current, &
473 : natoms_max, natoms_total, nkinds, nmols, nnames, nrep_list, nrep_molname
474 28 : INTEGER, ALLOCATABLE, DIMENSION(:) :: inds
475 28 : INTEGER, DIMENSION(:), POINTER :: iptr
476 : LOGICAL :: is_list, is_molname
477 28 : TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
478 : TYPE(molecule_kind_type), POINTER :: molecule_kind
479 28 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
480 : TYPE(molecule_type), POINTER :: molecule
481 28 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
482 :
483 28 : CALL timeset(routineN, handle)
484 :
485 : CALL cp_subsys_get(subsys, particle_set=particle_set, &
486 : molecule_set=molecule_set, &
487 28 : molecule_kind_set=molecule_kind_set)
488 28 : natoms_max = SIZE(particle_set)
489 28 : nkinds = SIZE(molecule_kind_set)
490 :
491 : CALL section_vals_val_get(input_section, "LIST", i_rep_section=i_rep_section, &
492 28 : n_rep_val=nrep_list, explicit=is_list)
493 : CALL section_vals_val_get(input_section, "MOLNAME", i_rep_section=i_rep_section, &
494 28 : n_rep_val=nrep_molname, explicit=is_molname)
495 :
496 : ! compute the number of atoms in the NEGF region, and check the validity of giben atomic indices
497 28 : natoms_total = 0
498 28 : IF (is_list .AND. nrep_list > 0) THEN
499 0 : DO irep = 1, nrep_list
500 0 : CALL section_vals_val_get(input_section, "LIST", i_rep_section=i_rep_section, i_rep_val=irep, i_vals=iptr)
501 :
502 0 : natoms_current = SIZE(iptr)
503 0 : DO iatom = 1, natoms_current
504 0 : IF (iptr(iatom) > natoms_max) THEN
505 0 : CALL integer_to_string(iptr(iatom), index_str)
506 0 : CALL integer_to_string(natoms_max, natoms_str)
507 : CALL cp_abort(__LOCATION__, &
508 : "NEGF: Atomic index "//TRIM(index_str)//" given in section "// &
509 : TRIM(input_section%section%name)//" exceeds the maximum number of atoms ("// &
510 0 : TRIM(natoms_str)//").")
511 : END IF
512 : END DO
513 :
514 0 : natoms_total = natoms_total + natoms_current
515 : END DO
516 : END IF
517 :
518 28 : IF (is_molname .AND. nrep_molname > 0) THEN
519 56 : DO irep = 1, nrep_molname
520 28 : CALL section_vals_val_get(input_section, "MOLNAME", i_rep_section=i_rep_section, i_rep_val=irep, c_vals=cptr)
521 28 : nnames = SIZE(cptr)
522 :
523 92 : DO iname = 1, nnames
524 144 : DO ikind = 1, nkinds
525 144 : IF (molecule_kind_set(ikind)%name .EQ. cptr(iname)) EXIT
526 : END DO
527 :
528 64 : IF (ikind <= nkinds) THEN
529 36 : molecule_kind => molecule_kind_set(ikind)
530 36 : CALL get_molecule_kind(molecule_kind, nmolecule=nmols, molecule_list=iptr)
531 :
532 180 : DO imol = 1, nmols
533 144 : molecule => molecule_set(iptr(imol))
534 144 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
535 144 : natoms_current = last_atom - first_atom + 1
536 180 : natoms_total = natoms_total + natoms_current
537 : END DO
538 : ELSE
539 : CALL cp_abort(__LOCATION__, &
540 : "NEGF: A molecule with the name '"//TRIM(cptr(iname))//"' mentioned in section "// &
541 0 : TRIM(input_section%section%name)//" has not been defined. Note that names are case sensitive.")
542 : END IF
543 : END DO
544 : END DO
545 : END IF
546 :
547 : ! create a list of atomic indices
548 28 : IF (natoms_total > 0) THEN
549 84 : ALLOCATE (atomlist(natoms_total))
550 :
551 28 : natoms_total = 0
552 :
553 28 : IF (is_list .AND. nrep_list > 0) THEN
554 0 : DO irep = 1, nrep_list
555 0 : CALL section_vals_val_get(input_section, "LIST", i_rep_section=i_rep_section, i_rep_val=irep, i_vals=iptr)
556 :
557 0 : natoms_current = SIZE(iptr)
558 0 : atomlist(natoms_total + 1:natoms_total + natoms_current) = iptr(1:natoms_current)
559 0 : natoms_total = natoms_total + natoms_current
560 : END DO
561 : END IF
562 :
563 28 : IF (is_molname .AND. nrep_molname > 0) THEN
564 56 : DO irep = 1, nrep_molname
565 28 : CALL section_vals_val_get(input_section, "MOLNAME", i_rep_section=i_rep_section, i_rep_val=irep, c_vals=cptr)
566 28 : nnames = SIZE(cptr)
567 :
568 92 : DO iname = 1, nnames
569 144 : DO ikind = 1, nkinds
570 144 : IF (molecule_kind_set(ikind)%name .EQ. cptr(iname)) EXIT
571 : END DO
572 :
573 64 : IF (ikind <= nkinds) THEN
574 36 : molecule_kind => molecule_kind_set(ikind)
575 36 : CALL get_molecule_kind(molecule_kind, nmolecule=nmols, molecule_list=iptr)
576 :
577 180 : DO imol = 1, nmols
578 144 : molecule => molecule_set(iptr(imol))
579 144 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
580 :
581 324 : DO natoms_current = first_atom, last_atom
582 144 : natoms_total = natoms_total + 1
583 288 : atomlist(natoms_total) = natoms_current
584 : END DO
585 : END DO
586 : END IF
587 : END DO
588 : END DO
589 : END IF
590 :
591 : ! remove duplicated atoms
592 84 : ALLOCATE (inds(natoms_total))
593 28 : CALL sort(atomlist, natoms_total, inds)
594 28 : DEALLOCATE (inds)
595 :
596 28 : natoms_current = 1
597 144 : DO iatom = natoms_current + 1, natoms_total
598 144 : IF (atomlist(iatom) /= atomlist(natoms_current)) THEN
599 116 : natoms_current = natoms_current + 1
600 116 : atomlist(natoms_current) = atomlist(iatom)
601 : END IF
602 : END DO
603 :
604 28 : IF (natoms_current < natoms_total) THEN
605 0 : CALL MOVE_ALLOC(atomlist, inds)
606 :
607 0 : ALLOCATE (atomlist(natoms_current))
608 0 : atomlist(1:natoms_current) = inds(1:natoms_current)
609 0 : DEALLOCATE (inds)
610 : END IF
611 : END IF
612 :
613 28 : CALL timestop(handle)
614 28 : END SUBROUTINE read_negf_atomlist
615 0 : END MODULE negf_control_types
|
