From 1b84f40cd1e5fca921fa993ecbba26c908d950af Mon Sep 17 00:00:00 2001 From: Juerg Hutter Date: Wed, 15 Jul 2026 17:30:31 +0200 Subject: [PATCH] xTB spin polarisation electronic hessian --- src/qs_linres_kernel.F | 5 +- src/qs_tddfpt2_forces.F | 3 +- src/response_solver.F | 2 +- src/xtb_ehess.F | 9 +- src/xtb_spinpol.F | 116 +++++++++++++++++++++- tests/xTB/regtest-spinpol/TEST_FILES.toml | 2 + tests/xTB/regtest-spinpol/ch2o_polar.inp | 63 ++++++++++++ 7 files changed, 194 insertions(+), 6 deletions(-) create mode 100644 tests/xTB/regtest-spinpol/ch2o_polar.inp diff --git a/src/qs_linres_kernel.F b/src/qs_linres_kernel.F index 72b5147d11..82abe2463a 100644 --- a/src/qs_linres_kernel.F +++ b/src/qs_linres_kernel.F @@ -572,7 +572,7 @@ CONTAINS REAL(dp), ALLOCATABLE, DIMENSION(:) :: mcharge, mcharge1 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, charges, charges1 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set - TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmat, rho_ao TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, matrix_p1, matrix_s TYPE(dft_control_type), POINTER :: dft_control TYPE(linres_control_type), POINTER :: linres_control @@ -651,7 +651,8 @@ CONTAINS mcharge1(iatom) = SUM(charges1(iatom, :)) END DO ! Coulomb Kernel - CALL xtb_coulomb_hessian(qs_env, p_env%kpp1, charges1, mcharge1, mcharge) + pmat => matrix_p1(:, 1) + CALL xtb_coulomb_hessian(qs_env, p_env%kpp1, charges1, mcharge1, mcharge, pmat) ! DEALLOCATE (charges, mcharge, charges1, mcharge1) END IF diff --git a/src/qs_tddfpt2_forces.F b/src/qs_tddfpt2_forces.F index 9e33c94673..537e329616 100644 --- a/src/qs_tddfpt2_forces.F +++ b/src/qs_tddfpt2_forces.F @@ -1436,7 +1436,8 @@ CONTAINS mcharge1(iatom) = SUM(charges1(iatom, :)) END DO ! Coulomb Kernel - CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge) + CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge, & + matrix_pe) ! DEALLOCATE (charges, mcharge, charges1, mcharge1) END IF diff --git a/src/response_solver.F b/src/response_solver.F index 6bc3b01205..b2e01757bc 100644 --- a/src/response_solver.F +++ b/src/response_solver.F @@ -2772,7 +2772,7 @@ CONTAINS mcharge1(iatom) = SUM(charges1(iatom, :)) END DO ! Coulomb Kernel - CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge) + CALL xtb_coulomb_hessian(qs_env, matrix_hz, charges1, mcharge1, mcharge, mpa) CALL calc_xtb_ehess_force(qs_env, p_matrix, mpa, charges, mcharge, charges1, & mcharge1, debug_forces) ! diff --git a/src/xtb_ehess.F b/src/xtb_ehess.F index a7c1f17cdf..8eb7a0a5af 100644 --- a/src/xtb_ehess.F +++ b/src/xtb_ehess.F @@ -51,6 +51,7 @@ MODULE xtb_ehess neighbor_list_set_p_type USE virial_types, ONLY: virial_type USE xtb_coulomb, ONLY: gamma_rab_sr + USE xtb_spinpol, ONLY: xtb_spinpol_hessian USE xtb_types, ONLY: get_xtb_atom_param,& xtb_atom_type #include "./base/base_uses.f90" @@ -72,13 +73,15 @@ CONTAINS !> \param charges1 ... !> \param mcharge1 ... !> \param mcharge ... +!> \param matrix_p1 ... ! ************************************************************************************************** - SUBROUTINE xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge) + SUBROUTINE xtb_coulomb_hessian(qs_env, ks_matrix, charges1, mcharge1, mcharge, matrix_p1) TYPE(qs_environment_type), POINTER :: qs_env TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix REAL(dp), DIMENSION(:, :) :: charges1 REAL(dp), DIMENSION(:) :: mcharge1, mcharge + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p1 CHARACTER(len=*), PARAMETER :: routineN = 'xtb_coulomb_hessian' @@ -284,6 +287,10 @@ CONTAINS DEALLOCATE (xgamma) END IF + IF (xtb_control%do_spinpol) THEN + CALL xtb_spinpol_hessian(qs_env, ks_matrix, matrix_p1) + END IF + IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN CPABORT("QMMM not available in xTB response calculations") END IF diff --git a/src/xtb_spinpol.F b/src/xtb_spinpol.F index da95cf1c7c..efdf0c78c1 100644 --- a/src/xtb_spinpol.F +++ b/src/xtb_spinpol.F @@ -58,7 +58,7 @@ MODULE xtb_spinpol CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_spinpol' - PUBLIC :: build_xtb_spinpol + PUBLIC :: build_xtb_spinpol, xtb_spinpol_hessian CONTAINS @@ -169,6 +169,7 @@ CONTAINS ! Calculate charges ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom)) aocg = 0.0_dp + bocg = 0.0_dp IF (nimg > 1) THEN matrix_s_kp => matrix_s(:, :) matrix_p_kp => matrix_p(1:1, :) @@ -449,6 +450,119 @@ CONTAINS END SUBROUTINE build_xtb_spinpol +! ************************************************************************************************** +!> \brief ... +!> \param qs_env ... +!> \param ks_matrix ... +!> \param matrix_p1 ... +! ************************************************************************************************** + SUBROUTINE xtb_spinpol_hessian(qs_env, ks_matrix, matrix_p1) + + TYPE(qs_environment_type), POINTER :: qs_env + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p1 + + CHARACTER(len=*), PARAMETER :: routineN = 'xtb_spinpol_hessian' + + INTEGER :: handle, ia, ib, icol, ikind, irow, & + jkind, la, lb, na, natom, natorb, nb, & + nimg, nkind, nsgf, nspins + INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of + INTEGER, DIMENSION(25) :: lao + LOGICAL :: found + REAL(KIND=dp) :: fval + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg1, bocg1 + REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk + REAL(KIND=dp), DIMENSION(3, 3) :: wall + REAL(KIND=dp), DIMENSION(:, :), POINTER :: aksb, bksb, sblock + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(dbcsr_iterator_type) :: iter + TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix + TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s + TYPE(dbcsr_type), POINTER :: s_matrix + TYPE(dft_control_type), POINTER :: dft_control + TYPE(kpoint_type), POINTER :: kpoints + TYPE(mp_para_env_type), POINTER :: para_env + TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set + TYPE(xtb_atom_type), POINTER :: xtb_kind + + CALL timeset(routineN, handle) + + CALL get_qs_env(qs_env, dft_control=dft_control) + nspins = dft_control%nspins + nimg = dft_control%nimages + + IF (nimg /= 1) THEN + CPABORT("No kpoints allowed in xTB response calculation") + END IF + + IF (nspins == 2) THEN + + CALL get_qs_env(qs_env, & + qs_kind_set=qs_kind_set, & + atomic_kind_set=atomic_kind_set) + + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of) + + CALL get_qs_env(qs_env, nkind=nkind, natom=natom) + CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf) + CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env) + + ! expand parameters + ALLOCATE (wabk(nsgf, nsgf, nkind)) + wabk = 0.0_dp + DO ikind = 1, nkind + CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind) + CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall) + DO ia = 1, natorb + la = lao(ia) + 1 + DO ib = 1, natorb + lb = lao(ib) + 1 + wabk(ia, ib, ikind) = wall(la, lb) + END DO + END DO + END DO + + ! Calculate response charges + ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom)) + aocg1 = 0.0_dp + bocg1 = 0.0_dp + s_matrix => matrix_s(1, 1)%matrix + p_matrix => matrix_p1(1:1) + CALL ao_charges(p_matrix, s_matrix, aocg1, para_env) + p_matrix => matrix_p1(2:2) + CALL ao_charges(p_matrix, s_matrix, bocg1, para_env) + aocg1 = 0.5_dp*aocg1 + bocg1 = 0.5_dp*bocg1 + + CALL dbcsr_iterator_start(iter, s_matrix) + DO WHILE (dbcsr_iterator_blocks_left(iter)) + CALL dbcsr_iterator_next_block(iter, irow, icol, sblock) + CALL dbcsr_get_block_p(matrix=ks_matrix(1)%matrix, & + row=irow, col=icol, block=aksb, found=found) + CPASSERT(found) + CALL dbcsr_get_block_p(matrix=ks_matrix(2)%matrix, & + row=irow, col=icol, block=bksb, found=found) + CPASSERT(found) + na = SIZE(aksb, 1) + nb = SIZE(aksb, 2) + ikind = kind_of(irow) + jkind = kind_of(icol) + fval = 1.00_dp + CALL ksupdate(aksb, bksb, sblock, na, nb, fval, & + wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), & + aocg1(1:na, irow), aocg1(1:nb, icol), & + bocg1(1:na, irow), bocg1(1:nb, icol)) + END DO + CALL dbcsr_iterator_stop(iter) + + DEALLOCATE (wabk) + DEALLOCATE (aocg1, bocg1) + END IF + + CALL timestop(handle) + + END SUBROUTINE xtb_spinpol_hessian + ! ************************************************************************************************** !> \brief ... !> \param aksb ... diff --git a/tests/xTB/regtest-spinpol/TEST_FILES.toml b/tests/xTB/regtest-spinpol/TEST_FILES.toml index d431d9102f..1a90adeb38 100644 --- a/tests/xTB/regtest-spinpol/TEST_FILES.toml +++ b/tests/xTB/regtest-spinpol/TEST_FILES.toml @@ -10,4 +10,6 @@ # "ch2o_kp_force.inp" = [] "ch2o_kp_stress.inp" = [] +# +"ch2o_polar.inp" = [] #EOF diff --git a/tests/xTB/regtest-spinpol/ch2o_polar.inp b/tests/xTB/regtest-spinpol/ch2o_polar.inp new file mode 100644 index 0000000000..01247c697d --- /dev/null +++ b/tests/xTB/regtest-spinpol/ch2o_polar.inp @@ -0,0 +1,63 @@ +&GLOBAL + PRINT_LEVEL LOW + PROJECT CH2O + RUN_TYPE DEBUG +&END GLOBAL + +&DEBUG + DE 0.0002 + DEBUG_DIPOLE F + DEBUG_FORCES F + DEBUG_POLARIZABILITY T + DEBUG_STRESS_TENSOR F + STOP_ON_MISMATCH T +&END DEBUG + +&FORCE_EVAL + &DFT + LSD + MULTIPLICITY 3 + &EFIELD + &END EFIELD + &PRINT + &MOMENTS ON + PERIODIC F + REFERENCE COM + &END MOMENTS + &END PRINT + &QS + METHOD xTB + &XTB + CHECK_ATOMIC_CHARGES F + SPIN_POLARISATION T + &END XTB + &END QS + &SCF + EPS_SCF 1.e-7 + MAX_SCF 100 + SCF_GUESS MOPAC + &END SCF + &END DFT + &PROPERTIES + &LINRES + EPS 1.e-10 + PRECONDITIONER FULL_ALL + &POLAR + DO_RAMAN T + PERIODIC_DIPOLE_OPERATOR F + &END POLAR + &END LINRES + &END PROPERTIES + &SUBSYS + &CELL + ABC 20.0 20.0 20.0 + PERIODIC NONE + &END CELL + &COORD + O 0.051368 0.000000 0.000000 + C 1.278612 0.000000 0.000000 + H 1.870460 0.939607 0.000000 + H 1.870460 -0.939607 0.000000 + &END COORD + &END SUBSYS +&END FORCE_EVAL