From a80b2e66393cf1cf841f561e07238061d72bca26 Mon Sep 17 00:00:00 2001 From: Sergey Chulkov Date: Thu, 10 May 2018 12:11:54 +0000 Subject: [PATCH] NEGF: contact-specific temperature, correct shift and scale factors svn-origin-rev: 18419 --- src/input_cp2k_negf.F | 45 ++- src/negf_control_types.F | 78 ++-- src/negf_env_types.F | 166 ++++++-- src/negf_green_methods.F | 9 +- src/negf_methods.F | 584 +++++++++++++++++++++++------ src/negf_qs_env_methods.F | 190 ---------- src/negf_qs_env_types.F | 151 -------- src/negf_vectors.F | 31 +- src/qs_environment_types.F | 23 +- src/qs_ks_methods.F | 18 +- tests/QS/regtest-negf/hchain.inp | 341 ----------------- tests/QS/regtest-negf/li_chain.inp | 376 +++++++++++++++++++ 12 files changed, 1105 insertions(+), 907 deletions(-) delete mode 100644 src/negf_qs_env_methods.F delete mode 100644 src/negf_qs_env_types.F delete mode 100644 tests/QS/regtest-negf/hchain.inp create mode 100644 tests/QS/regtest-negf/li_chain.inp diff --git a/src/input_cp2k_negf.F b/src/input_cp2k_negf.F index 6997e8eaf4..bac4e508cb 100644 --- a/src/input_cp2k_negf.F +++ b/src/input_cp2k_negf.F @@ -77,6 +77,12 @@ CONTAINS CALL section_add_subsection(section, subsection) CALL section_release(subsection) + CALL keyword_create(keyword, name="DISABLE_CACHE", & + description="Do not keep contact self-energy matrices for future reuse", & + default_l_val=.FALSE., lone_keyword_l_val=.TRUE.) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + ! convergence thresholds CALL keyword_create(keyword, name="EPS_DENSITY", & description="Target accuracy for electronic density.", & @@ -103,13 +109,6 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) - CALL keyword_create(keyword, name="TEMPERATURE", & - description="Electronic temperature.", & - n_var=1, type_of_var=real_t, unit_str="K", & - default_r_val=300.0_dp/kelvin) - CALL section_add_keyword(section, keyword) - CALL keyword_release(keyword) - CALL keyword_create(keyword, name="ENERGY_LBOUND", & description="Lower bound energy of the conductance band.", & n_var=1, type_of_var=real_t, unit_str="hartree", & @@ -125,8 +124,8 @@ CONTAINS CALL keyword_release(keyword) CALL keyword_create(keyword, name="HOMO_LUMO_GAP", & - description="The gap between the HOMO and some fictitions LUMO. This option is used as"// & - " an initial offset in determine the actual Fermi level of bulk contacts."// & + description="The gap between the HOMO and some fictitious LUMO. This option is used as"// & + " an initial offset to determine the actual Fermi level of bulk contacts."// & " It does not need to be exact HOMO-LUMO gap, just some value to start with.", & n_var=1, type_of_var=real_t, unit_str="hartree", & default_r_val=0.2_dp) @@ -188,6 +187,25 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, name="V_SHIFT", & + description="Initial value of the Hartree potential shift", & + n_var=1, type_of_var=real_t, unit_str="hartree", & + default_r_val=0.0_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="V_SHIFT_OFFSET", & + description="Initial offset to determine the optimal shift in Hartree potential.", & + n_var=1, type_of_var=real_t, default_r_val=0.10_dp) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + + CALL keyword_create(keyword, name="V_SHIFT_MAX_ITERS", & + description="Maximal number of iteration to determine the optimal shift in Hartree potential.", & + n_var=1, type_of_var=integer_t, default_i_val=30) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + ! PRINT subsection CALL section_create(subsection, "PRINT", "Print properties for the scattering region.", & repeats=.FALSE.) @@ -224,7 +242,7 @@ CONTAINS CALL section_create(section, "CONTACT", & description="Section defining the contact region of NEGF setup.", & - n_keywords=4, n_subsections=3, repeats=.TRUE.) + n_keywords=5, n_subsections=3, repeats=.TRUE.) NULLIFY (keyword, print_key, subsection) @@ -268,6 +286,13 @@ CONTAINS CALL section_add_keyword(section, keyword) CALL keyword_release(keyword) + CALL keyword_create(keyword, name="TEMPERATURE", & + description="Electronic temperature.", & + n_var=1, type_of_var=real_t, unit_str="K", & + default_r_val=300.0_dp/kelvin) + CALL section_add_keyword(section, keyword) + CALL keyword_release(keyword) + ! PRINT subsection CALL section_create(subsection, "PRINT", "Print properties for the given contact.", & repeats=.FALSE.) diff --git a/src/negf_control_types.F b/src/negf_control_types.F index dcf3eeaf48..16936698b3 100644 --- a/src/negf_control_types.F +++ b/src/negf_control_types.F @@ -53,6 +53,8 @@ MODULE negf_control_types LOGICAL :: refine_fermi_level !> Fermi level REAL(kind=dp) :: fermi_level + !> temperature [in a.u.] + REAL(kind=dp) :: temperature !> applied electric potential REAL(kind=dp) :: v_external END TYPE negf_control_contact_type @@ -70,6 +72,8 @@ MODULE negf_control_types !> atoms belonging to the scattering region as well as atoms belonging to !> screening regions of all the contacts INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist_S_screening + !> do not keep contact self-energy matrices + LOGICAL :: disable_cache !> convergence criteria for adaptive integration methods REAL(kind=dp) :: conv_density !> convergence criteria for iterative Lopez-Sancho algorithm @@ -78,15 +82,13 @@ MODULE negf_control_types REAL(kind=dp) :: conv_scf !> accuracy in mapping atoms between different force environments REAL(kind=dp) :: eps_geometry - !> temperature [in a.u.] - REAL(kind=dp) :: temperature !> applied bias [in a.u.] REAL(kind=dp) :: v_bias !> integration lower bound [in a.u.] REAL(kind=dp) :: energy_lbound !> infinitesimal offset along the imaginary axis [in a.u.] REAL(kind=dp) :: eta - !> initial guess in determine the actual Fermi level of bulk contacts [in a.u.] + !> initial guess to determine the actual Fermi level of bulk contacts [in a.u.] REAL(kind=dp) :: homo_lumo_gap !> number of residuals (poles of the Fermi function) INTEGER :: delta_npoles @@ -102,6 +104,12 @@ MODULE negf_control_types INTEGER :: max_scf !> minimal number of MPI processes to be used to compute Green's function per energy point INTEGER :: nprocs + !> shift in Hartree potential [in a.u.] + REAL(kind=dp) :: v_shift + !> initial offset to determine the correct shift in Hartree potential [in a.u.] + REAL(kind=dp) :: v_shift_offset + !> maximal number of iteration to determine the shift in Hartree potential + INTEGER :: v_shift_maxiters END TYPE negf_control_type PRIVATE :: read_negf_atomlist @@ -237,6 +245,13 @@ CONTAINS negf_control%contacts(i_rep)%compute_fermi_level = (.NOT. is_explicit) .OR. & negf_control%contacts(i_rep)%refine_fermi_level + CALL section_vals_val_get(subsection, "TEMPERATURE", & + r_val=negf_control%contacts(i_rep)%temperature, & + i_rep_section=i_rep) + IF (negf_control%contacts(i_rep)%temperature <= 0.0_dp) THEN + CALL cp_abort(__LOCATION__, "Electronic temperature must be > 0") + END IF + CALL section_vals_val_get(subsection, "ELECTRIC_POTENTIAL", & r_val=negf_control%contacts(i_rep)%v_external, & i_rep_section=i_rep) @@ -244,14 +259,14 @@ CONTAINS is_explicit = read_negf_atomlist(negf_control%atomlist_S, negf_section, 1, "SCATTERING_REGION", subsys) + CALL section_vals_val_get(negf_section, "DISABLE_CACHE", l_val=negf_control%disable_cache) + CALL section_vals_val_get(negf_section, "EPS_DENSITY", r_val=negf_control%conv_density) CALL section_vals_val_get(negf_section, "EPS_GREEN", r_val=negf_control%conv_green) CALL section_vals_val_get(negf_section, "EPS_SCF", r_val=negf_control%conv_scf) CALL section_vals_val_get(negf_section, "EPS_GEO", r_val=negf_control%eps_geometry) - CALL section_vals_val_get(negf_section, "TEMPERATURE", r_val=negf_control%temperature) - CALL section_vals_val_get(negf_section, "ENERGY_LBOUND", r_val=negf_control%energy_lbound) CALL section_vals_val_get(negf_section, "ETA", r_val=negf_control%eta) CALL section_vals_val_get(negf_section, "HOMO_LUMO_GAP", r_val=negf_control%homo_lumo_gap) @@ -269,38 +284,47 @@ CONTAINS CALL section_vals_val_get(negf_section, "NPROC_POINT", i_val=negf_control%nprocs) + CALL section_vals_val_get(negf_section, "V_SHIFT", r_val=negf_control%v_shift) + CALL section_vals_val_get(negf_section, "V_SHIFT_OFFSET", r_val=negf_control%v_shift_offset) + CALL section_vals_val_get(negf_section, "V_SHIFT_MAX_ITERS", i_val=negf_control%v_shift_maxiters) + ! check consistency IF (negf_control%eta < 0.0_dp) THEN CALL cp_abort(__LOCATION__, "ETA must be >= 0") END IF - IF (negf_control%temperature <= 0.0_dp) THEN - CALL cp_abort(__LOCATION__, "Electronic temperature must be > 0") + IF (n_rep > 0) THEN + delta_npoles_min = NINT(0.5_dp*(negf_control%eta/(pi*MAXVAL(negf_control%contacts(:)%temperature))+1.0_dp)) + ELSE + delta_npoles_min = 1 END IF - delta_npoles_min = NINT(0.5_dp*(negf_control%eta/(pi*negf_control%temperature)+1.0_dp)) - IF (negf_control%delta_npoles < delta_npoles_min) THEN - eta_max = REAL(2*negf_control%delta_npoles-1, kind=dp)*pi*negf_control%temperature - temp_current = negf_control%temperature*kelvin - temp_min = negf_control%eta/(pi*REAL(2*negf_control%delta_npoles-1, kind=dp))*kelvin + IF (n_rep > 0) THEN + eta_max = REAL(2*negf_control%delta_npoles-1, kind=dp)*pi*MAXVAL(negf_control%contacts(:)%temperature) + temp_current = MAXVAL(negf_control%contacts(:)%temperature)*kelvin + temp_min = negf_control%eta/(pi*REAL(2*negf_control%delta_npoles-1, kind=dp))*kelvin - WRITE (eta_current_str, '(ES11.4E2)') negf_control%eta - WRITE (eta_max_str, '(ES11.4E2)') eta_max - WRITE (npoles_current_str, '(I11)') negf_control%delta_npoles - WRITE (npoles_min_str, '(I11)') delta_npoles_min - WRITE (temp_current_str, '(F11.3)') temp_current - WRITE (temp_min_str, '(F11.3)') temp_min + WRITE (eta_current_str, '(ES11.4E2)') negf_control%eta + WRITE (eta_max_str, '(ES11.4E2)') eta_max + WRITE (npoles_current_str, '(I11)') negf_control%delta_npoles + WRITE (npoles_min_str, '(I11)') delta_npoles_min + WRITE (temp_current_str, '(F11.3)') temp_current + WRITE (temp_min_str, '(F11.3)') temp_min - CALL cp_abort(__LOCATION__, & - "Parameter DELTA_NPOLES must be at least "//TRIM(ADJUSTL(npoles_min_str))// & - " (instead of "//TRIM(ADJUSTL(npoles_current_str))// & - ") for given TEMPERATURE ("//TRIM(ADJUSTL(temp_current_str))// & - " K) and ETA ("//TRIM(ADJUSTL(eta_current_str))// & - "). Alternatively you can increase TEMPERATURE above "//TRIM(ADJUSTL(temp_min_str))// & - " K, or decrease ETA below "//TRIM(ADJUSTL(eta_max_str))// & - ". Please keep in mind that very tight ETA may result in dramatical precision loss"// & - " due to inversion of ill-conditioned matrices.") + CALL cp_abort(__LOCATION__, & + "Parameter DELTA_NPOLES must be at least "//TRIM(ADJUSTL(npoles_min_str))// & + " (instead of "//TRIM(ADJUSTL(npoles_current_str))// & + ") for given TEMPERATURE ("//TRIM(ADJUSTL(temp_current_str))// & + " K) and ETA ("//TRIM(ADJUSTL(eta_current_str))// & + "). Alternatively you can increase TEMPERATURE above "//TRIM(ADJUSTL(temp_min_str))// & + " K, or decrease ETA below "//TRIM(ADJUSTL(eta_max_str))// & + ". Please keep in mind that very tight ETA may result in dramatical precision loss"// & + " due to inversion of ill-conditioned matrices.") + ELSE + ! no leads have been defined, so calculation will abort anyway + negf_control%delta_npoles = delta_npoles_min + END IF END IF ! expand scattering region by adding atoms from contact screening regions diff --git a/src/negf_env_types.F b/src/negf_env_types.F index 50bf6fabd5..3ea33fc1ab 100644 --- a/src/negf_env_types.F +++ b/src/negf_env_types.F @@ -47,18 +47,18 @@ MODULE negf_env_types negf_copy_sym_dbcsr_to_fm_submat,& negf_reference_contact_matrix,& number_of_atomic_orbitals - USE negf_qs_env_types, ONLY: negf_qs_env_type USE negf_subgroup_types, ONLY: negf_subgroup_env_type USE negf_vectors, ONLY: contact_direction_vector,& projection_on_direction_vector USE particle_types, ONLY: particle_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type - USE pw_methods, ONLY: pw_copy USE pw_pool_types, ONLY: pw_pool_create_pw,& + pw_pool_give_back_pw,& pw_pool_type USE pw_types, ONLY: REALDATA3D,& REALSPACE,& + pw_p_type,& pw_type USE qs_density_mixing_types, ONLY: mixing_storage_create,& mixing_storage_release,& @@ -67,6 +67,7 @@ MODULE negf_env_types USE qs_energy_init, ONLY: qs_energies_init USE qs_environment_types, ONLY: get_qs_env,& qs_environment_type + USE qs_integrate_potential, ONLY: integrate_v_rspace USE qs_mo_types, ONLY: get_mo_set,& mo_set_p_type USE qs_rho_types, ONLY: qs_rho_get,& @@ -90,6 +91,7 @@ MODULE negf_env_types ! ************************************************************************************************** TYPE negf_env_contact_type REAL(kind=dp), DIMENSION(3) :: direction_vector, origin + REAL(kind=dp), DIMENSION(3) :: direction_vector_bias, origin_bias !> an axis towards the secondary contact unit cell which coincides with the transport direction !> 0 (undefined), 1 (+x), 2 (+y), 3 (+z), -1 (-x), -2 (-y), -3 (-z) INTEGER :: direction_axis @@ -125,6 +127,8 @@ MODULE negf_env_types TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:, :) :: h_sc !> overlap matrix of the scattering region TYPE(cp_fm_type), POINTER :: s_s => null() + !> an external Hartree potential in atomic basis set representation + TYPE(cp_fm_type), POINTER :: v_hartree_s => null() !> overlap matrix of the scattering region -- contact interface for every contact ([ncontacts]) TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: s_sc !> structure needed for density mixing @@ -176,9 +180,6 @@ CONTAINS TYPE(force_env_p_type), DIMENSION(:), POINTER :: sub_force_env TYPE(negf_atom_map_contact_type), ALLOCATABLE, & DIMENSION(:) :: map_contact - TYPE(negf_qs_env_type), POINTER :: negf_qs_env - TYPE(pw_env_type), POINTER :: pw_env - TYPE(pw_pool_type), POINTER :: pw_pool TYPE(pw_type), POINTER :: v_hartree_rspace TYPE(qs_environment_type), POINTER :: qs_env, qs_env_contact TYPE(qs_subsys_type), POINTER :: subsys, subsys_contact @@ -229,9 +230,7 @@ CONTAINS CALL qs_energies_init(qs_env, calc_forces=.FALSE.) CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, & matrix_s_kp=matrix_s_kp, matrix_ks_kp=matrix_ks_kp, & - negf_env=negf_qs_env, pw_env=pw_env, subsys=subsys, & - v_hartree_rspace=v_hartree_rspace) - CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool) + subsys=subsys, v_hartree_rspace=v_hartree_rspace) IF (do_kpoints) THEN CPABORT("k-points are currently not supported for device FORCE_EVAL") @@ -262,12 +261,6 @@ CONTAINS END IF END DO - ! obtain a Konh-Sham and overlap matrices for each contact - ALLOCATE (negf_qs_env%contacts(ncontacts)) - - ! The following code is kept commented until we find a way to obtain matrix elements in a consistent way - ! *********************************************************************************************************************** - ! ! stage 2: obtain relevant Kohn-Sham matrix blocks for each contact DO icontact = 1, ncontacts IF (log_unit > 0) & @@ -279,28 +272,12 @@ CONTAINS CALL negf_env_contact_init_matrices(contact_env=negf_env%contacts(icontact), sub_env=sub_env, & qs_env_contact=qs_env_contact, matrix_s_device=matrix_s_kp(1, 1)%matrix) - ! - ! CALL negf_get_bulk_contact_matrix(negf_control_aux%contacts(icontact)%ks_block_map, & - ! negf_control_aux%contacts(icontact)%ks_block_storage,& - ! matrix_ks_kp, negf_control%contacts(icontact)%atomlist_bulk, & - ! map_contact(icontact)%atom_map, & - ! matrix_ks_kp_contact, cell_to_index, blacs_env%para_env) - - negf_qs_env%contacts(icontact)%direction_vector = negf_env%contacts(icontact)%direction_vector - negf_qs_env%contacts(icontact)%origin_scattering = negf_env%contacts(icontact)%origin - negf_qs_env%contacts(icontact)%v_external = negf_control%contacts(icontact)%v_external END DO - !negf_control_aux%update_ks_matrix = .TRUE. - ! *********************************************************************************************************************** ! obtain an initial KS-matrix for the scattering region CALL qs_energies(qs_env, consistent_energies=.FALSE., calc_forces=.FALSE.) CALL negf_env_device_init_matrices(negf_env, negf_control, sub_env, qs_env) - ! reference Hartree potential - CALL pw_pool_create_pw(pw_pool, negf_qs_env%v_hartree_ref, use_data=REALDATA3D, in_space=REALSPACE) - CALL pw_copy(v_hartree_rspace, negf_qs_env%v_hartree_ref) - ! electron density mixing; ! the input section below should be consistent with the subroutine create_negf_section() NULLIFY (negf_env%mixing_storage) @@ -341,6 +318,8 @@ CONTAINS CALL contact_direction_vector(contact_env%origin, & contact_env%direction_vector, & + contact_env%origin_bias, & + contact_env%direction_vector_bias, & contact_control%atomlist_screening, & contact_control%atomlist_bulk, & subsys_device) @@ -563,8 +542,12 @@ CONTAINS LOGICAL :: do_kpoints TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type) :: hmat TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp, matrix_s_kp TYPE(dft_control_type), POINTER :: dft_control + TYPE(pw_env_type), POINTER :: pw_env + TYPE(pw_p_type) :: v_hartree + TYPE(pw_pool_type), POINTER :: pw_pool TYPE(qs_subsys_type), POINTER :: subsys CALL timeset(routineN, handle) @@ -576,7 +559,10 @@ CONTAINS matrix_ks_kp=matrix_ks_kp, & matrix_s_kp=matrix_s_kp, & para_env=para_env, & + pw_env=pw_env, & subsys=subsys) + CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool) + IF (do_kpoints) THEN CALL cp_abort(__LOCATION__, & "K-points in device region have not been implemented yet.") @@ -589,7 +575,7 @@ CONTAINS nao_s = number_of_atomic_orbitals(subsys, negf_control%atomlist_S_screening) ! ++ create matrices: h_s, s_s - NULLIFY (negf_env%s_s, fm_struct) + NULLIFY (negf_env%s_s, negf_env%v_hartree_s, fm_struct) ALLOCATE (negf_env%h_s(nspins)) CALL cp_fm_struct_create(fm_struct, nrow_global=nao_s, ncol_global=nao_s, context=sub_env%blacs_env) @@ -598,6 +584,7 @@ CONTAINS NULLIFY (negf_env%h_s(ispin)%matrix) CALL cp_fm_create(negf_env%h_s(ispin)%matrix, fm_struct) END DO + CALL cp_fm_create(negf_env%v_hartree_s, fm_struct) CALL cp_fm_struct_release(fm_struct) ! ++ create matrices: h_sc, s_sc @@ -634,6 +621,30 @@ CONTAINS subsys=subsys, mpi_comm_global=para_env%group, & do_upper_diag=.TRUE., do_lower=.TRUE.) + ! v_hartree_s + NULLIFY (hmat%matrix) + CALL dbcsr_init_p(hmat%matrix) + CALL dbcsr_copy(matrix_b=hmat%matrix, matrix_a=matrix_s_kp(1, 1)%matrix) + CALL dbcsr_set(hmat%matrix, 0.0_dp) + + NULLIFY (v_hartree%pw) + CALL pw_pool_create_pw(pw_pool, v_hartree%pw, use_data=REALDATA3D, in_space=REALSPACE) + CALL negf_env_init_v_hartree(v_hartree%pw, negf_env%contacts, negf_control%contacts) + + CALL integrate_v_rspace(v_rspace=v_hartree, hmat=hmat, qs_env=qs_env, & + calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE.) + + CALL pw_pool_give_back_pw(pw_pool, v_hartree%pw) + + CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=hmat%matrix, & + fm=negf_env%v_hartree_s, & + atomlist_row=negf_control%atomlist_S_screening, & + atomlist_col=negf_control%atomlist_S_screening, & + subsys=subsys, mpi_comm_global=para_env%group, & + do_upper_diag=.TRUE., do_lower=.TRUE.) + + CALL dbcsr_deallocate_matrix(hmat%matrix) + ! extract matrices: h_sc, s_sc DO icontact = 1, ncontacts DO ispin = 1, nspins @@ -657,6 +668,95 @@ CONTAINS CALL timestop(handle) END SUBROUTINE negf_env_device_init_matrices +! ************************************************************************************************** +!> \brief Contribution to the Hartree potential related to the external bias voltage. +!> \param v_hartree Hartree potential (modified on exit) +!> \param contact_env NEGF environment for every contact +!> \param contact_control NEGF control for every contact +!> \author Sergey Chulkov +! ************************************************************************************************** + SUBROUTINE negf_env_init_v_hartree(v_hartree, contact_env, contact_control) + TYPE(pw_type), POINTER :: v_hartree + TYPE(negf_env_contact_type), DIMENSION(:), & + INTENT(in) :: contact_env + TYPE(negf_control_contact_type), DIMENSION(:), & + INTENT(in) :: contact_control + + CHARACTER(len=*), PARAMETER :: routineN = 'negf_env_init_v_hartree', & + routineP = moduleN//':'//routineN + REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp) + + INTEGER :: dx, dy, dz, handle, icontact, ix, iy, & + iz, lx, ly, lz, ncontacts, ux, uy, uz + REAL(kind=dp) :: dvol, pot, proj, v1, v2 + REAL(kind=dp), DIMENSION(3) :: dirvector_bias, point_coord, & + point_indices, vector + + CALL timeset(routineN, handle) + + ncontacts = SIZE(contact_env) + CPASSERT(SIZE(contact_control) == ncontacts) + CPASSERT(ncontacts == 2) + + dirvector_bias = contact_env(2)%origin_bias-contact_env(1)%origin_bias + v1 = contact_control(1)%v_external + v2 = contact_control(2)%v_external + + lx = v_hartree%pw_grid%bounds_local(1, 1) + ux = v_hartree%pw_grid%bounds_local(2, 1) + ly = v_hartree%pw_grid%bounds_local(1, 2) + uy = v_hartree%pw_grid%bounds_local(2, 2) + lz = v_hartree%pw_grid%bounds_local(1, 3) + uz = v_hartree%pw_grid%bounds_local(2, 3) + + dx = v_hartree%pw_grid%npts(1)/2 + dy = v_hartree%pw_grid%npts(2)/2 + dz = v_hartree%pw_grid%npts(3)/2 + + dvol = v_hartree%pw_grid%dvol + + DO iz = lz, uz + point_indices(3) = REAL(iz+dz, kind=dp) + DO iy = ly, uy + point_indices(2) = REAL(iy+dy, kind=dp) + + DO ix = lx, ux + point_indices(1) = REAL(ix+dx, kind=dp) + point_coord(:) = MATMUL(v_hartree%pw_grid%dh, point_indices) + + vector = point_coord-contact_env(1)%origin_bias + proj = projection_on_direction_vector(vector, dirvector_bias) + IF (proj+threshold >= 0.0_dp .AND. proj-threshold <= 1.0_dp) THEN + ! scattering region + ! proj == 0 we are at the first contact boundary + ! proj == 1 we are at the second contact boundary + IF (proj < 0.0_dp) THEN + proj = 0.0_dp + ELSE IF (proj > 1.0_dp) THEN + proj = 1.0_dp + END IF + pot = v1+(v2-v1)*proj + ELSE + pot = 0.0_dp + DO icontact = 1, ncontacts + vector = point_coord-contact_env(icontact)%origin_bias + proj = projection_on_direction_vector(vector, contact_env(icontact)%direction_vector_bias) + + IF (proj+threshold >= 0.0_dp .AND. proj-threshold <= 1.0_dp) THEN + pot = contact_control(icontact)%v_external + EXIT + END IF + END DO + END IF + + v_hartree%cr3d(ix, iy, iz) = pot*dvol + END DO + END DO + END DO + + CALL timestop(handle) + END SUBROUTINE negf_env_init_v_hartree + ! ************************************************************************************************** !> \brief Detect the axis towards secondary unit cell. !> \param direction_vector direction vector @@ -1000,6 +1100,10 @@ CONTAINS DEALLOCATE (negf_env%s_sc) END IF + ! v_hartree_s + IF (ASSOCIATED(negf_env%v_hartree_s)) & + CALL cp_fm_release(negf_env%v_hartree_s) + ! density mixing CALL mixing_storage_release(negf_env%mixing_storage) diff --git a/src/negf_green_methods.F b/src/negf_green_methods.F index 8795ce3669..b1dcc64b16 100644 --- a/src/negf_green_methods.F +++ b/src/negf_green_methods.F @@ -310,13 +310,14 @@ CONTAINS !> \param self_energy_ret_sum sum of contact self-energies !> \param h_s Kohn-Sham matrix block of the scattering region !> \param s_s overlap matrix block of the scattering region +!> \param v_hartree_s contribution to the Kohn-Sham matrix from the external potential !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE negf_retarded_green_function(g_ret_s, omega, self_energy_ret_sum, h_s, s_s) + SUBROUTINE negf_retarded_green_function(g_ret_s, omega, self_energy_ret_sum, h_s, s_s, v_hartree_s) TYPE(cp_cfm_type), POINTER :: g_ret_s COMPLEX(kind=dp), INTENT(in) :: omega TYPE(cp_cfm_type), POINTER :: self_energy_ret_sum - TYPE(cp_fm_type), POINTER :: h_s, s_s + TYPE(cp_fm_type), POINTER :: h_s, s_s, v_hartree_s CHARACTER(len=*), PARAMETER :: routineN = 'negf_retarded_green_function', & routineP = moduleN//':'//routineN @@ -327,9 +328,11 @@ CONTAINS ! g_ret_s = [ omega * S_S - H_S - self_energy_left - self_energy_right]^{-1} ! - ! omega * S_S - H_S + ! omega * S_S - H_S - V_Hartree CALL cp_fm_to_cfm(msourcer=s_s, mtarget=g_ret_s) CALL cp_cfm_scale_and_add_fm(omega, g_ret_s, z_mone, h_s) + IF (ASSOCIATED(v_hartree_s)) & + CALL cp_cfm_scale_and_add_fm(z_one, g_ret_s, z_one, v_hartree_s) ! g_ret_s = [omega * S_S - H_S - \sum_{contact} self_energy_{contact}^{ret.} ]^-1 CALL cp_cfm_scale_and_add(z_one, g_ret_s, z_mone, self_energy_ret_sum) diff --git a/src/negf_methods.F b/src/negf_methods.F index 857983bbe5..c2a0cb1091 100644 --- a/src/negf_methods.F +++ b/src/negf_methods.F @@ -20,7 +20,7 @@ MODULE negf_methods copy_cfm_info_type, cp_cfm_cleanup_copy_general, cp_cfm_create, & cp_cfm_finish_copy_general, cp_cfm_get_info, cp_cfm_get_submatrix, cp_cfm_p_type, & cp_cfm_release, cp_cfm_retain, cp_cfm_set_submatrix, cp_cfm_start_copy_general, & - cp_cfm_to_fm + cp_cfm_to_fm, cp_cfm_type USE cp_control_types, ONLY: dft_control_type USE cp_fm_basic_linalg, ONLY: cp_fm_scale,& cp_fm_scale_and_add,& @@ -97,7 +97,6 @@ MODULE negf_methods USE negf_matrix_utils, ONLY: invert_cell_to_index,& negf_copy_fm_submat_to_dbcsr,& negf_copy_sym_dbcsr_to_fm_submat - USE negf_qs_env_types, ONLY: negf_qs_env_type USE negf_subgroup_types, ONLY: negf_sub_env_create,& negf_sub_env_release,& negf_subgroup_env_type @@ -199,7 +198,7 @@ CONTAINS ncontacts = SIZE(negf_control%contacts) DO icontact = 1, ncontacts NULLIFY (qs_env) - CALL force_env_get(sub_force_env(icontact)%force_env, qs_env=qs_env) + CALL force_env_get(sub_force_env(negf_control%contacts(icontact)%force_env_index)%force_env, qs_env=qs_env) CALL guess_fermi_level(icontact, negf_env, negf_control, sub_env, qs_env, log_unit) print_section => section_vals_get_subs_vals(negf_contact_section, "PRINT", i_rep_section=icontact) @@ -214,8 +213,9 @@ CONTAINS print_unit = cp_print_key_unit_nr(logger, print_section, "DOS", & middle_name=TRIM(ADJUSTL(contact_id_str)), extension=".dos") - CALL negf_print_dos(print_unit, energy_min, energy_max, npoints, negf_env, & - negf_control, sub_env, base_contact=icontact, just_contact=icontact) + CALL negf_print_dos(print_unit, energy_min, energy_max, npoints, & + v_shift=0.0_dp, negf_env=negf_env, negf_control=negf_control, & + sub_env=sub_env, base_contact=icontact, just_contact=icontact) CALL cp_print_key_finished_output(print_unit, logger, print_section, "DOS") END IF @@ -224,7 +224,9 @@ CONTAINS IF (ncontacts > 1) THEN NULLIFY (qs_env) CALL force_env_get(force_env, qs_env=qs_env) - CALL converge_density(negf_env, negf_control, sub_env, qs_env, base_contact=1, log_unit=log_unit) + CALL shift_potential(negf_env, negf_control, sub_env, qs_env, base_contact=1, log_unit=log_unit) + + CALL converge_density(negf_env, negf_control, sub_env, qs_env, negf_control%v_shift, base_contact=1, log_unit=log_unit) ! current CALL get_qs_env(qs_env, dft_control=dft_control) @@ -233,9 +235,13 @@ CONTAINS current = 0.0_dp DO ispin = 1, nspins - current = current+negf_compute_current(1, 2, negf_env, negf_control, sub_env, ispin, blacs_env) + current = current+ & + negf_compute_current(1, 2, v_shift=negf_control%v_shift, negf_env=negf_env, negf_control=negf_control, & + sub_env=sub_env, ispin=ispin, blacs_env_global=blacs_env) END DO + IF (nspins == 1) current = 2.0_dp*current + IF (log_unit > 0) THEN WRITE (log_unit, '(/,T2,A,T60,ES20.7E2)') "NEGF| Electric current (A)", current END IF @@ -253,8 +259,9 @@ CONTAINS print_unit = cp_print_key_unit_nr(logger, print_section, "DOS", & middle_name=TRIM(ADJUSTL(contact_id_str)), extension=".dos") - CALL negf_print_dos(print_unit, energy_min, energy_max, npoints, negf_env, & - negf_control, sub_env, base_contact=1) + CALL negf_print_dos(print_unit, energy_min, energy_max, npoints, negf_control%v_shift, & + negf_env=negf_env, negf_control=negf_control, & + sub_env=sub_env, base_contact=1) CALL cp_print_key_finished_output(print_unit, logger, print_section, "DOS") END IF @@ -271,8 +278,9 @@ CONTAINS print_unit = cp_print_key_unit_nr(logger, print_section, "TRANSMISSION", & middle_name=TRIM(ADJUSTL(contact_id_str)), extension=".transm") - CALL negf_print_transmission(print_unit, energy_min, energy_max, npoints, negf_env, & - negf_control, sub_env, contact_id1=1, contact_id2=2) + CALL negf_print_transmission(print_unit, energy_min, energy_max, npoints, negf_control%v_shift, & + negf_env=negf_env, negf_control=negf_control, & + sub_env=sub_env, contact_id1=1, contact_id2=2) CALL cp_print_key_finished_output(print_unit, logger, print_section, "TRANSMISSION") END IF @@ -410,7 +418,7 @@ CONTAINS DEALLOCATE (index_to_cell) IF (log_unit > 0) THEN - WRITE (temperature_str, '(F11.3)') negf_control%temperature*kelvin + WRITE (temperature_str, '(F11.3)') negf_control%contacts(contact_id)%temperature*kelvin WRITE (log_unit, '(/,T2,A,I0,A)') "COMPUTE FERMI LEVEL OF CONTACT ", & contact_id, " AT "//TRIM(ADJUSTL(temperature_str))//" KELVIN" WRITE (log_unit, '(/,T2,A,T60,F20.10,/)') "Electronic density of the isolated contact unit cell:", & @@ -427,8 +435,7 @@ CONTAINS END DO ! Use orbital energies of HOMO and LUMO as reference points and then - ! refine Fermi level by using a simple linear interpolation technique - + ! refine the Fermi level by using a simple linear interpolation technique IF (negf_control%homo_lumo_gap > 0.0_dp) THEN IF (negf_control%contacts(contact_id)%refine_fermi_level) THEN fermi_level_min = negf_control%contacts(contact_id)%fermi_level @@ -447,9 +454,9 @@ CONTAINS step = 0 lbound_cpath = CMPLX(negf_control%energy_lbound, negf_control%eta, kind=dp) - delta_au = REAL(negf_control%delta_npoles, kind=dp)*twopi*negf_control%temperature - offset_au = REAL(negf_control%gamma_kT, kind=dp)*negf_control%temperature - energy_ubound_minus_fermi = -2.0_dp*LOG(negf_control%conv_density)*negf_control%temperature + delta_au = REAL(negf_control%delta_npoles, kind=dp)*twopi*negf_control%contacts(contact_id)%temperature + offset_au = REAL(negf_control%gamma_kT, kind=dp)*negf_control%contacts(contact_id)%temperature + energy_ubound_minus_fermi = -2.0_dp*LOG(negf_control%conv_density)*negf_control%contacts(contact_id)%temperature t1 = m_walltime() DO @@ -475,6 +482,8 @@ CONTAINS DO ispin = 1, nspins CALL negf_init_rho_equiv_residuals(rho_ao_fm=rho_ao_fm, & + v_shift=0.0_dp, & + ignore_bias=.TRUE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -484,6 +493,8 @@ CONTAINS CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm, & stats=stats, & + v_shift=0.0_dp, & + ignore_bias=.TRUE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -499,6 +510,8 @@ CONTAINS CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm, & stats=stats, & + v_shift=0.0_dp, & + ignore_bias=.TRUE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -562,7 +575,7 @@ CONTAINS END IF IF (log_unit > 0) THEN - WRITE (temperature_str, '(F11.3)') negf_control%temperature*kelvin + WRITE (temperature_str, '(F11.3)') negf_control%contacts(contact_id)%temperature*kelvin WRITE (log_unit, '(/,T2,A,I0)') "NEGF| Contact No. ", contact_id WRITE (log_unit, '(T2,A,T62,F18.8)') "NEGF| Fermi level at "//TRIM(ADJUSTL(temperature_str))// & " Kelvin (a.u.):", negf_control%contacts(contact_id)%fermi_level @@ -574,49 +587,40 @@ CONTAINS END SUBROUTINE guess_fermi_level ! ************************************************************************************************** -!> \brief Converge electronic density of the scattering region. +!> \brief Compute shift in Hartree potential !> \param negf_env NEGF environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment !> \param qs_env QuickStep environment !> \param base_contact index of the reference contact !> \param log_unit output unit -!> \par History -!> * 06.2017 created [Sergey Chulkov] ! ************************************************************************************************** - SUBROUTINE converge_density(negf_env, negf_control, sub_env, qs_env, base_contact, log_unit) + SUBROUTINE shift_potential(negf_env, negf_control, sub_env, qs_env, base_contact, log_unit) TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env TYPE(qs_environment_type), POINTER :: qs_env INTEGER, INTENT(in) :: base_contact, log_unit - CHARACTER(LEN=*), PARAMETER :: routineN = 'converge_density', & + CHARACTER(LEN=*), PARAMETER :: routineN = 'shift_potential', & routineP = moduleN//':'//routineN - REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp) COMPLEX(kind=dp) :: lbound_cpath, ubound_cpath, ubound_lpath - INTEGER :: handle, icontact, image, ispin, & - iter_count, jcontact, nao, ncontacts, & - nimages, nspins + INTEGER :: handle, ispin, iter_count, nao, & + ncontacts, nspins LOGICAL :: do_kpoints - REAL(kind=dp) :: iter_delta, mu_base, nelectrons, & - nelectrons_diff, t1, t2, trace, & - v_base, v_contact + REAL(kind=dp) :: mu_base, nelectrons_guess, nelectrons_max, nelectrons_min, nelectrons_ref, & + t1, t2, temperature, trace, v_shift_guess, v_shift_max, v_shift_min TYPE(cp_blacs_env_type), POINTER :: blacs_env - TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: rho_ao_delta_fm, rho_ao_new_fm + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: rho_ao_fm TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_fm_type), POINTER :: matrix_s_fm, work_fm TYPE(cp_para_env_type), POINTER :: para_env - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_initial_kp, matrix_ks_qs_kp, & - rho_ao_initial_kp, rho_ao_new_kp, & - rho_ao_qs_kp + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_qs_kp TYPE(dft_control_type), POINTER :: dft_control TYPE(green_functions_cache_type), ALLOCATABLE, & - DIMENSION(:) :: g_surf_circular, g_surf_linear, & - g_surf_nonequiv + DIMENSION(:) :: g_surf_circular, g_surf_linear TYPE(integration_status_type) :: stats - TYPE(negf_qs_env_type), POINTER :: negf_qs_env TYPE(qs_rho_type), POINTER :: rho_struct TYPE(qs_subsys_type), POINTER :: subsys @@ -628,14 +632,282 @@ CONTAINS CALL timeset(routineN, handle) - CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, dft_control=dft_control, negf_env=negf_qs_env, & + CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, dft_control=dft_control, & + para_env=para_env, rho=rho_struct, subsys=subsys) + CPASSERT(.NOT. do_kpoints) + + ! apply external NEGF potential + t1 = m_walltime() + + ! need a globally distributed overlap matrix in order to compute integration errors + IF (sub_env%ngroups > 1) THEN + NULLIFY (matrix_s_fm, fm_struct) + + CALL cp_fm_get_info(negf_env%s_s, nrow_global=nao) + CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nao, context=blacs_env) + + CALL cp_fm_create(matrix_s_fm, fm_struct) + CALL cp_fm_struct_release(fm_struct) + + IF (sub_env%group_distribution(sub_env%mepos_global) == 0) THEN + work_fm => negf_env%s_s + ELSE + NULLIFY (work_fm) + END IF + + CALL cp_fm_copy_general(work_fm, matrix_s_fm, para_env) + ELSE + matrix_s_fm => negf_env%s_s + CALL cp_fm_retain(matrix_s_fm) + CALL cp_fm_get_info(matrix_s_fm, matrix_struct=fm_struct) + END IF + + CALL cp_fm_get_info(matrix_s_fm, matrix_struct=fm_struct) + + nspins = SIZE(negf_env%h_s) + + mu_base = negf_control%contacts(base_contact)%fermi_level + + ! keep the initial charge density matrix and Kohn-Sham matrix + CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_qs_kp) + + ! extract the reference density matrix blocks + nelectrons_ref = 0.0_dp + ALLOCATE (rho_ao_fm(nspins)) + DO ispin = 1, nspins + NULLIFY (rho_ao_fm(ispin)%matrix) + CALL cp_fm_create(rho_ao_fm(ispin)%matrix, fm_struct) + + CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=rho_ao_qs_kp(ispin, 1)%matrix, & + fm=rho_ao_fm(ispin)%matrix, & + atomlist_row=negf_control%atomlist_S_screening, & + atomlist_col=negf_control%atomlist_S_screening, & + subsys=subsys, mpi_comm_global=para_env%group, & + do_upper_diag=.TRUE., do_lower=.TRUE.) + + CALL cp_fm_trace(rho_ao_fm(ispin)%matrix, matrix_s_fm, trace) + nelectrons_ref = nelectrons_ref+trace + END DO + + IF (log_unit > 0) THEN + WRITE (log_unit, '(/,T2,A)') "COMPUTE SHIFT IN HARTREE POTENTIAL" + WRITE (log_unit, '(/,T2,A,T55,F25.14,/)') "Initial electronic density of the scattering region:", -1.0_dp*nelectrons_ref + WRITE (log_unit, '(T3,A)') "Step Integration method Time V shift Convergence (density)" + WRITE (log_unit, '(T3,78("-"))') + END IF + + temperature = negf_control%contacts(base_contact)%temperature + + ! integration limits: C-path (arch) + lbound_cpath = CMPLX(negf_control%energy_lbound, negf_control%eta, kind=dp) + ubound_cpath = CMPLX(mu_base-REAL(negf_control%gamma_kT, kind=dp)*temperature, & + REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp) + + ! integration limits: L-path (linear) + ubound_lpath = CMPLX(mu_base-LOG(negf_control%conv_density)*temperature, & + REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp) + + v_shift_min = negf_control%v_shift + v_shift_max = negf_control%v_shift+negf_control%v_shift_offset + + ALLOCATE (g_surf_circular(nspins), g_surf_linear(nspins)) + + DO iter_count = 1, negf_control%v_shift_maxiters + SELECT CASE (iter_count) + CASE (1) + v_shift_guess = v_shift_min + CASE (2) + v_shift_guess = v_shift_max + CASE DEFAULT + v_shift_guess = v_shift_min-(nelectrons_min-nelectrons_ref)* & + (v_shift_max-v_shift_min)/(nelectrons_max-nelectrons_min) + END SELECT + + ! compute an updated density matrix + CALL integration_status_reset(stats) + + DO ispin = 1, nspins + ! closed contour: residuals + CALL negf_init_rho_equiv_residuals(rho_ao_fm=rho_ao_fm(ispin)%matrix, & + v_shift=v_shift_guess, & + ignore_bias=.TRUE., & + negf_env=negf_env, & + negf_control=negf_control, & + sub_env=sub_env, & + ispin=ispin, & + base_contact=base_contact) + + ! closed contour: C-path + CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm(ispin)%matrix, & + stats=stats, & + v_shift=v_shift_guess, & + ignore_bias=.TRUE., & + negf_env=negf_env, & + negf_control=negf_control, & + sub_env=sub_env, & + ispin=ispin, & + base_contact=base_contact, & + integr_lbound=lbound_cpath, & + integr_ubound=ubound_cpath, & + matrix_s_global=matrix_s_fm, & + is_circular=.TRUE., & + g_surf_cache=g_surf_circular(ispin)) + IF (negf_control%disable_cache) & + CALL green_functions_cache_release(g_surf_circular(ispin)) + + ! closed contour: L-path + CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_fm(ispin)%matrix, & + stats=stats, & + v_shift=v_shift_guess, & + ignore_bias=.TRUE., & + negf_env=negf_env, & + negf_control=negf_control, & + sub_env=sub_env, & + ispin=ispin, & + base_contact=base_contact, & + integr_lbound=ubound_cpath, & + integr_ubound=ubound_lpath, & + matrix_s_global=matrix_s_fm, & + is_circular=.FALSE., & + g_surf_cache=g_surf_linear(ispin)) + IF (negf_control%disable_cache) & + CALL green_functions_cache_release(g_surf_linear(ispin)) + END DO + + IF (nspins > 1) THEN + DO ispin = 1, nspins + CALL cp_fm_scale_and_add(1.0_dp, rho_ao_fm(1)%matrix, 1.0_dp, rho_ao_fm(ispin)%matrix) + END DO + ELSE + CALL cp_fm_scale(2.0_dp, rho_ao_fm(1)%matrix) + END IF + + CALL cp_fm_trace(rho_ao_fm(1)%matrix, matrix_s_fm, nelectrons_guess) + + t2 = m_walltime() + + IF (log_unit > 0) THEN + WRITE (log_unit, '(T2,I5,T12,A,T32,F8.1,T42,F15.8,T60,ES20.5E2)') & + iter_count, get_method_description_string(stats, negf_control%integr_method), & + t2-t1, v_shift_guess, nelectrons_guess-nelectrons_ref + END IF + + IF (ABS(nelectrons_guess-nelectrons_ref) < negf_control%conv_scf) EXIT + + ! compute correction + SELECT CASE (iter_count) + CASE (1) + nelectrons_min = nelectrons_guess + CASE (2) + nelectrons_max = nelectrons_guess + CASE DEFAULT + IF (v_shift_guess < v_shift_min) THEN + v_shift_max = v_shift_min + nelectrons_max = nelectrons_min + v_shift_min = v_shift_guess + nelectrons_min = nelectrons_guess + ELSE IF (v_shift_guess > v_shift_max) THEN + v_shift_min = v_shift_max + nelectrons_min = nelectrons_max + v_shift_max = v_shift_guess + nelectrons_max = nelectrons_guess + ELSE IF (v_shift_max-v_shift_guess < v_shift_guess-v_shift_min) THEN + v_shift_max = v_shift_guess + nelectrons_max = nelectrons_guess + ELSE + v_shift_min = v_shift_guess + nelectrons_min = nelectrons_guess + END IF + END SELECT + + t1 = t2 + END DO + + negf_control%v_shift = v_shift_guess + + IF (log_unit > 0) THEN + WRITE (log_unit, '(T2,A,T62,F18.8)') "NEGF| Shift in Hartree potential", negf_control%v_shift + END IF + + DO ispin = nspins, 1, -1 + CALL green_functions_cache_release(g_surf_circular(ispin)) + CALL green_functions_cache_release(g_surf_linear(ispin)) + END DO + DEALLOCATE (g_surf_circular, g_surf_linear) + + DO ispin = nspins, 1, -1 + CALL cp_fm_release(rho_ao_fm(ispin)%matrix) + END DO + DEALLOCATE (rho_ao_fm) + + IF (ASSOCIATED(matrix_s_fm)) & + CALL cp_fm_release(matrix_s_fm) + + CALL timestop(handle) + END SUBROUTINE shift_potential + +! ************************************************************************************************** +!> \brief Converge electronic density of the scattering region. +!> \param negf_env NEGF environment +!> \param negf_control NEGF control +!> \param sub_env NEGF parallel (sub)group environment +!> \param qs_env QuickStep environment +!> \param v_shift shift in Hartree potential +!> \param base_contact index of the reference contact +!> \param log_unit output unit +!> \par History +!> * 06.2017 created [Sergey Chulkov] +! ************************************************************************************************** + SUBROUTINE converge_density(negf_env, negf_control, sub_env, qs_env, v_shift, base_contact, log_unit) + TYPE(negf_env_type), INTENT(in) :: negf_env + TYPE(negf_control_type), POINTER :: negf_control + TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env + TYPE(qs_environment_type), POINTER :: qs_env + REAL(kind=dp), INTENT(in) :: v_shift + INTEGER, INTENT(in) :: base_contact, log_unit + + CHARACTER(LEN=*), PARAMETER :: routineN = 'converge_density', & + routineP = moduleN//':'//routineN + REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp) + + COMPLEX(kind=dp) :: lbound_cpath, ubound_cpath, ubound_lpath + INTEGER :: handle, icontact, image, ispin, & + iter_count, nao, ncontacts, nimages, & + nspins + LOGICAL :: do_kpoints + REAL(kind=dp) :: iter_delta, mu_base, nelectrons, & + nelectrons_diff, t1, t2, temperature, & + trace, v_base, v_contact + TYPE(cp_blacs_env_type), POINTER :: blacs_env + TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: rho_ao_delta_fm, rho_ao_new_fm + TYPE(cp_fm_struct_type), POINTER :: fm_struct + TYPE(cp_fm_type), POINTER :: ao_ao_fm_global, matrix_s_fm, work_fm + TYPE(cp_para_env_type), POINTER :: para_env + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_initial_kp, matrix_ks_qs_kp, & + rho_ao_initial_kp, rho_ao_new_kp, & + rho_ao_qs_kp + TYPE(dft_control_type), POINTER :: dft_control + TYPE(green_functions_cache_type), ALLOCATABLE, & + DIMENSION(:) :: g_surf_circular, g_surf_linear, & + g_surf_nonequiv + TYPE(integration_status_type) :: stats + TYPE(qs_rho_type), POINTER :: rho_struct + TYPE(qs_subsys_type), POINTER :: subsys + + ncontacts = SIZE(negf_control%contacts) + ! nothing to do + IF (.NOT. (ALLOCATED(negf_env%h_s) .AND. ALLOCATED(negf_env%h_sc) .AND. & + ASSOCIATED(negf_env%s_s) .AND. ALLOCATED(negf_env%s_sc))) RETURN + IF (ncontacts < 2) RETURN + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, blacs_env=blacs_env, do_kpoints=do_kpoints, dft_control=dft_control, & matrix_ks_kp=matrix_ks_qs_kp, para_env=para_env, rho=rho_struct, subsys=subsys) CPASSERT(.NOT. do_kpoints) ! apply external NEGF potential t1 = m_walltime() - negf_qs_env%apply_bias = .TRUE. - CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.) ! need a globally distributed overlap matrix in order to compute integration errors IF (sub_env%ngroups > 1) THEN @@ -708,17 +980,21 @@ CONTAINS atomlist_col=negf_control%atomlist_S_screening, & subsys=subsys, mpi_comm_global=para_env%group, & do_upper_diag=.TRUE., do_lower=.TRUE.) + CALL cp_fm_trace(rho_ao_delta_fm(ispin)%matrix, matrix_s_fm, trace) nelectrons = nelectrons+trace END DO + NULLIFY (ao_ao_fm_global) + IF (sub_env%ngroups > 1) & + CALL cp_fm_create(ao_ao_fm_global, fm_struct) + ! mixing storage allocation IF (negf_env%mixing_method >= gspace_mixing_nr) THEN CALL mixing_allocate(qs_env, negf_env%mixing_method, nspins=nspins, mixing_store=negf_env%mixing_storage) CALL mixing_init(negf_env%mixing_method, rho_struct, negf_env%mixing_storage, para_env) END IF - ! TO DO: charge density form all device-contact off-diagonal blocks needs to be added as well IF (log_unit > 0) THEN WRITE (log_unit, '(/,T2,A)') "NEGF SELF-CONSISTENT PROCEDURE" WRITE (log_unit, '(/,T2,A,T55,F25.14,/)') "Initial electronic density of the scattering region:", -1.0_dp*nelectrons @@ -726,18 +1002,20 @@ CONTAINS WRITE (log_unit, '(T3,78("-"))') END IF + temperature = negf_control%contacts(base_contact)%temperature + DO icontact = 1, ncontacts IF (icontact /= base_contact) THEN v_contact = negf_control%contacts(icontact)%v_external ! integration limits: C-path (arch) lbound_cpath = CMPLX(negf_control%energy_lbound, negf_control%eta, kind=dp) - ubound_cpath = CMPLX(mu_base-REAL(negf_control%gamma_kT, kind=dp)*negf_control%temperature, & - REAL(negf_control%delta_npoles, kind=dp)*twopi*negf_control%temperature, kind=dp) + ubound_cpath = CMPLX(mu_base-REAL(negf_control%gamma_kT, kind=dp)*temperature, & + REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp) ! integration limits: L-path (linear) - ubound_lpath = CMPLX(mu_base-LOG(negf_control%conv_density)*negf_control%temperature, & - REAL(negf_control%delta_npoles, kind=dp)*twopi*negf_control%temperature, kind=dp) + ubound_lpath = CMPLX(mu_base-LOG(negf_control%conv_density)*temperature, & + REAL(negf_control%delta_npoles, kind=dp)*twopi*temperature, kind=dp) ALLOCATE (g_surf_circular(nspins), g_surf_linear(nspins), g_surf_nonequiv(nspins)) @@ -748,6 +1026,8 @@ CONTAINS DO ispin = 1, nspins ! closed contour: residuals CALL negf_init_rho_equiv_residuals(rho_ao_fm=rho_ao_new_fm(ispin)%matrix, & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -757,6 +1037,8 @@ CONTAINS ! closed contour: C-path CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_new_fm(ispin)%matrix, & stats=stats, & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -767,10 +1049,14 @@ CONTAINS matrix_s_global=matrix_s_fm, & is_circular=.TRUE., & g_surf_cache=g_surf_circular(ispin)) + IF (negf_control%disable_cache) & + CALL green_functions_cache_release(g_surf_circular(ispin)) ! closed contour: L-path CALL negf_add_rho_equiv_low(rho_ao_fm=rho_ao_new_fm(ispin)%matrix, & stats=stats, & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -781,12 +1067,15 @@ CONTAINS matrix_s_global=matrix_s_fm, & is_circular=.FALSE., & g_surf_cache=g_surf_linear(ispin)) + IF (negf_control%disable_cache) & + CALL green_functions_cache_release(g_surf_linear(ispin)) ! non-equilibrium part IF (ABS(negf_control%contacts(icontact)%v_external- & negf_control%contacts(base_contact)%v_external) >= threshold) THEN CALL negf_add_rho_nonequiv(rho_ao_fm=rho_ao_new_fm(ispin)%matrix, & stats=stats, & + v_shift=v_shift, & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -794,6 +1083,8 @@ CONTAINS base_contact=base_contact, & matrix_s_global=matrix_s_fm, & g_surf_cache=g_surf_nonequiv(ispin)) + IF (negf_control%disable_cache) & + CALL green_functions_cache_release(g_surf_nonequiv(ispin)) END IF END DO @@ -888,17 +1179,7 @@ CONTAINS atomlist_col=negf_control%atomlist_S_screening, & subsys=subsys, mpi_comm_global=para_env%group, & do_upper_diag=.TRUE., do_lower=.TRUE.) - END DO - DO jcontact = 1, ncontacts - DO ispin = 1, nspins - CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_qs_kp(ispin, 1)%matrix, & - fm=negf_env%h_sc(ispin, jcontact)%matrix, & - atomlist_row=negf_control%atomlist_S_screening, & - atomlist_col=negf_env%contacts(jcontact)%atomlist_cell0, & - subsys=subsys, mpi_comm_global=para_env%group, & - do_upper_diag=.TRUE., do_lower=.TRUE.) - END DO END DO END DO @@ -919,11 +1200,13 @@ CONTAINS END IF END DO + CALL cp_fm_release(ao_ao_fm_global) + DO ispin = nspins, 1, -1 CALL cp_fm_release(rho_ao_new_fm(ispin)%matrix) CALL cp_fm_release(rho_ao_delta_fm(ispin)%matrix) END DO - DEALLOCATE (rho_ao_new_fm, rho_ao_delta_fm) + DEALLOCATE (rho_ao_delta_fm, rho_ao_new_fm) DO image = 1, nimages DO ispin = 1, nspins @@ -937,9 +1220,9 @@ CONTAINS END DO DEALLOCATE (matrix_ks_initial_kp, rho_ao_new_kp, rho_ao_initial_kp) - CALL cp_fm_release(matrix_s_fm) + IF (ASSOCIATED(matrix_s_fm)) & + CALL cp_fm_release(matrix_s_fm) - negf_qs_env%apply_bias = .FALSE. CALL timestop(handle) END SUBROUTINE converge_density @@ -1003,6 +1286,8 @@ CONTAINS ! ************************************************************************************************** !> \brief Compute the retarded Green's function and related properties at a set of points in parallel. !> \param omega list of energy points +!> \param v_shift shift in Hartree potential +!> \param ignore_bias ignore v_external from negf_control !> \param negf_env NEGF environment !> \param negf_control NEGF control !> \param sub_env (sub)group environment @@ -1028,10 +1313,13 @@ CONTAINS !> \par History !> * 07.2017 created [Sergey Chulkov] ! ************************************************************************************************** - SUBROUTINE negf_retarded_green_function_batch(omega, negf_env, negf_control, sub_env, ispin, g_surf_contacts, & + SUBROUTINE negf_retarded_green_function_batch(omega, v_shift, ignore_bias, negf_env, negf_control, sub_env, ispin, & + g_surf_contacts, & g_ret_s, g_ret_scale, gamma_contacts, gret_gamma_gadv, dos, & transm_coeff, transm_contact1, transm_contact2, just_contact) COMPLEX(kind=dp), DIMENSION(:), INTENT(in) :: omega + REAL(kind=dp), INTENT(in) :: v_shift + LOGICAL, INTENT(in) :: ignore_bias TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -1054,6 +1342,7 @@ CONTAINS INTEGER :: handle, icontact, igroup, ipoint, & ncontacts, npoints, nrows + REAL(kind=dp) :: v_external TYPE(copy_cfm_info_type), ALLOCATABLE, & DIMENSION(:) :: info1 TYPE(copy_cfm_info_type), ALLOCATABLE, & @@ -1062,6 +1351,7 @@ CONTAINS zwork1_contacts, zwork2_contacts TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:, :) :: gamma_contacts_group, & gret_gamma_gadv_group + TYPE(cp_cfm_type), POINTER :: matrix_s_cfm TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_fm_type), POINTER :: g_ret_imag, matrix_s TYPE(cp_para_env_type), POINTER :: para_env @@ -1071,6 +1361,8 @@ CONTAINS ncontacts = SIZE(negf_env%contacts) CPASSERT(SIZE(negf_control%contacts) == ncontacts) + NULLIFY (matrix_s_cfm) + IF (PRESENT(just_contact)) THEN CPASSERT(just_contact <= ncontacts) ncontacts = 2 @@ -1078,6 +1370,8 @@ CONTAINS CPASSERT(ncontacts >= 2) + IF (ignore_bias) v_external = 0.0_dp + IF (PRESENT(transm_coeff) .OR. PRESENT(transm_contact1) .OR. PRESENT(transm_contact2)) THEN CPASSERT(PRESENT(transm_coeff)) CPASSERT(PRESENT(transm_contact1)) @@ -1115,7 +1409,8 @@ CONTAINS END DO END IF - IF (PRESENT(g_ret_s) .OR. PRESENT(gret_gamma_gadv) .OR. PRESENT(dos) .OR. PRESENT(transm_coeff)) THEN + IF (PRESENT(g_ret_s) .OR. PRESENT(gret_gamma_gadv) .OR. & + PRESENT(dos) .OR. PRESENT(transm_coeff)) THEN ALLOCATE (g_ret_s_group(npoints)) IF (sub_env%ngroups > 1 .OR. .NOT. PRESENT(g_ret_s)) THEN @@ -1219,8 +1514,10 @@ CONTAINS ELSE ! contact self energies DO icontact = 1, ncontacts + IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external + CALL negf_contact_self_energy(self_energy_c=self_energy_contacts(icontact)%matrix, & - omega=omega(ipoint)-negf_control%contacts(icontact)%v_external, & + omega=omega(ipoint)-v_external, & g_surf_c=g_surf_contacts(icontact, ipoint)%matrix, & h_sc0=negf_env%h_sc(ispin, icontact)%matrix, & s_sc0=negf_env%s_sc(icontact)%matrix, & @@ -1247,16 +1544,25 @@ CONTAINS ! retarded Green's function for the scattering region IF (PRESENT(just_contact)) THEN CALL negf_retarded_green_function(g_ret_s=g_ret_s_group(ipoint)%matrix, & - omega=omega(ipoint), & + omega=omega(ipoint)-v_shift, & self_energy_ret_sum=self_energy_contacts(1)%matrix, & h_s=negf_env%contacts(just_contact)%h_00(ispin)%matrix, & - s_s=negf_env%contacts(just_contact)%s_00) - ELSE + s_s=negf_env%contacts(just_contact)%s_00, & + v_hartree_s=null()) + ELSE IF (ignore_bias) THEN CALL negf_retarded_green_function(g_ret_s=g_ret_s_group(ipoint)%matrix, & - omega=omega(ipoint), & + omega=omega(ipoint)-v_shift, & self_energy_ret_sum=self_energy_contacts(1)%matrix, & h_s=negf_env%h_s(ispin)%matrix, & - s_s=negf_env%s_s) + s_s=negf_env%s_s, & + v_hartree_s=null()) + ELSE + CALL negf_retarded_green_function(g_ret_s=g_ret_s_group(ipoint)%matrix, & + omega=omega(ipoint)-v_shift, & + self_energy_ret_sum=self_energy_contacts(1)%matrix, & + h_s=negf_env%h_s(ispin)%matrix, & + s_s=negf_env%s_s, & + v_hartree_s=negf_env%v_hartree_s) END IF IF (PRESENT(g_ret_scale)) THEN @@ -1449,11 +1755,13 @@ CONTAINS CALL cp_cfm_trace(g_ret_s_group(ipoint)%matrix, & self_energy_contacts(transm_contact2)%matrix, & transm_coeff(ipoint)) + IF (sub_env%para_env%mepos /= 0) transm_coeff(ipoint) = 0.0_dp END IF END DO + ! transmission coefficients are scaled by 2/pi CALL mp_sum(transm_coeff, sub_env%mpi_comm_global) - transm_coeff(:) = 1.0_dp/pi*transm_coeff(:) + !transm_coeff(:) = 0.5_dp/pi*transm_coeff(:) END IF ! -- deallocate temporary matrices @@ -1464,6 +1772,8 @@ CONTAINS DEALLOCATE (g_ret_s_group) END IF + IF (ASSOCIATED(matrix_s_cfm)) CALL cp_cfm_release(matrix_s_cfm) + IF (ALLOCATED(gamma_contacts_group)) THEN DO ipoint = npoints, 1, -1 DO icontact = ncontacts, 1, -1 @@ -1537,6 +1847,8 @@ CONTAINS ! ************************************************************************************************** !> \brief Compute contribution to the density matrix from the poles of the Fermi function. !> \param rho_ao_fm density matrix (initialised on exit) +!> \param v_shift shift in Hartree potential +!> \param ignore_bias ignore v_external from negf_control !> \param negf_env NEGF environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment @@ -1545,9 +1857,11 @@ CONTAINS !> \param just_contact ... !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE negf_init_rho_equiv_residuals(rho_ao_fm, negf_env, negf_control, sub_env, ispin, & - base_contact, just_contact) + SUBROUTINE negf_init_rho_equiv_residuals(rho_ao_fm, v_shift, ignore_bias, negf_env, & + negf_control, sub_env, ispin, base_contact, just_contact) TYPE(cp_fm_type), POINTER :: rho_ao_fm + REAL(kind=dp), INTENT(in) :: v_shift + LOGICAL, INTENT(in) :: ignore_bias TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -1560,7 +1874,8 @@ CONTAINS COMPLEX(kind=dp), ALLOCATABLE, DIMENSION(:) :: omega INTEGER :: handle, icontact, ipole, ncontacts, & npoles - REAL(kind=dp) :: mu_base, pi_temperature, temperature + REAL(kind=dp) :: mu_base, pi_temperature, temperature, & + v_external TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: g_ret_s TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_para_env_type), POINTER :: para_env @@ -1568,8 +1883,14 @@ CONTAINS CALL timeset(routineN, handle) - temperature = negf_control%temperature - mu_base = negf_control%contacts(base_contact)%fermi_level+negf_control%contacts(base_contact)%v_external + temperature = negf_control%contacts(base_contact)%temperature + IF (ignore_bias) THEN + mu_base = negf_control%contacts(base_contact)%fermi_level + v_external = 0.0_dp + ELSE + mu_base = negf_control%contacts(base_contact)%fermi_level+negf_control%contacts(base_contact)%v_external + END IF + pi_temperature = pi*temperature npoles = negf_control%delta_npoles @@ -1611,6 +1932,8 @@ CONTAINS END DO ELSE DO icontact = 1, ncontacts + IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external + CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, :), & omega=omega(:), & h0=negf_env%contacts(icontact)%h_00(ispin)%matrix, & @@ -1618,12 +1941,14 @@ CONTAINS h1=negf_env%contacts(icontact)%h_01(ispin)%matrix, & s1=negf_env%contacts(icontact)%s_01, & sub_env=sub_env, & - v_external=negf_control%contacts(icontact)%v_external, & + v_external=v_external, & conv=negf_control%conv_green, transp=.FALSE.) END DO END IF CALL negf_retarded_green_function_batch(omega=omega(:), & + v_shift=v_shift, & + ignore_bias=ignore_bias, & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -1655,6 +1980,8 @@ CONTAINS !> \brief Compute equilibrium contribution to the density matrix. !> \param rho_ao_fm density matrix (initialised on exit) !> \param stats integration statistics (updated on exit) +!> \param v_shift shift in Hartree potential +!> \param ignore_bias ignore v_external from negf_control !> \param negf_env NEGF environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment @@ -1665,13 +1992,16 @@ CONTAINS !> \param matrix_s_global globally distributed overlap matrix !> \param is_circular compute the integral along the circular path !> \param g_surf_cache set of precomputed surface Green's functions (updated on exit) -!> \param just_contact ... +!> \param just_contact ... !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE negf_add_rho_equiv_low(rho_ao_fm, stats, negf_env, negf_control, sub_env, ispin, base_contact, & - integr_lbound, integr_ubound, matrix_s_global, is_circular, g_surf_cache, just_contact) + SUBROUTINE negf_add_rho_equiv_low(rho_ao_fm, stats, v_shift, ignore_bias, negf_env, negf_control, sub_env, & + ispin, base_contact, integr_lbound, integr_ubound, matrix_s_global, & + is_circular, g_surf_cache, just_contact) TYPE(cp_fm_type), POINTER :: rho_ao_fm TYPE(integration_status_type), INTENT(inout) :: stats + REAL(kind=dp), INTENT(in) :: v_shift + LOGICAL, INTENT(in) :: ignore_bias TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -1689,7 +2019,8 @@ CONTAINS INTEGER :: handle, icontact, interval_id, ipoint, max_points, min_points, ncontacts, & npoints, npoints_exist, npoints_tmp, npoints_total, shape_id LOGICAL :: do_surface_green - REAL(kind=dp) :: conv_integr, mu_base, temperature + REAL(kind=dp) :: conv_integr, mu_base, temperature, & + v_external TYPE(ccquad_type) :: cc_env TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: zdata, zdata_tmp TYPE(cp_fm_struct_type), POINTER :: fm_struct @@ -1702,10 +2033,17 @@ CONTAINS ! convergence criteria for the integral of the retarded Green's function. This integral needs to be ! computed for both spin-components and needs to be scaled by -1/pi to obtain the electron density. conv_integr = 0.5_dp*negf_control%conv_density*pi - mu_base = negf_control%contacts(base_contact)%fermi_level+negf_control%contacts(base_contact)%v_external + + IF (ignore_bias) THEN + mu_base = negf_control%contacts(base_contact)%fermi_level + v_external = 0.0_dp + ELSE + mu_base = negf_control%contacts(base_contact)%fermi_level+negf_control%contacts(base_contact)%v_external + END IF + min_points = negf_control%integr_min_points max_points = negf_control%integr_max_points - temperature = negf_control%temperature + temperature = negf_control%contacts(base_contact)%temperature ncontacts = SIZE(negf_env%contacts) CPASSERT(base_contact <= ncontacts) @@ -1772,6 +2110,8 @@ CONTAINS END DO ELSE DO icontact = 1, ncontacts + IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external + CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, npoints_total+1:), & omega=xnodes(1:npoints), & h0=negf_env%contacts(icontact)%h_00(ispin)%matrix, & @@ -1779,7 +2119,7 @@ CONTAINS h1=negf_env%contacts(icontact)%h_01(ispin)%matrix, & s1=negf_env%contacts(icontact)%s_01, & sub_env=sub_env, & - v_external=negf_control%contacts(icontact)%v_external, & + v_external=v_external, & conv=negf_control%conv_green, transp=.FALSE.) END DO END IF @@ -1796,6 +2136,8 @@ CONTAINS END IF CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), & + v_shift=v_shift, & + ignore_bias=ignore_bias, & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -1898,6 +2240,8 @@ CONTAINS END DO ELSE DO icontact = 1, ncontacts + IF (.NOT. ignore_bias) v_external = negf_control%contacts(icontact)%v_external + CALL negf_surface_green_function_batch(g_surf=g_surf_cache%g_surf_contacts(icontact, npoints_total+1:), & omega=xnodes(1:npoints), & h0=negf_env%contacts(icontact)%h_00(ispin)%matrix, & @@ -1905,7 +2249,7 @@ CONTAINS h1=negf_env%contacts(icontact)%h_01(ispin)%matrix, & s1=negf_env%contacts(icontact)%s_01, & sub_env=sub_env, & - v_external=negf_control%contacts(icontact)%v_external, & + v_external=v_external, & conv=negf_control%conv_green, transp=.FALSE.) END DO END IF @@ -1920,6 +2264,8 @@ CONTAINS END IF CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), & + v_shift=v_shift, & + ignore_bias=ignore_bias, & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -1976,6 +2322,7 @@ CONTAINS !> \brief Compute non-equilibrium contribution to the density matrix. !> \param rho_ao_fm density matrix (initialised on exit) !> \param stats integration statistics (updated on exit) +!> \param v_shift shift in Hartree potential !> \param negf_env NEGF environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment @@ -1985,10 +2332,11 @@ CONTAINS !> \param g_surf_cache set of precomputed surface Green's functions (updated on exit) !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE negf_add_rho_nonequiv(rho_ao_fm, stats, negf_env, negf_control, sub_env, ispin, & - base_contact, matrix_s_global, g_surf_cache) + SUBROUTINE negf_add_rho_nonequiv(rho_ao_fm, stats, v_shift, negf_env, negf_control, sub_env, & + ispin, base_contact, matrix_s_global, g_surf_cache) TYPE(cp_fm_type), POINTER :: rho_ao_fm TYPE(integration_status_type), INTENT(inout) :: stats + REAL(kind=dp), INTENT(in) :: v_shift TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -2006,8 +2354,9 @@ CONTAINS max_points, min_points, ncontacts, & npoints, npoints_total LOGICAL :: do_surface_green - REAL(kind=dp) :: conv_density, conv_integr, eta, mu_base, & - mu_contact, temperature + REAL(kind=dp) :: conv_density, conv_integr, eta, & + ln_conv_density, mu_base, mu_contact, & + temperature_base, temperature_contact TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:, :) :: zdata TYPE(cp_fm_struct_type), POINTER :: fm_struct TYPE(cp_fm_type), POINTER :: integral_real @@ -2026,9 +2375,10 @@ CONTAINS mu_base = negf_control%contacts(base_contact)%fermi_level+negf_control%contacts(base_contact)%v_external min_points = negf_control%integr_min_points max_points = negf_control%integr_max_points - temperature = negf_control%temperature + temperature_base = negf_control%contacts(base_contact)%temperature eta = negf_control%eta conv_density = negf_control%conv_density + ln_conv_density = LOG(conv_density) ! convergence criteria for the integral. This integral needs to be computed for both ! spin-components and needs to be scaled by -1/pi to obtain the electron density. @@ -2037,9 +2387,12 @@ CONTAINS DO icontact = 1, ncontacts IF (icontact /= base_contact) THEN mu_contact = negf_control%contacts(icontact)%fermi_level+negf_control%contacts(icontact)%v_external + temperature_contact = negf_control%contacts(icontact)%temperature - integr_lbound = CMPLX(MIN(mu_base, mu_contact)+LOG(conv_density)*temperature, eta, kind=dp) - integr_ubound = CMPLX(MAX(mu_base, mu_contact)-LOG(conv_density)*temperature, eta, kind=dp) + integr_lbound = CMPLX(MIN(mu_base+ln_conv_density*temperature_base, & + mu_contact+ln_conv_density*temperature_contact), eta, kind=dp) + integr_ubound = CMPLX(MAX(mu_base-ln_conv_density*temperature_base, & + mu_contact-ln_conv_density*temperature_contact), eta, kind=dp) do_surface_green = .NOT. ALLOCATED(g_surf_cache%tnodes) @@ -2089,6 +2442,8 @@ CONTAINS END DO CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -2098,9 +2453,9 @@ CONTAINS DO ipoint = 1, npoints fermi_base = fermi_function(CMPLX(REAL(xnodes(ipoint), kind=dp)-mu_base, 0.0_dp, kind=dp), & - temperature) + temperature_base) fermi_contact = fermi_function(CMPLX(REAL(xnodes(ipoint), kind=dp)-mu_contact, 0.0_dp, kind=dp), & - temperature) + temperature_contact) CALL cp_cfm_scale(fermi_contact-fermi_base, zdata(icontact, ipoint)%matrix) END DO @@ -2191,17 +2546,19 @@ CONTAINS !> \brief Compute electric current for one spin-channel through the scattering region. !> \param contact_id1 reference contact !> \param contact_id2 another contact +!> \param v_shift shift in Hartree potential !> \param negf_env NEFG environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment !> \param ispin spin conponent to proceed !> \param blacs_env_global global BLACS environment -!> \return electric current +!> \return electric current in Amper !> \author Sergey Chulkov ! ************************************************************************************************** - FUNCTION negf_compute_current(contact_id1, contact_id2, negf_env, negf_control, sub_env, ispin, & + FUNCTION negf_compute_current(contact_id1, contact_id2, v_shift, negf_env, negf_control, sub_env, ispin, & blacs_env_global) RESULT(current) INTEGER, INTENT(in) :: contact_id1, contact_id2 + REAL(kind=dp), INTENT(in) :: v_shift TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -2220,9 +2577,8 @@ CONTAINS INTEGER :: handle, icontact, ipoint, max_points, & min_points, ncontacts, npoints, & npoints_total - REAL(kind=dp) :: conv_density, energy, eta, mu_contact1, & - mu_contact2, temperature, v_contact1, & - v_contact2 + REAL(kind=dp) :: conv_density, energy, eta, ln_conv_density, mu_contact1, mu_contact2, & + temperature_contact1, temperature_contact2, v_contact1, v_contact2 TYPE(cp_cfm_p_type), ALLOCATABLE, DIMENSION(:) :: zdata TYPE(cp_fm_struct_type), POINTER :: fm_struct_single TYPE(cp_fm_type), POINTER :: weights @@ -2252,12 +2608,16 @@ CONTAINS min_points = negf_control%integr_min_points max_points = negf_control%integr_max_points - temperature = negf_control%temperature + temperature_contact1 = negf_control%contacts(contact_id1)%temperature + temperature_contact2 = negf_control%contacts(contact_id2)%temperature eta = negf_control%eta conv_density = negf_control%conv_density + ln_conv_density = LOG(conv_density) - integr_lbound = CMPLX(MIN(mu_contact1, mu_contact2)+LOG(conv_density)*temperature, eta, kind=dp) - integr_ubound = CMPLX(MAX(mu_contact1, mu_contact2)-LOG(conv_density)*temperature, eta, kind=dp) + integr_lbound = CMPLX(MIN(mu_contact1+ln_conv_density*temperature_contact1, & + mu_contact2+ln_conv_density*temperature_contact2), eta, kind=dp) + integr_lbound = CMPLX(MAX(mu_contact1-ln_conv_density*temperature_contact1, & + mu_contact2-ln_conv_density*temperature_contact2), eta, kind=dp) npoints_total = 0 npoints = min_points @@ -2288,6 +2648,8 @@ CONTAINS END DO CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints), & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -2301,8 +2663,8 @@ CONTAINS CALL cp_cfm_create(zdata(ipoint)%matrix, fm_struct_single) energy = REAL(xnodes(ipoint), kind=dp) - fermi_contact1 = fermi_function(CMPLX(energy-mu_contact1, 0.0_dp, kind=dp), temperature) - fermi_contact2 = fermi_function(CMPLX(energy-mu_contact2, 0.0_dp, kind=dp), temperature) + fermi_contact1 = fermi_function(CMPLX(energy-mu_contact1, 0.0_dp, kind=dp), temperature_contact1) + fermi_contact2 = fermi_function(CMPLX(energy-mu_contact2, 0.0_dp, kind=dp), temperature_contact2) transmission(1, 1) = transm_coeff(ipoint)*(fermi_contact1-fermi_contact2) CALL cp_cfm_set_submatrix(zdata(ipoint)%matrix, transmission) @@ -2325,8 +2687,8 @@ CONTAINS CALL cp_cfm_get_submatrix(sr_env%integral, transmission) - ! an additional scale factor (1/pi) was already applied to transmission coefficients - current = 0.5_dp*REAL(transmission(1, 1), kind=dp)*e_charge/seconds + ! minus sign is due to negative charge of electrons + current = -0.5_dp/pi*REAL(transmission(1, 1), kind=dp)*e_charge/seconds CALL cp_fm_release(weights) CALL cp_fm_struct_release(fm_struct_single) @@ -2343,6 +2705,7 @@ CONTAINS !> \param energy_min energy point to start with !> \param energy_max energy point to end with !> \param npoints number of points to compute +!> \param v_shift shift in Hartree potential !> \param negf_env NEFG environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment @@ -2351,11 +2714,12 @@ CONTAINS !> \param volume unit cell volume !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE negf_print_dos(log_unit, energy_min, energy_max, npoints, negf_env, & + SUBROUTINE negf_print_dos(log_unit, energy_min, energy_max, npoints, v_shift, negf_env, & negf_control, sub_env, base_contact, just_contact, volume) INTEGER, INTENT(in) :: log_unit REAL(kind=dp), INTENT(in) :: energy_min, energy_max INTEGER, INTENT(in) :: npoints + REAL(kind=dp), INTENT(in) :: v_shift TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -2457,6 +2821,8 @@ CONTAINS END IF CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints_bundle), & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & @@ -2493,6 +2859,7 @@ CONTAINS !> \param energy_min energy point to start with !> \param energy_max energy point to end with !> \param npoints number of points to compute +!> \param v_shift shift in Hartree potential !> \param negf_env NEFG environment !> \param negf_control NEGF control !> \param sub_env NEGF parallel (sub)group environment @@ -2500,11 +2867,12 @@ CONTAINS !> \param contact_id2 index of another contact !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE negf_print_transmission(log_unit, energy_min, energy_max, npoints, negf_env, & + SUBROUTINE negf_print_transmission(log_unit, energy_min, energy_max, npoints, v_shift, negf_env, & negf_control, sub_env, contact_id1, contact_id2) INTEGER, INTENT(in) :: log_unit REAL(kind=dp), INTENT(in) :: energy_min, energy_max INTEGER, INTENT(in) :: npoints + REAL(kind=dp), INTENT(in) :: v_shift TYPE(negf_env_type), INTENT(in) :: negf_env TYPE(negf_control_type), POINTER :: negf_control TYPE(negf_subgroup_env_type), INTENT(in) :: sub_env @@ -2527,7 +2895,7 @@ CONTAINS IF (log_unit > 0) THEN WRITE (log_unit, '(A)') "# Transmission coefficient for the scattering region" - WRITE (log_unit, '(A,T10,A,T40,A)') "#", "Energy (a.u.)", "Transmission coefficient (a.u.)" + WRITE (log_unit, '(A,T10,A,T40,A)') "#", "Energy (a.u.)", "Transmission coefficient (G0 = 2 e^2/h)" WRITE (log_unit, '("#", T3,78("-"))') END IF @@ -2542,6 +2910,10 @@ CONTAINS rscale = 1.0_dp END IF + ! print transmission coefficients in terms of G0 = 2 * e^2 / h = 1 / pi ; + ! transmission coefficients returned by negf_retarded_green_function_batch() are already multiplied by 2 / pi + rscale = 0.5_dp*rscale + npoints_bundle = 4*sub_env%ngroups IF (npoints_bundle > npoints) npoints_bundle = npoints @@ -2576,6 +2948,8 @@ CONTAINS END DO CALL negf_retarded_green_function_batch(omega=xnodes(1:npoints_bundle), & + v_shift=v_shift, & + ignore_bias=.FALSE., & negf_env=negf_env, & negf_control=negf_control, & sub_env=sub_env, & diff --git a/src/negf_qs_env_methods.F b/src/negf_qs_env_methods.F deleted file mode 100644 index d696fc1556..0000000000 --- a/src/negf_qs_env_methods.F +++ /dev/null @@ -1,190 +0,0 @@ -!--------------------------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 2000 - 2018 CP2K developers group ! -!--------------------------------------------------------------------------------------------------! - -! ************************************************************************************************** -!> \brief Derived types for NEGF based quantum transport calculations -! ************************************************************************************************** -MODULE negf_qs_env_methods - USE dbcsr_api, ONLY: dbcsr_get_block_p,& - dbcsr_p_type - USE kinds, ONLY: dp - USE negf_qs_env_types, ONLY: negf_qs_env_contact_type,& - negf_qs_env_type - USE negf_vectors, ONLY: projection_on_direction_vector - USE pw_types, ONLY: pw_p_type,& - pw_type -#include "./base/base_uses.f90" - - IMPLICIT NONE - PRIVATE - - CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_qs_env_methods' - LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE. - - PUBLIC :: negf_update_ks_matrix, negf_modify_hartree_pot - -CONTAINS - -! ************************************************************************************************** -!> \brief Ensure that Kohn-Sham matrix blocks which corresponds to a bulk contact region are -!> identical to the actual Kohn-Sham matrix blocks from bulk calculation. -!> \param matrix_ks_kp Kohn-Sham matrix to update -!> \param negf_qs_env NEGF_QS environment -!> \author Sergey Chulkov -! ************************************************************************************************** - SUBROUTINE negf_update_ks_matrix(matrix_ks_kp, negf_qs_env) - TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp - TYPE(negf_qs_env_type), POINTER :: negf_qs_env - - CHARACTER(len=*), PARAMETER :: routineN = 'negf_update_ks_matrix', & - routineP = moduleN//':'//routineN - - INTEGER :: block_index, handle, iatom_col, iatom_row, iblock, icontact, image, ispin, & - ispin_contact, nblocks, ncontacts, nimages, nspins, nspins_contacts - LOGICAL :: found, is_transpose - REAL(kind=dp), DIMENSION(:, :), POINTER :: sm_block - TYPE(negf_qs_env_contact_type), POINTER :: contact - - CALL timeset(routineN, handle) - - IF (ALLOCATED(negf_qs_env%contacts)) THEN - nspins = SIZE(matrix_ks_kp, 1) - nimages = SIZE(matrix_ks_kp, 2) - - ncontacts = SIZE(negf_qs_env%contacts) - - DO icontact = 1, ncontacts - contact => negf_qs_env%contacts(icontact) - nblocks = SIZE(contact%ks_block_map) - nspins_contacts = SIZE(contact%ks_block_storage, 2) - - ispin_contact = 0 - DO ispin = 1, nspins - IF (ispin_contact < nspins_contacts) & - ispin_contact = ispin_contact+1 - - DO iblock = 1, nblocks - iatom_row = contact%ks_block_map(iblock)%iatom_row - iatom_col = contact%ks_block_map(iblock)%iatom_col - image = contact%ks_block_map(iblock)%image - - CPASSERT(image <= nimages) - - block_index = contact%ks_block_map(iblock)%block_index - is_transpose = block_index < 0 - block_index = ABS(block_index) - - CALL dbcsr_get_block_p(matrix=matrix_ks_kp(ispin, image)%matrix, & - row=iatom_row, col=iatom_col, block=sm_block, found=found) - CPASSERT(found) - - IF (block_index > 0) THEN - IF (is_transpose) THEN - sm_block(:, :) = & - TRANSPOSE(contact%ks_block_storage(block_index, ispin_contact)%matrix) - ELSE - sm_block(:, :) = & - contact%ks_block_storage(block_index, ispin_contact)%matrix(:, :) - END IF - ELSE - sm_block(:, :) = 0.0_dp - END IF - END DO - END DO - END DO - END IF - - CALL timestop(handle) - END SUBROUTINE negf_update_ks_matrix - -! ************************************************************************************************** -!> \brief Modify the Hartree potential by applying an external potential. -!> \param v_hartree Hartree potential to update -!> \param negf_qs_env NEGF_QS environment -!> \author Sergey Chulkov -! ************************************************************************************************** - SUBROUTINE negf_modify_hartree_pot(v_hartree, negf_qs_env) - TYPE(pw_p_type), INTENT(inout) :: v_hartree - TYPE(negf_qs_env_type), POINTER :: negf_qs_env - - CHARACTER(len=*), PARAMETER :: routineN = 'negf_modify_hartree_pot', & - routineP = moduleN//':'//routineN - REAL(kind=dp), PARAMETER :: threshold = 16.0_dp*EPSILON(0.0_dp) - - INTEGER :: dx, dy, dz, handle, icontact, ix, iy, & - iz, lx, ly, lz, ncontacts, ux, uy, uz - REAL(kind=dp) :: dvol, pot, proj, v1, v2 - REAL(kind=dp), DIMENSION(3) :: dirvector_bias, point_coord, & - point_indices, vector - TYPE(pw_type), POINTER :: v_hartree_ref - - CALL timeset(routineN, handle) - - IF (ALLOCATED(negf_qs_env%contacts)) THEN - ncontacts = SIZE(negf_qs_env%contacts) - CPASSERT(ncontacts == 2) - - dirvector_bias = negf_qs_env%contacts(2)%origin_scattering-negf_qs_env%contacts(1)%origin_scattering - v1 = negf_qs_env%contacts(1)%v_external - v2 = negf_qs_env%contacts(2)%v_external - - lx = v_hartree%pw%pw_grid%bounds_local(1, 1) - ux = v_hartree%pw%pw_grid%bounds_local(2, 1) - ly = v_hartree%pw%pw_grid%bounds_local(1, 2) - uy = v_hartree%pw%pw_grid%bounds_local(2, 2) - lz = v_hartree%pw%pw_grid%bounds_local(1, 3) - uz = v_hartree%pw%pw_grid%bounds_local(2, 3) - - dx = v_hartree%pw%pw_grid%npts(1)/2 - dy = v_hartree%pw%pw_grid%npts(2)/2 - dz = v_hartree%pw%pw_grid%npts(3)/2 - - dvol = v_hartree%pw%pw_grid%dvol - v_hartree_ref => negf_qs_env%v_hartree_ref - - DO iz = lz, uz - point_indices(3) = REAL(iz+dz, kind=dp) - DO iy = ly, uy - point_indices(2) = REAL(iy+dy, kind=dp) - - DO ix = lx, ux - point_indices(1) = REAL(ix+dx, kind=dp) - point_coord(:) = MATMUL(v_hartree%pw%pw_grid%dh, point_indices) - - vector = point_coord-negf_qs_env%contacts(1)%origin_scattering - proj = projection_on_direction_vector(vector, dirvector_bias) - IF (proj+threshold >= 0.0_dp .AND. proj-threshold <= 1.0_dp) THEN - ! scattering region - ! proj == 0 we are at the first contact boundary - ! proj == 1 we are at the second contact boundary - IF (proj < 0.0_dp) THEN - proj = 0.0_dp - ELSE IF (proj > 1.0_dp) THEN - proj = 1.0_dp - END IF - pot = v1+(v2-v1)*proj - v_hartree%pw%cr3d(ix, iy, iz) = v_hartree%pw%cr3d(ix, iy, iz)+pot*dvol - ELSE - pot = 0.0_dp - DO icontact = 1, ncontacts - vector = point_coord-negf_qs_env%contacts(icontact)%origin_scattering - proj = projection_on_direction_vector(vector, negf_qs_env%contacts(icontact)%direction_vector) - - IF (proj+threshold >= 0.0_dp .AND. proj-threshold <= 1.0_dp) THEN - pot = negf_qs_env%contacts(icontact)%v_external - EXIT - END IF - END DO - END IF - - v_hartree%pw%cr3d(ix, iy, iz) = v_hartree_ref%cr3d(ix, iy, iz)+pot*dvol - END DO - END DO - END DO - END IF - - CALL timestop(handle) - END SUBROUTINE negf_modify_hartree_pot -END MODULE negf_qs_env_methods diff --git a/src/negf_qs_env_types.F b/src/negf_qs_env_types.F deleted file mode 100644 index e904976043..0000000000 --- a/src/negf_qs_env_types.F +++ /dev/null @@ -1,151 +0,0 @@ -!--------------------------------------------------------------------------------------------------! -! CP2K: A general program to perform molecular dynamics simulations ! -! Copyright (C) 2000 - 2018 CP2K developers group ! -!--------------------------------------------------------------------------------------------------! - -! ************************************************************************************************** -!> \brief Part of the NEGF environment that needs to be stored within qs_environment. -!> It is used to update the Hartree potential and the Kohn-Sham matrix within qs_scf(). -! ************************************************************************************************** -MODULE negf_qs_env_types - USE kinds, ONLY: dp - USE negf_alloc_types, ONLY: negf_allocatable_rmatrix - USE pw_env_types, ONLY: pw_env_get,& - pw_env_type - USE pw_pool_types, ONLY: pw_pool_give_back_pw,& - pw_pool_type - USE pw_types, ONLY: pw_type -#include "./base/base_uses.f90" - - IMPLICIT NONE - PRIVATE - - CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_qs_env_types' - LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE. - -! ************************************************************************************************** -!> \brief ... -!> \author Sergey Chulkov -! ************************************************************************************************** - TYPE negf_matrix_block_map_type - !> column and row atomic indices - INTEGER :: iatom_col, iatom_row - !> cell replica - INTEGER :: image - !> absolute value is the index of the block within the block data storage; - !> 0 all matrix elements of the given block are equal to zero, - !> >0 block is stored simply as block_storage(block_index), - !> <0 block is stored as TRANSPOSE(block_storage(-block_index)) - INTEGER :: block_index - END TYPE negf_matrix_block_map_type - -! ************************************************************************************************** -!> \brief Part of the NEGF environment stored within the QuickStep environment. -!> \author Sergey Chulkov -! ************************************************************************************************** - TYPE negf_qs_env_contact_type - REAL(kind=dp), DIMENSION(3) :: direction_vector, origin_scattering - REAL(kind=dp) :: v_external - !> atom_row, atom_col, image, and the index of the actual data block in 'block_storage' - !> all indices are identical for both spins - TYPE(negf_matrix_block_map_type), ALLOCATABLE, & - DIMENSION(:) :: ks_block_map - !> collection of DBSCR Kohn-Sham matrix blocks ([index, spin]) - TYPE(negf_allocatable_rmatrix), ALLOCATABLE, & - DIMENSION(:, :) :: ks_block_storage - END TYPE negf_qs_env_contact_type - -! ************************************************************************************************** -!> \brief data to update the Hartree potential and the Kohn-Sham matrix within qs_scf(). -!> \author Sergey Chulkov -! ************************************************************************************************** - TYPE negf_qs_env_type - LOGICAL :: apply_bias - LOGICAL :: update_ks_matrix - TYPE(negf_qs_env_contact_type), ALLOCATABLE, & - DIMENSION(:) :: contacts - !> reference Hartree potential (with no bias) - TYPE(pw_type), POINTER :: v_hartree_ref - END TYPE negf_qs_env_type - - PUBLIC :: negf_qs_env_type, negf_qs_env_contact_type, negf_matrix_block_map_type - PUBLIC :: negf_qs_env_create, negf_qs_env_release - -CONTAINS - -! ************************************************************************************************** -!> \brief Allocate control options for Non-equilibrium Green's Function calculation. -!> \param negf_qs_env an object to create -!> \par History -!> * 10.2017 created [Sergey Chulkov] -! ************************************************************************************************** - SUBROUTINE negf_qs_env_create(negf_qs_env) - TYPE(negf_qs_env_type), POINTER :: negf_qs_env - - CHARACTER(len=*), PARAMETER :: routineN = 'negf_qs_env_create', & - routineP = moduleN//':'//routineN - - INTEGER :: handle - - CPASSERT(.NOT. ASSOCIATED(negf_qs_env)) - CALL timeset(routineN, handle) - - ALLOCATE (negf_qs_env) - negf_qs_env%apply_bias = .FALSE. - negf_qs_env%update_ks_matrix = .FALSE. - NULLIFY (negf_qs_env%v_hartree_ref) - - CALL timestop(handle) - END SUBROUTINE negf_qs_env_create - -! ************************************************************************************************** -!> \brief Release memory allocated for NEGF control options. -!> \param negf_qs_env an object to release -!> \param pw_env plane-wave environment used to create the reference Hartree potential -!> \par History -!> * 10.2017 created [Sergey Chulkov] -! ************************************************************************************************** - SUBROUTINE negf_qs_env_release(negf_qs_env, pw_env) - TYPE(negf_qs_env_type), POINTER :: negf_qs_env - TYPE(pw_env_type), POINTER :: pw_env - - CHARACTER(len=*), PARAMETER :: routineN = 'negf_qs_env_release', & - routineP = moduleN//':'//routineN - - INTEGER :: handle, iblock, icontact, ispin - TYPE(pw_pool_type), POINTER :: pw_pool - - CALL timeset(routineN, handle) - - IF (ASSOCIATED(negf_qs_env)) THEN - IF (ALLOCATED(negf_qs_env%contacts)) THEN - DO icontact = SIZE(negf_qs_env%contacts), 1, -1 - IF (ALLOCATED(negf_qs_env%contacts(icontact)%ks_block_map)) & - DEALLOCATE (negf_qs_env%contacts(icontact)%ks_block_map) - - IF (ALLOCATED(negf_qs_env%contacts(icontact)%ks_block_storage)) THEN - DO ispin = SIZE(negf_qs_env%contacts(icontact)%ks_block_storage, 2), 1, -1 - DO iblock = SIZE(negf_qs_env%contacts(icontact)%ks_block_storage, 1), 1, -1 - IF (ALLOCATED(negf_qs_env%contacts(icontact)%ks_block_storage(iblock, ispin)%matrix)) & - DEALLOCATE (negf_qs_env%contacts(icontact)%ks_block_storage(iblock, ispin)%matrix) - END DO - END DO - - DEALLOCATE (negf_qs_env%contacts(icontact)%ks_block_storage) - END IF - END DO - - DEALLOCATE (negf_qs_env%contacts) - END IF - - IF (ASSOCIATED(negf_qs_env%v_hartree_ref)) THEN - CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool) - CALL pw_pool_give_back_pw(pw_pool, negf_qs_env%v_hartree_ref) - END IF - - DEALLOCATE (negf_qs_env) - END IF - - CALL timestop(handle) - END SUBROUTINE negf_qs_env_release -END MODULE negf_qs_env_types diff --git a/src/negf_vectors.F b/src/negf_vectors.F index 30b65055f6..dcd74126c8 100644 --- a/src/negf_vectors.F +++ b/src/negf_vectors.F @@ -25,15 +25,22 @@ CONTAINS ! ************************************************************************************************** !> \brief compute direction vector of the given contact -!> \param origin origin -!> \param direction_vector direction vector -!> \param atomlist_screening atoms belonging to the contact's screening region -!> \param atomlist_bulk atoms belonging to the contact's bulk region -!> \param subsys QuickStep subsystem +!> \param origin origin +!> \param direction_vector direction vector +!> \param origin_bias origin which will be used to apply the external bias +!> (in contrast with 'origin' it does not include screening region) +!> \param direction_vector_bias direction vector which will be used to apply the external bias +!> (together with 'origin_bias' it defines a contact region where +!> the external potential is kept constant) +!> \param atomlist_screening atoms belonging to the contact's screening region +!> \param atomlist_bulk atoms belonging to the contact's bulk region +!> \param subsys QuickStep subsystem !> \author Sergey Chulkov ! ************************************************************************************************** - SUBROUTINE contact_direction_vector(origin, direction_vector, atomlist_screening, atomlist_bulk, subsys) - REAL(kind=dp), DIMENSION(3), INTENT(out) :: origin, direction_vector + SUBROUTINE contact_direction_vector(origin, direction_vector, origin_bias, direction_vector_bias, & + atomlist_screening, atomlist_bulk, subsys) + REAL(kind=dp), DIMENSION(3), INTENT(out) :: origin, direction_vector, origin_bias, & + direction_vector_bias INTEGER, DIMENSION(:), INTENT(in) :: atomlist_screening, atomlist_bulk TYPE(qs_subsys_type), POINTER :: subsys @@ -42,7 +49,7 @@ CONTAINS INTEGER :: handle, iatom, natoms_bulk, & natoms_screening, nparticles - REAL(kind=dp) :: proj, proj_max, proj_min + REAL(kind=dp) :: proj, proj_max, proj_min, proj_min_bias REAL(kind=dp), DIMENSION(3) :: vector TYPE(particle_type), DIMENSION(:), POINTER :: particle_set @@ -77,11 +84,13 @@ CONTAINS ! the point 'center_of_coords0' belongs to the first unit cell, so the lowest projection of any point ! from the first unit cell on the direction vector 'center_of_coords1 - center_of_coords0' should be <= 0 proj_min = 0.0_dp + proj_min_bias = 0.0_dp DO iatom = 1, natoms_screening vector = particle_set(atomlist_screening(iatom))%r-origin proj = projection_on_direction_vector(vector, direction_vector) IF (proj < proj_min) proj_min = proj + IF (proj > proj_min_bias) proj_min_bias = proj END DO ! the point 'center_of_coords1' belongs to the given contact, so the highest projection should be >= 1 @@ -93,11 +102,13 @@ CONTAINS IF (proj > proj_max) proj_max = proj END DO - ! adjust the origin, so it lies on a plane between the given contact and the scaterring region + ! adjust the origin, so it lies on a plane between the given contact and the scattering region + origin_bias = origin+proj_min_bias*direction_vector origin = origin+proj_min*direction_vector ! rescale the vector, so the last atom of the given contact and the point 'origin + direction_vector' lie on - ! the same plane parallel to the 'origin' plane -- which separates the contact from the scaterring region. + ! the same plane parallel to the 'origin' plane -- which separates the contact from the scattering region. + direction_vector_bias = (proj_max-proj_min_bias)*direction_vector direction_vector = (proj_max-proj_min)*direction_vector CALL timestop(handle) diff --git a/src/qs_environment_types.F b/src/qs_environment_types.F index 27ee2085f5..3b695d7943 100644 --- a/src/qs_environment_types.F +++ b/src/qs_environment_types.F @@ -77,9 +77,6 @@ MODULE qs_environment_types mp2_type USE mscfg_types, ONLY: molecular_scf_guess_env_destroy,& molecular_scf_guess_env_type - USE negf_qs_env_types, ONLY: negf_qs_env_create,& - negf_qs_env_release,& - negf_qs_env_type USE particle_types, ONLY: particle_type USE pw_env_types, ONLY: pw_env_type USE pw_types, ONLY: pw_p_type,& @@ -297,8 +294,6 @@ MODULE qs_environment_types ! Embedding potential TYPE(pw_p_type), POINTER :: embed_pot TYPE(pw_p_type), POINTER :: spin_embed_pot - ! NEGF external potential - TYPE(negf_qs_env_type), POINTER :: negf_env END TYPE qs_environment_type ! ************************************************************************************************** @@ -474,7 +469,6 @@ CONTAINS !> \param variable_history ... !> \param embed_pot ... !> \param spin_embed_pot ... -!> \param negf_env data needed for NEGF calculation !> \date 23.01.2002 !> \author MK !> \version 1.0 @@ -505,7 +499,7 @@ CONTAINS lri_env, lri_density, dispersion_env, vee, rho_external, external_vxc, mask, & mp2_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, & s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, & - gradient_history, variable_history, embed_pot, spin_embed_pot, negf_env) + gradient_history, variable_history, embed_pot, spin_embed_pot) TYPE(qs_environment_type), POINTER :: qs_env TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, & POINTER :: atomic_kind_set @@ -636,7 +630,6 @@ CONTAINS TYPE(almo_scf_env_type), OPTIONAL, POINTER :: almo_scf_env REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: gradient_history, variable_history TYPE(pw_p_type), OPTIONAL, POINTER :: embed_pot, spin_embed_pot - TYPE(negf_qs_env_type), OPTIONAL, POINTER :: negf_env CHARACTER(len=*), PARAMETER :: routineN = 'get_qs_env', routineP = moduleN//':'//routineN @@ -702,9 +695,6 @@ CONTAINS IF (PRESENT(embed_pot)) embed_pot => qs_env%embed_pot IF (PRESENT(spin_embed_pot)) spin_embed_pot => qs_env%spin_embed_pot - ! NEGF - IF (PRESENT(negf_env)) negf_env => qs_env%negf_env - IF (PRESENT(rho_atom_set)) & CALL get_local_rho(qs_env%local_rho_set, rho_atom_set=rho_atom_set) IF (PRESENT(rho0_atom_set)) & @@ -869,10 +859,6 @@ CONTAINS NULLIFY (qs_env%embed_pot) NULLIFY (qs_env%spin_embed_pot) - ! NEGF - NULLIFY (qs_env%negf_env) - CALL negf_qs_env_create(qs_env%negf_env) - NULLIFY (qs_env%ks_env) NULLIFY (qs_env%ks_qmmm_env) NULLIFY (qs_env%wf_history) @@ -1318,7 +1304,6 @@ CONTAINS CHARACTER(len=*), PARAMETER :: routineN = 'qs_env_release', routineP = moduleN//':'//routineN INTEGER :: i - TYPE(pw_env_type), POINTER :: pw_env IF (ASSOCIATED(qs_env)) THEN CPASSERT(qs_env%ref_count > 0) @@ -1428,12 +1413,6 @@ CONTAINS CALL release_active_space_type(qs_env%active_space) ENDIF - ! NEGF - IF (ASSOCIATED(qs_env%negf_env)) THEN - CALL get_ks_env(qs_env%ks_env, pw_env=pw_env) - CALL negf_qs_env_release(qs_env%negf_env, pw_env) - END IF - CALL qs_charges_release(qs_env%qs_charges) CALL qs_ks_release(qs_env%ks_env) CALL qs_ks_qmmm_release(qs_env%ks_qmmm_env) diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index b3fbeff888..3928633402 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -79,9 +79,6 @@ MODULE qs_ks_methods lri_kind_type USE mathlib, ONLY: abnormal_value USE message_passing, ONLY: mp_sum - USE negf_qs_env_methods, ONLY: negf_modify_hartree_pot,& - negf_update_ks_matrix - USE negf_qs_env_types, ONLY: negf_qs_env_type USE pw_env_types, ONLY: pw_env_get,& pw_env_type USE pw_methods, ONLY: pw_axpy,& @@ -235,7 +232,6 @@ CONTAINS TYPE(lri_density_type), POINTER :: lri_density TYPE(lri_environment_type), POINTER :: lri_env TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int - TYPE(negf_qs_env_type), POINTER :: negf_qs_env TYPE(pw_env_type), POINTER :: pw_env TYPE(pw_p_type) :: rho_tot_gspace, v_efield_rspace, v_hartree_gspace, v_hartree_rspace, & v_minus_veps, v_sccs_rspace, v_sic_rspace, v_spin_ddapc_rest_r @@ -289,8 +285,7 @@ CONTAINS rho=rho, & rho_core=rho_core, & rho_xc=rho_xc, & - energy=energy, & - negf_env=negf_qs_env) + energy=energy) CALL qs_rho_get(rho, rho_r=rho_r, rho_ao_kp=rho_ao) @@ -534,11 +529,6 @@ CONTAINS END IF CALL pw_pool_give_back_pw(auxbas_pw_pool, rho_tot_gspace%pw) - ! external NEGF potential - IF (negf_qs_env%apply_bias) THEN - CALL negf_modify_hartree_pot(v_hartree_rspace, negf_qs_env) - END IF - ! calculate the density matrix for the fitted mo_coeffs IF (dft_control%do_admm) THEN CALL hfx_admm_init(qs_env) @@ -860,12 +850,6 @@ CONTAINS END IF END IF - IF (.NOT. just_energy .AND. negf_qs_env%update_ks_matrix) THEN - ! ensure that Konh-Sham matrix elements which correspond to the bulk - ! contact region are identical to the actual bulk calculation - CALL negf_update_ks_matrix(ks_matrix, negf_qs_env) - END IF - ! Sum all energy terms to obtain the total energy energy%total = energy%core_overlap+energy%core_self+energy%core+energy%hartree+ & energy%hartree_1c+energy%exc+energy%exc1+energy%ex+ & diff --git a/tests/QS/regtest-negf/hchain.inp b/tests/QS/regtest-negf/hchain.inp deleted file mode 100644 index a363af32db..0000000000 --- a/tests/QS/regtest-negf/hchain.inp +++ /dev/null @@ -1,341 +0,0 @@ -&GLOBAL - PROJECT hchain - RUN_TYPE negf - PRINT_LEVEL low - EXTENDED_FFT_LENGTHS .true. -&END GLOBAL - -&MULTIPLE_FORCE_EVALS - FORCE_EVAL_ORDER 1 2 3 - MULTIPLE_SUBSYS .true. -&END - -&NEGF - &CONTACT - &BULK_CONTACT - MOLNAME L1 L2 - &END - &SCREENING_REGION - MOLNAME L0 - &END - - # indices starting from 0 -- primary force_env - BULK_FORCE_EVAL_SECTION 1 - - #FERMI_LEVEL 0.09790647 - ELECTRIC_POTENTIAL [eV] 0 - - &PRINT - &DOS - FILENAME lead1 - &END - &END - &END CONTACT - - &CONTACT - &BULK_CONTACT - MOLNAME R1 R2 - &END - &SCREENING_REGION - MOLNAME R0 - &END - - BULK_FORCE_EVAL_SECTION 2 - #FERMI_LEVEL 0.09790647 - ELECTRIC_POTENTIAL [eV] 0 - &END CONTACT - - &SCATTERING_REGION - MOLNAME S - &END SCATTERING_REGION - - &MIXING - ALPHA 0.2 - METHOD BROYDEN_MIXING - NBUFFER 8 - &END MIXING - - INTEGRATION_METHOD simpson - INTEGRATION_MIN_POINTS 16 - INTEGRATION_MAX_POINTS 128 - - ENERGY_LBOUND -2.0 - EPS_DENSITY 1e-3 - EPS_SCF 1e-3 - MAX_SCF 100 - - NPROC_POINT 1 -&END NEGF - - -&FORCE_EVAL - METHOD Quickstep - - &DFT - BASIS_SET_FILE_NAME BASIS_MOLOPT - POTENTIAL_FILE_NAME GTH_POTENTIALS - - &MGRID - CUTOFF 300 - &END MGRID - - &QS - METHOD gpw - EPS_DEFAULT 1e-10 - &END QS - - &SCF - MAX_SCF 10 - EPS_SCF 1e-6 - SCF_GUESS atomic - ADDED_MOS 100 - - &SMEAR - ELECTRONIC_TEMPERATURE 293 - METHOD FERMI_DIRAC - &END SMEAR - - &MIXING - METHOD BROYDEN_MIXING - ALPHA 0.3 - NBUFFER 8 - &END MIXING - - &PRINT - &RESTART - FILENAME device - &END RESTART - &END PRINT - &END SCF - - &POISSON - PERIODIC xyz - &END POISSON - - &XC - &XC_FUNCTIONAL pbe - &END XC_FUNCTIONAL - &END XC - &END DFT - - &SUBSYS - &CELL - ABC 8.00000 8.00000 40.7000 - PERIODIC xyz - &END CELL - - &COORD -H 0.00000000 0.00000000 0.00000000 -H 0.00000000 0.00000000 0.74000000 -H 0.00000000 0.00000000 1.48000000 -H 0.00000000 0.00000000 2.22000000 -H 0.00000000 0.00000000 2.96000000 -H 0.00000000 0.00000000 3.70000000 -H 0.00000000 0.00000000 4.44000000 -H 0.00000000 0.00000000 5.18000000 -H 0.00000000 0.00000000 5.92000000 -H 0.00000000 0.00000000 6.66000000 -H 0.00000000 0.00000000 7.40000000 L2 -H 0.00000000 0.00000000 8.14000000 L2 -H 0.00000000 0.00000000 8.88000000 L2 -H 0.00000000 0.00000000 9.62000000 L2 -H 0.00000000 0.00000000 10.36000000 L2 -H 0.00000000 0.00000000 11.10000000 L1 -H 0.00000000 0.00000000 11.84000000 L1 -H 0.00000000 0.00000000 12.58000000 L1 -H 0.00000000 0.00000000 13.32000000 L1 -H 0.00000000 0.00000000 14.06000000 L1 -H 0.00000000 0.00000000 14.80000000 L0 -H 0.00000000 0.00000000 15.54000000 L0 -H 0.00000000 0.00000000 16.28000000 L0 -H 0.00000000 0.00000000 17.02000000 L0 -H 0.00000000 0.00000000 17.76000000 L0 -H 0.00000000 0.00000000 18.50000000 S -H 0.00000000 0.00000000 19.24000000 S -H 0.00000000 0.00000000 19.98000000 S -H 0.00000000 0.00000000 20.72000000 S -H 0.00000000 0.00000000 21.46000000 S -H 0.00000000 0.00000000 22.20000000 R0 -H 0.00000000 0.00000000 22.94000000 R0 -H 0.00000000 0.00000000 23.68000000 R0 -H 0.00000000 0.00000000 24.42000000 R0 -H 0.00000000 0.00000000 25.16000000 R0 -H 0.00000000 0.00000000 25.90000000 R1 -H 0.00000000 0.00000000 26.64000000 R1 -H 0.00000000 0.00000000 27.38000000 R1 -H 0.00000000 0.00000000 28.12000000 R1 -H 0.00000000 0.00000000 28.86000000 R1 -H 0.00000000 0.00000000 29.60000000 R2 -H 0.00000000 0.00000000 30.34000000 R2 -H 0.00000000 0.00000000 31.08000000 R2 -H 0.00000000 0.00000000 31.82000000 R2 -H 0.00000000 0.00000000 32.56000000 R2 -H 0.00000000 0.00000000 33.30000000 -H 0.00000000 0.00000000 34.04000000 -H 0.00000000 0.00000000 34.78000000 -H 0.00000000 0.00000000 35.52000000 -H 0.00000000 0.00000000 36.26000000 -H 0.00000000 0.00000000 37.00000000 -H 0.00000000 0.00000000 37.74000000 -H 0.00000000 0.00000000 38.48000000 -H 0.00000000 0.00000000 39.22000000 -H 0.00000000 0.00000000 39.96000000 - &END COORD - - &KIND H - BASIS_SET SZV-MOLOPT-GTH - POTENTIAL GTH-PBE-q1 - &END KIND - &END SUBSYS -&END FORCE_EVAL - -&FORCE_EVAL - METHOD Quickstep - - &DFT - BASIS_SET_FILE_NAME BASIS_MOLOPT - POTENTIAL_FILE_NAME GTH_POTENTIALS - WFN_RESTART_FILE_NAME hchain-lead1.kp - - &MGRID - CUTOFF 300 - &END MGRID - - &QS - METHOD gpw - EPS_DEFAULT 1e-10 - &END QS - - &SCF - MAX_SCF 10 - EPS_SCF 1e-6 - SCF_GUESS atomic - ADDED_MOS 100 - - &SMEAR - ELECTRONIC_TEMPERATURE 293 - METHOD FERMI_DIRAC - &END SMEAR - - &MIXING - METHOD BROYDEN_MIXING - ALPHA 0.3 - NBUFFER 8 - &END MIXING - - &PRINT - &RESTART - FILENAME lead1 - &END RESTART - &END PRINT - &END SCF - - &POISSON - PERIODIC xyz - &END POISSON - - &XC - &XC_FUNCTIONAL pbe - &END XC_FUNCTIONAL - &END XC - - &KPOINTS - SCHEME MONKHORST-PACK 4 4 4 - &END KPOINTS - &END DFT - - &SUBSYS - &CELL - ABC 8.00000 8.00000 3.70000 - PERIODIC xyz - &END CELL - - &COORD -H 0.00000000 0.00000000 11.10000000 L1 -H 0.00000000 0.00000000 11.84000000 L1 -H 0.00000000 0.00000000 12.58000000 L1 -H 0.00000000 0.00000000 13.32000000 L1 -H 0.00000000 0.00000000 14.06000000 L1 - &END COORD - - &KIND H - BASIS_SET SZV-MOLOPT-GTH - POTENTIAL GTH-PBE-q1 - &END KIND - &END SUBSYS -&END FORCE_EVAL - -&FORCE_EVAL - METHOD Quickstep - - &DFT - BASIS_SET_FILE_NAME BASIS_MOLOPT - POTENTIAL_FILE_NAME GTH_POTENTIALS - - &MGRID - CUTOFF 300 - &END MGRID - - &QS - METHOD gpw - EPS_DEFAULT 1e-10 - &END QS - - &SCF - MAX_SCF 10 - EPS_SCF 1e-6 - SCF_GUESS atomic - ADDED_MOS 100 - - &SMEAR - ELECTRONIC_TEMPERATURE 293 - METHOD FERMI_DIRAC - &END SMEAR - - &MIXING - METHOD BROYDEN_MIXING - ALPHA 0.3 - NBUFFER 8 - &END MIXING - - &PRINT - &RESTART - FILENAME lead2 - &END RESTART - &END PRINT - &END SCF - - &POISSON - PERIODIC xyz - &END POISSON - - &XC - &XC_FUNCTIONAL pbe - &END XC_FUNCTIONAL - &END XC - - &KPOINTS - SCHEME MONKHORST-PACK 4 4 4 - &END KPOINTS - &END DFT - - &SUBSYS - &CELL - ABC 8.00000 8.00000 3.70000 - PERIODIC xyz - &END CELL - - &COORD -H 0.00000000 0.00000000 25.90000000 R1 -H 0.00000000 0.00000000 26.64000000 R1 -H 0.00000000 0.00000000 27.38000000 R1 -H 0.00000000 0.00000000 28.12000000 R1 -H 0.00000000 0.00000000 28.86000000 R1 - &END COORD - - &KIND H - BASIS_SET SZV-MOLOPT-GTH - POTENTIAL GTH-PBE-q1 - &END KIND - &END SUBSYS -&END FORCE_EVAL diff --git a/tests/QS/regtest-negf/li_chain.inp b/tests/QS/regtest-negf/li_chain.inp new file mode 100644 index 0000000000..b5eb2eddb2 --- /dev/null +++ b/tests/QS/regtest-negf/li_chain.inp @@ -0,0 +1,376 @@ +@SET project li_chain + +@SET v_bias_volt 0.0 +@SET temperature 293.0 + +# The following two parameters are computed automatically during the NEGF run. +# You can uncomment keywords 'FERMI_LEVEL' and 'V_SHIFT' to skip these stages. +@SET fermi_level -0.08987575 +@SET v_shift 0.00280782 + +@SET basis_set_path BASIS_MOLOPT +@SET potential_path GTH_POTENTIALS +@SET basis_set_li SZV-MOLOPT-SR-GTH +@SET potential_li GTH-PADE-q3 + +@SET xc_functional pade + +@SET cutoff 300 +@SET added_mos 100 +@SET max_scf 50 +@SET eps_scf 1e-7 + +@SET lx 8.0 +@SET ly 8.0 + +&GLOBAL + PROJECT ${project} + RUN_TYPE negf + PRINT_LEVEL low + EXTENDED_FFT_LENGTHS .true. +&END GLOBAL + +&MULTIPLE_FORCE_EVALS + # Order in which the input file parser internally stores different force environments. + # + # 'FORCE_EVAL_ORDER 1 2 3' means that the 1st force environment -- + # which holds the entire system -- will have the internal index '0', + # while the 2nd and 3rd force environments -- which correspond to the + # 1st and 2nd semi-infinite electrodes respectively -- will have + # the internal indices '1' and '2'. Note that internal indices start from 0. + # + FORCE_EVAL_ORDER 1 2 3 + + # Different force environments contain non-identical molecular structures. + # By default, the input parser ignores 'subsys' sections from all but the force + # environment with the internal index '0'. + MULTIPLE_SUBSYS .true. +&END + +&NEGF + &CONTACT + # the force environment with the internal index '1' should be used + # to compute 'bulk' properties of this (1st) semi-infinite electrode. + FORCE_EVAL_SECTION 1 + + &BULK_REGION + # This section is used to map matrix elements between the primary + # (with index 0) and the 'bulk' (with index 1) force environments. + # + # Atomic indices (and/or names of molecular fragments) are relative + # to the primary force environment. + MOLNAME L1 L2 + &END + &SCREENING_REGION + MOLNAME L0 + &END + + # Uncomment the next line to skip refining the Fermi level for this contact + #FERMI_LEVEL ${fermi_level} + + ELECTRIC_POTENTIAL [eV] ${v_bias_volt}/2 + TEMPERATURE [K] ${temperature} + &END CONTACT + + &CONTACT + # the second semi-infinite electrode + FORCE_EVAL_SECTION 2 + + &BULK_REGION + MOLNAME R1 R2 + &END + &SCREENING_REGION + MOLNAME R0 + &END + + #FERMI_LEVEL ${fermi_level} + + ELECTRIC_POTENTIAL [eV] -${v_bias_volt}/2 + TEMPERATURE [K] ${temperature} + &END CONTACT + + &SCATTERING_REGION + MOLNAME S + &END SCATTERING_REGION + + # Density mixing is not needed in zero-bias case, as NEGF self-consistent + # procedure will converge in one iteration. However, when an external bias + # is applied, a density mixing scheme in inverse space is typically + # much stable then the default direct mixing of new and old density matrices. + &MIXING + ALPHA 0.2 + METHOD broyden_mixing + NBUFFER 8 + &END MIXING + + # Integration lower bound (in Hartree) + ENERGY_LBOUND -2.0 + + # Maximum error in each matrix element of the density matrix + EPS_DENSITY 1e-5 + + MAX_SCF 16 + + # Boundary conditions define solution of the Poisson equation up to + # an additive constant, which affects Hamiltonian matrix elements within + # the NEGF procedure. The following parameters control evaluation of + # the correct shift in Hartree potential: + # * V_SHIFT -- initial guess for the shift + # * V_SHIFT_MAX_ITERS -- maximal number of iterations + #V_SHIFT ${v_shift} + V_SHIFT_MAX_ITERS 20 + + &PRINT + &DOS + FILENAME device + FROM_ENERGY -0.2 + TILL_ENERGY 0.1 + N_GRIDPOINTS 301 + &END + &TRANSMISSION + FILENAME transm + FROM_ENERGY -0.2 + TILL_ENERGY 0.1 + N_GRIDPOINTS 301 + &END + &END + + # Number of parallel processes per energy point. + # 'NPROC_POINT 1' means that if you run this job using 8 MPI processes, + # than Green's functions at up to 8 (8 / NPROC_POINT) energy points + # will be computed simultaneously. + # + NPROC_POINT 1 +&END NEGF + +# *** Primary force environment *** +&FORCE_EVAL + METHOD Quickstep + + &DFT + BASIS_SET_FILE_NAME ${basis_set_path} + POTENTIAL_FILE_NAME ${potential_path} + + &MGRID + CUTOFF ${cutoff} + &END MGRID + + &SCF + MAX_SCF ${max_scf} + EPS_SCF ${eps_scf} + SCF_GUESS atomic + ADDED_MOS ${added_mos} + + &SMEAR + ELECTRONIC_TEMPERATURE [K] ${temperature} + METHOD fermi_dirac + &END SMEAR + + &MIXING + ALPHA 0.3 + METHOD broyden_mixing + NBUFFER 8 + &END MIXING + + &PRINT + &RESTART + FILENAME RESTART_0 + &END RESTART + &END PRINT + &END SCF + + &POISSON + PERIODIC xyz + &END POISSON + + &XC + &XC_FUNCTIONAL ${xc_functional} + &END XC_FUNCTIONAL + &END XC + &END DFT + + &SUBSYS + &CELL + ABC ${lx} ${ly} 83.720 + PERIODIC xyz + &END CELL + + &COORD + Li 0.00000000 0.00000000 0.00000000 L2 + Li 0.00000000 0.00000000 2.99000000 L2 + Li 0.00000000 0.00000000 5.98000000 L2 + Li 0.00000000 0.00000000 8.97000000 L2 + Li 0.00000000 0.00000000 11.96000000 L1 + Li 0.00000000 0.00000000 14.95000000 L1 + Li 0.00000000 0.00000000 17.94000000 L1 + Li 0.00000000 0.00000000 20.93000000 L1 + Li 0.00000000 0.00000000 23.92000000 L0 + Li 0.00000000 0.00000000 26.91000000 L0 + Li 0.00000000 0.00000000 29.90000000 L0 + Li 0.00000000 0.00000000 32.89000000 L0 + Li 0.00000000 0.00000000 35.88000000 S + Li 0.00000000 0.00000000 38.87000000 S + Li 0.00000000 0.00000000 41.86000000 S + Li 0.00000000 0.00000000 44.85000000 S + Li 0.00000000 0.00000000 47.84000000 R0 + Li 0.00000000 0.00000000 50.83000000 R0 + Li 0.00000000 0.00000000 53.82000000 R0 + Li 0.00000000 0.00000000 56.81000000 R0 + Li 0.00000000 0.00000000 59.80000000 R1 + Li 0.00000000 0.00000000 62.79000000 R1 + Li 0.00000000 0.00000000 65.78000000 R1 + Li 0.00000000 0.00000000 68.77000000 R1 + Li 0.00000000 0.00000000 71.76000000 R2 + Li 0.00000000 0.00000000 74.75000000 R2 + Li 0.00000000 0.00000000 77.74000000 R2 + Li 0.00000000 0.00000000 80.73000000 R2 + &END COORD + + &KIND Li + BASIS_SET ${basis_set_li} + POTENTIAL ${potential_li} + &END KIND + &END SUBSYS +&END FORCE_EVAL + +# *** The 1st (left) semi-infinite electrode *** +&FORCE_EVAL + METHOD Quickstep + + &DFT + BASIS_SET_FILE_NAME ${basis_set_path} + POTENTIAL_FILE_NAME ${potential_path} + + &MGRID + CUTOFF ${cutoff} + &END MGRID + + &SCF + MAX_SCF ${max_scf} + EPS_SCF ${eps_scf} + SCF_GUESS atomic + ADDED_MOS ${added_mos} + + &SMEAR + ELECTRONIC_TEMPERATURE [K] ${temperature} + METHOD fermi_dirac + &END SMEAR + + &MIXING + ALPHA 0.3 + METHOD broyden_mixing + NBUFFER 8 + &END MIXING + + &PRINT + &RESTART + FILENAME RESTART_1 + &END RESTART + &END PRINT + &END SCF + + &POISSON + PERIODIC xyz + &END POISSON + + &XC + &XC_FUNCTIONAL ${xc_functional} + &END XC_FUNCTIONAL + &END XC + + &KPOINTS + SCHEME MONKHORST-PACK 4 4 8 + &END KPOINTS + &END DFT + + &SUBSYS + # a primary unit cell of the 1st semi-infinite electrode + &CELL + ABC ${lx} ${ly} 11.960 + PERIODIC xyz + &END CELL + + &COORD + Li 0.00000000 0.00000000 11.96000000 L1 + Li 0.00000000 0.00000000 14.95000000 L1 + Li 0.00000000 0.00000000 17.94000000 L1 + Li 0.00000000 0.00000000 20.93000000 L1 + &END COORD + + &KIND Li + BASIS_SET ${basis_set_li} + POTENTIAL ${potential_li} + &END KIND + &END SUBSYS +&END FORCE_EVAL + +# *** The 2nd (right) semi-infinite electrode *** +&FORCE_EVAL + METHOD Quickstep + + &DFT + BASIS_SET_FILE_NAME ${basis_set_path} + POTENTIAL_FILE_NAME ${potential_path} + + &MGRID + CUTOFF ${cutoff} + &END MGRID + + &SCF + MAX_SCF ${max_scf} + EPS_SCF ${eps_scf} + SCF_GUESS atomic + ADDED_MOS ${added_mos} + + &SMEAR + ELECTRONIC_TEMPERATURE [K] ${temperature} + METHOD fermi_dirac + &END SMEAR + + &MIXING + ALPHA 0.3 + METHOD broyden_mixing + NBUFFER 8 + &END MIXING + + &PRINT + &RESTART + FILENAME RESTART_2 + &END RESTART + &END PRINT + &END SCF + + &POISSON + PERIODIC xyz + &END POISSON + + &XC + &XC_FUNCTIONAL ${xc_functional} + &END XC_FUNCTIONAL + &END XC + + &KPOINTS + SCHEME MONKHORST-PACK 4 4 8 + &END KPOINTS + &END DFT + + &SUBSYS + # a primary unit cell of the 2nd semi-infinite electrode + &CELL + ABC ${lx} ${ly} 11.960 + PERIODIC xyz + &END CELL + + &COORD + Li 0.00000000 0.00000000 59.80000000 R1 + Li 0.00000000 0.00000000 62.79000000 R1 + Li 0.00000000 0.00000000 65.78000000 R1 + Li 0.00000000 0.00000000 68.77000000 R1 + &END COORD + + &KIND Li + BASIS_SET ${basis_set_li} + POTENTIAL ${potential_li} + &END KIND + &END SUBSYS +&END FORCE_EVAL