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