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