From 48db7ef82b9abb09d4866a50036cdab8ac993d68 Mon Sep 17 00:00:00 2001 From: Dmitry Ryndyk Date: Mon, 20 Apr 2026 13:35:38 +0200 Subject: [PATCH] NEGF: xTB added --- src/negf_env_types.F | 2 + src/negf_methods.F | 32 ++- src/qs_ks_methods.F | 2 +- src/qs_scf_post_tb.F | 2 +- tests/QS/regtest-negf-xtb-1/H28_k.inp | 285 ++++++++++++++++++++ tests/QS/regtest-negf-xtb-1/TEST_FILES.toml | 7 + tests/TEST_DIRS | 1 + 7 files changed, 322 insertions(+), 9 deletions(-) create mode 100644 tests/QS/regtest-negf-xtb-1/H28_k.inp create mode 100644 tests/QS/regtest-negf-xtb-1/TEST_FILES.toml diff --git a/src/negf_env_types.F b/src/negf_env_types.F index b9fb915dee..7ddea351b0 100644 --- a/src/negf_env_types.F +++ b/src/negf_env_types.F @@ -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 diff --git a/src/negf_methods.F b/src/negf_methods.F index 8d8c55a554..0328a25233 100644 --- a/src/negf_methods.F +++ b/src/negf_methods.F @@ -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, & diff --git a/src/qs_ks_methods.F b/src/qs_ks_methods.F index f97872e617..88ce5c86d8 100644 --- a/src/qs_ks_methods.F +++ b/src/qs_ks_methods.F @@ -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 diff --git a/src/qs_scf_post_tb.F b/src/qs_scf_post_tb.F index 509e85bc24..548dd177e6 100644 --- a/src/qs_scf_post_tb.F +++ b/src/qs_scf_post_tb.F @@ -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 ! ************************************************************************************************** diff --git a/tests/QS/regtest-negf-xtb-1/H28_k.inp b/tests/QS/regtest-negf-xtb-1/H28_k.inp new file mode 100644 index 0000000000..056beee564 --- /dev/null +++ b/tests/QS/regtest-negf-xtb-1/H28_k.inp @@ -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 diff --git a/tests/QS/regtest-negf-xtb-1/TEST_FILES.toml b/tests/QS/regtest-negf-xtb-1/TEST_FILES.toml new file mode 100644 index 0000000000..8910e5e613 --- /dev/null +++ b/tests/QS/regtest-negf-xtb-1/TEST_FILES.toml @@ -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}] diff --git a/tests/TEST_DIRS b/tests/TEST_DIRS index 1b21d81c4b..1b0b7a36df 100644 --- a/tests/TEST_DIRS +++ b/tests/TEST_DIRS @@ -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