NEGF: xTB added

This commit is contained in:
Dmitry Ryndyk 2026-04-20 13:35:38 +02:00 committed by Matthias Krack
parent ff7b8ebec0
commit 48db7ef82b
7 changed files with 322 additions and 9 deletions

View file

@ -77,6 +77,7 @@ MODULE negf_env_types
mo_set_type
USE qs_rho_types, ONLY: qs_rho_get,&
qs_rho_type
USE qs_scf_post_tb, ONLY: rebuild_pw_env
USE qs_subsys_types, ONLY: qs_subsys_get,&
qs_subsys_type
#include "./base/base_uses.f90"
@ -1313,6 +1314,7 @@ CONTAINS
para_env=para_env, &
pw_env=pw_env, &
subsys=subsys)
IF (dft_control%qs_control%xtb) CALL rebuild_pw_env(qs_env)
CALL pw_env_get(pw_env, auxbas_pw_pool=pw_pool)
IF (do_kpoints) THEN

View file

@ -109,15 +109,19 @@ MODULE negf_methods
evolt,&
kelvin,&
seconds
USE qs_density_mixing_types, ONLY: direct_mixing_nr,&
gspace_mixing_nr
USE qs_density_mixing_types, ONLY: broyden_mixing_nr,&
direct_mixing_nr,&
gspace_mixing_nr,&
multisecant_mixing_nr,&
pulay_mixing_nr
USE qs_energy, ONLY: qs_energies
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gspace_mixing, ONLY: gspace_mixing
USE qs_ks_methods, ONLY: qs_ks_build_kohn_sham_matrix
USE qs_mixing_utils, ONLY: mixing_allocate,&
USE qs_ks_methods, ONLY: rebuild_ks_matrix
USE qs_mixing_utils, ONLY: charge_mixing_init,&
mixing_allocate,&
mixing_init
USE qs_rho_methods, ONLY: qs_rho_update_rho
USE qs_rho_types, ONLY: qs_rho_get,&
@ -1127,6 +1131,18 @@ CONTAINS
IF (log_unit > 0) THEN
WRITE (log_unit, '(/,T2,A)') "NEGF SELF-CONSISTENT PROCEDURE"
WRITE (log_unit, "( ' ------------------------------')")
IF (negf_env%mixing_method == direct_mixing_nr) THEN
WRITE (log_unit, '(T3,A)') "Mixing method: Direct mixing of new and old density matrices"
END IF
IF (negf_env%mixing_method == broyden_mixing_nr) THEN
WRITE (log_unit, '(T3,A)') "Mixing method: Broyden mixing"
END IF
IF (negf_env%mixing_method == pulay_mixing_nr) THEN
WRITE (log_unit, '(T3,A)') "Mixing method: Pulay mixing"
END IF
IF (negf_env%mixing_method == multisecant_mixing_nr) THEN
WRITE (log_unit, '(T3,A)') "Mixing method: Multisecant scheme for mixing"
END IF
END IF
IF (negf_control%update_HS .AND. (.NOT. negf_control%is_dft_entire)) &
@ -1312,8 +1328,10 @@ CONTAINS
! 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)
IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
CPABORT('TB Code not available')
IF (dft_control%qs_control%dftb) THEN
CPABORT('DFTB Code not available')
ELSE IF (dft_control%qs_control%xtb) THEN
CALL charge_mixing_init(negf_env%mixing_storage)
ELSE IF (dft_control%qs_control%semi_empirical) THEN
CPABORT('SE Code not possible')
ELSE
@ -1504,7 +1522,7 @@ CONTAINS
! update KS-matrix
IF (negf_control%update_HS) THEN
CALL qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
CALL rebuild_ks_matrix(qs_env, calculate_forces=.FALSE., just_energy=.FALSE.)
! extract blocks from the updated Kohn-Sham matrix
DO ispin = 1, nspins
CALL negf_copy_sym_dbcsr_to_fm_submat(matrix=matrix_ks_qs_kp(ispin, 1)%matrix, &

View file

@ -145,7 +145,7 @@ MODULE qs_ks_methods
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods'
PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
qs_ks_allocate_basics, evaluate_core_matrix_traces
qs_ks_allocate_basics, evaluate_core_matrix_traces, rebuild_ks_matrix
CONTAINS

View file

@ -149,7 +149,7 @@ MODULE qs_scf_post_tb
! Global parameters
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_scf_post_tb'
PUBLIC :: scf_post_calculation_tb, make_lumo_tb
PUBLIC :: scf_post_calculation_tb, make_lumo_tb, rebuild_pw_env
! **************************************************************************************************

View file

@ -0,0 +1,285 @@
@SET Project H28_k
@SET OT F
@SET LCD F
@SET UPDATE_HS T
@SET RESTART_SCF F
@SET V_bias_volt 0.0 #1.0
@SET Temperature 100.
# The following two parameters are computed automatically during the NEGF run.
# You can use keywords 'FERMI_LEVEL' to give this value and skip calculation
# or 'FERMI_LEVEL'+'REFINE_FERMI_LEVEL' to compute starting from this value.
# If 'FERMI_LEVEL' is not given, it is calculated automatically.
# 'V_SHIFT' is used to give initial value.
@SET Fermi_level -0.38560885
@SET V_shift 0.00005408
@SET cutoff 300
@SET added_mos 10
@SET max_scf 1000
@SET eps_scf 1e-8
@SET lx 6.0
@SET ly 6.0
!--------------------------------------------------------------------------------------------------!
&GLOBAL
EXTENDED_FFT_LENGTHS .true.
PREFERRED_DIAG_LIBRARY SL
PRINT_LEVEL LOW #SILENT #LOW #MEDIUM #HIGH #DEBUG
PROJECT ${Project}
RUN_TYPE NEGF
&END GLOBAL
&NEGF
# Integration lower bound (in Hartree)
ENERGY_LBOUND -2.0
!clenshaw-curtis simpson
! INTEGRATION_MIN_POINTS 64
! INTEGRATION_MAX_POINTS 768
# Maximum error in each matrix element of the density matrix
EPS_DENSITY 1e-7
EPS_GREEN 1e-7
EPS_SCF 1e-8
ETA 1e-7
HOMO_LUMO_GAP 0.5
# 'simpson' is slower, but works without FFTW3 library
# clenshaw-curtis is faster, but requires FFTW3 library
INTEGRATION_METHOD simpson
MAX_SCF 0
# 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 2
# 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 100
&CONTACT
ELECTRIC_POTENTIAL [eV] ${V_bias_volt}/2
FERMI_LEVEL ${Fermi_level}
!FERMI_LEVEL_SHIFTED -0.12306442
!REFINE_FERMI_LEVEL
TEMPERATURE [K] ${Temperature}
# 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 'electrode' (with index 1) force environments.
#
# Atomic indices (and/or names of molecular fragments) are relative
# to the primary force environment.
MOLNAME L1 L2
! LIST 1..8
&CELL
MOLNAME L1
! LIST 5..8
&END CELL
&CELL
MOLNAME L2
! LIST 1..4
&END CELL
&END BULK_REGION
&RESTART
! FILENAME Test
!READ_WRITE_HS
! READ_WRITE_SGF
&END RESTART
&SCREENING_REGION
MOLNAME L0
! LIST 9..12
&END SCREENING_REGION
!&PRINT
! &DOS
! FILENAME left
! FROM_ENERGY [eV] -15
! TILL_ENERGY [eV] 15
! N_GRIDPOINTS 401
! &END
!&END
&END CONTACT
&CONTACT
ELECTRIC_POTENTIAL [eV] -${V_bias_volt}/2
FERMI_LEVEL ${Fermi_level}
!REFINE_FERMI_LEVEL
TEMPERATURE [K] ${Temperature}
# the second semi-infinite electrode
! FORCE_EVAL_SECTION 2
&BULK_REGION
MOLNAME R1 R2
&CELL
MOLNAME R1
&END CELL
&CELL
MOLNAME R2
&END CELL
&END BULK_REGION
&RESTART
! FILENAME Test
!READ_WRITE_HS
! READ_WRITE_SGF
&END RESTART
&SCREENING_REGION
MOLNAME R0
&END SCREENING_REGION
&END CONTACT
# 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
! METHOD pulay_mixing
!direct_p_mixing broyden_mixing pulay_mixing
ALPHA 0.05
NBUFFER 8
PULAY_ALPHA 0.2
! PULAY_BETA 0.2
&END MIXING
&PRINT
&DOS
FILENAME device
FROM_ENERGY [eV] -15
N_GRIDPOINTS 401
TILL_ENERGY [eV] 15
&END DOS
&PROGRAM_RUN_INFO low
PRINT_LEVEL high
!high debug
&END PROGRAM_RUN_INFO
&RESTART OFF
BACKUP_COPIES 1
!FILENAME Test
&EACH
NEGF_SCF 10
&END EACH
&END RESTART
&TRANSMISSION
FILENAME transm
FROM_ENERGY [eV] -15
N_GRIDPOINTS 401
TILL_ENERGY [eV] 15
&END TRANSMISSION
&END PRINT
&SCATTERING_REGION
MOLNAME S
!LIST 13..16
&RESTART
! FILENAME Test
!READ_WRITE_HS
! READ_WRITE_SGF
&END RESTART
&END SCATTERING_REGION
&SCF
RESTART_SCF $RESTART_SCF
UPDATE_HS $UPDATE_HS
&END SCF
! DISABLE_CACHE
&END NEGF
# *** Primary force environment ***
&FORCE_EVAL
METHOD Quickstep
&DFT
LSD $LCD
MULTIPLICITY 1
&MGRID
CUTOFF ${cutoff}
&END MGRID
&POISSON
PERIODIC xyz
&EWALD
ALPHA 1.0
&END EWALD
&END POISSON
&QS
METHOD XTB
&XTB
CHECK_ATOMIC_CHARGES T
DO_EWALD T
GFN_TYPE 1
! &TBLITE
! METHOD GFN2
! &END
&END XTB
&END QS
&SCF
MAX_SCF ${max_scf}
EPS_SCF ${eps_scf}
SCF_GUESS restart
!atomic restart
@IF ($OT == T)
&OT
MINIMIZER DIIS
N_DIIS 12
PRECONDITIONER FULL_SINGLE_INVERSE
! PRECONDITIONER FULL_ALL
! PRECONDITIONER FULL_KINETIC
! PRECOND_SOLVER INVERSE_CHOLESKY
! ALGORITHM IRAC
! ENERGY_GAP 0.001
! LINESEARCH 2PNT
! STEPSIZE 0.1
&END OT
@ENDIF
@IF ($OT == F)
ADDED_MOS ${added_mos}
&SMEAR
ELECTRONIC_TEMPERATURE [K] 300
METHOD fermi_dirac
&END SMEAR
&MIXING
ALPHA 0.2
METHOD broyden_mixing
!broyden_mixing DIRECT_P_MIXING
NBUFFER 8
&END MIXING
&DIAGONALIZATION
&END DIAGONALIZATION
@ENDIF
&PRINT
&RESTART
FILENAME RESTART_0
&END RESTART
&END PRINT
&END SCF
&END DFT
&SUBSYS
&CELL
ABC ${lx} ${ly} 42.
PERIODIC xyz
&END CELL
&COORD
H 0.00000000 0.00000000 0.00000000 L2
H 0.00000000 0.00000000 1.50000000 L2
H 0.00000000 0.00000000 3.00000000 L2
H 0.00000000 0.00000000 4.50000000 L2
H 0.00000000 0.00000000 6.00000000 L1
H 0.00000000 0.00000000 7.50000000 L1
H 0.00000000 0.00000000 9.00000000 L1
H 0.00000000 0.00000000 10.50000000 L1
H 0.00000000 0.00000000 12.00000000 L0
H 0.00000000 0.00000000 13.50000000 L0
H 0.00000000 0.00000000 15.00000000 L0
H 0.00000000 0.00000000 16.50000000 L0
H 0.00000000 0.00000000 18.00000000 S
H 0.00000000 0.00000000 19.50000000 S
H 0.00000000 0.00000000 21.00000000 S
H 0.00000000 0.00000000 22.50000000 S
H 0.00000000 0.00000000 24.00000000 R0
H 0.00000000 0.00000000 25.50000000 R0
H 0.00000000 0.00000000 27.00000000 R0
H 0.00000000 0.00000000 28.50000000 R0
H 0.00000000 0.00000000 30.00000000 R1
H 0.00000000 0.00000000 31.50000000 R1
H 0.00000000 0.00000000 33.00000000 R1
H 0.00000000 0.00000000 34.50000000 R1
H 0.00000000 0.00000000 36.00000000 R2
H 0.00000000 0.00000000 37.50000000 R2
H 0.00000000 0.00000000 39.00000000 R2
H 0.00000000 0.00000000 40.50000000 R2
&END COORD
&END SUBSYS
&END FORCE_EVAL

View file

@ -0,0 +1,7 @@
# runs are executed in the same order as in this file
# the second field tells which test should be run in order to compare with the last available output
# e.g. 0 means do not compare anything, running is enough
# 1 compares the last total energy in the file
# for details see cp2k/tools/do_regtest
#
"H28_k.inp" = [{matcher="E_total", tol=1e-9, ref=-13.31730669072749}]

View file

@ -370,6 +370,7 @@ FE/regtest-3
QS/regtest-properties/resp
QS/regtest-negf-1 fftw3
QS/regtest-negf-2 fftw3
QS/regtest-negf-xtb-1 fftw3
NEB/regtest-1
QS/regtest-double-hybrid-stress-numer-meta libxc !ifx
QS/regtest-elpa-qr elpa mpiranks==1||mpiranks%2==0