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

xTB spin polarisation electronic hessian

Browse source
This commit is contained in:
Juerg Hutter 2026-07-15 17:30:31 +02:00
parent 285b52e4fb
commit 1b84f40cd1
7 changed files with 194 additions and 6 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/qs_linres_kernel.F
View file

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

3
src/qs_tddfpt2_forces.F
View file

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

2
src/response_solver.F
View file

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

9
src/xtb_ehess.F
View file

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

116
src/xtb_spinpol.F
View file

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

2
tests/xTB/regtest-spinpol/TEST_FILES.toml
View file

@ -10,4 +10,6 @@
#
"ch2o_kp_force.inp" = []
"ch2o_kp_stress.inp" = []
#
"ch2o_polar.inp" = []
#EOF

63
tests/xTB/regtest-spinpol/ch2o_polar.inp Normal file
View file

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