Line data Source code
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 : ! **************************************************************************************************
12 : MODULE xtb_spinpol
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind,&
15 : get_atomic_kind_set
16 : USE atprop_types, ONLY: atprop_type
17 : USE bibliography, ONLY: Neugebauer2023,&
18 : cite_reference
19 : USE cell_types, ONLY: cell_type
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
22 : dbcsr_iterator_blocks_left,&
23 : dbcsr_iterator_next_block,&
24 : dbcsr_iterator_start,&
25 : dbcsr_iterator_stop,&
26 : dbcsr_iterator_type,&
27 : dbcsr_p_type,&
28 : dbcsr_type
29 : USE kinds, ONLY: dp
30 : USE kpoint_types, ONLY: get_kpoint_info,&
31 : kpoint_type
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE mulliken, ONLY: ao_charges
34 : USE particle_types, ONLY: particle_type
35 : USE qs_energy_types, ONLY: qs_energy_type
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type
38 : USE qs_force_types, ONLY: qs_force_type
39 : USE qs_kind_types, ONLY: get_qs_kind,&
40 : get_qs_kind_set,&
41 : qs_kind_type
42 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
43 : neighbor_list_iterate,&
44 : neighbor_list_iterator_create,&
45 : neighbor_list_iterator_p_type,&
46 : neighbor_list_iterator_release,&
47 : neighbor_list_set_p_type
48 : USE sap_kind_types, ONLY: sap_int_type
49 : USE virial_methods, ONLY: virial_pair_force
50 : USE virial_types, ONLY: virial_type
51 : USE xtb_types, ONLY: get_xtb_atom_param,&
52 : xtb_atom_type
53 : #include "./base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 :
57 : PRIVATE
58 :
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_spinpol'
60 :
61 : PUBLIC :: build_xtb_spinpol, xtb_spinpol_hessian, xtb_spinpol_hforce
62 :
63 : CONTAINS
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 1678 : 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 1678 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
89 : INTEGER, DIMENSION(25) :: lao
90 : INTEGER, DIMENSION(3) :: cellind
91 1678 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
92 : LOGICAL :: defined, found, use_virial
93 : REAL(KIND=dp) :: dr, espin, fi, fval
94 1678 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: docg
95 1678 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, bocg, pam, pbm, wab
96 1678 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
97 : REAL(KIND=dp), DIMENSION(3) :: fij, rij
98 : REAL(KIND=dp), DIMENSION(3, 3) :: wall
99 1678 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: aksb, bksb, dsblock, pamat, pbmat, sblock
100 1678 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dsint
101 1678 : 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 1678 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
106 1678 : 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
111 : TYPE(neighbor_list_iterator_p_type), &
112 1678 : DIMENSION(:), POINTER :: nl_iterator
113 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
114 1678 : POINTER :: n_list
115 1678 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
116 1678 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
117 1678 : 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 1678 : CALL timeset(routineN, handle)
122 :
123 1678 : energy%xtb_spinpol = 0.0_dp
124 :
125 1678 : CALL get_qs_env(qs_env, dft_control=dft_control)
126 1678 : nspins = dft_control%nspins
127 1678 : nimg = dft_control%nimages
128 :
129 1678 : IF (nspins == 2) THEN
130 :
131 1678 : 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 1678 : atprop=atprop)
140 :
141 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
142 : kind_of=kind_of, &
143 1678 : atom_of_kind=atom_of_kind)
144 :
145 1678 : use_virial = .FALSE.
146 1678 : IF (calculate_forces) THEN
147 12 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
148 : END IF
149 :
150 1678 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
151 1678 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
152 1678 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
153 :
154 : ! expand parameters
155 8390 : ALLOCATE (wabk(nsgf, nsgf, nkind))
156 1678 : wabk = 0.0_dp
157 6648 : DO ikind = 1, nkind
158 4970 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
159 4970 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall)
160 24602 : DO ia = 1, natorb
161 17954 : la = lao(ia) + 1
162 100898 : DO ib = 1, natorb
163 77974 : lb = lao(ib) + 1
164 95928 : wabk(ia, ib, ikind) = wall(la, lb)
165 : END DO
166 : END DO
167 : END DO
168 :
169 : ! Calculate charges
170 10068 : ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
171 1678 : aocg = 0.0_dp
172 1678 : bocg = 0.0_dp
173 1678 : IF (nimg > 1) THEN
174 786 : matrix_s_kp => matrix_s(:, :)
175 786 : matrix_p_kp => matrix_p(1:1, :)
176 786 : CALL ao_charges(matrix_p_kp, matrix_s_kp, aocg, para_env)
177 786 : matrix_p_kp => matrix_p(2:2, :)
178 786 : CALL ao_charges(matrix_p_kp, matrix_s_kp, bocg, para_env)
179 : ELSE
180 892 : s_matrix => matrix_s(1, 1)%matrix
181 892 : p_matrix => matrix_p(1:1, 1)
182 892 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
183 892 : p_matrix => matrix_p(2:2, 1)
184 892 : CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
185 : END IF
186 :
187 : ! calculate energy
188 6648 : DO ikind = 1, nkind
189 4970 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=na)
190 4970 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
191 4970 : CALL get_xtb_atom_param(xtb_kind, defined=defined, natorb=natorb)
192 4970 : IF (.NOT. defined .OR. natorb < 1) CYCLE
193 29820 : ALLOCATE (docg(natorb), wab(natorb, natorb))
194 100898 : wab(1:natorb, 1:natorb) = wabk(1:natorb, 1:natorb, ikind)
195 16480 : DO iatom = 1, na
196 11510 : atom_a = atomic_kind_set(ikind)%atom_list(iatom)
197 11510 : docg = 0.0_dp
198 47692 : docg(1:natorb) = aocg(1:natorb, atom_a) - bocg(1:natorb, atom_a)
199 241916 : espin = 0.5_dp*DOT_PRODUCT(docg, MATMUL(wab, docg))
200 11510 : energy%xtb_spinpol = energy%xtb_spinpol + espin
201 16480 : IF (atprop%energy) THEN
202 0 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + espin
203 : END IF
204 : END DO
205 16588 : DEALLOCATE (docg, wab)
206 : END DO
207 :
208 : ! Forces and Virial
209 1678 : IF (calculate_forces) THEN
210 12 : CALL get_qs_env(qs_env=qs_env, force=force)
211 12 : NULLIFY (cell_to_index)
212 12 : IF (nimg > 1) THEN
213 4 : NULLIFY (kpoints)
214 4 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
215 4 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
216 : END IF
217 12 : IF (nimg == 1) THEN
218 : ! no k-points; all matrices have been transformed to periodic bsf
219 8 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
220 44 : DO WHILE (dbcsr_iterator_blocks_left(iter))
221 36 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
222 36 : ikind = kind_of(irow)
223 36 : atom_i = atom_of_kind(irow)
224 36 : jkind = kind_of(icol)
225 36 : atom_j = atom_of_kind(icol)
226 :
227 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
228 36 : row=irow, col=icol, block=pamat, found=found)
229 36 : CPASSERT(found)
230 : CALL dbcsr_get_block_p(matrix=matrix_p(2, 1)%matrix, &
231 36 : row=irow, col=icol, block=pbmat, found=found)
232 36 : CPASSERT(found)
233 :
234 36 : na = SIZE(pamat, 1)
235 36 : nb = SIZE(pamat, 2)
236 :
237 152 : DO i = 1, 3
238 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
239 108 : row=irow, col=icol, block=dsblock, found=found)
240 108 : 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 108 : aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
245 :
246 108 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
247 252 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
248 : END DO
249 :
250 : END DO
251 8 : CALL dbcsr_iterator_stop(iter)
252 : ! use dsint list
253 8 : IF (use_virial .AND. 0 == 0) THEN
254 4 : CPASSERT(ASSOCIATED(sap_int))
255 16 : DO ikind = 1, nkind
256 52 : DO jkind = 1, nkind
257 36 : iac = ikind + nkind*(jkind - 1)
258 36 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist)) CYCLE
259 48 : DO ia = 1, sap_int(iac)%nalist
260 18 : IF (.NOT. ASSOCIATED(sap_int(iac)%alist(ia)%clist)) CYCLE
261 18 : iatom = sap_int(iac)%alist(ia)%aatom
262 265 : DO ic = 1, sap_int(iac)%alist(ia)%nclist
263 211 : jatom = sap_int(iac)%alist(ia)%clist(ic)%catom
264 844 : rij = sap_int(iac)%alist(ia)%clist(ic)%rac
265 844 : dr = SQRT(SUM(rij(:)**2))
266 229 : IF (dr > 1.e-6_dp) THEN
267 203 : dsint => sap_int(iac)%alist(ia)%clist(ic)%acint
268 203 : icol = MAX(iatom, jatom)
269 203 : irow = MIN(iatom, jatom)
270 : CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, &
271 203 : row=irow, col=icol, block=pamat, found=found)
272 203 : CPASSERT(found)
273 : CALL dbcsr_get_block_p(matrix=matrix_p(2, 1)%matrix, &
274 203 : row=irow, col=icol, block=pbmat, found=found)
275 203 : CPASSERT(found)
276 203 : IF (irow == iatom) THEN
277 119 : na = SIZE(pamat, 1)
278 119 : nb = SIZE(pamat, 2)
279 714 : ALLOCATE (pam(na, nb), pbm(na, nb))
280 1341 : pam(1:na, 1:nb) = pamat(1:na, 1:nb)
281 1341 : pbm(1:na, 1:nb) = pbmat(1:na, 1:nb)
282 : ELSE
283 84 : na = SIZE(pamat, 2)
284 84 : nb = SIZE(pamat, 1)
285 504 : ALLOCATE (pam(na, nb), pbm(na, nb))
286 1070 : pam(1:na, 1:nb) = TRANSPOSE(pamat(1:nb, 1:na))
287 1070 : pbm(1:na, 1:nb) = TRANSPOSE(pbmat(1:nb, 1:na))
288 : END IF
289 :
290 812 : 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 609 : bocg(1:na, iatom), bocg(1:nb, jatom))
295 812 : fij(i) = fi
296 : END DO
297 203 : fi = 1.0_dp
298 203 : IF (iatom == jatom) fi = 0.5_dp
299 203 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
300 203 : 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 4 : NULLIFY (n_list)
310 4 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
311 4 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
312 792 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
313 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
314 788 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
315 :
316 3152 : dr = SQRT(SUM(rij**2))
317 788 : IF (iatom == jatom .AND. dr < 1.0e-6_dp) CYCLE
318 :
319 780 : icol = MAX(iatom, jatom)
320 780 : irow = MIN(iatom, jatom)
321 :
322 780 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
323 780 : CPASSERT(ic > 0)
324 :
325 780 : IF (irow == iatom) THEN
326 462 : iknd = ikind
327 462 : jknd = jkind
328 462 : atom_i = atom_of_kind(iatom)
329 462 : atom_j = atom_of_kind(jatom)
330 : rij = rij
331 : ELSE
332 318 : iknd = jkind
333 318 : jknd = ikind
334 318 : atom_i = atom_of_kind(jatom)
335 318 : atom_j = atom_of_kind(iatom)
336 1272 : rij = -rij
337 : END IF
338 : !
339 : CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
340 780 : row=irow, col=icol, block=pamat, found=found)
341 780 : CPASSERT(found)
342 : CALL dbcsr_get_block_p(matrix=matrix_p(2, ic)%matrix, &
343 780 : row=irow, col=icol, block=pbmat, found=found)
344 780 : CPASSERT(found)
345 :
346 780 : na = SIZE(pamat, 1)
347 780 : nb = SIZE(pamat, 2)
348 :
349 780 : fij = 0.0_dp
350 3120 : DO i = 1, 3
351 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, ic)%matrix, &
352 2340 : row=irow, col=icol, block=dsblock, found=found)
353 2340 : 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 2340 : aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
358 :
359 2340 : force(iknd)%rho_elec(i, atom_i) = force(iknd)%rho_elec(i, atom_i) + fi
360 2340 : force(jknd)%rho_elec(i, atom_j) = force(jknd)%rho_elec(i, atom_j) - fi
361 5460 : fij(i) = fi
362 : END DO
363 784 : IF (use_virial) THEN
364 390 : fi = 1.0_dp
365 390 : IF (iatom == jatom) fi = 0.5_dp
366 390 : CALL virial_pair_force(virial%pv_virial, fi, fij, rij)
367 : END IF
368 :
369 : END DO
370 4 : CALL neighbor_list_iterator_release(nl_iterator)
371 :
372 : END IF
373 : END IF
374 :
375 : ! KS matrix
376 1678 : IF (.NOT. just_energy) THEN
377 1678 : IF (nimg > 1) THEN
378 786 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
379 786 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
380 : END IF
381 1678 : IF (nimg == 1) THEN
382 : ! no k-points; all matrices have been transformed to periodic bsf
383 892 : CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
384 36827 : DO WHILE (dbcsr_iterator_blocks_left(iter))
385 35935 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
386 : CALL dbcsr_get_block_p(matrix=ks_matrix(1, 1)%matrix, &
387 35935 : row=irow, col=icol, block=aksb, found=found)
388 35935 : CPASSERT(found)
389 : CALL dbcsr_get_block_p(matrix=ks_matrix(2, 1)%matrix, &
390 35935 : row=irow, col=icol, block=bksb, found=found)
391 35935 : CPASSERT(found)
392 35935 : na = SIZE(aksb, 1)
393 35935 : nb = SIZE(aksb, 2)
394 35935 : ikind = kind_of(irow)
395 35935 : jkind = kind_of(icol)
396 35935 : 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 36827 : aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
400 : END DO
401 892 : CALL dbcsr_iterator_stop(iter)
402 : ELSE
403 786 : CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
404 786 : CALL neighbor_list_iterator_create(nl_iterator, n_list)
405 155628 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
406 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
407 154842 : iatom=iatom, jatom=jatom, r=rij, cell=cellind)
408 :
409 154842 : icol = MAX(iatom, jatom)
410 154842 : irow = MIN(iatom, jatom)
411 :
412 154842 : ic = cell_to_index(cellind(1), cellind(2), cellind(3))
413 154842 : CPASSERT(ic > 0)
414 :
415 154842 : ikind = kind_of(irow)
416 154842 : jkind = kind_of(icol)
417 :
418 : CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
419 154842 : row=irow, col=icol, block=sblock, found=found)
420 154842 : CPASSERT(found)
421 : CALL dbcsr_get_block_p(matrix=ks_matrix(1, ic)%matrix, &
422 154842 : row=irow, col=icol, block=aksb, found=found)
423 154842 : CPASSERT(found)
424 : CALL dbcsr_get_block_p(matrix=ks_matrix(2, ic)%matrix, &
425 154842 : row=irow, col=icol, block=bksb, found=found)
426 154842 : CPASSERT(found)
427 :
428 154842 : na = SIZE(aksb, 1)
429 154842 : nb = SIZE(aksb, 2)
430 154842 : 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 155628 : aocg(1:na, irow), aocg(1:nb, icol), bocg(1:na, irow), bocg(1:nb, icol))
434 : END DO
435 786 : CALL neighbor_list_iterator_release(nl_iterator)
436 : END IF
437 :
438 : END IF
439 :
440 1678 : DEALLOCATE (wabk)
441 3356 : DEALLOCATE (aocg, bocg)
442 : END IF
443 :
444 1678 : CALL timestop(handle)
445 :
446 3356 : END SUBROUTINE build_xtb_spinpol
447 :
448 : ! **************************************************************************************************
449 : !> \brief ...
450 : !> \param qs_env ...
451 : !> \param ks_matrix ...
452 : !> \param matrix_p1 ...
453 : ! **************************************************************************************************
454 34 : 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 34 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
465 : INTEGER, DIMENSION(25) :: lao
466 : LOGICAL :: found
467 : REAL(KIND=dp) :: fval
468 34 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg1, bocg1
469 34 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
470 : REAL(KIND=dp), DIMENSION(3, 3) :: wall
471 34 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: aksb, bksb, sblock
472 34 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
473 : TYPE(dbcsr_iterator_type) :: iter
474 34 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p_matrix
475 34 : 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 34 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
480 : TYPE(xtb_atom_type), POINTER :: xtb_kind
481 :
482 34 : CALL timeset(routineN, handle)
483 :
484 34 : CALL get_qs_env(qs_env, dft_control=dft_control)
485 34 : nspins = dft_control%nspins
486 34 : nimg = dft_control%nimages
487 :
488 34 : IF (nimg /= 1) THEN
489 0 : CPABORT("No kpoints allowed in xTB response calculation")
490 : END IF
491 :
492 34 : IF (nspins == 2) THEN
493 :
494 : CALL get_qs_env(qs_env, &
495 : qs_kind_set=qs_kind_set, &
496 34 : atomic_kind_set=atomic_kind_set)
497 :
498 34 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
499 :
500 34 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
501 34 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
502 34 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s, para_env=para_env)
503 :
504 : ! expand parameters
505 170 : ALLOCATE (wabk(nsgf, nsgf, nkind))
506 34 : wabk = 0.0_dp
507 120 : DO ikind = 1, nkind
508 86 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
509 86 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall)
510 396 : DO ia = 1, natorb
511 276 : la = lao(ia) + 1
512 1330 : DO ib = 1, natorb
513 968 : lb = lao(ib) + 1
514 1244 : wabk(ia, ib, ikind) = wall(la, lb)
515 : END DO
516 : END DO
517 : END DO
518 :
519 : ! Calculate response charges
520 204 : ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
521 34 : aocg1 = 0.0_dp
522 34 : bocg1 = 0.0_dp
523 34 : s_matrix => matrix_s(1, 1)%matrix
524 34 : p_matrix => matrix_p1(1:1)
525 34 : CALL ao_charges(p_matrix, s_matrix, aocg1, para_env)
526 34 : p_matrix => matrix_p1(2:2)
527 34 : CALL ao_charges(p_matrix, s_matrix, bocg1, para_env)
528 634 : aocg1 = 0.5_dp*aocg1
529 634 : bocg1 = 0.5_dp*bocg1
530 :
531 34 : CALL dbcsr_iterator_start(iter, s_matrix)
532 172 : DO WHILE (dbcsr_iterator_blocks_left(iter))
533 138 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
534 : CALL dbcsr_get_block_p(matrix=ks_matrix(1)%matrix, &
535 138 : row=irow, col=icol, block=aksb, found=found)
536 138 : CPASSERT(found)
537 : CALL dbcsr_get_block_p(matrix=ks_matrix(2)%matrix, &
538 138 : row=irow, col=icol, block=bksb, found=found)
539 138 : CPASSERT(found)
540 138 : na = SIZE(aksb, 1)
541 138 : nb = SIZE(aksb, 2)
542 138 : ikind = kind_of(irow)
543 138 : jkind = kind_of(icol)
544 138 : 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 172 : bocg1(1:na, irow), bocg1(1:nb, icol))
549 : END DO
550 34 : CALL dbcsr_iterator_stop(iter)
551 :
552 34 : DEALLOCATE (wabk)
553 102 : DEALLOCATE (aocg1, bocg1)
554 : END IF
555 :
556 34 : CALL timestop(handle)
557 :
558 68 : END SUBROUTINE xtb_spinpol_hessian
559 :
560 : ! **************************************************************************************************
561 : !> \brief ...
562 : !> \param qs_env ...
563 : !> \param matrix_p0 ...
564 : !> \param matrix_p1 ...
565 : ! **************************************************************************************************
566 2 : 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 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
576 : INTEGER, DIMENSION(25) :: lao
577 : LOGICAL :: found
578 : REAL(KIND=dp) :: fi
579 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aocg, aocg1, bocg, bocg1
580 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: wabk
581 : REAL(KIND=dp), DIMENSION(3, 3) :: wall
582 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dsblock, p0amat, p0bmat, p1amat, p1bmat, &
583 2 : sblock
584 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
585 : TYPE(dbcsr_iterator_type) :: iter
586 2 : 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 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
591 2 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
592 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
593 : TYPE(xtb_atom_type), POINTER :: xtb_kind
594 :
595 2 : CALL timeset(routineN, handle)
596 :
597 2 : CALL get_qs_env(qs_env, dft_control=dft_control)
598 2 : nspins = dft_control%nspins
599 2 : nimg = dft_control%nimages
600 2 : IF (nimg /= 1) THEN
601 0 : CPABORT("xTB response forces for spin polarisation Hamiltonian not available")
602 : END IF
603 :
604 2 : 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 2 : 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 2 : atom_of_kind=atom_of_kind)
614 :
615 2 : CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
616 2 : CALL get_qs_kind_set(qs_kind_set, maxsgf=nsgf)
617 :
618 2 : CALL get_qs_env(qs_env, matrix_s=matrix_s, para_env=para_env)
619 :
620 : ! expand parameters
621 10 : ALLOCATE (wabk(nsgf, nsgf, nkind))
622 2 : wabk = 0.0_dp
623 6 : DO ikind = 1, nkind
624 4 : CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
625 4 : CALL get_xtb_atom_param(xtb_kind, natorb=natorb, lao=lao, wall=wall)
626 18 : DO ia = 1, natorb
627 12 : la = lao(ia) + 1
628 56 : DO ib = 1, natorb
629 40 : lb = lao(ib) + 1
630 52 : wabk(ia, ib, ikind) = wall(la, lb)
631 : END DO
632 : END DO
633 : END DO
634 :
635 : ! Calculate charges
636 2 : s_matrix => matrix_s(1)%matrix
637 12 : ALLOCATE (aocg(nsgf, natom), bocg(nsgf, natom))
638 2 : aocg = 0.0_dp
639 2 : bocg = 0.0_dp
640 2 : p_matrix => matrix_p0(1:1)
641 2 : CALL ao_charges(p_matrix, s_matrix, aocg, para_env)
642 2 : p_matrix => matrix_p0(2:2)
643 2 : CALL ao_charges(p_matrix, s_matrix, bocg, para_env)
644 : ! Calculate response charges
645 12 : ALLOCATE (aocg1(nsgf, natom), bocg1(nsgf, natom))
646 2 : aocg1 = 0.0_dp
647 2 : bocg1 = 0.0_dp
648 2 : p_matrix => matrix_p1(1:1)
649 2 : CALL ao_charges(p_matrix, s_matrix, aocg1, para_env)
650 2 : p_matrix => matrix_p1(2:2)
651 2 : CALL ao_charges(p_matrix, s_matrix, bocg1, para_env)
652 32 : aocg1 = 0.5_dp*aocg1
653 32 : bocg1 = 0.5_dp*bocg1
654 :
655 : ! calculate forces
656 2 : CALL get_qs_env(qs_env=qs_env, force=force)
657 : ! no k-points; all matrices have been transformed to periodic bsf
658 2 : CALL dbcsr_iterator_start(iter, s_matrix)
659 8 : DO WHILE (dbcsr_iterator_blocks_left(iter))
660 6 : CALL dbcsr_iterator_next_block(iter, irow, icol, sblock)
661 6 : ikind = kind_of(irow)
662 6 : atom_i = atom_of_kind(irow)
663 6 : jkind = kind_of(icol)
664 6 : atom_j = atom_of_kind(icol)
665 :
666 : CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
667 6 : row=irow, col=icol, block=p0amat, found=found)
668 6 : CPASSERT(found)
669 : CALL dbcsr_get_block_p(matrix=matrix_p0(2)%matrix, &
670 6 : row=irow, col=icol, block=p0bmat, found=found)
671 6 : CPASSERT(found)
672 : CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
673 6 : row=irow, col=icol, block=p1amat, found=found)
674 6 : CPASSERT(found)
675 : CALL dbcsr_get_block_p(matrix=matrix_p1(2)%matrix, &
676 6 : row=irow, col=icol, block=p1bmat, found=found)
677 6 : CPASSERT(found)
678 :
679 6 : na = SIZE(p0amat, 1)
680 6 : nb = SIZE(p0amat, 2)
681 :
682 26 : DO i = 1, 3
683 : CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
684 18 : row=irow, col=icol, block=dsblock, found=found)
685 18 : 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 18 : bocg1(1:na, irow), bocg1(1:nb, icol))
694 :
695 18 : force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
696 42 : force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
697 : END DO
698 :
699 : END DO
700 6 : CALL dbcsr_iterator_stop(iter)
701 :
702 : END IF
703 :
704 2 : CALL timestop(handle)
705 :
706 4 : 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 190915 : SUBROUTINE ksupdate(aksb, bksb, sb, na, nb, fval, &
724 190915 : 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 381830 : REAL(KIND=dp), DIMENSION(na) :: dqa, wa
735 381830 : REAL(KIND=dp), DIMENSION(na, nb) :: wqab
736 190915 : REAL(KIND=dp), DIMENSION(nb) :: dqb, wb
737 :
738 807168 : dqa = qai - qbi
739 657124 : dqb = qaj - qbj
740 3889639 : wa = MATMUL(wabi, dqa)
741 2589187 : wb = MATMUL(wabj, dqb)
742 657124 : DO ib = 1, nb
743 2230435 : DO ia = 1, na
744 2039520 : wqab(ia, ib) = fval*sb(ia, ib)*(wa(ia) + wb(ib))
745 : END DO
746 : END DO
747 :
748 2230435 : aksb = aksb + wqab
749 2230435 : bksb = bksb - wqab
750 :
751 190915 : 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 3057 : SUBROUTINE fupdate(fij, pa, pb, ds, na, nb, &
769 3057 : 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 6114 : REAL(KIND=dp), DIMENSION(na) :: dpsa, dqa, wa
779 6114 : REAL(KIND=dp), DIMENSION(na, nb) :: dpab
780 3057 : REAL(KIND=dp), DIMENSION(nb) :: dpsb, dqb, wb
781 :
782 12429 : dqa = qai - qbi
783 10545 : dqb = qaj - qbj
784 56634 : wa = MATMUL(wabi, dqa)
785 41562 : wb = MATMUL(wabj, dqb)
786 34269 : dpab = pa - pb
787 12429 : dpsa = 0.0_dp
788 10545 : dpsb = 0.0_dp
789 10545 : DO ib = 1, nb
790 34269 : DO ia = 1, na
791 23724 : dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
792 31212 : dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
793 : END DO
794 : END DO
795 :
796 22974 : fij = SUM(wa*dpsa) + SUM(wb*dpsb)
797 :
798 3057 : 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 36 : SUBROUTINE f2update(fij, p0a, p0b, p1a, p1b, ds, na, nb, wabi, wabj, &
822 36 : 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 36 : REAL(KIND=dp), DIMENSION(na) :: dpsa, dqa, wa
832 36 : REAL(KIND=dp), DIMENSION(na, nb) :: dpab
833 18 : REAL(KIND=dp), DIMENSION(nb) :: dpsb, dqb, wb
834 :
835 : fij = 0.0_dp
836 :
837 72 : dqa = qai - qbi
838 60 : dqb = qaj - qbj
839 324 : wa = MATMUL(wabi, dqa)
840 228 : wb = MATMUL(wabj, dqb)
841 192 : dpab = p1a - p1b
842 72 : dpsa = 0.0_dp
843 60 : dpsb = 0.0_dp
844 60 : DO ib = 1, nb
845 192 : DO ia = 1, na
846 132 : dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
847 174 : dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
848 : END DO
849 : END DO
850 :
851 132 : fij = fij + SUM(wa*dpsa) + SUM(wb*dpsb)
852 :
853 72 : dqa = rai - rbi
854 60 : dqb = raj - rbj
855 324 : wa = MATMUL(wabi, dqa)
856 228 : wb = MATMUL(wabj, dqb)
857 192 : dpab = p0a - p0b
858 72 : dpsa = 0.0_dp
859 60 : dpsb = 0.0_dp
860 60 : DO ib = 1, nb
861 192 : DO ia = 1, na
862 132 : dpsa(ia) = dpsa(ia) + dpab(ia, ib)*ds(ia, ib)
863 174 : dpsb(ib) = dpsb(ib) + dpab(ia, ib)*ds(ia, ib)
864 : END DO
865 : END DO
866 :
867 132 : fij = fij + SUM(wa*dpsa) + SUM(wb*dpsb)
868 :
869 18 : END SUBROUTINE f2update
870 :
871 11510 : END MODULE xtb_spinpol
|
