(git:71c3ab0)
Loading...
Searching...
No Matches
xtb_spinpol.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculation of Spin Polarisation contributions in xTB
10!> \author JGH
11! **************************************************************************************************
16 USE atprop_types, ONLY: atprop_type
17 USE bibliography, ONLY: neugebauer2023,&
18 cite_reference
19 USE cell_types, ONLY: cell_type
29 USE kinds, ONLY: dp
30 USE kpoint_types, ONLY: get_kpoint_info,&
33 USE mulliken, ONLY: ao_charges
39 USE qs_kind_types, ONLY: get_qs_kind,&
50 USE virial_types, ONLY: virial_type
51 USE xtb_types, ONLY: get_xtb_atom_param,&
53#include "./base/base_uses.f90"
54
55 IMPLICIT NONE
56
57 PRIVATE
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_spinpol'
60
62
63CONTAINS
64
65! **************************************************************************************************
66!> \brief ...
67!> \param qs_env ...
68!> \param ks_matrix ...
69!> \param matrix_p ...
70!> \param energy ...
71!> \param sap_int ...
72!> \param calculate_forces ...
73!> \param just_energy ...
74! **************************************************************************************************
75 SUBROUTINE build_xtb_spinpol(qs_env, ks_matrix, matrix_p, energy, &
76 sap_int, calculate_forces, just_energy)
77
78 TYPE(qs_environment_type), POINTER :: qs_env
79 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, matrix_p
80 TYPE(qs_energy_type), POINTER :: energy
81 TYPE(sap_int_type), DIMENSION(:), POINTER :: sap_int
82 LOGICAL, INTENT(in) :: calculate_forces, just_energy
83
84 CHARACTER(len=*), PARAMETER :: routinen = 'build_xtb_spinpol'
85
86 INTEGER :: atom_a, atom_i, atom_j, handle, i, ia, iac, iatom, ib, ic, icol, ikind, iknd, &
87 irow, jatom, jkind, jknd, la, lb, na, natom, natorb, nb, nimg, nkind, nsgf, nspins
88 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
89 INTEGER, DIMENSION(25) :: lao
90 INTEGER, DIMENSION(3) :: cellind
91 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
92 LOGICAL :: defined, found, use_virial
93 REAL(kind=dp) :: dr, espin, fi, fval
94 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: docg
95 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, bocg, pam, pbm, wab
96 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
97 REAL(kind=dp), DIMENSION(3) :: fij, rij
98 REAL(kind=dp), DIMENSION(3, 3) :: wall
99 REAL(kind=dp), DIMENSION(:, :), POINTER :: aksb, bksb, dsblock, pamat, pbmat, sblock
100 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dsint
101 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
102 TYPE(atprop_type), POINTER :: atprop
103 TYPE(cell_type), POINTER :: cell
104 TYPE(dbcsr_iterator_type) :: iter
105 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
106 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p_kp, matrix_s, matrix_s_kp
107 TYPE(dbcsr_type), POINTER :: s_matrix
108 TYPE(dft_control_type), POINTER :: dft_control
109 TYPE(kpoint_type), POINTER :: kpoints
110 TYPE(mp_para_env_type), POINTER :: para_env
112 DIMENSION(:), POINTER :: nl_iterator
113 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
114 POINTER :: n_list
115 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
116 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
117 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118 TYPE(virial_type), POINTER :: virial
119 TYPE(xtb_atom_type), POINTER :: xtb_kind
120
121 CALL timeset(routinen, handle)
122
123 energy%xtb_spinpol = 0.0_dp
124
125 CALL get_qs_env(qs_env, dft_control=dft_control)
126 nspins = dft_control%nspins
127 nimg = dft_control%nimages
128
129 IF (nspins == 2) THEN
130
131 CALL cite_reference(neugebauer2023)
132
133 CALL get_qs_env(qs_env, &
134 qs_kind_set=qs_kind_set, &
135 particle_set=particle_set, &
136 atomic_kind_set=atomic_kind_set, &
137 cell=cell, &
138 virial=virial, &
139 atprop=atprop)
140
141 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
142 kind_of=kind_of, &
143 atom_of_kind=atom_of_kind)
144
145 use_virial = .false.
146 IF (calculate_forces) THEN
147 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
148 END IF
149
150 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
151 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
152 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
153
154 ! expand parameters
155 ALLOCATE (wabk(nsgf, nsgf, nkind))
156 wabk = 0.0_dp
157 DO ikind = 1, nkind
158 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
159 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall)
160 DO ia = 1, natorb
161 la = lao(ia) + 1
162 DO ib = 1, natorb
163 lb = lao(ib) + 1
164 wabk(ia, ib, ikind) = wall(la, lb)
165 END DO
166 END DO
167 END DO
168
169 ! Calculate charges
170 ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
171 aocg = 0.0_dp
172 bocg = 0.0_dp
173 IF (nimg > 1) THEN
174 matrix_s_kp => matrix_s(:, :)
175 matrix_p_kp => matrix_p(1:1, :)
176 CALL ao_charges(matrix_p_kp, matrix_s_kp, aocg, para_env)
177 matrix_p_kp => matrix_p(2:2, :)
178 CALL ao_charges(matrix_p_kp, matrix_s_kp, bocg, para_env)
179 ELSE
180 s_matrix => matrix_s(1, 1)%matrix
181 p_matrix => matrix_p(1:1, 1)
182 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
183 p_matrix => matrix_p(2:2, 1)
184 CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
185 END IF
186
187 ! calculate energy
188 DO ikind = 1, nkind
189 CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
190 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
191 CALL get_xtb_atom_param(xtb_kind, defined=defined, natorb=natorb)
192 IF (.NOT. defined .OR. natorb < 1) cycle
193 ALLOCATE (docg(natorb), wab(natorb, natorb))
194 wab(1:natorb, 1:natorb) = wabk(1:natorb, 1:natorb, ikind)
195 DO iatom = 1, na
196 atom_a = atomic_kind_set(ikind)%atom_list(iatom)
197 docg = 0.0_dp
198 docg(1:natorb) = aocg(1:natorb, atom_a) - bocg(1:natorb, atom_a)
199 espin = 0.5_dp*dot_product(docg, matmul(wab, docg))
200 energy%xtb_spinpol = energy%xtb_spinpol + espin
201 IF (atprop%energy) THEN
202 atprop%atecoul(iatom) = atprop%atecoul(iatom) + espin
203 END IF
204 END DO
205 DEALLOCATE (docg, wab)
206 END DO
207
208 ! Forces and Virial
209 IF (calculate_forces) THEN
210 CALL get_qs_env(qs_env=qs_env, force=force)
211 NULLIFY (cell_to_index)
212 IF (nimg > 1) THEN
213 NULLIFY (kpoints)
214 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
215 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
216 END IF
217 IF (nimg == 1) THEN
218 ! no k-points; all matrices have been transformed to periodic bsf
219 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
220 DO WHILE (dbcsr_iterator_blocks_left(iter))
221 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
222 ikind = kind_of(irow)
223 atom_i = atom_of_kind(irow)
224 jkind = kind_of(icol)
225 atom_j = atom_of_kind(icol)
226
227 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
228 row=irow, col=icol, block=pamat, found=found)
229 cpassert(found)
230 CALL dbcsr_get_block_p(matrix=matrix_p(2, 1)%matrix, &
231 row=irow, col=icol, block=pbmat, found=found)
232 cpassert(found)
233
234 na = SIZE(pamat, 1)
235 nb = SIZE(pamat, 2)
236
237 DO i = 1, 3
238 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
239 row=irow, col=icol, block=dsblock, found=found)
240 cpassert(found)
241
242 CALL fupdate(fi, pamat, pbmat, dsblock, na, nb, &
243 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
244 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
245
246 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
247 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
248 END DO
249
250 END DO
251 CALL dbcsr_iterator_stop(iter)
252 ! use dsint list
253 IF (use_virial .AND. 0 == 0) THEN
254 cpassert(ASSOCIATED(sap_int))
255 DO ikind = 1, nkind
256 DO jkind = 1, nkind
257 iac = ikind + nkind*(jkind - 1)
258 IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) cycle
259 DO ia = 1, sap_int(iac)%nalist
260 IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) cycle
261 iatom = sap_int(iac)%alist(ia)%aatom
262 DO ic = 1, sap_int(iac)%alist(ia)%nclist
263 jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
264 rij = sap_int(iac)%alist(ia)%clist(ic)%rac
265 dr = sqrt(sum(rij(:)**2))
266 IF (dr > 1.e-6_dp) THEN
267 dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
268 icol = max(iatom, jatom)
269 irow = min(iatom, jatom)
270 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
271 row=irow, col=icol, block=pamat, found=found)
272 cpassert(found)
273 CALL dbcsr_get_block_p(matrix=matrix_p(2, 1)%matrix, &
274 row=irow, col=icol, block=pbmat, found=found)
275 cpassert(found)
276 IF (irow == iatom) THEN
277 na = SIZE(pamat, 1)
278 nb = SIZE(pamat, 2)
279 ALLOCATE (pam(na, nb), pbm(na, nb))
280 pam(1:na, 1:nb) = pamat(1:na, 1:nb)
281 pbm(1:na, 1:nb) = pbmat(1:na, 1:nb)
282 ELSE
283 na = SIZE(pamat, 2)
284 nb = SIZE(pamat, 1)
285 ALLOCATE (pam(na, nb), pbm(na, nb))
286 pam(1:na, 1:nb) = transpose(pamat(1:nb, 1:na))
287 pbm(1:na, 1:nb) = transpose(pbmat(1:nb, 1:na))
288 END IF
289
290 DO i = 1, 3
291 CALL fupdate(fi, pam, pbm, dsint(:, :, i), na, nb, &
292 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
293 aocg(1:na, iatom), aocg(1:nb, jatom), &
294 bocg(1:na, iatom), bocg(1:nb, jatom))
295 fij(i) = fi
296 END DO
297 fi = 1.0_dp
298 IF (iatom == jatom) fi = 0.5_dp
299 CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
300 DEALLOCATE (pam, pbm)
301
302 END IF
303 END DO
304 END DO
305 END DO
306 END DO
307 END IF
308 ELSE
309 NULLIFY (n_list)
310 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
311 CALL neighbor_list_iterator_create(nl_iterator, n_list)
312 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
313 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
314 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
315
316 dr = sqrt(sum(rij**2))
317 IF (iatom == jatom .AND. dr < 1.0e-6_dp) cycle
318
319 icol = max(iatom, jatom)
320 irow = min(iatom, jatom)
321
322 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
323 cpassert(ic > 0)
324
325 IF (irow == iatom) THEN
326 iknd = ikind
327 jknd = jkind
328 atom_i = atom_of_kind(iatom)
329 atom_j = atom_of_kind(jatom)
330 rij = rij
331 ELSE
332 iknd = jkind
333 jknd = ikind
334 atom_i = atom_of_kind(jatom)
335 atom_j = atom_of_kind(iatom)
336 rij = -rij
337 END IF
338 !
339 CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
340 row=irow, col=icol, block=pamat, found=found)
341 cpassert(found)
342 CALL dbcsr_get_block_p(matrix=matrix_p(2, ic)%matrix, &
343 row=irow, col=icol, block=pbmat, found=found)
344 cpassert(found)
345
346 na = SIZE(pamat, 1)
347 nb = SIZE(pamat, 2)
348
349 fij = 0.0_dp
350 DO i = 1, 3
351 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
352 row=irow, col=icol, block=dsblock, found=found)
353 cpassert(found)
354
355 CALL fupdate(fi, pamat, pbmat, dsblock, na, nb, &
356 wabk(1:na, 1:na, iknd), wabk(1:nb, 1:nb, jknd), &
357 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
358
359 force(iknd)%rho_elec(i, atom_i) = force(iknd)%rho_elec(i, atom_i) + fi
360 force(jknd)%rho_elec(i, atom_j) = force(jknd)%rho_elec(i, atom_j) - fi
361 fij(i) = fi
362 END DO
363 IF (use_virial) THEN
364 fi = 1.0_dp
365 IF (iatom == jatom) fi = 0.5_dp
366 CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
367 END IF
368
369 END DO
370 CALL neighbor_list_iterator_release(nl_iterator)
371
372 END IF
373 END IF
374
375 ! KS matrix
376 IF (.NOT. just_energy) THEN
377 IF (nimg > 1) THEN
378 CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
379 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
380 END IF
381 IF (nimg == 1) THEN
382 ! no k-points; all matrices have been transformed to periodic bsf
383 CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
384 DO WHILE (dbcsr_iterator_blocks_left(iter))
385 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
386 CALL dbcsr_get_block_p(matrix=ks_matrix(1, 1)%matrix, &
387 row=irow, col=icol, block=aksb, found=found)
388 cpassert(found)
389 CALL dbcsr_get_block_p(matrix=ks_matrix(2, 1)%matrix, &
390 row=irow, col=icol, block=bksb, found=found)
391 cpassert(found)
392 na = SIZE(aksb, 1)
393 nb = SIZE(aksb, 2)
394 ikind = kind_of(irow)
395 jkind = kind_of(icol)
396 fval = 0.5_dp
397 CALL ksupdate(aksb, bksb, sblock, na, nb, fval, &
398 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
399 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
400 END DO
401 CALL dbcsr_iterator_stop(iter)
402 ELSE
403 CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
404 CALL neighbor_list_iterator_create(nl_iterator, n_list)
405 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
406 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
407 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
408
409 icol = max(iatom, jatom)
410 irow = min(iatom, jatom)
411
412 ic = cell_to_index(cellind(1), cellind(2), cellind(3))
413 cpassert(ic > 0)
414
415 ikind = kind_of(irow)
416 jkind = kind_of(icol)
417
418 CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
419 row=irow, col=icol, block=sblock, found=found)
420 cpassert(found)
421 CALL dbcsr_get_block_p(matrix=ks_matrix(1, ic)%matrix, &
422 row=irow, col=icol, block=aksb, found=found)
423 cpassert(found)
424 CALL dbcsr_get_block_p(matrix=ks_matrix(2, ic)%matrix, &
425 row=irow, col=icol, block=bksb, found=found)
426 cpassert(found)
427
428 na = SIZE(aksb, 1)
429 nb = SIZE(aksb, 2)
430 fval = 0.5_dp
431 CALL ksupdate(aksb, bksb, sblock, na, nb, fval, &
432 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
433 aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
434 END DO
435 CALL neighbor_list_iterator_release(nl_iterator)
436 END IF
437
438 END IF
439
440 DEALLOCATE (wabk)
441 DEALLOCATE (aocg, bocg)
442 END IF
443
444 CALL timestop(handle)
445
446 END SUBROUTINE build_xtb_spinpol
447
448! **************************************************************************************************
449!> \brief ...
450!> \param qs_env ...
451!> \param ks_matrix ...
452!> \param matrix_p1 ...
453! **************************************************************************************************
454 SUBROUTINE xtb_spinpol_hessian(qs_env, ks_matrix, matrix_p1)
455
456 TYPE(qs_environment_type), POINTER :: qs_env
457 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p1
458
459 CHARACTER(len=*), PARAMETER :: routinen = 'xtb_spinpol_hessian'
460
461 INTEGER :: handle, ia, ib, icol, ikind, irow, &
462 jkind, la, lb, na, natom, natorb, nb, &
463 nimg, nkind, nsgf, nspins
464 INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
465 INTEGER, DIMENSION(25) :: lao
466 LOGICAL :: found
467 REAL(kind=dp) :: fval
468 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg1, bocg1
469 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
470 REAL(kind=dp), DIMENSION(3, 3) :: wall
471 REAL(kind=dp), DIMENSION(:, :), POINTER :: aksb, bksb, sblock
472 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
473 TYPE(dbcsr_iterator_type) :: iter
474 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
475 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
476 TYPE(dbcsr_type), POINTER :: s_matrix
477 TYPE(dft_control_type), POINTER :: dft_control
478 TYPE(mp_para_env_type), POINTER :: para_env
479 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
480 TYPE(xtb_atom_type), POINTER :: xtb_kind
481
482 CALL timeset(routinen, handle)
483
484 CALL get_qs_env(qs_env, dft_control=dft_control)
485 nspins = dft_control%nspins
486 nimg = dft_control%nimages
487
488 IF (nimg /= 1) THEN
489 cpabort("No kpoints allowed in xTB response calculation")
490 END IF
491
492 IF (nspins == 2) THEN
493
494 CALL get_qs_env(qs_env, &
495 qs_kind_set=qs_kind_set, &
496 atomic_kind_set=atomic_kind_set)
497
498 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
499
500 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
501 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
502 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
503
504 ! expand parameters
505 ALLOCATE (wabk(nsgf, nsgf, nkind))
506 wabk = 0.0_dp
507 DO ikind = 1, nkind
508 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
509 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall)
510 DO ia = 1, natorb
511 la = lao(ia) + 1
512 DO ib = 1, natorb
513 lb = lao(ib) + 1
514 wabk(ia, ib, ikind) = wall(la, lb)
515 END DO
516 END DO
517 END DO
518
519 ! Calculate response charges
520 ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
521 aocg1 = 0.0_dp
522 bocg1 = 0.0_dp
523 s_matrix => matrix_s(1, 1)%matrix
524 p_matrix => matrix_p1(1:1)
525 CALL ao_charges(p_matrix, s_matrix, aocg1, para_env)
526 p_matrix => matrix_p1(2:2)
527 CALL ao_charges(p_matrix, s_matrix, bocg1, para_env)
528 aocg1 = 0.5_dp*aocg1
529 bocg1 = 0.5_dp*bocg1
530
531 CALL dbcsr_iterator_start(iter, s_matrix)
532 DO WHILE (dbcsr_iterator_blocks_left(iter))
533 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
534 CALL dbcsr_get_block_p(matrix=ks_matrix(1)%matrix, &
535 row=irow, col=icol, block=aksb, found=found)
536 cpassert(found)
537 CALL dbcsr_get_block_p(matrix=ks_matrix(2)%matrix, &
538 row=irow, col=icol, block=bksb, found=found)
539 cpassert(found)
540 na = SIZE(aksb, 1)
541 nb = SIZE(aksb, 2)
542 ikind = kind_of(irow)
543 jkind = kind_of(icol)
544 fval = 1.00_dp
545 CALL ksupdate(aksb, bksb, sblock, na, nb, fval, &
546 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
547 aocg1(1:na, irow), aocg1(1:nb, icol), &
548 bocg1(1:na, irow), bocg1(1:nb, icol))
549 END DO
550 CALL dbcsr_iterator_stop(iter)
551
552 DEALLOCATE (wabk)
553 DEALLOCATE (aocg1, bocg1)
554 END IF
555
556 CALL timestop(handle)
557
558 END SUBROUTINE xtb_spinpol_hessian
559
560! **************************************************************************************************
561!> \brief ...
562!> \param qs_env ...
563!> \param matrix_p0 ...
564!> \param matrix_p1 ...
565! **************************************************************************************************
566 SUBROUTINE xtb_spinpol_hforce(qs_env, matrix_p0, matrix_p1)
567 TYPE(qs_environment_type), POINTER :: qs_env
568 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p0, matrix_p1
569
570 CHARACTER(len=*), PARAMETER :: routinen = 'xtb_spinpol_hforce'
571
572 INTEGER :: atom_i, atom_j, handle, i, ia, ib, icol, &
573 ikind, irow, jkind, la, lb, na, natom, &
574 natorb, nb, nimg, nkind, nsgf, nspins
575 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
576 INTEGER, DIMENSION(25) :: lao
577 LOGICAL :: found
578 REAL(kind=dp) :: fi
579 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, bocg, bocg1
580 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
581 REAL(kind=dp), DIMENSION(3, 3) :: wall
582 REAL(kind=dp), DIMENSION(:, :), POINTER :: dsblock, p0amat, p0bmat, p1amat, p1bmat, &
583 sblock
584 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
585 TYPE(dbcsr_iterator_type) :: iter
586 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, p_matrix
587 TYPE(dbcsr_type), POINTER :: s_matrix
588 TYPE(dft_control_type), POINTER :: dft_control
589 TYPE(mp_para_env_type), POINTER :: para_env
590 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
591 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
592 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
593 TYPE(xtb_atom_type), POINTER :: xtb_kind
594
595 CALL timeset(routinen, handle)
596
597 CALL get_qs_env(qs_env, dft_control=dft_control)
598 nspins = dft_control%nspins
599 nimg = dft_control%nimages
600 IF (nimg /= 1) THEN
601 cpabort("xTB response forces for spin polarisation Hamiltonian not available")
602 END IF
603
604 IF (nspins == 2) THEN
605
606 CALL get_qs_env(qs_env, &
607 qs_kind_set=qs_kind_set, &
608 particle_set=particle_set, &
609 atomic_kind_set=atomic_kind_set)
610
611 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
612 kind_of=kind_of, &
613 atom_of_kind=atom_of_kind)
614
615 CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
616 CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
617
618 CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env)
619
620 ! expand parameters
621 ALLOCATE (wabk(nsgf, nsgf, nkind))
622 wabk = 0.0_dp
623 DO ikind = 1, nkind
624 CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
625 CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall)
626 DO ia = 1, natorb
627 la = lao(ia) + 1
628 DO ib = 1, natorb
629 lb = lao(ib) + 1
630 wabk(ia, ib, ikind) = wall(la, lb)
631 END DO
632 END DO
633 END DO
634
635 ! Calculate charges
636 s_matrix => matrix_s(1)%matrix
637 ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
638 aocg = 0.0_dp
639 bocg = 0.0_dp
640 p_matrix => matrix_p0(1:1)
641 CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
642 p_matrix => matrix_p0(2:2)
643 CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
644 ! Calculate response charges
645 ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
646 aocg1 = 0.0_dp
647 bocg1 = 0.0_dp
648 p_matrix => matrix_p1(1:1)
649 CALL ao_charges(p_matrix, s_matrix, aocg1, para_env)
650 p_matrix => matrix_p1(2:2)
651 CALL ao_charges(p_matrix, s_matrix, bocg1, para_env)
652 aocg1 = 0.5_dp*aocg1
653 bocg1 = 0.5_dp*bocg1
654
655 ! calculate forces
656 CALL get_qs_env(qs_env=qs_env, force=force)
657 ! no k-points; all matrices have been transformed to periodic bsf
658 CALL dbcsr_iterator_start(iter, s_matrix)
659 DO WHILE (dbcsr_iterator_blocks_left(iter))
660 CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
661 ikind = kind_of(irow)
662 atom_i = atom_of_kind(irow)
663 jkind = kind_of(icol)
664 atom_j = atom_of_kind(icol)
665
666 CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
667 row=irow, col=icol, block=p0amat, found=found)
668 cpassert(found)
669 CALL dbcsr_get_block_p(matrix=matrix_p0(2)%matrix, &
670 row=irow, col=icol, block=p0bmat, found=found)
671 cpassert(found)
672 CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
673 row=irow, col=icol, block=p1amat, found=found)
674 cpassert(found)
675 CALL dbcsr_get_block_p(matrix=matrix_p1(2)%matrix, &
676 row=irow, col=icol, block=p1bmat, found=found)
677 cpassert(found)
678
679 na = SIZE(p0amat, 1)
680 nb = SIZE(p0amat, 2)
681
682 DO i = 1, 3
683 CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
684 row=irow, col=icol, block=dsblock, found=found)
685 cpassert(found)
686
687 fi = 0.0_dp
688 CALL f2update(fi, p0amat, p0bmat, p1amat, p1bmat, dsblock, na, nb, &
689 wabk(1:na, 1:na, ikind), wabk(1:nb, 1:nb, jkind), &
690 aocg(1:na, irow), aocg(1:nb, icol), &
691 bocg(1:na, irow), bocg(1:nb, icol), &
692 aocg1(1:na, irow), aocg1(1:nb, icol), &
693 bocg1(1:na, irow), bocg1(1:nb, icol))
694
695 force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
696 force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
697 END DO
698
699 END DO
700 CALL dbcsr_iterator_stop(iter)
701
702 END IF
703
704 CALL timestop(handle)
705
706 END SUBROUTINE xtb_spinpol_hforce
707
708! **************************************************************************************************
709!> \brief ...
710!> \param aksb ...
711!> \param bksb ...
712!> \param sb ...
713!> \param na ...
714!> \param nb ...
715!> \param fval ...
716!> \param wabi ...
717!> \param wabj ...
718!> \param qai ...
719!> \param qaj ...
720!> \param qbi ...
721!> \param qbj ...
722! **************************************************************************************************
723 SUBROUTINE ksupdate(aksb, bksb, sb, na, nb, fval, &
724 wabi, wabj, qai, qaj, qbi, qbj)
725
726 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: aksb, bksb
727 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: sb
728 INTEGER, INTENT(IN) :: na, nb
729 REAL(kind=dp), INTENT(IN) :: fval
730 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: wabi, wabj
731 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: qai, qaj, qbi, qbj
732
733 INTEGER :: ia, ib
734 REAL(kind=dp), DIMENSION(na) :: dqa, wa
735 REAL(kind=dp), DIMENSION(na, nb) :: wqab
736 REAL(kind=dp), DIMENSION(nb) :: dqb, wb
737
738 dqa = qai - qbi
739 dqb = qaj - qbj
740 wa = matmul(wabi, dqa)
741 wb = matmul(wabj, dqb)
742 DO ib = 1, nb
743 DO ia = 1, na
744 wqab(ia, ib) = fval*sb(ia, ib)*(wa(ia) + wb(ib))
745 END DO
746 END DO
747
748 aksb = aksb + wqab
749 bksb = bksb - wqab
750
751 END SUBROUTINE ksupdate
752
753! **************************************************************************************************
754!> \brief ...
755!> \param fij ...
756!> \param pa ...
757!> \param pb ...
758!> \param ds ...
759!> \param na ...
760!> \param nb ...
761!> \param wabi ...
762!> \param wabj ...
763!> \param qai ...
764!> \param qaj ...
765!> \param qbi ...
766!> \param qbj ...
767! **************************************************************************************************
768 SUBROUTINE fupdate(fij, pa, pb, ds, na, nb, &
769 wabi, wabj, qai, qaj, qbi, qbj)
770
771 REAL(kind=dp), INTENT(OUT) :: fij
772 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: pa, pb, ds
773 INTEGER, INTENT(IN) :: na, nb
774 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: wabi, wabj
775 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: qai, qaj, qbi, qbj
776
777 INTEGER :: ia, ib
778 REAL(kind=dp), DIMENSION(na) :: dpsa, dqa, wa
779 REAL(kind=dp), DIMENSION(na, nb) :: dpab
780 REAL(kind=dp), DIMENSION(nb) :: dpsb, dqb, wb
781
782 dqa = qai - qbi
783 dqb = qaj - qbj
784 wa = matmul(wabi, dqa)
785 wb = matmul(wabj, dqb)
786 dpab = pa - pb
787 dpsa = 0.0_dp
788 dpsb = 0.0_dp
789 DO ib = 1, nb
790 DO ia = 1, na
791 dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
792 dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
793 END DO
794 END DO
795
796 fij = sum(wa*dpsa) + sum(wb*dpsb)
797
798 END SUBROUTINE fupdate
799
800! **************************************************************************************************
801!> \brief ...
802!> \param fij ...
803!> \param p0a ...
804!> \param p0b ...
805!> \param p1a ...
806!> \param p1b ...
807!> \param ds ...
808!> \param na ...
809!> \param nb ...
810!> \param wabi ...
811!> \param wabj ...
812!> \param qai ...
813!> \param qaj ...
814!> \param qbi ...
815!> \param qbj ...
816!> \param rai ...
817!> \param raj ...
818!> \param rbi ...
819!> \param rbj ...
820! **************************************************************************************************
821 SUBROUTINE f2update(fij, p0a, p0b, p1a, p1b, ds, na, nb, wabi, wabj, &
822 qai, qaj, qbi, qbj, rai, raj, rbi, rbj)
823
824 REAL(kind=dp), INTENT(OUT) :: fij
825 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: p0a, p0b, p1a, p1b, ds
826 INTEGER, INTENT(IN) :: na, nb
827 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: wabi, wabj
828 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: qai, qaj, qbi, qbj, rai, raj, rbi, rbj
829
830 INTEGER :: ia, ib
831 REAL(kind=dp), DIMENSION(na) :: dpsa, dqa, wa
832 REAL(kind=dp), DIMENSION(na, nb) :: dpab
833 REAL(kind=dp), DIMENSION(nb) :: dpsb, dqb, wb
834
835 fij = 0.0_dp
836
837 dqa = qai - qbi
838 dqb = qaj - qbj
839 wa = matmul(wabi, dqa)
840 wb = matmul(wabj, dqb)
841 dpab = p1a - p1b
842 dpsa = 0.0_dp
843 dpsb = 0.0_dp
844 DO ib = 1, nb
845 DO ia = 1, na
846 dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
847 dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
848 END DO
849 END DO
850
851 fij = fij + sum(wa*dpsa) + sum(wb*dpsb)
852
853 dqa = rai - rbi
854 dqb = raj - rbj
855 wa = matmul(wabi, dqa)
856 wb = matmul(wabj, dqb)
857 dpab = p0a - p0b
858 dpsa = 0.0_dp
859 dpsb = 0.0_dp
860 DO ib = 1, nb
861 DO ia = 1, na
862 dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
863 dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
864 END DO
865 END DO
866
867 fij = fij + sum(wa*dpsa) + sum(wb*dpsb)
868
869 END SUBROUTINE f2update
870
871END MODULE xtb_spinpol
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Holds information on atomic properties.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public neugebauer2023
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
logical function, public dbcsr_iterator_blocks_left(iterator)
...
subroutine, public dbcsr_iterator_stop(iterator)
...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered)
Retrieve information from a kpoint environment.
Interface to the message passing library MPI.
compute mulliken charges we (currently) define them as c_i = 1/2 [ (PS)_{ii} + (SP)_{ii} ]
Definition mulliken.F:13
Define the data structure for the particle information.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
General overlap type integrals containers.
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Calculation of Spin Polarisation contributions in xTB.
Definition xtb_spinpol.F:12
subroutine, public xtb_spinpol_hessian(qs_env, ks_matrix, matrix_p1)
...
subroutine, public build_xtb_spinpol(qs_env, ks_matrix, matrix_p, energy, sap_int, calculate_forces, just_energy)
...
Definition xtb_spinpol.F:77
subroutine, public xtb_spinpol_hforce(qs_env, matrix_p0, matrix_p1)
...
Definition of the xTB parameter types.
Definition xtb_types.F:20
subroutine, public get_xtb_atom_param(xtb_parameter, symbol, aname, typ, defined, z, zeff, natorb, lmax, nao, lao, rcut, rcov, kx, eta, xgamma, alpha, zneff, nshell, nval, lval, kpoly, kappa, wall, hen, zeta, xi, kappa0, alpg, occupation, electronegativity, chmax, en, kqat2, kcn, kq)
...
Definition xtb_types.F:203
Provides all information about an atomic kind.
type for the atomic properties
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
Contains information about kpoints.
stores all the informations relevant to an mpi environment
Provides all information about a quickstep kind.