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

xTB spinpol response forces

Browse source
This commit is contained in:
Juerg Hutter 2026-07-17 12:08:40 +02:00
parent 559f9a7b6f
commit c036843911
4 changed files with 306 additions and 3 deletions
Download patch file Download diff file Expand all files Collapse all files

13
src/xtb_ehess_force.F
View file

@ -54,6 +54,7 @@ MODULE xtb_ehess_force
USE virial_types, ONLY: virial_type
USE xtb_coulomb, ONLY: dgamma_rab_sr,&
gamma_rab_sr
USE xtb_spinpol, ONLY: xtb_spinpol_hforce
USE xtb_types, ONLY: get_xtb_atom_param,&
xtb_atom_type
#include "./base/base_uses.f90"
@ -436,6 +437,18 @@ CONTAINS
alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
END IF
IF (xtb_control%do_spinpol) THEN
IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
!
CALL xtb_spinpol_hforce(qs_env, matrix_p0, matrix_p1)
!
IF (debug_forces) THEN
fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
CALL para_env%sum(fodeb)
IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*Hspin[P] ", fodeb
END IF
END IF
! QMMM
IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
CPABORT("Not Available")

228
src/xtb_spinpol.F
View file

@ -58,7 +58,7 @@ MODULE xtb_spinpol
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_spinpol'
PUBLIC :: build_xtb_spinpol, xtb_spinpol_hessian
PUBLIC :: build_xtb_spinpol, xtb_spinpol_hessian, xtb_spinpol_hforce
CONTAINS
@ -563,6 +563,156 @@ CONTAINS
END SUBROUTINE xtb_spinpol_hessian
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param matrix_p0 ...
!> \param matrix_p1 ...
! **************************************************************************************************
SUBROUTINE xtb_spinpol_hforce(qs_env, matrix_p0, matrix_p1)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p0, matrix_p1
CHARACTER(len=*), PARAMETER :: routineN = 'xtb_spinpol_hforce'
INTEGER :: atom_i, atom_j, handle, i, ia, ib, icol, &
ikind, irow, jkind, la, lb, na, natom, &
natorb, nb, nimg, nkind, nsgf, nspins
INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
INTEGER, DIMENSION(25) :: lao
LOGICAL :: found
REAL(KIND=dp) :: fi, fval
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, bocg, bocg1
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
REAL(KIND=dp), DIMENSION(3, 3) :: wall
REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, p0amat, p0bmat, p1amat, p1bmat, &
sblock
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dbcsr_iterator_type) :: iter
TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, p_matrix
TYPE(dbcsr_type), POINTER :: s_matrix
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
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("xTB response forces for spin polarisation Hamiltonian not available")
END IF
IF (nspins == 2) THEN
CALL get_qs_env(qs_env, &
qs_kind_set=qs_kind_set, &
particle_set=particle_set, &
atomic_kind_set=atomic_kind_set)
CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
kind_of=kind_of, &
atom_of_kind=atom_of_kind)
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=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 charges
s_matrix => matrix_s(1)%matrix
ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
aocg = 0.0_dp
bocg = 0.0_dp
p_matrix => matrix_p0(1:1)
CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
p_matrix => matrix_p0(2:2)
CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
! Calculate response charges
ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
aocg1 = 0.0_dp
bocg1 = 0.0_dp
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
! calculate forces
CALL get_qs_env(qs_env=qs_env, force=force)
! no k-points; all matrices have been transformed to periodic bsf
CALL dbcsr_iterator_start(iter, s_matrix)
DO WHILE (dbcsr_iterator_blocks_left(iter))
CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
ikind = kind_of(irow)
atom_i = atom_of_kind(irow)
jkind = kind_of(icol)
atom_j = atom_of_kind(icol)
CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
row=irow, col=icol, block=p0amat, found=found)
CPASSERT(found)
CALL dbcsr_get_block_p(matrix=matrix_p0(2)%matrix, &
row=irow, col=icol, block=p0bmat, found=found)
CPASSERT(found)
CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
row=irow, col=icol, block=p1amat, found=found)
CPASSERT(found)
CALL dbcsr_get_block_p(matrix=matrix_p1(2)%matrix, &
row=irow, col=icol, block=p1bmat, found=found)
CPASSERT(found)
na = SIZE(p0amat, 1)
nb = SIZE(p0amat, 2)
fval = 1.0_dp
DO i = 1, 3
CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
row=irow, col=icol, block=dsblock, found=found)
CPASSERT(found)
fi = 0.0_dp
CALL f2update(fi, p0amat, p0bmat, p1amat, p1bmat, dsblock, na, nb, fval, &
wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
aocg(1:na, irow), aocg(1:nb, icol), &
bocg(1:na, irow), bocg(1:nb, icol), &
aocg1(1:na, irow), aocg1(1:nb, icol), &
bocg1(1:na, irow), bocg1(1:nb, icol))
force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
END DO
END DO
CALL dbcsr_iterator_stop(iter)
END IF
CALL timestop(handle)
END SUBROUTINE xtb_spinpol_hforce
! **************************************************************************************************
!> \brief ...
!> \param aksb ...
@ -628,8 +778,7 @@ CONTAINS
wabi, wabj, qai, qaj, qbi, qbj)
REAL(KIND=dp), INTENT(OUT) :: fij
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pa, pb
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ds
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: pa, pb, ds
INTEGER, INTENT(IN) :: na, nb
REAL(KIND=dp), INTENT(IN) :: fval
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: wabi, wabj
@ -658,4 +807,77 @@ CONTAINS
END SUBROUTINE fupdate
! **************************************************************************************************
!> \brief ...
!> \param fij ...
!> \param p0a ...
!> \param p0b ...
!> \param p1a ...
!> \param p1b ...
!> \param ds ...
!> \param na ...
!> \param nb ...
!> \param fval ...
!> \param wabi ...
!> \param wabj ...
!> \param qai ...
!> \param qaj ...
!> \param qbi ...
!> \param qbj ...
!> \param rai ...
!> \param raj ...
!> \param rbi ...
!> \param rbj ...
! **************************************************************************************************
SUBROUTINE f2update(fij, p0a, p0b, p1a, p1b, ds, na, nb, fval, wabi, wabj, &
qai, qaj, qbi, qbj, rai, raj, rbi, rbj)
REAL(KIND=dp), INTENT(OUT) :: fij
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: p0a, p0b, p1a, p1b, ds
INTEGER, INTENT(IN) :: na, nb
REAL(KIND=dp), INTENT(IN) :: fval
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: wabi, wabj
REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: qai, qaj, qbi, qbj, rai, raj, rbi, rbj
INTEGER :: ia, ib
REAL(KIND=dp), DIMENSION(na) :: dpsa, dqa, wa
REAL(KIND=dp), DIMENSION(na, nb) :: dpab
REAL(KIND=dp), DIMENSION(nb) :: dpsb, dqb, wb
fij = 0.0_dp
dqa = qai - qbi
dqb = qaj - qbj
wa = MATMUL(wabi, dqa)
wb = MATMUL(wabj, dqb)
dpab = p1a - p1b
dpsa = 0.0_dp
dpsb = 0.0_dp
DO ib = 1, nb
DO ia = 1, na
dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
END DO
END DO
fij = fij + SUM(wa*dpsa) + SUM(wb*dpsb)
dqa = rai - rbi
dqb = raj - rbj
wa = MATMUL(wabi, dqa)
wb = MATMUL(wabj, dqb)
dpab = p0a - p0b
dpsa = 0.0_dp
dpsb = 0.0_dp
DO ib = 1, nb
DO ia = 1, na
dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
END DO
END DO
fij = fij + SUM(wa*dpsa) + SUM(wb*dpsb)
END SUBROUTINE f2update
END MODULE xtb_spinpol

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

