comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

NEGF: contact-specific temperature, correct shift and scale factors

Browse source 
svn-origin-rev: 18419
This commit is contained in:
Sergey Chulkov 2018-05-10 12:11:54 +00:00
parent e2eae28d33
commit a80b2e6639
12 changed files with 1105 additions and 907 deletions
Download patch file Download diff file Expand all files Collapse all files

45
src/input_cp2k_negf.F
View file

@ -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.)

78
src/negf_control_types.F
View file

@ -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

166
src/negf_env_types.F
View file

@ -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)

9
src/negf_green_methods.F
View file

@ -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)

584
src/negf_methods.F
View file

File diff suppressed because it is too large Load diff

190
src/negf_qs_env_methods.F
View file

@ -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

151
src/negf_qs_env_types.F
View file

@ -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

31
src/negf_vectors.F
View file

@ -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)

23
src/qs_environment_types.F
View file

@ -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)

18
src/qs_ks_methods.F
View file

@ -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+ &

341
tests/QS/regtest-negf/hchain.inp
View file

@ -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

376
tests/QS/regtest-negf/li_chain.inp Normal file
View file

@ -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