32#include "./base/base_uses.f90"
37 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'negf_control_types'
38 LOGICAL,
PARAMETER,
PRIVATE :: debug_this_module = .true.
49 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atomlist_bulk, atomlist_screening
52 DIMENSION(:) :: atomlist_cell
54 INTEGER :: force_env_index = -1
56 LOGICAL :: compute_fermi_level = .false.
58 LOGICAL :: refine_fermi_level = .false.
60 LOGICAL :: shift_fermi_level = .false.
62 REAL(kind=
dp) :: fermi_level = -1.0_dp
64 REAL(kind=
dp) :: fermi_level_shifted = -1.0_dp
66 REAL(kind=
dp) :: temperature = -1.0_dp
68 REAL(kind=
dp) :: v_external = 0.0_dp
78 DIMENSION(:) :: contacts
80 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atomlist_s
83 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atomlist_s_screening
85 LOGICAL :: disable_cache = .false.
87 REAL(kind=
dp) :: conv_density = -1.0_dp
89 REAL(kind=
dp) :: conv_green = -1.0_dp
91 REAL(kind=
dp) :: conv_scf = -1.0_dp
93 REAL(kind=
dp) :: eps_geometry = -1.0_dp
95 REAL(kind=
dp) :: v_bias = -1.0_dp
97 REAL(kind=
dp) :: energy_lbound = -1.0_dp
99 REAL(kind=
dp) :: eta = -1.0_dp
101 REAL(kind=
dp) :: homo_lumo_gap = -1.0_dp
103 INTEGER :: delta_npoles = -1
105 INTEGER :: gamma_kt = -1
107 INTEGER :: integr_method = -1
109 INTEGER :: integr_min_points = -1
111 INTEGER :: integr_max_points = -1
113 INTEGER :: max_scf = -1
115 INTEGER :: nprocs = -1
117 REAL(kind=
dp) :: v_shift = -1.0_dp
119 REAL(kind=
dp) :: v_shift_offset = -1.0_dp
121 INTEGER :: v_shift_maxiters = -1
124 PRIVATE :: read_negf_atomlist
137 CHARACTER(len=*),
PARAMETER :: routinen =
'negf_control_create'
141 cpassert(.NOT.
ASSOCIATED(negf_control))
142 CALL timeset(routinen, handle)
144 ALLOCATE (negf_control)
146 CALL timestop(handle)
158 CHARACTER(len=*),
PARAMETER :: routinen =
'negf_control_release'
160 INTEGER :: handle, i, j
162 CALL timeset(routinen, handle)
164 IF (
ASSOCIATED(negf_control))
THEN
165 IF (
ALLOCATED(negf_control%atomlist_S))
DEALLOCATE (negf_control%atomlist_S)
166 IF (
ALLOCATED(negf_control%atomlist_S_screening))
DEALLOCATE (negf_control%atomlist_S_screening)
168 IF (
ALLOCATED(negf_control%contacts))
THEN
169 DO i =
SIZE(negf_control%contacts), 1, -1
170 IF (
ALLOCATED(negf_control%contacts(i)%atomlist_bulk)) &
171 DEALLOCATE (negf_control%contacts(i)%atomlist_bulk)
173 IF (
ALLOCATED(negf_control%contacts(i)%atomlist_screening)) &
174 DEALLOCATE (negf_control%contacts(i)%atomlist_screening)
176 IF (
ALLOCATED(negf_control%contacts(i)%atomlist_cell))
THEN
177 DO j =
SIZE(negf_control%contacts(i)%atomlist_cell), 1, -1
178 IF (
ALLOCATED(negf_control%contacts(i)%atomlist_cell(j)%vector)) &
179 DEALLOCATE (negf_control%contacts(i)%atomlist_cell(j)%vector)
181 DEALLOCATE (negf_control%contacts(i)%atomlist_cell)
185 DEALLOCATE (negf_control%contacts)
188 DEALLOCATE (negf_control)
191 CALL timestop(handle)
205 CHARACTER(len=*),
PARAMETER :: routinen =
'read_negf_control'
207 CHARACTER(len=default_string_length) :: contact_id_str, eta_current_str, eta_max_str, &
208 npoles_current_str, npoles_min_str, temp_current_str, temp_min_str
209 INTEGER :: delta_npoles_min, handle, i2_rep, i_rep, &
210 n2_rep, n_rep, natoms_current, &
211 natoms_total, run_type
212 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: inds
213 LOGICAL :: do_negf, is_explicit
214 REAL(kind=
dp) :: eta_max, temp_current, temp_min
216 negf_section, region_section
218 CALL timeset(routinen, handle)
226 CALL section_vals_get(contact_section, n_repetition=n_rep, explicit=is_explicit)
227 IF ((.NOT. is_explicit) .AND. do_negf)
THEN
228 CALL cp_abort(__location__, &
229 "At least one contact is needed for NEGF calculation.")
232 ALLOCATE (negf_control%contacts(n_rep))
237 IF ((.NOT. is_explicit) .AND. do_negf)
THEN
238 WRITE (contact_id_str,
'(I11)') i_rep
239 CALL cp_abort(__location__, &
240 "The screening region must be defined for the contact "//trim(adjustl(contact_id_str))//
".")
243 IF (is_explicit)
THEN
244 CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_screening, region_section, 1, subsys)
251 IF ((.NOT. is_explicit) .AND. do_negf)
THEN
252 WRITE (contact_id_str,
'(I11)') i_rep
253 CALL cp_abort(__location__, &
254 "The bulk region must be defined for the contact "//trim(adjustl(contact_id_str))//
".")
257 IF (is_explicit)
THEN
258 CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_bulk, region_section, 1, subsys)
262 i_val=negf_control%contacts(i_rep)%force_env_index, &
266 CALL section_vals_get(cell_section, n_repetition=n2_rep, explicit=is_explicit)
268 IF (((.NOT. is_explicit) .OR. n2_rep /= 2) .AND. negf_control%contacts(i_rep)%force_env_index <= 0 .AND. do_negf)
THEN
269 WRITE (contact_id_str,
'(I11)') i_rep
270 CALL cp_abort(__location__, &
271 "You must either provide indices of atoms belonging to two adjacent bulk unit cells "// &
272 "(BULK_REGION/CELL) for the contact, or the index of the FORCE_EVAL section (FORCE_EVAL_SECTION) "// &
273 "which will be used to construct Kohn-Sham matrix for the bulk contact "// &
274 trim(adjustl(contact_id_str))//
".")
277 IF (is_explicit .AND. n2_rep > 0)
THEN
278 ALLOCATE (negf_control%contacts(i_rep)%atomlist_cell(n2_rep))
280 DO i2_rep = 1, n2_rep
281 CALL read_negf_atomlist(negf_control%contacts(i_rep)%atomlist_cell(i2_rep)%vector, cell_section, i2_rep, subsys)
286 l_val=negf_control%contacts(i_rep)%refine_fermi_level, &
290 r_val=negf_control%contacts(i_rep)%fermi_level, &
291 i_rep_section=i_rep, explicit=is_explicit)
292 IF (.NOT. is_explicit) negf_control%contacts(i_rep)%refine_fermi_level = .false.
293 negf_control%contacts(i_rep)%compute_fermi_level = (.NOT. is_explicit) .OR. &
294 negf_control%contacts(i_rep)%refine_fermi_level
297 r_val=negf_control%contacts(i_rep)%fermi_level_shifted, &
298 i_rep_section=i_rep, explicit=is_explicit)
299 IF (is_explicit) negf_control%contacts(i_rep)%shift_fermi_level = .true.
302 r_val=negf_control%contacts(i_rep)%temperature, &
304 IF (negf_control%contacts(i_rep)%temperature <= 0.0_dp)
THEN
305 CALL cp_abort(__location__,
"Electronic temperature must be > 0")
309 r_val=negf_control%contacts(i_rep)%v_external, &
315 IF (is_explicit)
THEN
316 CALL read_negf_atomlist(negf_control%atomlist_S, region_section, 1, subsys)
334 CALL section_vals_val_get(negf_section,
"INTEGRATION_MIN_POINTS", i_val=negf_control%integr_min_points)
335 CALL section_vals_val_get(negf_section,
"INTEGRATION_MAX_POINTS", i_val=negf_control%integr_max_points)
337 IF (negf_control%integr_max_points < negf_control%integr_min_points) &
338 negf_control%integr_max_points = negf_control%integr_min_points
346 CALL section_vals_val_get(negf_section,
"V_SHIFT_MAX_ITERS", i_val=negf_control%v_shift_maxiters)
349 IF (negf_control%eta < 0.0_dp)
THEN
350 CALL cp_abort(__location__,
"ETA must be >= 0")
354 delta_npoles_min = nint(0.5_dp*(negf_control%eta/(
pi*maxval(negf_control%contacts(:)%temperature)) + 1.0_dp))
359 IF (negf_control%delta_npoles < delta_npoles_min)
THEN
361 eta_max = real(2*negf_control%delta_npoles - 1, kind=
dp)*
pi*maxval(negf_control%contacts(:)%temperature)
362 temp_current = maxval(negf_control%contacts(:)%temperature)*
kelvin
363 temp_min = negf_control%eta/(
pi*real(2*negf_control%delta_npoles - 1, kind=
dp))*
kelvin
365 WRITE (eta_current_str,
'(ES11.4E2)') negf_control%eta
366 WRITE (eta_max_str,
'(ES11.4E2)') eta_max
367 WRITE (npoles_current_str,
'(I11)') negf_control%delta_npoles
368 WRITE (npoles_min_str,
'(I11)') delta_npoles_min
369 WRITE (temp_current_str,
'(F11.3)') temp_current
370 WRITE (temp_min_str,
'(F11.3)') temp_min
372 CALL cp_abort(__location__, &
373 "Parameter DELTA_NPOLES must be at least "//trim(adjustl(npoles_min_str))// &
374 " (instead of "//trim(adjustl(npoles_current_str))// &
375 ") for given TEMPERATURE ("//trim(adjustl(temp_current_str))// &
376 " K) and ETA ("//trim(adjustl(eta_current_str))// &
377 "). Alternatively you can increase TEMPERATURE above "//trim(adjustl(temp_min_str))// &
378 " K, or decrease ETA below "//trim(adjustl(eta_max_str))// &
379 ". Please keep in mind that very tight ETA may result in dramatical precision loss"// &
380 " due to inversion of ill-conditioned matrices.")
383 negf_control%delta_npoles = delta_npoles_min
388 n_rep =
SIZE(negf_control%contacts)
389 IF (
ALLOCATED(negf_control%atomlist_S))
THEN
390 natoms_total =
SIZE(negf_control%atomlist_S)
396 IF (
ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening))
THEN
397 IF (
ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening)) &
398 natoms_total = natoms_total +
SIZE(negf_control%contacts(i_rep)%atomlist_screening)
402 IF (natoms_total > 0)
THEN
403 ALLOCATE (negf_control%atomlist_S_screening(natoms_total))
404 IF (
ALLOCATED(negf_control%atomlist_S))
THEN
405 natoms_total =
SIZE(negf_control%atomlist_S)
406 negf_control%atomlist_S_screening(1:natoms_total) = negf_control%atomlist_S(1:natoms_total)
412 IF (
ALLOCATED(negf_control%contacts(i_rep)%atomlist_screening))
THEN
413 natoms_current =
SIZE(negf_control%contacts(i_rep)%atomlist_screening)
415 negf_control%atomlist_S_screening(natoms_total + 1:natoms_total + natoms_current) = &
416 negf_control%contacts(i_rep)%atomlist_screening(1:natoms_current)
418 natoms_total = natoms_total + natoms_current
423 ALLOCATE (inds(natoms_total))
424 CALL sort(negf_control%atomlist_S_screening, natoms_total, inds)
428 DO i_rep = natoms_current + 1, natoms_total
429 IF (negf_control%atomlist_S_screening(i_rep) /= negf_control%atomlist_S_screening(natoms_current))
THEN
430 natoms_current = natoms_current + 1
431 negf_control%atomlist_S_screening(natoms_current) = negf_control%atomlist_S_screening(i_rep)
435 IF (natoms_current < natoms_total)
THEN
436 CALL move_alloc(negf_control%atomlist_S_screening, inds)
438 ALLOCATE (negf_control%atomlist_S_screening(natoms_current))
439 negf_control%atomlist_S_screening(1:natoms_current) = inds(1:natoms_current)
444 IF (do_negf .AND.
SIZE(negf_control%contacts) > 2)
THEN
445 CALL cp_abort(__location__, &
446 "General case (> 2 contacts) has not been implemented yet")
449 CALL timestop(handle)
459 SUBROUTINE read_negf_atomlist(atomlist, input_section, i_rep_section, subsys)
460 INTEGER,
ALLOCATABLE,
DIMENSION(:),
INTENT(out) :: atomlist
462 INTEGER,
INTENT(in) :: i_rep_section
465 CHARACTER(len=*),
PARAMETER :: routinen =
'read_negf_atomlist'
467 CHARACTER(len=default_string_length) :: index_str, natoms_str
468 CHARACTER(len=default_string_length), &
469 DIMENSION(:),
POINTER :: cptr
470 INTEGER :: first_atom, handle, iatom, ikind, imol, iname, irep, last_atom, natoms_current, &
471 natoms_max, natoms_total, nkinds, nmols, nnames, nrep_list, nrep_molname
472 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: inds
473 INTEGER,
DIMENSION(:),
POINTER :: iptr
474 LOGICAL :: is_list, is_molname
481 CALL timeset(routinen, handle)
484 molecule_set=molecule_set, &
485 molecule_kind_set=molecule_kind_set)
486 natoms_max =
SIZE(particle_set)
487 nkinds =
SIZE(molecule_kind_set)
490 n_rep_val=nrep_list, explicit=is_list)
492 n_rep_val=nrep_molname, explicit=is_molname)
496 IF (is_list .AND. nrep_list > 0)
THEN
497 DO irep = 1, nrep_list
498 CALL section_vals_val_get(input_section,
"LIST", i_rep_section=i_rep_section, i_rep_val=irep, i_vals=iptr)
500 natoms_current =
SIZE(iptr)
501 DO iatom = 1, natoms_current
502 IF (iptr(iatom) > natoms_max)
THEN
505 CALL cp_abort(__location__, &
506 "NEGF: Atomic index "//trim(index_str)//
" given in section "// &
507 trim(input_section%section%name)//
" exceeds the maximum number of atoms ("// &
508 trim(natoms_str)//
").")
512 natoms_total = natoms_total + natoms_current
516 IF (is_molname .AND. nrep_molname > 0)
THEN
517 DO irep = 1, nrep_molname
518 CALL section_vals_val_get(input_section,
"MOLNAME", i_rep_section=i_rep_section, i_rep_val=irep, c_vals=cptr)
523 IF (molecule_kind_set(ikind)%name == cptr(iname))
EXIT
526 IF (ikind <= nkinds)
THEN
527 molecule_kind => molecule_kind_set(ikind)
531 molecule => molecule_set(iptr(imol))
532 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
533 natoms_current = last_atom - first_atom + 1
534 natoms_total = natoms_total + natoms_current
537 CALL cp_abort(__location__, &
538 "NEGF: A molecule with the name '"//trim(cptr(iname))//
"' mentioned in section "// &
539 trim(input_section%section%name)//
" has not been defined. Note that names are case sensitive.")
546 IF (natoms_total > 0)
THEN
547 ALLOCATE (atomlist(natoms_total))
551 IF (is_list .AND. nrep_list > 0)
THEN
552 DO irep = 1, nrep_list
553 CALL section_vals_val_get(input_section,
"LIST", i_rep_section=i_rep_section, i_rep_val=irep, i_vals=iptr)
555 natoms_current =
SIZE(iptr)
556 atomlist(natoms_total + 1:natoms_total + natoms_current) = iptr(1:natoms_current)
557 natoms_total = natoms_total + natoms_current
561 IF (is_molname .AND. nrep_molname > 0)
THEN
562 DO irep = 1, nrep_molname
563 CALL section_vals_val_get(input_section,
"MOLNAME", i_rep_section=i_rep_section, i_rep_val=irep, c_vals=cptr)
568 IF (molecule_kind_set(ikind)%name == cptr(iname))
EXIT
571 IF (ikind <= nkinds)
THEN
572 molecule_kind => molecule_kind_set(ikind)
576 molecule => molecule_set(iptr(imol))
577 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
579 DO natoms_current = first_atom, last_atom
580 natoms_total = natoms_total + 1
581 atomlist(natoms_total) = natoms_current
590 ALLOCATE (inds(natoms_total))
591 CALL sort(atomlist, natoms_total, inds)
595 DO iatom = natoms_current + 1, natoms_total
596 IF (atomlist(iatom) /= atomlist(natoms_current))
THEN
597 natoms_current = natoms_current + 1
598 atomlist(natoms_current) = atomlist(iatom)
602 IF (natoms_current < natoms_total)
THEN
603 CALL move_alloc(atomlist, inds)
605 ALLOCATE (atomlist(natoms_current))
606 atomlist(1:natoms_current) = inds(1:natoms_current)
611 CALL timestop(handle)
612 END SUBROUTINE read_negf_atomlist
types that represent a subsys, i.e. a part of the system
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
Define the molecule kind structure types and the corresponding functionality.
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Define the data structure for the molecule information.
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Allocatable vectors for NEGF based quantum transport calculations.
Input control types for NEGF based quantum transport calculations.
subroutine, public negf_control_create(negf_control)
allocate control options for Non-equilibrium Green's Function calculation
subroutine, public read_negf_control(negf_control, input, subsys)
Read NEGF input parameters.
subroutine, public negf_control_release(negf_control)
release memory allocated for NEGF control options
Define the data structure for the particle information.
Definition of physical constants:
real(kind=dp), parameter, public kelvin
Utilities for string manipulations.
subroutine, public integer_to_string(inumber, string)
Converts an integer number to a string. The WRITE statement will return an error message,...
All kind of helpful little routines.
represents a system: atoms, molecules, their pos,vel,...
Allocatable 1-D integer vector.
Input parameters related to the NEGF run.