@ -12,4 +12,5 @@
"ch2o_kp_stress.inp" = []
#
"ch2o_polar.inp" = []
"h2o_sTDA_spin.inp" = []
#EOF

67
tests/xTB/regtest-spinpol/h2o_sTDA_spin.inp Normal file
View file

@ -0,0 +1,67 @@
&GLOBAL
PRINT_LEVEL LOW
PROJECT H2O
RUN_TYPE DEBUG
&END GLOBAL
&DEBUG
CHECK_ATOM_FORCE 1 Z
DE 0.0002
DEBUG_DIPOLE .FALSE.
DEBUG_FORCES .TRUE.
DEBUG_POLARIZABILITY .FALSE.
DEBUG_STRESS_TENSOR .FALSE.
STOP_ON_MISMATCH F
&END DEBUG
&FORCE_EVAL
METHOD Quickstep
&DFT
LSD
MULTIPLICITY 3
&EXCITED_STATES T
STATE 1
&END EXCITED_STATES
&QS
METHOD xTB
&XTB
CHECK_ATOMIC_CHARGES F
SPIN_POLARISATION T
&END XTB
&END QS
&SCF
EPS_SCF 1.0E-8
MAX_SCF 50
SCF_GUESS MOPAC
&END SCF
&END DFT
&PRINT
&FORCES
&END FORCES
&END PRINT
&PROPERTIES
&TDDFPT
CONVERGENCE [eV] 1.0e-7
KERNEL sTDA
MAX_ITER 50
NSTATES 5
&STDA
FRACTION 0.50
&END STDA
&END TDDFPT
&END PROPERTIES
&SUBSYS
&CELL
ABC [angstrom] 6.0 6.0 6.0
&END CELL
&COORD
O 0.000000 0.000000 0.000000
H 0.000000 -0.757136 0.580545
H 0.000000 0.757136 0.580545
&END COORD
&TOPOLOGY
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&END SUBSYS
&END FORCE_EVAL