(git:d43dca4)
Loading...
Searching...
No Matches
qs_moments.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 Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
10!> \par History
11!> added angular moments (JGH 11.2012)
12!> \author JGH (20.07.2006)
13! **************************************************************************************************
14MODULE qs_moments
15 USE ai_angmom, ONLY: angmom
16 USE ai_moments, ONLY: contract_cossin, &
17 cossin, &
18 diff_momop, &
19 diff_momop2, &
20 diff_momop_velocity, &
21 moment
22 USE atomic_kind_types, ONLY: atomic_kind_type, &
23 get_atomic_kind
24 USE basis_set_types, ONLY: gto_basis_set_p_type, &
25 gto_basis_set_type
26 USE bibliography, ONLY: mattiat2019, &
27 cite_reference
28 USE block_p_types, ONLY: block_p_type
29 USE cell_types, ONLY: cell_type, &
30 pbc, &
31 get_cell
32 USE commutator_rpnl, ONLY: build_com_mom_nl
33 USE cp_cfm_basic_linalg, ONLY: cp_cfm_det
34 USE cp_cfm_types, ONLY: cp_cfm_create, &
35 cp_cfm_get_info, &
36 cp_cfm_release, &
37 cp_cfm_type
38 USE cp_control_types, ONLY: dft_control_type
39 USE cp_dbcsr_api, ONLY: &
40 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, &
41 dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_type, &
42 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
43 USE cp_dbcsr_contrib, ONLY: dbcsr_dot, &
44 dbcsr_trace
45 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
46 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply, &
47 dbcsr_allocate_matrix_set, &
48 dbcsr_deallocate_matrix_set
49 USE cp_fm_struct, ONLY: cp_fm_struct_create, &
50 cp_fm_struct_double, &
51 cp_fm_struct_release, &
52 cp_fm_struct_type
53 USE cp_fm_types, ONLY: cp_fm_create, &
54 cp_fm_get_info, &
55 cp_fm_release, &
56 cp_fm_set_all, &
57 cp_fm_type
58 USE cp_result_methods, ONLY: cp_results_erase, &
59 put_results
60 USE cp_result_types, ONLY: cp_result_type
61 USE distribution_1d_types, ONLY: distribution_1d_type
62 USE kinds, ONLY: default_string_length, &
63 dp, &
64 max_line_length
65 USE kpoint_k_r_trafo_simple, ONLY: replicate_rs_matrices, &
66 rs_to_kp
67 USE kpoint_methods, ONLY: kpoint_init_cell_index, &
68 kpoint_initialize
69 USE kpoint_types, ONLY: get_kpoint_info, &
70 kpoint_type, &
71 kpoint_create, &
72 kpoint_release, &
73 read_kpoint_section
74 USE mathconstants, ONLY: pi, &
75 twopi, &
76 gaussi, &
77 z_zero
78 USE message_passing, ONLY: mp_para_env_type
79 USE moments_utils, ONLY: get_reference_point
80 USE orbital_pointers, ONLY: current_maxl, &
81 indco, &
82 ncoset
83 USE parallel_gemm_api, ONLY: parallel_gemm
84 USE particle_methods, ONLY: get_particle_set
85 USE particle_types, ONLY: particle_type
86 USE physcon, ONLY: bohr, &
87 debye, &
88 angstrom
89 USE qs_environment_types, ONLY: get_qs_env, &
90 qs_environment_type
91 USE qs_kind_types, ONLY: get_qs_kind, &
92 get_qs_kind_set, &
93 qs_kind_type
94 USE qs_ks_types, ONLY: get_ks_env, &
95 qs_ks_env_type
96 USE qs_mo_types, ONLY: get_mo_set, &
97 mo_set_type
98 USE qs_neighbor_list_types, ONLY: get_iterator_info, &
99 neighbor_list_iterate, &
100 neighbor_list_iterator_create, &
101 neighbor_list_iterator_p_type, &
102 neighbor_list_iterator_release, &
103 neighbor_list_set_p_type
104 USE qs_operators_ao, ONLY: build_lin_mom_matrix
105 USE qs_rho_types, ONLY: qs_rho_get, &
106 qs_rho_type
107 USE rt_propagation_types, ONLY: get_rtp, &
108 rt_prop_type
109 USE cp_parser_methods, ONLY: read_float_object
110 USE input_section_types, ONLY: section_vals_get, &
111 section_vals_get_subs_vals, &
112 section_vals_type, &
113 section_vals_val_get
114 USE mathlib, ONLY: geeig_right, &
115 gemm_square
116 USE string_utilities, ONLY: uppercase
117
118#include "./base/base_uses.f90"
119
120 IMPLICIT NONE
121
122 PRIVATE
123
124 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
125
126 ! Public subroutines
127 PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
128 PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
129 PUBLIC :: qs_moment_berry_phase, qs_moment_locop
130 PUBLIC :: qs_moment_kpoints, build_local_moment_matrix_rs_img
131 PUBLIC :: qs_moment_kpoints_deep
132 PUBLIC :: dipole_deriv_ao
133 PUBLIC :: build_local_moments_der_matrix
134 PUBLIC :: build_dsdv_moments
135 PUBLIC :: dipole_velocity_deriv
136 PUBLIC :: calculate_commutator_nl_terms, op_orbbas, op_orbbas_rtp, &
137 print_moments, print_moments_nl, set_label
138
139CONTAINS
140
141! *****************************************************************************
142!> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
143!> to the basis function on the right
144!> difdip(beta, alpha) = < mu | r_beta | ∂_alpha nu > * (mu - nu)
145!> \param qs_env ...
146!> \param difdip ...
147!> \param order The order of the derivative (1 for dipole moment)
148!> \param lambda The atom on which we take the derivative
149!> \param rc ...
150!> \author Edward Ditler
151! **************************************************************************************************
152 SUBROUTINE dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
153 TYPE(qs_environment_type), POINTER :: qs_env
154 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: difdip
155 INTEGER, INTENT(IN) :: order, lambda
156 REAL(kind=dp), DIMENSION(3) :: rc
157
158 CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_velocity_deriv'
159
160 INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
161 last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
162 sgfb
163 LOGICAL :: found
164 REAL(dp) :: dab
165 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
166 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
167 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
168 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab, difmab2
169 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
170 TYPE(cell_type), POINTER :: cell
171 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
172 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
173 TYPE(neighbor_list_iterator_p_type), &
174 DIMENSION(:), POINTER :: nl_iterator
175 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
176 POINTER :: sab_all
177 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
178 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
179 TYPE(qs_kind_type), POINTER :: qs_kind
180
181 CALL timeset(routinen, handle)
182
183 NULLIFY (cell, particle_set, qs_kind_set, sab_all)
184 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
185 qs_kind_set=qs_kind_set, sab_all=sab_all)
186 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
187 maxco=ldab, maxsgf=maxsgf)
188
189 nkind = SIZE(qs_kind_set)
190 natom = SIZE(particle_set)
191
192 m_dim = ncoset(order) - 1
193
194 ALLOCATE (basis_set_list(nkind))
195
196 ALLOCATE (mab(ldab, ldab, m_dim))
197 ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
198 ALLOCATE (work(ldab, maxsgf))
199 ALLOCATE (mint(3, 3))
200 ALLOCATE (mint2(3, 3))
201
202 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
203 difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
204 work(1:ldab, 1:maxsgf) = 0.0_dp
205
206 DO i = 1, 3
207 DO j = 1, 3
208 NULLIFY (mint(i, j)%block)
209 NULLIFY (mint2(i, j)%block)
210 END DO
211 END DO
212
213 ! Set the basis_set_list(nkind) to point to the corresponding basis sets
214 DO ikind = 1, nkind
215 qs_kind => qs_kind_set(ikind)
216 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
217 IF (ASSOCIATED(basis_set_a)) THEN
218 basis_set_list(ikind)%gto_basis_set => basis_set_a
219 ELSE
220 NULLIFY (basis_set_list(ikind)%gto_basis_set)
221 END IF
222 END DO
223
224 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
225 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
226 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
227 iatom=iatom, jatom=jatom, r=rab)
228
229 basis_set_a => basis_set_list(ikind)%gto_basis_set
230 basis_set_b => basis_set_list(jkind)%gto_basis_set
231 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
232 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
233
234 associate( &
235 ! basis ikind
236 first_sgfa => basis_set_a%first_sgf, &
237 la_max => basis_set_a%lmax, &
238 la_min => basis_set_a%lmin, &
239 npgfa => basis_set_a%npgf, &
240 nsgfa => basis_set_a%nsgf_set, &
241 rpgfa => basis_set_a%pgf_radius, &
242 set_radius_a => basis_set_a%set_radius, &
243 sphi_a => basis_set_a%sphi, &
244 zeta => basis_set_a%zet, &
245 ! basis jkind, &
246 first_sgfb => basis_set_b%first_sgf, &
247 lb_max => basis_set_b%lmax, &
248 lb_min => basis_set_b%lmin, &
249 npgfb => basis_set_b%npgf, &
250 nsgfb => basis_set_b%nsgf_set, &
251 rpgfb => basis_set_b%pgf_radius, &
252 set_radius_b => basis_set_b%set_radius, &
253 sphi_b => basis_set_b%sphi, &
254 zetb => basis_set_b%zet)
255
256 nseta = basis_set_a%nset
257 nsetb = basis_set_b%nset
258
259 IF (inode == 1) last_jatom = 0
260
261 ! this guarantees minimum image convention
262 ! anything else would not make sense
263 IF (jatom == last_jatom) THEN
264 cycle
265 END IF
266
267 last_jatom = jatom
268
269 irow = iatom
270 icol = jatom
271
272 DO i = 1, 3
273 DO j = 1, 3
274 NULLIFY (mint(i, j)%block)
275 CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
276 row=irow, col=icol, block=mint(i, j)%block, &
277 found=found)
278 cpassert(found)
279 mint(i, j)%block = 0._dp
280 END DO
281 END DO
282
283 ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
284 ra = pbc(particle_set(iatom)%r(:), cell)
285 rb(:) = ra(:) + rab(:)
286 rac = pbc(rc, ra, cell)
287 rbc = rac + rab
288 dab = norm2(rab)
289
290 DO iset = 1, nseta
291 ncoa = npgfa(iset)*ncoset(la_max(iset))
292 sgfa = first_sgfa(1, iset)
293 DO jset = 1, nsetb
294 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
295 ncob = npgfb(jset)*ncoset(lb_max(jset))
296 sgfb = first_sgfb(1, jset)
297 ldab = max(ncoa, ncob)
298 lda = ncoset(la_max(iset))*npgfa(iset)
299 ldb = ncoset(lb_max(jset))*npgfb(jset)
300 ALLOCATE (difmab(lda, ldb, m_dim, 3))
301
302 ! Calculate integral difmab(beta, alpha) = (a|r_beta|db_alpha)
303 ! difmab(beta, alpha) = < a | r_beta | ∂_alpha b >
304 ! difmab(j, idir) = < a | r_j | ∂_idir b >
305 CALL diff_momop_velocity(la_max(iset), npgfa(iset), zeta(:, iset), &
306 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
307 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
308 difmab, lambda=lambda, iatom=iatom, jatom=jatom)
309
310 ! *** Contraction step ***
311
312 DO idir = 1, 3 ! derivative of AO function
313 DO j = 1, 3 ! position operator r_j
314 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
315 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
316 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
317 0.0_dp, work(1, 1), SIZE(work, 1))
318
319 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
320 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
321 work(1, 1), SIZE(work, 1), &
322 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
323 SIZE(mint(j, idir)%block, 1))
324 END DO !j
325 END DO !idir
326 DEALLOCATE (difmab)
327 END DO !jset
328 END DO !iset
329 END associate
330 END do!iterator
331
332 CALL neighbor_list_iterator_release(nl_iterator)
333
334 DO i = 1, 3
335 DO j = 1, 3
336 NULLIFY (mint(i, j)%block)
337 END DO
338 END DO
339
340 DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
341
342 CALL timestop(handle)
343 END SUBROUTINE dipole_velocity_deriv
344
345! **************************************************************************************************
346!> \brief Builds the moments for the derivative of the overlap with respect to nuclear velocities
347!> \param qs_env ...
348!> \param moments ...
349!> \param nmoments ...
350!> \param ref_point ...
351!> \param ref_points ...
352!> \param basis_type ...
353!> \author Edward Ditler
354! **************************************************************************************************
355 SUBROUTINE build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
356
357 TYPE(qs_environment_type), POINTER :: qs_env
358 TYPE(dbcsr_p_type), DIMENSION(:) :: moments
359 INTEGER, INTENT(IN) :: nmoments
360 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
361 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
362 OPTIONAL :: ref_points
363 CHARACTER(len=*), OPTIONAL :: basis_type
364
365 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_dsdv_moments'
366
367 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
368 maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
369 INTEGER, DIMENSION(3) :: image_cell
370 LOGICAL :: found
371 REAL(kind=dp) :: dab
372 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
373 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
374 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
375 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
376 TYPE(cell_type), POINTER :: cell
377 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
378 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
379 TYPE(neighbor_list_iterator_p_type), &
380 DIMENSION(:), POINTER :: nl_iterator
381 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
382 POINTER :: sab_orb
383 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
384 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
385 TYPE(qs_kind_type), POINTER :: qs_kind
386
387 IF (nmoments < 1) RETURN
388
389 CALL timeset(routinen, handle)
390
391 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
392
393 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
394 cpassert(SIZE(moments) >= nm)
395
396 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
397 CALL get_qs_env(qs_env=qs_env, &
398 qs_kind_set=qs_kind_set, &
399 particle_set=particle_set, cell=cell, &
400 sab_orb=sab_orb)
401
402 nkind = SIZE(qs_kind_set)
403 natom = SIZE(particle_set)
404
405 ! Allocate work storage
406 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
407 maxco=maxco, maxsgf=maxsgf, &
408 basis_type=basis_type)
409
410 ALLOCATE (mab(maxco, maxco, nm))
411 mab(:, :, :) = 0.0_dp
412
413 ALLOCATE (work(maxco, maxsgf))
414 work(:, :) = 0.0_dp
415
416 ALLOCATE (mint(nm))
417 DO i = 1, nm
418 NULLIFY (mint(i)%block)
419 END DO
420
421 ALLOCATE (basis_set_list(nkind))
422 DO ikind = 1, nkind
423 qs_kind => qs_kind_set(ikind)
424 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
425 IF (ASSOCIATED(basis_set_a)) THEN
426 basis_set_list(ikind)%gto_basis_set => basis_set_a
427 ELSE
428 NULLIFY (basis_set_list(ikind)%gto_basis_set)
429 END IF
430 END DO
431 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
432 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
433 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
434 iatom=iatom, jatom=jatom, r=rab, cell=image_cell)
435 basis_set_a => basis_set_list(ikind)%gto_basis_set
436 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
437 basis_set_b => basis_set_list(jkind)%gto_basis_set
438 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
439 ! basis ikind
440 associate( &
441 first_sgfa => basis_set_a%first_sgf, &
442 la_max => basis_set_a%lmax, &
443 la_min => basis_set_a%lmin, &
444 npgfa => basis_set_a%npgf, &
445 nsgfa => basis_set_a%nsgf_set, &
446 rpgfa => basis_set_a%pgf_radius, &
447 set_radius_a => basis_set_a%set_radius, &
448 sphi_a => basis_set_a%sphi, &
449 zeta => basis_set_a%zet, &
450 ! basis jkind, &
451 first_sgfb => basis_set_b%first_sgf, &
452 lb_max => basis_set_b%lmax, &
453 lb_min => basis_set_b%lmin, &
454 npgfb => basis_set_b%npgf, &
455 nsgfb => basis_set_b%nsgf_set, &
456 rpgfb => basis_set_b%pgf_radius, &
457 set_radius_b => basis_set_b%set_radius, &
458 sphi_b => basis_set_b%sphi, &
459 zetb => basis_set_b%zet)
460
461 nseta = basis_set_a%nset
462 nsetb = basis_set_b%nset
463
464 IF (inode == 1) last_jatom = 0
465
466 ! this guarantees minimum image convention
467 ! anything else would not make sense
468 IF (jatom == last_jatom) THEN
469 cycle
470 END IF
471
472 last_jatom = jatom
473
474 IF (iatom <= jatom) THEN
475 irow = iatom
476 icol = jatom
477 ELSE
478 irow = jatom
479 icol = iatom
480 END IF
481
482 DO i = 1, nm
483 NULLIFY (mint(i)%block)
484 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
485 row=irow, col=icol, block=mint(i)%block, found=found)
486 mint(i)%block = 0._dp
487 END DO
488
489 ! fold atomic position back into unit cell
490 IF (PRESENT(ref_points)) THEN
491 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
492 ELSE IF (PRESENT(ref_point)) THEN
493 rc(:) = ref_point(:)
494 ELSE
495 rc(:) = 0._dp
496 END IF
497 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
498 ! by folding around the center, such screwing can be avoided for a proper choice of center.
499 ! we dont use PBC at this point
500
501 ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
502 ra(:) = particle_set(iatom)%r(:)
503 rb(:) = ra(:) + rab(:)
504 rac = pbc(rc, ra, cell)
505 rbc = rac + rab
506
507 dab = norm2(rab)
508
509 DO iset = 1, nseta
510
511 ncoa = npgfa(iset)*ncoset(la_max(iset))
512 sgfa = first_sgfa(1, iset)
513
514 DO jset = 1, nsetb
515
516 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
517
518 ncob = npgfb(jset)*ncoset(lb_max(jset))
519 sgfb = first_sgfb(1, jset)
520
521 ! Calculate the primitive integrals
522 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
523 rpgfa(:, iset), la_min(iset), &
524 lb_max(jset), npgfb(jset), zetb(:, jset), &
525 rpgfb(:, jset), nmoments, rac, rbc, mab)
526
527 ! Contraction step
528 DO i = 1, nm
529
530 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
531 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
532 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
533 0.0_dp, work(1, 1), SIZE(work, 1))
534
535 IF (iatom <= jatom) THEN
536
537 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
538 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
539 work(1, 1), SIZE(work, 1), &
540 1.0_dp, mint(i)%block(sgfa, sgfb), &
541 SIZE(mint(i)%block, 1))
542
543 ELSE
544
545 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
546 1.0_dp, work(1, 1), SIZE(work, 1), &
547 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
548 1.0_dp, mint(i)%block(sgfb, sgfa), &
549 SIZE(mint(i)%block, 1))
550
551 END IF
552
553 END DO
554
555 END DO
556 END DO
557 END associate
558
559 END DO ! iterator
560
561 CALL neighbor_list_iterator_release(nl_iterator)
562
563 ! Release work storage
564 DEALLOCATE (mab, basis_set_list)
565 DEALLOCATE (work)
566 DO i = 1, nm
567 NULLIFY (mint(i)%block)
568 END DO
569 DEALLOCATE (mint)
570
571 CALL timestop(handle)
572
573 END SUBROUTINE build_dsdv_moments
574
575! **************************************************************************************************
576!> \brief ...
577!> \param qs_env ...
578!> \param moments ...
579!> \param nmoments ...
580!> \param ref_point ...
581!> \param ref_points ...
582!> \param basis_type ...
583! **************************************************************************************************
584 SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
585
586 TYPE(qs_environment_type), POINTER :: qs_env
587 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
588 INTEGER, INTENT(IN) :: nmoments
589 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
590 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
591 OPTIONAL :: ref_points
592 CHARACTER(len=*), OPTIONAL :: basis_type
593
594 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix'
595
596 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
597 maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
598 LOGICAL :: found
599 REAL(kind=dp) :: dab
600 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
601 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
602 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
603 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
604 TYPE(cell_type), POINTER :: cell
605 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
606 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
607 TYPE(neighbor_list_iterator_p_type), &
608 DIMENSION(:), POINTER :: nl_iterator
609 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
610 POINTER :: sab_orb
611 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
612 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
613 TYPE(qs_kind_type), POINTER :: qs_kind
614
615 IF (nmoments < 1) RETURN
616
617 CALL timeset(routinen, handle)
618
619 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
620
621 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
622 cpassert(SIZE(moments) >= nm)
623
624 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
625 CALL get_qs_env(qs_env=qs_env, &
626 qs_kind_set=qs_kind_set, &
627 particle_set=particle_set, cell=cell, &
628 sab_orb=sab_orb)
629
630 nkind = SIZE(qs_kind_set)
631 natom = SIZE(particle_set)
632
633 ! Allocate work storage
634 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
635 maxco=maxco, maxsgf=maxsgf, &
636 basis_type=basis_type)
637
638 ALLOCATE (mab(maxco, maxco, nm))
639 mab(:, :, :) = 0.0_dp
640
641 ALLOCATE (work(maxco, maxsgf))
642 work(:, :) = 0.0_dp
643
644 ALLOCATE (mint(nm))
645 DO i = 1, nm
646 NULLIFY (mint(i)%block)
647 END DO
648
649 ALLOCATE (basis_set_list(nkind))
650 DO ikind = 1, nkind
651 qs_kind => qs_kind_set(ikind)
652 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
653 IF (ASSOCIATED(basis_set_a)) THEN
654 basis_set_list(ikind)%gto_basis_set => basis_set_a
655 ELSE
656 NULLIFY (basis_set_list(ikind)%gto_basis_set)
657 END IF
658 END DO
659 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
660 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
661 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
662 iatom=iatom, jatom=jatom, r=rab)
663 basis_set_a => basis_set_list(ikind)%gto_basis_set
664 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
665 basis_set_b => basis_set_list(jkind)%gto_basis_set
666 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
667 associate( &
668 ! basis ikind
669 first_sgfa => basis_set_a%first_sgf, &
670 la_max => basis_set_a%lmax, &
671 la_min => basis_set_a%lmin, &
672 npgfa => basis_set_a%npgf, &
673 nsgfa => basis_set_a%nsgf_set, &
674 rpgfa => basis_set_a%pgf_radius, &
675 set_radius_a => basis_set_a%set_radius, &
676 sphi_a => basis_set_a%sphi, &
677 zeta => basis_set_a%zet, &
678 ! basis jkind, &
679 first_sgfb => basis_set_b%first_sgf, &
680 lb_max => basis_set_b%lmax, &
681 lb_min => basis_set_b%lmin, &
682 npgfb => basis_set_b%npgf, &
683 nsgfb => basis_set_b%nsgf_set, &
684 rpgfb => basis_set_b%pgf_radius, &
685 set_radius_b => basis_set_b%set_radius, &
686 sphi_b => basis_set_b%sphi, &
687 zetb => basis_set_b%zet)
688
689 nseta = basis_set_a%nset
690 nsetb = basis_set_b%nset
691
692 IF (inode == 1) last_jatom = 0
693
694 ! this guarantees minimum image convention
695 ! anything else would not make sense
696 IF (jatom == last_jatom) THEN
697 cycle
698 END IF
699
700 last_jatom = jatom
701
702 IF (iatom <= jatom) THEN
703 irow = iatom
704 icol = jatom
705 ELSE
706 irow = jatom
707 icol = iatom
708 END IF
709
710 DO i = 1, nm
711 NULLIFY (mint(i)%block)
712 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
713 row=irow, col=icol, block=mint(i)%block, found=found)
714 mint(i)%block = 0._dp
715 END DO
716
717 ! fold atomic position back into unit cell
718 IF (PRESENT(ref_points)) THEN
719 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
720 ELSE IF (PRESENT(ref_point)) THEN
721 rc(:) = ref_point(:)
722 ELSE
723 rc(:) = 0._dp
724 END IF
725 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
726 ! by folding around the center, such screwing can be avoided for a proper choice of center.
727 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
728 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
729 ! we dont use PBC at this point
730 rab(:) = ra(:) - rb(:)
731 rac(:) = ra(:) - rc(:)
732 rbc(:) = rb(:) - rc(:)
733 dab = norm2(rab)
734
735 DO iset = 1, nseta
736
737 ncoa = npgfa(iset)*ncoset(la_max(iset))
738 sgfa = first_sgfa(1, iset)
739
740 DO jset = 1, nsetb
741
742 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
743
744 ncob = npgfb(jset)*ncoset(lb_max(jset))
745 sgfb = first_sgfb(1, jset)
746
747 ! Calculate the primitive integrals
748 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
749 rpgfa(:, iset), la_min(iset), &
750 lb_max(jset), npgfb(jset), zetb(:, jset), &
751 rpgfb(:, jset), nmoments, rac, rbc, mab)
752
753 ! Contraction step
754 DO i = 1, nm
755
756 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
757 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
758 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
759 0.0_dp, work(1, 1), SIZE(work, 1))
760
761 IF (iatom <= jatom) THEN
762
763 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
764 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
765 work(1, 1), SIZE(work, 1), &
766 1.0_dp, mint(i)%block(sgfa, sgfb), &
767 SIZE(mint(i)%block, 1))
768
769 ELSE
770
771 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
772 1.0_dp, work(1, 1), SIZE(work, 1), &
773 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
774 1.0_dp, mint(i)%block(sgfb, sgfa), &
775 SIZE(mint(i)%block, 1))
776
777 END IF
778
779 END DO
780
781 END DO
782 END DO
783 END associate
784
785 END DO
786 CALL neighbor_list_iterator_release(nl_iterator)
787
788 ! Release work storage
789 DEALLOCATE (mab, basis_set_list)
790 DEALLOCATE (work)
791 DO i = 1, nm
792 NULLIFY (mint(i)%block)
793 END DO
794 DEALLOCATE (mint)
795
796 CALL timestop(handle)
797
798 END SUBROUTINE build_local_moment_matrix
799
800! **************************************************************************************************
801!> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
802!> Optionally stores the multipole moments themselves for free.
803!> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
804!> Only first derivatives are performed, e. g. x d/dy
805!> \param qs_env ...
806!> \param moments_der will contain the derivatives of the multipole moments
807!> \param nmoments_der order of the moments with derivatives
808!> \param nmoments order of the multipole moments (no derivatives, same output as
809!> build_local_moment_matrix, needs moments as arguments to store results)
810!> \param ref_point ...
811!> \param moments contains the multipole moments, optionally for free, up to order nmoments
812!> \note
813!> Adapted from rRc_xyz_der_ao in qs_operators_ao
814! **************************************************************************************************
815 SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
816 ref_point, moments)
817 TYPE(qs_environment_type), POINTER :: qs_env
818 TYPE(dbcsr_p_type), DIMENSION(:, :), &
819 INTENT(INOUT), POINTER :: moments_der
820 INTEGER, INTENT(IN) :: nmoments_der, nmoments
821 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
822 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
823 OPTIONAL, POINTER :: moments
824
825 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moments_der_matrix'
826
827 INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
828 jatom, jkind, jpgf, jset, last_jatom, m_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
829 nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
830 LOGICAL :: found
831 REAL(kind=dp) :: dab
832 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
833 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
834 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
835 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
836 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
837 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
838 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
839 TYPE(cell_type), POINTER :: cell
840 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
841 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
842 TYPE(neighbor_list_iterator_p_type), &
843 DIMENSION(:), POINTER :: nl_iterator
844 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
845 POINTER :: sab_orb
846 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
847 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
848 TYPE(qs_kind_type), POINTER :: qs_kind
849
850 nmom_build = max(nmoments, nmoments_der) ! build moments up to order nmom_buiod
851 IF (nmom_build < 1) RETURN
852
853 CALL timeset(routinen, handle)
854
855 nders = 1 ! only first order derivatives
856 dimders = ncoset(nders) - 1
857
858 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
859 CALL get_qs_env(qs_env=qs_env, &
860 qs_kind_set=qs_kind_set, &
861 particle_set=particle_set, &
862 cell=cell, &
863 sab_orb=sab_orb)
864
865 nkind = SIZE(qs_kind_set)
866
867 ! Work storage
868 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
869 maxco=maxco, maxsgf=maxsgf)
870
871 IF (nmoments > 0) THEN
872 cpassert(PRESENT(moments))
873 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
874 cpassert(SIZE(moments) == nm)
875 ! storage for integrals
876 ALLOCATE (mab(maxco, maxco, nm))
877 ! blocks
878 mab(:, :, :) = 0.0_dp
879 ALLOCATE (mom_block(nm))
880 DO i = 1, nm
881 NULLIFY (mom_block(i)%block)
882 END DO
883 END IF
884
885 IF (nmoments_der > 0) THEN
886 m_dim = ncoset(nmoments_der) - 1
887 cpassert(SIZE(moments_der, dim=1) == m_dim)
888 cpassert(SIZE(moments_der, dim=2) == dimders)
889 ! storage for integrals
890 ALLOCATE (difmab(maxco, maxco, m_dim, dimders))
891 difmab(:, :, :, :) = 0.0_dp
892 ! blocks
893 ALLOCATE (mom_block_der(m_dim, dimders))
894 DO i = 1, m_dim
895 DO ider = 1, dimders
896 NULLIFY (mom_block_der(i, ider)%block)
897 END DO
898 END DO
899 END IF
900
901 ALLOCATE (work(maxco, maxsgf))
902 work(:, :) = 0.0_dp
903
904 NULLIFY (basis_set_a, basis_set_b, basis_set_list)
905 NULLIFY (qs_kind)
906 ALLOCATE (basis_set_list(nkind))
907 DO ikind = 1, nkind
908 qs_kind => qs_kind_set(ikind)
909 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
910 IF (ASSOCIATED(basis_set_a)) THEN
911 basis_set_list(ikind)%gto_basis_set => basis_set_a
912 ELSE
913 NULLIFY (basis_set_list(ikind)%gto_basis_set)
914 END IF
915 END DO
916
917 ! Calculate derivatives looping over neighbour list
918 NULLIFY (nl_iterator)
919 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
920 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
921 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
922 iatom=iatom, jatom=jatom, r=rab)
923 basis_set_a => basis_set_list(ikind)%gto_basis_set
924 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
925 basis_set_b => basis_set_list(jkind)%gto_basis_set
926 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
927 associate( &
928 ! basis ikind
929 first_sgfa => basis_set_a%first_sgf, &
930 la_max => basis_set_a%lmax, &
931 la_min => basis_set_a%lmin, &
932 npgfa => basis_set_a%npgf, &
933 nsgfa => basis_set_a%nsgf_set, &
934 rpgfa => basis_set_a%pgf_radius, &
935 set_radius_a => basis_set_a%set_radius, &
936 sphi_a => basis_set_a%sphi, &
937 zeta => basis_set_a%zet, &
938 ! basis jkind, &
939 first_sgfb => basis_set_b%first_sgf, &
940 lb_max => basis_set_b%lmax, &
941 lb_min => basis_set_b%lmin, &
942 npgfb => basis_set_b%npgf, &
943 nsgfb => basis_set_b%nsgf_set, &
944 rpgfb => basis_set_b%pgf_radius, &
945 set_radius_b => basis_set_b%set_radius, &
946 sphi_b => basis_set_b%sphi, &
947 zetb => basis_set_b%zet)
948
949 nseta = basis_set_a%nset
950 nsetb = basis_set_b%nset
951
952 ! reference point
953 IF (PRESENT(ref_point)) THEN
954 rc(:) = ref_point(:)
955 ELSE
956 rc(:) = 0._dp
957 END IF
958 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
959 ! by folding around the center, such screwing can be avoided for a proper choice of center.
960 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
961 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
962 ! we dont use PBC at this point
963 rab(:) = ra(:) - rb(:)
964 rac(:) = ra(:) - rc(:)
965 rbc(:) = rb(:) - rc(:)
966 dab = norm2(rab)
967
968 ! get blocks
969 IF (inode == 1) last_jatom = 0
970
971 IF (jatom == last_jatom) THEN
972 cycle
973 END IF
974
975 last_jatom = jatom
976
977 IF (iatom <= jatom) THEN
978 irow = iatom
979 icol = jatom
980 ELSE
981 irow = jatom
982 icol = iatom
983 END IF
984
985 IF (nmoments > 0) THEN
986 DO i = 1, nm
987 NULLIFY (mom_block(i)%block)
988 ! get block from pre calculated overlap matrix
989 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
990 row=irow, col=icol, block=mom_block(i)%block, found=found)
991 cpassert(found .AND. ASSOCIATED(mom_block(i)%block))
992 mom_block(i)%block = 0._dp
993 END DO
994 END IF
995 IF (nmoments_der > 0) THEN
996 DO i = 1, m_dim
997 DO ider = 1, dimders
998 NULLIFY (mom_block_der(i, ider)%block)
999 CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
1000 row=irow, col=icol, &
1001 block=mom_block_der(i, ider)%block, &
1002 found=found)
1003 cpassert(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
1004 mom_block_der(i, ider)%block = 0._dp
1005 END DO
1006 END DO
1007 END IF
1008
1009 DO iset = 1, nseta
1010
1011 ncoa = npgfa(iset)*ncoset(la_max(iset))
1012 sgfa = first_sgfa(1, iset)
1013
1014 DO jset = 1, nsetb
1015
1016 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1017
1018 ncob = npgfb(jset)*ncoset(lb_max(jset))
1019 sgfb = first_sgfb(1, jset)
1020
1021 NULLIFY (mab_tmp)
1022 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
1023 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
1024
1025 ! Calculate the primitive integrals (need l+1 for derivatives)
1026 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1027 rpgfa(:, iset), la_min(iset), &
1028 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
1029 rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
1030
1031 IF (nmoments_der > 0) THEN
1032 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1033 rpgfa(:, iset), la_min(iset), &
1034 lb_max(jset), npgfb(jset), zetb(:, jset), &
1035 rpgfb(:, jset), lb_min(jset), &
1036 nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp)
1037 END IF
1038
1039 IF (nmoments > 0) THEN
1040 ! copy subset of mab_tmp (l+1) to mab (l)
1041 mab = 0.0_dp
1042 DO ii = 1, nm
1043 na = 0
1044 nda = 0
1045 DO ipgf = 1, npgfa(iset)
1046 nb = 0
1047 ndb = 0
1048 DO jpgf = 1, npgfb(jset)
1049 DO j = 1, ncoset(lb_max(jset))
1050 DO i = 1, ncoset(la_max(iset))
1051 mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
1052 END DO ! i
1053 END DO ! j
1054 nb = nb + ncoset(lb_max(jset))
1055 ndb = ndb + ncoset(lb_max(jset) + 1)
1056 END DO ! jpgf
1057 na = na + ncoset(la_max(iset))
1058 nda = nda + ncoset(la_max(iset) + 1)
1059 END DO ! ipgf
1060 END DO
1061 ! Contraction step
1062 DO i = 1, nm
1063
1064 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1065 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1066 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1067 0.0_dp, work(1, 1), SIZE(work, 1))
1068
1069 IF (iatom <= jatom) THEN
1070 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1071 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1072 work(1, 1), SIZE(work, 1), &
1073 1.0_dp, mom_block(i)%block(sgfa, sgfb), &
1074 SIZE(mom_block(i)%block, 1))
1075 ELSE
1076 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1077 1.0_dp, work(1, 1), SIZE(work, 1), &
1078 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1079 1.0_dp, mom_block(i)%block(sgfb, sgfa), &
1080 SIZE(mom_block(i)%block, 1))
1081 END IF
1082 END DO
1083 END IF
1084
1085 IF (nmoments_der > 0) THEN
1086 DO i = 1, m_dim
1087 DO ider = 1, dimders
1088 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1089 1.0_dp, difmab(1, 1, i, ider), SIZE(difmab, 1), &
1090 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1091 0._dp, work(1, 1), SIZE(work, 1))
1092
1093 IF (iatom <= jatom) THEN
1094 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1095 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1096 work(1, 1), SIZE(work, 1), &
1097 1._dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
1098 SIZE(mom_block_der(i, ider)%block, 1))
1099 ELSE
1100 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1101 -1.0_dp, work(1, 1), SIZE(work, 1), &
1102 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1103 1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
1104 SIZE(mom_block_der(i, ider)%block, 1))
1105 END IF
1106 END DO
1107 END DO
1108 END IF
1109 DEALLOCATE (mab_tmp)
1110 END DO
1111 END DO
1112 END associate
1113 END DO
1114 CALL neighbor_list_iterator_release(nl_iterator)
1115
1116 ! deallocations
1117 DEALLOCATE (basis_set_list)
1118 DEALLOCATE (work)
1119 IF (nmoments > 0) THEN
1120 DEALLOCATE (mab)
1121 DO i = 1, nm
1122 NULLIFY (mom_block(i)%block)
1123 END DO
1124 DEALLOCATE (mom_block)
1125 END IF
1126 IF (nmoments_der > 0) THEN
1127 DEALLOCATE (difmab)
1128 DO i = 1, m_dim
1129 DO ider = 1, dimders
1130 NULLIFY (mom_block_der(i, ider)%block)
1131 END DO
1132 END DO
1133 DEALLOCATE (mom_block_der)
1134 END IF
1135
1136 CALL timestop(handle)
1137
1138 END SUBROUTINE build_local_moments_der_matrix
1139
1140! **************************************************************************************************
1141!> \brief ...
1142!> \param qs_env ...
1143!> \param magmom ...
1144!> \param nmoments ...
1145!> \param ref_point ...
1146!> \param ref_points ...
1147!> \param basis_type ...
1148! **************************************************************************************************
1149 SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
1150
1151 TYPE(qs_environment_type), POINTER :: qs_env
1152 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
1153 INTEGER, INTENT(IN) :: nmoments
1154 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
1155 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
1156 OPTIONAL :: ref_points
1157 CHARACTER(len=*), OPTIONAL :: basis_type
1158
1159 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_magmom_matrix'
1160
1161 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
1162 maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
1163 LOGICAL :: found
1164 REAL(kind=dp) :: dab
1165 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1166 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
1167 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
1168 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
1169 TYPE(cell_type), POINTER :: cell
1170 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1171 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1172 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1173 TYPE(neighbor_list_iterator_p_type), &
1174 DIMENSION(:), POINTER :: nl_iterator
1175 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1176 POINTER :: sab_orb
1177 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1178 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1179 TYPE(qs_kind_type), POINTER :: qs_kind
1180
1181 IF (nmoments < 1) RETURN
1182
1183 CALL timeset(routinen, handle)
1184
1185 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
1186 NULLIFY (matrix_s)
1187
1188 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1189
1190 ! magnetic dipoles/angular moments only
1191 nm = 3
1192
1193 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1194 CALL get_qs_env(qs_env=qs_env, &
1195 qs_kind_set=qs_kind_set, &
1196 particle_set=particle_set, cell=cell, &
1197 sab_orb=sab_orb)
1198
1199 nkind = SIZE(qs_kind_set)
1200 natom = SIZE(particle_set)
1201
1202 ! Allocate work storage
1203 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1204 maxco=maxco, maxsgf=maxsgf)
1205
1206 ALLOCATE (mab(maxco, maxco, nm))
1207 mab(:, :, :) = 0.0_dp
1208
1209 ALLOCATE (work(maxco, maxsgf))
1210 work(:, :) = 0.0_dp
1211
1212 ALLOCATE (mint(nm))
1213 DO i = 1, nm
1214 NULLIFY (mint(i)%block)
1215 END DO
1216
1217 ALLOCATE (basis_set_list(nkind))
1218 DO ikind = 1, nkind
1219 qs_kind => qs_kind_set(ikind)
1220 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1221 IF (ASSOCIATED(basis_set_a)) THEN
1222 basis_set_list(ikind)%gto_basis_set => basis_set_a
1223 ELSE
1224 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1225 END IF
1226 END DO
1227 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1228 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1229 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1230 iatom=iatom, jatom=jatom, r=rab)
1231 basis_set_a => basis_set_list(ikind)%gto_basis_set
1232 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1233 basis_set_b => basis_set_list(jkind)%gto_basis_set
1234 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1235 associate( &
1236 ! basis ikind
1237 first_sgfa => basis_set_a%first_sgf, &
1238 la_max => basis_set_a%lmax, &
1239 la_min => basis_set_a%lmin, &
1240 npgfa => basis_set_a%npgf, &
1241 nsgfa => basis_set_a%nsgf_set, &
1242 rpgfa => basis_set_a%pgf_radius, &
1243 set_radius_a => basis_set_a%set_radius, &
1244 sphi_a => basis_set_a%sphi, &
1245 zeta => basis_set_a%zet, &
1246 ! basis jkind, &
1247 first_sgfb => basis_set_b%first_sgf, &
1248 lb_max => basis_set_b%lmax, &
1249 lb_min => basis_set_b%lmin, &
1250 npgfb => basis_set_b%npgf, &
1251 nsgfb => basis_set_b%nsgf_set, &
1252 rpgfb => basis_set_b%pgf_radius, &
1253 set_radius_b => basis_set_b%set_radius, &
1254 sphi_b => basis_set_b%sphi, &
1255 zetb => basis_set_b%zet)
1256
1257 nseta = basis_set_a%nset
1258 nsetb = basis_set_b%nset
1259
1260 IF (iatom <= jatom) THEN
1261 irow = iatom
1262 icol = jatom
1263 ELSE
1264 irow = jatom
1265 icol = iatom
1266 END IF
1267
1268 DO i = 1, nm
1269 NULLIFY (mint(i)%block)
1270 CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
1271 row=irow, col=icol, block=mint(i)%block, found=found)
1272 mint(i)%block = 0._dp
1273 cpassert(ASSOCIATED(mint(i)%block))
1274 END DO
1275
1276 ! fold atomic position back into unit cell
1277 IF (PRESENT(ref_points)) THEN
1278 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
1279 ELSE IF (PRESENT(ref_point)) THEN
1280 rc(:) = ref_point(:)
1281 ELSE
1282 rc(:) = 0._dp
1283 END IF
1284 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
1285 ! by folding around the center, such screwing can be avoided for a proper choice of center.
1286 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
1287 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
1288 ! we dont use PBC at this point
1289 rab(:) = ra(:) - rb(:)
1290 rac(:) = ra(:) - rc(:)
1291 rbc(:) = rb(:) - rc(:)
1292 dab = norm2(rab)
1293
1294 DO iset = 1, nseta
1295
1296 ncoa = npgfa(iset)*ncoset(la_max(iset))
1297 sgfa = first_sgfa(1, iset)
1298
1299 DO jset = 1, nsetb
1300
1301 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1302
1303 ncob = npgfb(jset)*ncoset(lb_max(jset))
1304 sgfb = first_sgfb(1, jset)
1305
1306 ! Calculate the primitive integrals
1307 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
1308 rpgfa(:, iset), la_min(iset), &
1309 lb_max(jset), npgfb(jset), zetb(:, jset), &
1310 rpgfb(:, jset), rac, rbc, mab)
1311
1312 ! Contraction step
1313 DO i = 1, nm
1314 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1315 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1316 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1317 0.0_dp, work(1, 1), SIZE(work, 1))
1318
1319 IF (iatom <= jatom) THEN
1320 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1321 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1322 work(1, 1), SIZE(work, 1), &
1323 1.0_dp, mint(i)%block(sgfa, sgfb), &
1324 SIZE(mint(i)%block, 1))
1325 ELSE
1326 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1327 -1.0_dp, work(1, 1), SIZE(work, 1), &
1328 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1329 1.0_dp, mint(i)%block(sgfb, sgfa), &
1330 SIZE(mint(i)%block, 1))
1331 END IF
1332
1333 END DO
1334
1335 END DO
1336 END DO
1337 END associate
1338 END DO
1339 CALL neighbor_list_iterator_release(nl_iterator)
1340
1341 ! Release work storage
1342 DEALLOCATE (mab, basis_set_list)
1343 DEALLOCATE (work)
1344 DO i = 1, nm
1345 NULLIFY (mint(i)%block)
1346 END DO
1347 DEALLOCATE (mint)
1348
1349 CALL timestop(handle)
1350
1351 END SUBROUTINE build_local_magmom_matrix
1352
1353! **************************************************************************************************
1354!> \brief ...
1355!> \param qs_env ...
1356!> \param cosmat ...
1357!> \param sinmat ...
1358!> \param kvec ...
1359!> \param sab_orb_external ...
1360!> \param basis_type ...
1361! **************************************************************************************************
1362 SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1363
1364 TYPE(qs_environment_type), POINTER :: qs_env
1365 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1366 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1367 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1368 OPTIONAL, POINTER :: sab_orb_external
1369 CHARACTER(len=*), OPTIONAL :: basis_type
1370
1371 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_moment_matrix'
1372
1373 INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, ldsa, &
1374 ldsb, ldwork, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
1375 LOGICAL :: found
1376 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1377 REAL(kind=dp) :: dab
1378 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1379 TYPE(cell_type), POINTER :: cell
1380 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1381 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1382 TYPE(neighbor_list_iterator_p_type), &
1383 DIMENSION(:), POINTER :: nl_iterator
1384 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1385 POINTER :: sab_orb
1386 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1387 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1388 TYPE(qs_kind_type), POINTER :: qs_kind
1389
1390 CALL timeset(routinen, handle)
1391
1392 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1393 CALL get_qs_env(qs_env=qs_env, &
1394 qs_kind_set=qs_kind_set, &
1395 particle_set=particle_set, cell=cell, &
1396 sab_orb=sab_orb)
1397
1398 IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
1399
1400 CALL dbcsr_set(sinmat, 0.0_dp)
1401 CALL dbcsr_set(cosmat, 0.0_dp)
1402
1403 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
1404 ldab = ldwork
1405 ALLOCATE (cosab(ldab, ldab))
1406 ALLOCATE (sinab(ldab, ldab))
1407 ALLOCATE (work(ldwork, ldwork))
1408
1409 nkind = SIZE(qs_kind_set)
1410 natom = SIZE(particle_set)
1411
1412 ALLOCATE (basis_set_list(nkind))
1413 DO ikind = 1, nkind
1414 qs_kind => qs_kind_set(ikind)
1415 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1416 IF (ASSOCIATED(basis_set_a)) THEN
1417 basis_set_list(ikind)%gto_basis_set => basis_set_a
1418 ELSE
1419 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1420 END IF
1421 END DO
1422 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1423 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1424 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1425 iatom=iatom, jatom=jatom, r=rab)
1426 basis_set_a => basis_set_list(ikind)%gto_basis_set
1427 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1428 basis_set_b => basis_set_list(jkind)%gto_basis_set
1429 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1430 associate( &
1431 ! basis ikind
1432 first_sgfa => basis_set_a%first_sgf, &
1433 la_max => basis_set_a%lmax, &
1434 la_min => basis_set_a%lmin, &
1435 npgfa => basis_set_a%npgf, &
1436 nsgfa => basis_set_a%nsgf_set, &
1437 rpgfa => basis_set_a%pgf_radius, &
1438 set_radius_a => basis_set_a%set_radius, &
1439 sphi_a => basis_set_a%sphi, &
1440 zeta => basis_set_a%zet, &
1441 ! basis jkind, &
1442 first_sgfb => basis_set_b%first_sgf, &
1443 lb_max => basis_set_b%lmax, &
1444 lb_min => basis_set_b%lmin, &
1445 npgfb => basis_set_b%npgf, &
1446 nsgfb => basis_set_b%nsgf_set, &
1447 rpgfb => basis_set_b%pgf_radius, &
1448 set_radius_b => basis_set_b%set_radius, &
1449 sphi_b => basis_set_b%sphi, &
1450 zetb => basis_set_b%zet)
1451
1452 nseta = basis_set_a%nset
1453 nsetb = basis_set_b%nset
1454
1455 ldsa = SIZE(sphi_a, 1)
1456 ldsb = SIZE(sphi_b, 1)
1457
1458 IF (iatom <= jatom) THEN
1459 irow = iatom
1460 icol = jatom
1461 ELSE
1462 irow = jatom
1463 icol = iatom
1464 END IF
1465
1466 NULLIFY (cblock)
1467 CALL dbcsr_get_block_p(matrix=cosmat, &
1468 row=irow, col=icol, block=cblock, found=found)
1469 NULLIFY (sblock)
1470 CALL dbcsr_get_block_p(matrix=sinmat, &
1471 row=irow, col=icol, block=sblock, found=found)
1472 IF (ASSOCIATED(cblock) .AND. .NOT. ASSOCIATED(sblock) .OR. &
1473 .NOT. ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1474 cpabort("")
1475 END IF
1476
1477 IF (ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1478
1479 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1480 rb(:) = ra + rab
1481 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1482
1483 DO iset = 1, nseta
1484
1485 ncoa = npgfa(iset)*ncoset(la_max(iset))
1486 sgfa = first_sgfa(1, iset)
1487
1488 DO jset = 1, nsetb
1489
1490 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1491
1492 ncob = npgfb(jset)*ncoset(lb_max(jset))
1493 sgfb = first_sgfb(1, jset)
1494
1495 ! Calculate the primitive integrals
1496 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1497 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1498 ra, rb, kvec, cosab, sinab)
1499 CALL contract_cossin(cblock, sblock, &
1500 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1501 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1502 cosab, sinab, ldab, work, ldwork)
1503
1504 END DO
1505 END DO
1506
1507 END IF
1508 END associate
1509 END DO
1510 CALL neighbor_list_iterator_release(nl_iterator)
1511
1512 DEALLOCATE (cosab)
1513 DEALLOCATE (sinab)
1514 DEALLOCATE (work)
1515 DEALLOCATE (basis_set_list)
1516
1517 CALL timestop(handle)
1518
1519 END SUBROUTINE build_berry_moment_matrix
1520
1521! **************************************************************************************************
1522!> \brief ...
1523!> \param qs_env ...
1524!> \param cosmat ...
1525!> \param sinmat ...
1526!> \param kvec ...
1527!> \param basis_type ...
1528! **************************************************************************************************
1529 SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
1530
1531 TYPE(qs_environment_type), POINTER :: qs_env
1532 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1533 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1534 CHARACTER(len=*), OPTIONAL :: basis_type
1535
1536 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_kpoint_matrix'
1537
1538 INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1539 ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1540 INTEGER, DIMENSION(3) :: icell
1541 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1542 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1543 LOGICAL :: found, use_cell_mapping
1544 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1545 REAL(kind=dp) :: dab
1546 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1547 TYPE(cell_type), POINTER :: cell
1548 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1549 TYPE(dft_control_type), POINTER :: dft_control
1550 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1551 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1552 TYPE(kpoint_type), POINTER :: kpoints
1553 TYPE(neighbor_list_iterator_p_type), &
1554 DIMENSION(:), POINTER :: nl_iterator
1555 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1556 POINTER :: sab_orb
1557 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1558 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1559 TYPE(qs_kind_type), POINTER :: qs_kind
1560 TYPE(qs_ks_env_type), POINTER :: ks_env
1561
1562 CALL timeset(routinen, handle)
1563
1564 CALL get_qs_env(qs_env, &
1565 ks_env=ks_env, &
1566 dft_control=dft_control)
1567 nimg = dft_control%nimages
1568 IF (nimg > 1) THEN
1569 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1570 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1571 use_cell_mapping = .true.
1572 ELSE
1573 use_cell_mapping = .false.
1574 END IF
1575
1576 CALL get_qs_env(qs_env=qs_env, &
1577 qs_kind_set=qs_kind_set, &
1578 particle_set=particle_set, cell=cell, &
1579 sab_orb=sab_orb)
1580
1581 nkind = SIZE(qs_kind_set)
1582 natom = SIZE(particle_set)
1583 ALLOCATE (basis_set_list(nkind))
1584 DO ikind = 1, nkind
1585 qs_kind => qs_kind_set(ikind)
1586 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1587 IF (ASSOCIATED(basis_set)) THEN
1588 basis_set_list(ikind)%gto_basis_set => basis_set
1589 ELSE
1590 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1591 END IF
1592 END DO
1593
1594 ALLOCATE (row_blk_sizes(natom))
1595 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1596 basis=basis_set_list)
1597 CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1598 ! (re)allocate matrix sets
1599 CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1600 CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1601 DO i = 1, nimg
1602 ! sin
1603 ALLOCATE (sinmat(1, i)%matrix)
1604 CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1605 name="SINMAT", &
1606 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1607 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1608 CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1609 CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1610 ! cos
1611 ALLOCATE (cosmat(1, i)%matrix)
1612 CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1613 name="COSMAT", &
1614 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1615 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1616 CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1617 CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1618 END DO
1619
1620 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1621 ldab = ldwork
1622 ALLOCATE (cosab(ldab, ldab))
1623 ALLOCATE (sinab(ldab, ldab))
1624 ALLOCATE (work(ldwork, ldwork))
1625
1626 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1627 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1628 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1629 iatom=iatom, jatom=jatom, r=rab, cell=icell)
1630 basis_set_a => basis_set_list(ikind)%gto_basis_set
1631 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1632 basis_set_b => basis_set_list(jkind)%gto_basis_set
1633 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1634 associate( &
1635 ! basis ikind
1636 first_sgfa => basis_set_a%first_sgf, &
1637 la_max => basis_set_a%lmax, &
1638 la_min => basis_set_a%lmin, &
1639 npgfa => basis_set_a%npgf, &
1640 nsgfa => basis_set_a%nsgf_set, &
1641 rpgfa => basis_set_a%pgf_radius, &
1642 set_radius_a => basis_set_a%set_radius, &
1643 sphi_a => basis_set_a%sphi, &
1644 zeta => basis_set_a%zet, &
1645 ! basis jkind, &
1646 first_sgfb => basis_set_b%first_sgf, &
1647 lb_max => basis_set_b%lmax, &
1648 lb_min => basis_set_b%lmin, &
1649 npgfb => basis_set_b%npgf, &
1650 nsgfb => basis_set_b%nsgf_set, &
1651 rpgfb => basis_set_b%pgf_radius, &
1652 set_radius_b => basis_set_b%set_radius, &
1653 sphi_b => basis_set_b%sphi, &
1654 zetb => basis_set_b%zet)
1655
1656 nseta = basis_set_a%nset
1657 nsetb = basis_set_b%nset
1658
1659 ldsa = SIZE(sphi_a, 1)
1660 ldsb = SIZE(sphi_b, 1)
1661
1662 IF (iatom <= jatom) THEN
1663 irow = iatom
1664 icol = jatom
1665 ELSE
1666 irow = jatom
1667 icol = iatom
1668 END IF
1669
1670 IF (use_cell_mapping) THEN
1671 ic = cell_to_index(icell(1), icell(2), icell(3))
1672 cpassert(ic > 0)
1673 ELSE
1674 ic = 1
1675 END IF
1676
1677 NULLIFY (sblock)
1678 CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1679 row=irow, col=icol, block=sblock, found=found)
1680 cpassert(found)
1681 NULLIFY (cblock)
1682 CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1683 row=irow, col=icol, block=cblock, found=found)
1684 cpassert(found)
1685
1686 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1687 rb(:) = ra + rab
1688 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1689
1690 DO iset = 1, nseta
1691
1692 ncoa = npgfa(iset)*ncoset(la_max(iset))
1693 sgfa = first_sgfa(1, iset)
1694
1695 DO jset = 1, nsetb
1696
1697 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1698
1699 ncob = npgfb(jset)*ncoset(lb_max(jset))
1700 sgfb = first_sgfb(1, jset)
1701
1702 ! Calculate the primitive integrals
1703 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1704 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1705 ra, rb, kvec, cosab, sinab)
1706 CALL contract_cossin(cblock, sblock, &
1707 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1708 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1709 cosab, sinab, ldab, work, ldwork)
1710
1711 END DO
1712 END DO
1713 END associate
1714 END DO
1715 CALL neighbor_list_iterator_release(nl_iterator)
1716
1717 DEALLOCATE (cosab)
1718 DEALLOCATE (sinab)
1719 DEALLOCATE (work)
1720 DEALLOCATE (basis_set_list)
1721 DEALLOCATE (row_blk_sizes)
1722
1723 CALL timestop(handle)
1724
1725 END SUBROUTINE build_berry_kpoint_matrix
1726
1727! **************************************************************************************************
1728!> \brief ...
1729!> \param qs_env ...
1730!> \param magnetic ...
1731!> \param nmoments ...
1732!> \param reference ...
1733!> \param ref_point ...
1734!> \param unit_number ...
1735! **************************************************************************************************
1736 SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
1737
1738 TYPE(qs_environment_type), POINTER :: qs_env
1739 LOGICAL, INTENT(IN) :: magnetic
1740 INTEGER, INTENT(IN) :: nmoments, reference
1741 REAL(dp), DIMENSION(:), POINTER :: ref_point
1742 INTEGER, INTENT(IN) :: unit_number
1743
1744 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_berry_phase'
1745
1746 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1747 CHARACTER(LEN=default_string_length) :: description
1748 COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1749 zij(3, 3), zijk(3, 3, 3), &
1750 zijkl(3, 3, 3, 3), zphase(3), zz
1751 INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1752 iy, iz, j, k, l, nao, nm, nmo, nmom, &
1753 nmotot, tmp_dim
1754 LOGICAL :: floating, ghost, uniform
1755 REAL(dp) :: charge, ci(3), cij(3, 3), dd, occ, trace
1756 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1757 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1758 REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1759 TYPE(atomic_kind_type), POINTER :: atomic_kind
1760 TYPE(cell_type), POINTER :: cell
1761 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1762 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1763 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1764 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1765 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1766 TYPE(cp_fm_type), POINTER :: mo_coeff
1767 TYPE(cp_result_type), POINTER :: results
1768 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1769 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1770 TYPE(dft_control_type), POINTER :: dft_control
1771 TYPE(distribution_1d_type), POINTER :: local_particles
1772 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1773 TYPE(mp_para_env_type), POINTER :: para_env
1774 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1775 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1776 TYPE(qs_rho_type), POINTER :: rho
1777 TYPE(rt_prop_type), POINTER :: rtp
1778
1779 cpassert(ASSOCIATED(qs_env))
1780
1781 IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1782 IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1783 RETURN
1784 END IF
1785
1786 CALL timeset(routinen, handle)
1787
1788 ! restrict maximum moment available
1789 nmom = min(nmoments, 2)
1790
1791 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1792 ! rmom(:,1)=electronic
1793 ! rmom(:,2)=nuclear
1794 ! rmom(:,1)=total
1795 ALLOCATE (rmom(nm + 1, 3))
1796 ALLOCATE (rlab(nm + 1))
1797 rmom = 0.0_dp
1798 rlab = ""
1799 IF (magnetic) THEN
1800 nm = 3
1801 ALLOCATE (mmom(nm))
1802 mmom = 0._dp
1803 END IF
1804
1805 NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1806 local_particles, matrix_s, mos, rho_ao)
1807
1808 CALL get_qs_env(qs_env, &
1809 dft_control=dft_control, &
1810 rho=rho, &
1811 cell=cell, &
1812 results=results, &
1813 particle_set=particle_set, &
1814 qs_kind_set=qs_kind_set, &
1815 para_env=para_env, &
1816 local_particles=local_particles, &
1817 matrix_s=matrix_s, &
1818 mos=mos)
1819
1820 CALL qs_rho_get(rho, rho_ao=rho_ao)
1821
1822 NULLIFY (cosmat, sinmat)
1823 ALLOCATE (cosmat, sinmat)
1824 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1825 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1826 CALL dbcsr_set(cosmat, 0.0_dp)
1827 CALL dbcsr_set(sinmat, 0.0_dp)
1828
1829 ALLOCATE (op_fm_set(2, dft_control%nspins))
1830 ALLOCATE (opvec(dft_control%nspins))
1831 ALLOCATE (eigrmat(dft_control%nspins))
1832 nmotot = 0
1833 DO ispin = 1, dft_control%nspins
1834 NULLIFY (tmp_fm_struct, mo_coeff)
1835 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1836 nmotot = nmotot + nmo
1837 CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1838 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1839 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1840 DO i = 1, SIZE(op_fm_set, 1)
1841 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1842 END DO
1843 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1844 CALL cp_fm_struct_release(tmp_fm_struct)
1845 END DO
1846
1847 ! occupation
1848 DO ispin = 1, dft_control%nspins
1849 CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
1850 IF (.NOT. uniform) THEN
1851 cpwarn("Berry phase moments for non uniform MOs' occupation numbers not implemented")
1852 END IF
1853 END DO
1854
1855 ! reference point
1856 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1857 rcc = pbc(rcc, cell)
1858
1859 ! label
1860 DO l = 1, nm
1861 ix = indco(1, l + 1)
1862 iy = indco(2, l + 1)
1863 iz = indco(3, l + 1)
1864 CALL set_label(rlab(l + 1), ix, iy, iz)
1865 END DO
1866
1867 ! nuclear contribution
1868 DO ia = 1, SIZE(particle_set)
1869 atomic_kind => particle_set(ia)%atomic_kind
1870 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1871 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1872 IF (.NOT. ghost .AND. .NOT. floating) THEN
1873 rmom(1, 2) = rmom(1, 2) - charge
1874 END IF
1875 END DO
1876 ria = twopi*matmul(cell%h_inv, rcc)
1877 zphase = cmplx(cos(ria), sin(ria), dp)**rmom(1, 2)
1878
1879 zi = 0._dp
1880 zij = 0._dp
1881 zijk = 0._dp
1882 zijkl = 0._dp
1883
1884 DO l = 1, nmom
1885 SELECT CASE (l)
1886 CASE (1)
1887 ! Dipole
1888 zi(:) = cmplx(1._dp, 0._dp, dp)
1889 DO ia = 1, SIZE(particle_set)
1890 atomic_kind => particle_set(ia)%atomic_kind
1891 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1892 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1893 IF (.NOT. ghost .AND. .NOT. floating) THEN
1894 ria = particle_set(ia)%r
1895 ria = pbc(ria, cell)
1896 DO i = 1, 3
1897 kvec(:) = twopi*cell%h_inv(i, :)
1898 dd = sum(kvec(:)*ria(:))
1899 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1900 zi(i) = zi(i)*zdeta
1901 END DO
1902 END IF
1903 END DO
1904 zi = zi*zphase
1905 ci = aimag(log(zi))/twopi
1906 qq = aimag(log(zi))
1907 rmom(2:4, 2) = matmul(cell%hmat, ci)
1908 CASE (2)
1909 ! Quadrupole
1910 cpabort("Berry phase moments bigger than 1 not implemented")
1911 zij(:, :) = cmplx(1._dp, 0._dp, dp)
1912 DO ia = 1, SIZE(particle_set)
1913 atomic_kind => particle_set(ia)%atomic_kind
1914 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1915 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge)
1916 ria = particle_set(ia)%r
1917 ria = pbc(ria, cell)
1918 DO i = 1, 3
1919 DO j = i, 3
1920 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1921 dd = sum(kvec(:)*ria(:))
1922 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1923 zij(i, j) = zij(i, j)*zdeta
1924 zij(j, i) = zij(i, j)
1925 END DO
1926 END DO
1927 END DO
1928 DO i = 1, 3
1929 DO j = 1, 3
1930 zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1931 zz = zij(i, j)/zi(i)/zi(j)
1932 cij(i, j) = aimag(log(zz))/twopi
1933 END DO
1934 END DO
1935 cij = 0.5_dp*cij/twopi/twopi
1936 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1937 DO k = 4, 9
1938 ix = indco(1, k + 1)
1939 iy = indco(2, k + 1)
1940 iz = indco(3, k + 1)
1941 IF (ix == 0) THEN
1942 rmom(k + 1, 2) = cij(iy, iz)
1943 ELSE IF (iy == 0) THEN
1944 rmom(k + 1, 2) = cij(ix, iz)
1945 ELSE IF (iz == 0) THEN
1946 rmom(k + 1, 2) = cij(ix, iy)
1947 END IF
1948 END DO
1949 CASE (3)
1950 ! Octapole
1951 cpabort("Berry phase moments bigger than 2 not implemented")
1952 CASE (4)
1953 ! Hexadecapole
1954 cpabort("Berry phase moments bigger than 3 not implemented")
1955 CASE DEFAULT
1956 cpabort("Berry phase moments bigger than 4 not implemented")
1957 END SELECT
1958 END DO
1959
1960 ! electronic contribution
1961
1962 ria = twopi*real(nmotot, dp)*occ*matmul(cell%h_inv, rcc)
1963 xphase = cmplx(cos(ria), sin(ria), dp)
1964
1965 ! charge
1966 trace = 0.0_dp
1967 DO ispin = 1, dft_control%nspins
1968 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1969 rmom(1, 1) = rmom(1, 1) + trace
1970 END DO
1971
1972 zi = 0._dp
1973 zij = 0._dp
1974 zijk = 0._dp
1975 zijkl = 0._dp
1976
1977 DO l = 1, nmom
1978 SELECT CASE (l)
1979 CASE (1)
1980 ! Dipole
1981 DO i = 1, 3
1982 kvec(:) = twopi*cell%h_inv(i, :)
1983 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1984 IF (qs_env%run_rtp) THEN
1985 CALL get_qs_env(qs_env, rtp=rtp)
1986 CALL get_rtp(rtp, mos_new=mos_new)
1987 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1988 ELSE
1989 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1990 END IF
1991 zdet = cmplx(1._dp, 0._dp, dp)
1992 DO ispin = 1, dft_control%nspins
1993 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1994 DO idim = 1, tmp_dim
1995 eigrmat(ispin)%local_data(:, idim) = &
1996 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1997 -op_fm_set(2, ispin)%local_data(:, idim), dp)
1998 END DO
1999 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
2000 CALL cp_cfm_det(eigrmat(ispin), zdeta)
2001 zdet = zdet*zdeta
2002 IF (dft_control%nspins == 1) THEN
2003 zdet = zdet*zdeta
2004 END IF
2005 END DO
2006 zi(i) = zdet
2007 END DO
2008 zi = zi*xphase
2009 CASE (2)
2010 ! Quadrupole
2011 cpabort("Berry phase moments bigger than 1 not implemented")
2012 DO i = 1, 3
2013 DO j = i, 3
2014 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
2015 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
2016 IF (qs_env%run_rtp) THEN
2017 CALL get_qs_env(qs_env, rtp=rtp)
2018 CALL get_rtp(rtp, mos_new=mos_new)
2019 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
2020 ELSE
2021 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
2022 END IF
2023 zdet = cmplx(1._dp, 0._dp, dp)
2024 DO ispin = 1, dft_control%nspins
2025 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
2026 DO idim = 1, tmp_dim
2027 eigrmat(ispin)%local_data(:, idim) = &
2028 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
2029 -op_fm_set(2, ispin)%local_data(:, idim), dp)
2030 END DO
2031 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
2032 CALL cp_cfm_det(eigrmat(ispin), zdeta)
2033 zdet = zdet*zdeta
2034 IF (dft_control%nspins == 1) THEN
2035 zdet = zdet*zdeta
2036 END IF
2037 END DO
2038 zij(i, j) = zdet*xphase(i)*xphase(j)
2039 zij(j, i) = zdet*xphase(i)*xphase(j)
2040 END DO
2041 END DO
2042 CASE (3)
2043 ! Octapole
2044 cpabort("Berry phase moments bigger than 2 not implemented")
2045 CASE (4)
2046 ! Hexadecapole
2047 cpabort("Berry phase moments bigger than 3 not implemented")
2048 CASE DEFAULT
2049 cpabort("Berry phase moments bigger than 4 not implemented")
2050 END SELECT
2051 END DO
2052 DO l = 1, nmom
2053 SELECT CASE (l)
2054 CASE (1)
2055 ! Dipole (apply periodic (2 Pi) boundary conditions)
2056 ci = aimag(log(zi))
2057 DO i = 1, 3
2058 IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
2059 IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
2060 END DO
2061 rmom(2:4, 1) = matmul(cell%hmat, ci)/twopi
2062 CASE (2)
2063 ! Quadrupole
2064 cpabort("Berry phase moments bigger than 1 not implemented")
2065 DO i = 1, 3
2066 DO j = 1, 3
2067 zz = zij(i, j)/zi(i)/zi(j)
2068 cij(i, j) = aimag(log(zz))/twopi
2069 END DO
2070 END DO
2071 cij = 0.5_dp*cij/twopi/twopi
2072 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
2073 DO k = 4, 9
2074 ix = indco(1, k + 1)
2075 iy = indco(2, k + 1)
2076 iz = indco(3, k + 1)
2077 IF (ix == 0) THEN
2078 rmom(k + 1, 1) = cij(iy, iz)
2079 ELSE IF (iy == 0) THEN
2080 rmom(k + 1, 1) = cij(ix, iz)
2081 ELSE IF (iz == 0) THEN
2082 rmom(k + 1, 1) = cij(ix, iy)
2083 END IF
2084 END DO
2085 CASE (3)
2086 ! Octapole
2087 cpabort("Berry phase moments bigger than 2 not implemented")
2088 CASE (4)
2089 ! Hexadecapole
2090 cpabort("Berry phase moments bigger than 3 not implemented")
2091 CASE DEFAULT
2092 cpabort("Berry phase moments bigger than 4 not implemented")
2093 END SELECT
2094 END DO
2095
2096 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2097 description = "[DIPOLE]"
2098 CALL cp_results_erase(results=results, description=description)
2099 CALL put_results(results=results, description=description, &
2100 values=rmom(2:4, 3))
2101 IF (magnetic) THEN
2102 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true., mmom=mmom)
2103 ELSE
2104 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true.)
2105 END IF
2106
2107 DEALLOCATE (rmom)
2108 DEALLOCATE (rlab)
2109 IF (magnetic) THEN
2110 DEALLOCATE (mmom)
2111 END IF
2112
2113 CALL dbcsr_deallocate_matrix(cosmat)
2114 CALL dbcsr_deallocate_matrix(sinmat)
2115
2116 CALL cp_fm_release(opvec)
2117 CALL cp_fm_release(op_fm_set)
2118 DO ispin = 1, dft_control%nspins
2119 CALL cp_cfm_release(eigrmat(ispin))
2120 END DO
2121 DEALLOCATE (eigrmat)
2122
2123 CALL timestop(handle)
2124
2125 END SUBROUTINE qs_moment_berry_phase
2126
2127! **************************************************************************************************
2128!> \brief ...
2129!> \param cosmat ...
2130!> \param sinmat ...
2131!> \param mos ...
2132!> \param op_fm_set ...
2133!> \param opvec ...
2134! **************************************************************************************************
2135 SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
2136
2137 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2138 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2139 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2140 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: opvec
2141
2142 INTEGER :: i, nao, nmo
2143 TYPE(cp_fm_type), POINTER :: mo_coeff
2144
2145 DO i = 1, SIZE(op_fm_set, 2) ! spin
2146 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2147 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
2148 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2149 op_fm_set(1, i))
2150 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
2151 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2152 op_fm_set(2, i))
2153 END DO
2154
2155 END SUBROUTINE op_orbbas
2156
2157! **************************************************************************************************
2158!> \brief ...
2159!> \param cosmat ...
2160!> \param sinmat ...
2161!> \param mos ...
2162!> \param op_fm_set ...
2163!> \param mos_new ...
2164! **************************************************************************************************
2165 SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
2166
2167 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2168 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2169 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2170 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
2171
2172 INTEGER :: i, icol, lcol, nao, newdim, nmo
2173 LOGICAL :: double_col, double_row
2174 TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
2175 TYPE(cp_fm_type) :: work, work1, work2
2176 TYPE(cp_fm_type), POINTER :: mo_coeff
2177
2178 DO i = 1, SIZE(op_fm_set, 2) ! spin
2179 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2180 CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
2181 double_col = .true.
2182 double_row = .false.
2183 CALL cp_fm_struct_double(newstruct, &
2184 mos_new(2*i)%matrix_struct, &
2185 mos_new(2*i)%matrix_struct%context, &
2186 double_col, &
2187 double_row)
2188
2189 CALL cp_fm_create(work, matrix_struct=newstruct)
2190 CALL cp_fm_create(work1, matrix_struct=newstruct)
2191 CALL cp_fm_create(work2, matrix_struct=newstruct)
2192 CALL cp_fm_get_info(work, ncol_global=newdim)
2193
2194 CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
2195 DO icol = 1, lcol
2196 work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
2197 work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
2198 END DO
2199
2200 CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
2201 CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
2202
2203 DO icol = 1, lcol
2204 work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
2205 work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
2206 END DO
2207
2208 CALL cp_fm_release(work1)
2209 CALL cp_fm_release(work2)
2210
2211 CALL cp_fm_struct_double(newstruct1, &
2212 op_fm_set(1, i)%matrix_struct, &
2213 op_fm_set(1, i)%matrix_struct%context, &
2214 double_col, &
2215 double_row)
2216
2217 CALL cp_fm_create(work1, matrix_struct=newstruct1)
2218
2219 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
2220 work, 0.0_dp, work1)
2221
2222 DO icol = 1, lcol
2223 op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
2224 op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
2225 END DO
2226
2227 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
2228 work, 0.0_dp, work1)
2229
2230 DO icol = 1, lcol
2231 op_fm_set(1, i)%local_data(:, icol) = &
2232 op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
2233 op_fm_set(2, i)%local_data(:, icol) = &
2234 op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
2235 END DO
2236
2237 CALL cp_fm_release(work)
2238 CALL cp_fm_release(work1)
2239 CALL cp_fm_struct_release(newstruct)
2240 CALL cp_fm_struct_release(newstruct1)
2241
2242 END DO
2243
2244 END SUBROUTINE op_orbbas_rtp
2245
2246! **************************************************************************************************
2247!> \brief ...
2248!> \param qs_env ...
2249!> \param magnetic ...
2250!> \param nmoments ...
2251!> \param reference ...
2252!> \param ref_point ...
2253!> \param unit_number ...
2254!> \param vel_reprs ...
2255!> \param com_nl ...
2256! **************************************************************************************************
2257 SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
2258
2259 TYPE(qs_environment_type), POINTER :: qs_env
2260 LOGICAL, INTENT(IN) :: magnetic
2261 INTEGER, INTENT(IN) :: nmoments, reference
2262 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
2263 INTEGER, INTENT(IN) :: unit_number
2264 LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
2265
2266 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_locop'
2267
2268 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
2269 CHARACTER(LEN=default_string_length) :: description
2270 INTEGER :: akind, handle, i, ia, iatom, idir, &
2271 ikind, ispin, ix, iy, iz, l, nm, nmom, &
2272 order
2273 LOGICAL :: my_com_nl, my_velreprs
2274 REAL(dp) :: charge, dd, strace, trace
2275 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
2276 nlcom_rv, nlcom_rvr, nlcom_rxrv, &
2277 qupole_der, rmom_vel
2278 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
2279 REAL(dp), DIMENSION(3) :: rcc, ria
2280 TYPE(atomic_kind_type), POINTER :: atomic_kind
2281 TYPE(cell_type), POINTER :: cell
2282 TYPE(cp_result_type), POINTER :: results
2283 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
2284 rho_ao
2285 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
2286 TYPE(dbcsr_type), POINTER :: tmp_ao
2287 TYPE(dft_control_type), POINTER :: dft_control
2288 TYPE(distribution_1d_type), POINTER :: local_particles
2289 TYPE(mp_para_env_type), POINTER :: para_env
2290 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2291 POINTER :: sab_all, sab_orb
2292 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2293 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2294 TYPE(qs_rho_type), POINTER :: rho
2295
2296 cpassert(ASSOCIATED(qs_env))
2297
2298 CALL timeset(routinen, handle)
2299
2300 my_velreprs = .false.
2301 IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
2302 IF (PRESENT(com_nl)) my_com_nl = com_nl
2303 IF (my_velreprs) CALL cite_reference(mattiat2019)
2304
2305 ! reference point
2306 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
2307
2308 ! only allow for moments up to maxl set by basis
2309 nmom = min(nmoments, current_maxl)
2310 ! electronic contribution
2311 NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
2312 CALL get_qs_env(qs_env, &
2313 dft_control=dft_control, &
2314 rho=rho, &
2315 cell=cell, &
2316 results=results, &
2317 particle_set=particle_set, &
2318 qs_kind_set=qs_kind_set, &
2319 para_env=para_env, &
2320 matrix_s=matrix_s, &
2321 sab_all=sab_all, &
2322 sab_orb=sab_orb)
2323
2324 IF (my_com_nl) THEN
2325 IF ((nmom >= 1) .AND. my_velreprs) THEN
2326 ALLOCATE (nlcom_rv(3))
2327 nlcom_rv(:) = 0._dp
2328 END IF
2329 IF ((nmom >= 2) .AND. my_velreprs) THEN
2330 ALLOCATE (nlcom_rrv(6))
2331 nlcom_rrv(:) = 0._dp
2332 ALLOCATE (nlcom_rvr(6))
2333 nlcom_rvr(:) = 0._dp
2334 ALLOCATE (nlcom_rrv_vrr(6))
2335 nlcom_rrv_vrr(:) = 0._dp
2336 END IF
2337 IF (magnetic) THEN
2338 ALLOCATE (nlcom_rxrv(3))
2339 nlcom_rxrv = 0._dp
2340 END IF
2341 ! Calculate non local correction terms
2342 CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
2343 END IF
2344
2345 NULLIFY (moments)
2346 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
2347 CALL dbcsr_allocate_matrix_set(moments, nm)
2348 DO i = 1, nm
2349 ALLOCATE (moments(i)%matrix)
2350 IF (my_velreprs .AND. (nmom >= 2)) THEN
2351 CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
2352 matrix_type=dbcsr_type_symmetric)
2353 CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
2354 ELSE
2355 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
2356 END IF
2357 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
2358 END DO
2359
2360 ! calculate derivatives if quadrupole in vel. reprs. is requested
2361 IF (my_velreprs .AND. (nmom >= 2)) THEN
2362 NULLIFY (moments_der)
2363 CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
2364 DO i = 1, 3 ! x, y, z
2365 DO idir = 1, 3 ! d/dx, d/dy, d/dz
2366 CALL dbcsr_init_p(moments_der(i, idir)%matrix)
2367 CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
2368 matrix_type=dbcsr_type_antisymmetric)
2369 CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
2370 CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
2371 END DO
2372 END DO
2373 CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
2374 ELSE
2375 CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
2376 END IF
2377
2378 CALL qs_rho_get(rho, rho_ao=rho_ao)
2379
2380 ALLOCATE (rmom(nm + 1, 3))
2381 ALLOCATE (rlab(nm + 1))
2382 rmom = 0.0_dp
2383 rlab = ""
2384
2385 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2386 ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
2387 ! symmetric matrices)
2388 NULLIFY (tmp_ao)
2389 CALL dbcsr_init_p(tmp_ao)
2390 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2391 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2392 CALL dbcsr_set(tmp_ao, 0.0_dp)
2393 END IF
2394
2395 trace = 0.0_dp
2396 DO ispin = 1, dft_control%nspins
2397 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
2398 rmom(1, 1) = rmom(1, 1) + trace
2399 END DO
2400
2401 DO i = 1, SIZE(moments)
2402 strace = 0._dp
2403 DO ispin = 1, dft_control%nspins
2404 IF (my_velreprs .AND. nmoments >= 2) THEN
2405 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
2406 0.0_dp, tmp_ao)
2407 CALL dbcsr_trace(tmp_ao, trace)
2408 ELSE
2409 CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
2410 END IF
2411 strace = strace + trace
2412 END DO
2413 rmom(i + 1, 1) = strace
2414 END DO
2415
2416 CALL dbcsr_deallocate_matrix_set(moments)
2417
2418 ! nuclear contribution
2419 CALL get_qs_env(qs_env=qs_env, &
2420 local_particles=local_particles)
2421 DO ikind = 1, SIZE(local_particles%n_el)
2422 DO ia = 1, local_particles%n_el(ikind)
2423 iatom = local_particles%list(ikind)%array(ia)
2424 ! fold atomic positions back into unit cell
2425 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
2426 ria = ria - rcc
2427 atomic_kind => particle_set(iatom)%atomic_kind
2428 CALL get_atomic_kind(atomic_kind, kind_number=akind)
2429 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
2430 rmom(1, 2) = rmom(1, 2) - charge
2431 DO l = 1, nm
2432 ix = indco(1, l + 1)
2433 iy = indco(2, l + 1)
2434 iz = indco(3, l + 1)
2435 dd = 1._dp
2436 IF (ix > 0) dd = dd*ria(1)**ix
2437 IF (iy > 0) dd = dd*ria(2)**iy
2438 IF (iz > 0) dd = dd*ria(3)**iz
2439 rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
2440 CALL set_label(rlab(l + 1), ix, iy, iz)
2441 END DO
2442 END DO
2443 END DO
2444 CALL para_env%sum(rmom(:, 2))
2445 rmom(:, :) = -rmom(:, :)
2446 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2447
2448 ! magnetic moments
2449 IF (magnetic) THEN
2450 NULLIFY (magmom)
2451 CALL dbcsr_allocate_matrix_set(magmom, 3)
2452 DO i = 1, SIZE(magmom)
2453 CALL dbcsr_init_p(magmom(i)%matrix)
2454 CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2455 matrix_type=dbcsr_type_antisymmetric)
2456 CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2457 CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2458 END DO
2459
2460 CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2461
2462 ALLOCATE (mmom(SIZE(magmom)))
2463 mmom(:) = 0.0_dp
2464 IF (qs_env%run_rtp) THEN
2465 ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2466 ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2467 ! There may be other cases, where the imaginary part of the density is relevant
2468 NULLIFY (rho_ao)
2469 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2470 END IF
2471 ! if the density is purely real this is an expensive way to calculate zero
2472 DO i = 1, SIZE(magmom)
2473 strace = 0._dp
2474 DO ispin = 1, dft_control%nspins
2475 CALL dbcsr_set(tmp_ao, 0.0_dp)
2476 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2477 0.0_dp, tmp_ao)
2478 CALL dbcsr_trace(tmp_ao, trace)
2479 strace = strace + trace
2480 END DO
2481 mmom(i) = strace
2482 END DO
2483
2484 CALL dbcsr_deallocate_matrix_set(magmom)
2485 END IF
2486
2487 ! velocity representations
2488 IF (my_velreprs) THEN
2489 ALLOCATE (rmom_vel(nm))
2490 rmom_vel = 0.0_dp
2491
2492 DO order = 1, nmom
2493 SELECT CASE (order)
2494
2495 CASE (1) ! expectation value of momentum
2496 NULLIFY (momentum)
2497 CALL dbcsr_allocate_matrix_set(momentum, 3)
2498 DO i = 1, 3
2499 CALL dbcsr_init_p(momentum(i)%matrix)
2500 CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2501 matrix_type=dbcsr_type_antisymmetric)
2502 CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2503 CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2504 END DO
2505 CALL build_lin_mom_matrix(qs_env, momentum)
2506
2507 ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2508 IF (qs_env%run_rtp) THEN
2509 NULLIFY (rho_ao)
2510 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2511 DO idir = 1, SIZE(momentum)
2512 strace = 0._dp
2513 DO ispin = 1, dft_control%nspins
2514 CALL dbcsr_set(tmp_ao, 0.0_dp)
2515 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2516 0.0_dp, tmp_ao)
2517 CALL dbcsr_trace(tmp_ao, trace)
2518 strace = strace + trace
2519 END DO
2520 rmom_vel(idir) = rmom_vel(idir) + strace
2521 END DO
2522 END IF
2523
2524 CALL dbcsr_deallocate_matrix_set(momentum)
2525
2526 CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2527 ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2528 qupole_der = 0._dp
2529
2530 NULLIFY (rho_ao)
2531 CALL qs_rho_get(rho, rho_ao=rho_ao)
2532
2533 ! Calculate expectation value over real part
2534 trace = 0._dp
2535 DO i = 1, 3
2536 DO idir = 1, 3
2537 strace = 0._dp
2538 DO ispin = 1, dft_control%nspins
2539 CALL dbcsr_set(tmp_ao, 0._dp)
2540 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2541 CALL dbcsr_trace(tmp_ao, trace)
2542 strace = strace + trace
2543 END DO
2544 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2545 END DO
2546 END DO
2547
2548 IF (qs_env%run_rtp) THEN
2549 NULLIFY (rho_ao)
2550 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2551
2552 ! Calculate expectation value over imaginary part
2553 trace = 0._dp
2554 DO i = 1, 3
2555 DO idir = 1, 3
2556 strace = 0._dp
2557 DO ispin = 1, dft_control%nspins
2558 CALL dbcsr_set(tmp_ao, 0._dp)
2559 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2560 CALL dbcsr_trace(tmp_ao, trace)
2561 strace = strace + trace
2562 END DO
2563 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2564 END DO
2565 END DO
2566 END IF
2567
2568 ! calculate vel. reprs. of quadrupole moment from derivatives
2569 rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2570 rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2571 rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2572 rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2573 rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2574 rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2575
2576 DEALLOCATE (qupole_der)
2577 CASE DEFAULT
2578 END SELECT
2579 END DO
2580 END IF
2581
2582 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2583 CALL dbcsr_deallocate_matrix(tmp_ao)
2584 END IF
2585 IF (my_velreprs .AND. (nmoments >= 2)) THEN
2586 CALL dbcsr_deallocate_matrix_set(moments_der)
2587 END IF
2588
2589 description = "[DIPOLE]"
2590 CALL cp_results_erase(results=results, description=description)
2591 CALL put_results(results=results, description=description, &
2592 values=rmom(2:4, 3))
2593
2594 IF (magnetic .AND. my_velreprs) THEN
2595 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom, rmom_vel=rmom_vel)
2596 ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2597 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom)
2598 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2599 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., rmom_vel=rmom_vel)
2600 ELSE
2601 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false.)
2602 END IF
2603
2604 IF (my_com_nl) THEN
2605 IF (magnetic) THEN
2606 mmom(:) = nlcom_rxrv(:)
2607 END IF
2608 IF (my_velreprs .AND. (nmom >= 1)) THEN
2609 DEALLOCATE (rmom_vel)
2610 ALLOCATE (rmom_vel(21))
2611 rmom_vel(1:3) = nlcom_rv
2612 END IF
2613 IF (my_velreprs .AND. (nmom >= 2)) THEN
2614 rmom_vel(4:9) = nlcom_rrv
2615 rmom_vel(10:15) = nlcom_rvr
2616 rmom_vel(16:21) = nlcom_rrv_vrr
2617 END IF
2618 IF (magnetic .AND. .NOT. my_velreprs) THEN
2619 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2620 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2621 CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2622 ELSE IF (my_velreprs .AND. magnetic) THEN
2623 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2624 END IF
2625
2626 END IF
2627
2628 IF (my_com_nl) THEN
2629 IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2630 IF (nmom >= 2 .AND. my_velreprs) THEN
2631 DEALLOCATE (nlcom_rrv)
2632 DEALLOCATE (nlcom_rvr)
2633 DEALLOCATE (nlcom_rrv_vrr)
2634 END IF
2635 IF (magnetic) DEALLOCATE (nlcom_rxrv)
2636 END IF
2637
2638 DEALLOCATE (rmom)
2639 DEALLOCATE (rlab)
2640 IF (magnetic) THEN
2641 DEALLOCATE (mmom)
2642 END IF
2643 IF (my_velreprs) THEN
2644 DEALLOCATE (rmom_vel)
2645 END IF
2646
2647 CALL timestop(handle)
2648
2649 END SUBROUTINE qs_moment_locop
2650
2651! **************************************************************************************************
2652!> \brief ...
2653!> \param label ...
2654!> \param ix ...
2655!> \param iy ...
2656!> \param iz ...
2657! **************************************************************************************************
2658 SUBROUTINE set_label(label, ix, iy, iz)
2659 CHARACTER(LEN=*), INTENT(OUT) :: label
2660 INTEGER, INTENT(IN) :: ix, iy, iz
2661
2662 INTEGER :: i
2663
2664 label = ""
2665 DO i = 1, ix
2666 WRITE (label(i:), "(A1)") "X"
2667 END DO
2668 DO i = ix + 1, ix + iy
2669 WRITE (label(i:), "(A1)") "Y"
2670 END DO
2671 DO i = ix + iy + 1, ix + iy + iz
2672 WRITE (label(i:), "(A1)") "Z"
2673 END DO
2674
2675 END SUBROUTINE set_label
2676
2677! **************************************************************************************************
2678!> \brief ...
2679!> \param unit_number ...
2680!> \param nmom ...
2681!> \param rmom ...
2682!> \param rlab ...
2683!> \param rcc ...
2684!> \param cell ...
2685!> \param periodic ...
2686!> \param mmom ...
2687!> \param rmom_vel ...
2688! **************************************************************************************************
2689 SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
2690 INTEGER, INTENT(IN) :: unit_number, nmom
2691 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2692 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2693 REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2694 TYPE(cell_type), POINTER :: cell
2695 LOGICAL :: periodic
2696 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2697
2698 INTEGER :: i, i0, i1, j, l
2699 REAL(dp) :: dd
2700
2701 IF (unit_number > 0) THEN
2702 DO l = 0, nmom
2703 SELECT CASE (l)
2704 CASE (0)
2705 WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2706 WRITE (unit_number, "(T3,A)") "Charges"
2707 WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2708 "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2709 CASE (1)
2710 IF (periodic) THEN
2711 WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2712 WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
2713 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
2714 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
2715 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
2716 ELSE
2717 WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2718 END IF
2719 dd = sqrt(sum(rmom(2:4, 3)**2))*debye
2720 WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
2721 WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
2722 (trim(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
2723 CASE (2)
2724 WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
2725 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2726 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
2727 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2728 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
2729 CASE (3)
2730 WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
2731 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2732 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
2733 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2734 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
2735 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2736 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
2737 CASE (4)
2738 WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
2739 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2740 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
2741 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2742 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
2743 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2744 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
2745 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2746 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
2747 CASE DEFAULT
2748 WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
2749 " L=", l
2750 i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2751 i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2752 dd = debye/(bohr)**(l - 1)
2753 DO i = i0, i1, 3
2754 WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
2755 (trim(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, min(i1, i + 2))
2756 END DO
2757 END SELECT
2758 END DO
2759 IF (PRESENT(mmom)) THEN
2760 IF (nmom >= 1) THEN
2761 dd = sqrt(sum(mmom(1:3)**2))
2762 WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2763 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2764 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2765 END IF
2766 END IF
2767 IF (PRESENT(rmom_vel)) THEN
2768 DO l = 1, nmom
2769 SELECT CASE (l)
2770 CASE (1)
2771 dd = sqrt(sum(rmom_vel(1:3)**2))
2772 WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2773 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2774 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2775 CASE (2)
2776 WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2777 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2778 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2779 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2780 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2781 CASE DEFAULT
2782 END SELECT
2783 END DO
2784 END IF
2785 END IF
2786
2787 END SUBROUTINE print_moments
2788
2789! **************************************************************************************************
2790!> \brief ...
2791!> \param unit_number ...
2792!> \param nmom ...
2793!> \param rlab ...
2794!> \param mmom ...
2795!> \param rmom_vel ...
2796! **************************************************************************************************
2797 SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2798 INTEGER, INTENT(IN) :: unit_number, nmom
2799 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2800 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2801
2802 INTEGER :: i, l
2803 REAL(dp) :: dd
2804
2805 IF (unit_number > 0) THEN
2806 IF (PRESENT(mmom)) THEN
2807 IF (nmom >= 1) THEN
2808 dd = sqrt(sum(mmom(1:3)**2))
2809 WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2810 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2811 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2812 END IF
2813 END IF
2814 IF (PRESENT(rmom_vel)) THEN
2815 DO l = 1, nmom
2816 SELECT CASE (l)
2817 CASE (1)
2818 dd = sqrt(sum(rmom_vel(1:3)**2))
2819 WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2820 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2821 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2822 CASE (2)
2823 WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2824 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2825 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2826 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2827 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2828 WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2829 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2830 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2831 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2832 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2833 WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2834 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2835 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2836 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2837 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2838 CASE DEFAULT
2839 END SELECT
2840 END DO
2841 END IF
2842 END IF
2843
2844 END SUBROUTINE print_moments_nl
2845
2846! **************************************************************************************************
2847!> \brief Calculate the expectation value of operators related to non-local potential:
2848!> [r, Vnl], noted rv
2849!> r x [r,Vnl], noted rxrv
2850!> [rr,Vnl], noted rrv
2851!> r x Vnl x r, noted rvr
2852!> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2853!> Note that the 3 first operator are commutator while the 2 last
2854!> are not. For reading clarity the same notation is used for all 5
2855!> operators.
2856!> \param qs_env ...
2857!> \param nlcom_rv ...
2858!> \param nlcom_rxrv ...
2859!> \param nlcom_rrv ...
2860!> \param nlcom_rvr ...
2861!> \param nlcom_rrv_vrr ...
2862!> \param ref_point ...
2863! **************************************************************************************************
2864 SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2865 nlcom_rrv_vrr, ref_point)
2866
2867 TYPE(qs_environment_type), POINTER :: qs_env
2868 REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2869 nlcom_rvr, nlcom_rrv_vrr
2870 REAL(dp), DIMENSION(3) :: ref_point
2871
2872 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_commutator_nl_terms'
2873
2874 INTEGER :: handle, ind, ispin
2875 LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2876 calc_rvr, calc_rxrv
2877 REAL(dp) :: eps_ppnl, strace, trace
2878 TYPE(cell_type), POINTER :: cell
2879 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2880 matrix_rvr, matrix_rxrv, matrix_s, &
2881 rho_ao
2882 TYPE(dbcsr_type), POINTER :: tmp_ao
2883 TYPE(dft_control_type), POINTER :: dft_control
2884 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2885 POINTER :: sab_all, sab_orb, sap_ppnl
2886 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2887 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2888 TYPE(qs_rho_type), POINTER :: rho
2889
2890 CALL timeset(routinen, handle)
2891
2892 calc_rv = .false.
2893 calc_rxrv = .false.
2894 calc_rrv = .false.
2895 calc_rvr = .false.
2896 calc_rrv_vrr = .false.
2897
2898 ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2899 ! The real part of the density matrix rho_ao is symmetric so that
2900 ! the expectation value of real density matrix is zero. Hence, if
2901 ! the density matrix is real, no need to compute these quantities.
2902 ! This is not the case for rvr and rrv_vrr which are symmetric.
2903
2904 IF (ALLOCATED(nlcom_rv)) THEN
2905 nlcom_rv(:) = 0._dp
2906 IF (qs_env%run_rtp) calc_rv = .true.
2907 END IF
2908 IF (ALLOCATED(nlcom_rxrv)) THEN
2909 nlcom_rxrv(:) = 0._dp
2910 IF (qs_env%run_rtp) calc_rxrv = .true.
2911 END IF
2912 IF (ALLOCATED(nlcom_rrv)) THEN
2913 nlcom_rrv(:) = 0._dp
2914 IF (qs_env%run_rtp) calc_rrv = .true.
2915 END IF
2916 IF (ALLOCATED(nlcom_rvr)) THEN
2917 nlcom_rvr(:) = 0._dp
2918 calc_rvr = .true.
2919 END IF
2920 IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2921 nlcom_rrv_vrr(:) = 0._dp
2922 calc_rrv_vrr = .true.
2923 END IF
2924
2925 IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv .OR. calc_rvr .OR. calc_rrv_vrr)) THEN
2926 CALL timestop(handle)
2927 RETURN
2928 END IF
2929
2930 NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2931 CALL get_qs_env(qs_env, &
2932 cell=cell, &
2933 dft_control=dft_control, &
2934 matrix_s=matrix_s, &
2935 particle_set=particle_set, &
2936 qs_kind_set=qs_kind_set, &
2937 rho=rho, &
2938 sab_orb=sab_orb, &
2939 sab_all=sab_all, &
2940 sap_ppnl=sap_ppnl)
2941
2942 eps_ppnl = dft_control%qs_control%eps_ppnl
2943
2944 ! Allocate storage
2945 NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2946 IF (calc_rv) THEN
2947 CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2948 DO ind = 1, 3
2949 CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2950 CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2951 matrix_type=dbcsr_type_antisymmetric)
2952 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2953 CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2954 END DO
2955 END IF
2956
2957 IF (calc_rxrv) THEN
2958 CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2959 DO ind = 1, 3
2960 CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2961 CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2962 matrix_type=dbcsr_type_antisymmetric)
2963 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2964 CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2965 END DO
2966 END IF
2967
2968 IF (calc_rrv) THEN
2969 CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2970 DO ind = 1, 6
2971 CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2972 CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2973 matrix_type=dbcsr_type_antisymmetric)
2974 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2975 CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2976 END DO
2977 END IF
2978
2979 IF (calc_rvr) THEN
2980 CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2981 DO ind = 1, 6
2982 CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2983 CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2984 matrix_type=dbcsr_type_symmetric)
2985 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2986 CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2987 END DO
2988 END IF
2989 IF (calc_rrv_vrr) THEN
2990 CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2991 DO ind = 1, 6
2992 CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2993 CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2994 matrix_type=dbcsr_type_symmetric)
2995 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2996 CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2997 END DO
2998 END IF
2999
3000 ! calculate evaluation of operators in AO basis set
3001 CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
3002 matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
3003 matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
3004
3005 ! Calculate expectation values
3006 ! Real part
3007 NULLIFY (tmp_ao)
3008 CALL dbcsr_init_p(tmp_ao)
3009 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
3010 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
3011 CALL dbcsr_set(tmp_ao, 0.0_dp)
3012
3013 IF (calc_rvr .OR. calc_rrv_vrr) THEN
3014 NULLIFY (rho_ao)
3015 CALL qs_rho_get(rho, rho_ao=rho_ao)
3016
3017 IF (calc_rvr) THEN
3018 trace = 0._dp
3019 DO ind = 1, SIZE(matrix_rvr)
3020 strace = 0._dp
3021 DO ispin = 1, dft_control%nspins
3022 CALL dbcsr_set(tmp_ao, 0.0_dp)
3023 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
3024 0.0_dp, tmp_ao)
3025 CALL dbcsr_trace(tmp_ao, trace)
3026 strace = strace + trace
3027 END DO
3028 nlcom_rvr(ind) = nlcom_rvr(ind) + strace
3029 END DO
3030 END IF
3031
3032 IF (calc_rrv_vrr) THEN
3033 trace = 0._dp
3034 DO ind = 1, SIZE(matrix_rrv_vrr)
3035 strace = 0._dp
3036 DO ispin = 1, dft_control%nspins
3037 CALL dbcsr_set(tmp_ao, 0.0_dp)
3038 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
3039 0.0_dp, tmp_ao)
3040 CALL dbcsr_trace(tmp_ao, trace)
3041 strace = strace + trace
3042 END DO
3043 nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
3044 END DO
3045 END IF
3046 END IF
3047
3048 ! imagninary part of the density matrix
3049 NULLIFY (rho_ao)
3050 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
3051
3052 IF (calc_rv) THEN
3053 trace = 0._dp
3054 DO ind = 1, SIZE(matrix_rv)
3055 strace = 0._dp
3056 DO ispin = 1, dft_control%nspins
3057 CALL dbcsr_set(tmp_ao, 0.0_dp)
3058 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
3059 0.0_dp, tmp_ao)
3060 CALL dbcsr_trace(tmp_ao, trace)
3061 strace = strace + trace
3062 END DO
3063 nlcom_rv(ind) = nlcom_rv(ind) + strace
3064 END DO
3065 END IF
3066
3067 IF (calc_rrv) THEN
3068 trace = 0._dp
3069 DO ind = 1, SIZE(matrix_rrv)
3070 strace = 0._dp
3071 DO ispin = 1, dft_control%nspins
3072 CALL dbcsr_set(tmp_ao, 0.0_dp)
3073 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
3074 0.0_dp, tmp_ao)
3075 CALL dbcsr_trace(tmp_ao, trace)
3076 strace = strace + trace
3077 END DO
3078 nlcom_rrv(ind) = nlcom_rrv(ind) + strace
3079 END DO
3080 END IF
3081
3082 IF (calc_rxrv) THEN
3083 trace = 0._dp
3084 DO ind = 1, SIZE(matrix_rxrv)
3085 strace = 0._dp
3086 DO ispin = 1, dft_control%nspins
3087 CALL dbcsr_set(tmp_ao, 0.0_dp)
3088 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
3089 0.0_dp, tmp_ao)
3090 CALL dbcsr_trace(tmp_ao, trace)
3091 strace = strace + trace
3092 END DO
3093 nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
3094 END DO
3095 END IF
3096 CALL dbcsr_deallocate_matrix(tmp_ao)
3097 IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
3098 IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
3099 IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
3100 IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
3101 IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
3102
3103 CALL timestop(handle)
3104 END SUBROUTINE calculate_commutator_nl_terms
3105
3106! *****************************************************************************
3107!> \brief ...
3108!> \param qs_env ...
3109!> \param difdip ...
3110!> \param deltaR ...
3111!> \param order ...
3112!> \param rcc ...
3113!> \note calculate matrix elements <a|r_beta |db/dR_alpha > + <da/dR_alpha | r_beta | b >
3114!> be aware: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > only valid
3115!> if alpha .neq.beta
3116!> if alpha=beta: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > - < a | b >
3117!> modified from qs_efield_mo_derivatives
3118!> SL July 2015
3119! **************************************************************************************************
3120 SUBROUTINE dipole_deriv_ao(qs_env, difdip, deltaR, order, rcc)
3121 TYPE(qs_environment_type), POINTER :: qs_env
3122 TYPE(dbcsr_p_type), DIMENSION(:, :), &
3123 INTENT(INOUT), POINTER :: difdip
3124 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
3125 POINTER :: deltar
3126 INTEGER, INTENT(IN) :: order
3127 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3128
3129 CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_deriv_ao'
3130
3131 INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
3132 last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
3133 sgfb
3134 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3135 npgfb, nsgfa, nsgfb
3136 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3137 LOGICAL :: found
3138 REAL(dp) :: dab
3139 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
3140 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3141 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
3142 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
3143 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
3144 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
3145 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: difmab2
3146 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
3147 TYPE(cell_type), POINTER :: cell
3148 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3149 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3150 TYPE(neighbor_list_iterator_p_type), &
3151 DIMENSION(:), POINTER :: nl_iterator
3152 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3153 POINTER :: sab_all, sab_orb
3154 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3155 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3156 TYPE(qs_kind_type), POINTER :: qs_kind
3157
3158 CALL timeset(routinen, handle)
3159
3160 NULLIFY (cell, particle_set, qs_kind_set, sab_orb, sab_all)
3161 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
3162 qs_kind_set=qs_kind_set, sab_orb=sab_orb, sab_all=sab_all)
3163 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
3164 maxco=ldab, maxsgf=maxsgf)
3165
3166 nkind = SIZE(qs_kind_set)
3167 natom = SIZE(particle_set)
3168
3169 m_dim = ncoset(order) - 1
3170
3171 IF (PRESENT(rcc)) THEN
3172 rc = rcc
3173 ELSE
3174 rc = 0._dp
3175 END IF
3176
3177 ALLOCATE (basis_set_list(nkind))
3178
3179 ALLOCATE (mab(ldab, ldab, m_dim))
3180 ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
3181 ALLOCATE (work(ldab, maxsgf))
3182 ALLOCATE (mint(3, 3))
3183 ALLOCATE (mint2(3, 3))
3184
3185 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
3186 difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
3187 work(1:ldab, 1:maxsgf) = 0.0_dp
3188
3189 DO i = 1, 3
3190 DO j = 1, 3
3191 NULLIFY (mint(i, j)%block)
3192 NULLIFY (mint2(i, j)%block)
3193 END DO
3194 END DO
3195
3196 ! Set the basis_set_list(nkind) to point to the corresponding basis sets
3197 DO ikind = 1, nkind
3198 qs_kind => qs_kind_set(ikind)
3199 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
3200 IF (ASSOCIATED(basis_set_a)) THEN
3201 basis_set_list(ikind)%gto_basis_set => basis_set_a
3202 ELSE
3203 NULLIFY (basis_set_list(ikind)%gto_basis_set)
3204 END IF
3205 END DO
3206
3207 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
3208 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3209 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
3210 iatom=iatom, jatom=jatom, r=rab)
3211
3212 basis_set_a => basis_set_list(ikind)%gto_basis_set
3213 basis_set_b => basis_set_list(jkind)%gto_basis_set
3214 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
3215 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
3216
3217 ! basis ikind
3218 first_sgfa => basis_set_a%first_sgf
3219 la_max => basis_set_a%lmax
3220 la_min => basis_set_a%lmin
3221 npgfa => basis_set_a%npgf
3222 nseta = basis_set_a%nset
3223 nsgfa => basis_set_a%nsgf_set
3224 rpgfa => basis_set_a%pgf_radius
3225 set_radius_a => basis_set_a%set_radius
3226 sphi_a => basis_set_a%sphi
3227 zeta => basis_set_a%zet
3228 ! basis jkind
3229 first_sgfb => basis_set_b%first_sgf
3230 lb_max => basis_set_b%lmax
3231 lb_min => basis_set_b%lmin
3232 npgfb => basis_set_b%npgf
3233 nsetb = basis_set_b%nset
3234 nsgfb => basis_set_b%nsgf_set
3235 rpgfb => basis_set_b%pgf_radius
3236 set_radius_b => basis_set_b%set_radius
3237 sphi_b => basis_set_b%sphi
3238 zetb => basis_set_b%zet
3239
3240 IF (inode == 1) last_jatom = 0
3241
3242 ! this guarentees minimum image convention
3243 ! anything else would not make sense
3244 IF (jatom == last_jatom) THEN
3245 cycle
3246 END IF
3247
3248 last_jatom = jatom
3249
3250 irow = iatom
3251 icol = jatom
3252
3253 DO i = 1, 3
3254 DO j = 1, 3
3255 NULLIFY (mint(i, j)%block)
3256 CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
3257 row=irow, col=icol, block=mint(i, j)%block, &
3258 found=found)
3259 cpassert(found)
3260 END DO
3261 END DO
3262
3263 ra(:) = particle_set(iatom)%r(:)
3264 rb(:) = particle_set(jatom)%r(:)
3265 rab(:) = pbc(rb, ra, cell)
3266 rac(:) = pbc(ra - rc, cell)
3267 rbc(:) = pbc(rb - rc, cell)
3268 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
3269
3270 DO iset = 1, nseta
3271 ncoa = npgfa(iset)*ncoset(la_max(iset))
3272 sgfa = first_sgfa(1, iset)
3273 DO jset = 1, nsetb
3274 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
3275 ncob = npgfb(jset)*ncoset(lb_max(jset))
3276 sgfb = first_sgfb(1, jset)
3277 ldab = max(ncoa, ncob)
3278 lda = ncoset(la_max(iset))*npgfa(iset)
3279 ldb = ncoset(lb_max(jset))*npgfb(jset)
3280 ALLOCATE (difmab(lda, ldb, m_dim, 3))
3281
3282 ! Calculate integral (da|r|b)
3283 CALL diff_momop2(la_max(iset), npgfa(iset), zeta(:, iset), &
3284 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
3285 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
3286 difmab, deltar=deltar, iatom=iatom, jatom=jatom)
3287
3288! *** Contraction step ***
3289
3290 DO idir = 1, 3 ! derivative of AO function
3291 DO j = 1, 3 ! position operator r_j
3292 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
3293 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
3294 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
3295 0.0_dp, work(1, 1), SIZE(work, 1))
3296
3297 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
3298 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
3299 work(1, 1), SIZE(work, 1), &
3300 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
3301 SIZE(mint(j, idir)%block, 1))
3302 END DO !j
3303 END DO !idir
3304 DEALLOCATE (difmab)
3305 END DO !jset
3306 END DO !iset
3307 END do!iterator
3308
3309 CALL neighbor_list_iterator_release(nl_iterator)
3310
3311 DO i = 1, 3
3312 DO j = 1, 3
3313 NULLIFY (mint(i, j)%block)
3314 END DO
3315 END DO
3316
3317 DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
3318
3319 CALL timestop(handle)
3320 END SUBROUTINE dipole_deriv_ao
3321
3322! **************************************************************************************************
3323!> \brief Get list of kpoints from input to compute dipole moment elements
3324!> \param qs_env ...
3325!> \param xkp ...
3326!> \param special_pnts ...
3327!> \author Shridhar Shanbhag
3328! **************************************************************************************************
3329 SUBROUTINE get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
3330 TYPE(qs_environment_type), POINTER :: qs_env
3331 TYPE(section_vals_type), POINTER :: kpnts, kpset
3332 REAL(kind=dp), DIMENSION(3, 3) :: hmat
3333
3334 CHARACTER(LEN=*), PARAMETER :: routinen = 'get_xkp_for_dipole_calc'
3335
3336 CHARACTER(LEN=default_string_length) :: ustr
3337 TYPE(kpoint_type), POINTER :: kpoint_work
3338 TYPE(cell_type), POINTER :: cell
3339 CHARACTER(LEN=default_string_length), &
3340 DIMENSION(:), POINTER :: strptr
3341 CHARACTER(LEN=default_string_length), &
3342 DIMENSION(:), POINTER :: special_pnts, spname
3343 CHARACTER(LEN=max_line_length) :: error_message
3344 INTEGER :: handle, i, ik, ikk, ip, &
3345 n_ptr, npline, nkp
3346 LOGICAL :: explicit_kpnts, explicit_kpset
3347 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: kspecial, xkp
3348 REAL(kind=dp), DIMENSION(3) :: kpptr
3349 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3350
3351 CALL timeset(routinen, handle)
3352 kpset => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINT_SET")
3353 kpnts => section_vals_get_subs_vals(qs_env%input, "DFT%PRINT%MOMENTS%KPOINTS")
3354 CALL section_vals_get(kpset, explicit=explicit_kpset)
3355 CALL section_vals_get(kpnts, explicit=explicit_kpnts)
3356 IF (explicit_kpset .AND. explicit_kpnts) &
3357 cpabort("Both KPOINT_SET and KPOINTS present in MOMENTS section")
3358
3359 IF (explicit_kpset) THEN
3360 CALL get_qs_env(qs_env, cell=cell)
3361 CALL get_cell(cell, h=hmat)
3362 CALL section_vals_val_get(kpset, "NPOINTS", i_val=npline)
3363 CALL section_vals_val_get(kpset, "UNITS", c_val=ustr)
3364 CALL uppercase(ustr)
3365 CALL section_vals_val_get(kpset, "SPECIAL_POINT", n_rep_val=n_ptr)
3366 cpassert(n_ptr > 0)
3367 ALLOCATE (kspecial(3, n_ptr))
3368 ALLOCATE (spname(n_ptr))
3369 DO ip = 1, n_ptr
3370 CALL section_vals_val_get(kpset, "SPECIAL_POINT", i_rep_val=ip, c_vals=strptr)
3371 IF (SIZE(strptr(:), 1) == 4) THEN
3372 spname(ip) = strptr(1)
3373 DO i = 1, 3
3374 CALL read_float_object(strptr(i + 1), kpptr(i), error_message)
3375 IF (len_trim(error_message) > 0) cpabort(trim(error_message))
3376 END DO
3377 ELSE IF (SIZE(strptr(:), 1) == 3) THEN
3378 spname(ip) = "not specified"
3379 DO i = 1, 3
3380 CALL read_float_object(strptr(i), kpptr(i), error_message)
3381 IF (len_trim(error_message) > 0) cpabort(trim(error_message))
3382 END DO
3383 ELSE
3384 cpabort("Input SPECIAL_POINT invalid")
3385 END IF
3386 SELECT CASE (ustr)
3387 CASE ("B_VECTOR")
3388 kspecial(1:3, ip) = kpptr(1:3)
3389 CASE ("CART_ANGSTROM")
3390 kspecial(1:3, ip) = (kpptr(1)*hmat(1, 1:3) + &
3391 kpptr(2)*hmat(2, 1:3) + &
3392 kpptr(3)*hmat(3, 1:3))/twopi*angstrom
3393 CASE ("CART_BOHR")
3394 kspecial(1:3, ip) = (kpptr(1)*hmat(1, 1:3) + &
3395 kpptr(2)*hmat(2, 1:3) + &
3396 kpptr(3)*hmat(3, 1:3))/twopi
3397 CASE DEFAULT
3398 cpabort("Unknown unit <"//trim(ustr)//"> specified for k-point definition")
3399 END SELECT
3400 END DO
3401 nkp = (n_ptr - 1)*npline + 1
3402 cpassert(nkp >= 1)
3403
3404 ! Initialize environment and calculate MOs
3405 ALLOCATE (xkp(3, nkp))
3406 ALLOCATE (special_pnts(nkp))
3407 special_pnts(:) = ""
3408 xkp(1:3, 1) = kspecial(1:3, 1)
3409 ikk = 1
3410 special_pnts(ikk) = spname(1)
3411 DO ik = 2, n_ptr
3412 DO ip = 1, npline
3413 ikk = ikk + 1
3414 xkp(1:3, ikk) = kspecial(1:3, ik - 1) + &
3415 REAL(ip, kind=dp)/real(npline, kind=dp)* &
3416 (kspecial(1:3, ik) - kspecial(1:3, ik - 1))
3417 END DO
3418 special_pnts(ikk) = spname(ik)
3419 END DO
3420 DEALLOCATE (spname, kspecial)
3421 ELSE IF (explicit_kpnts) THEN
3422 CALL get_qs_env(qs_env, particle_set=particle_set, cell=cell)
3423 CALL get_cell(cell, h=hmat)
3424 NULLIFY (kpoint_work)
3425 CALL kpoint_create(kpoint_work)
3426 CALL read_kpoint_section(kpoint_work, kpnts, hmat)
3427 CALL kpoint_initialize(kpoint_work, particle_set, cell)
3428 nkp = kpoint_work%nkp
3429 ALLOCATE (xkp(3, nkp))
3430 ALLOCATE (special_pnts(nkp))
3431 special_pnts(:) = ""
3432 xkp(1:3, :) = kpoint_work%xkp(1:3, :)
3433 CALL kpoint_release(kpoint_work)
3434 ELSE
3435 ! use k-point mesh from DFT calculation
3436 CALL get_qs_env(qs_env, kpoints=kpoint_work)
3437 nkp = kpoint_work%nkp
3438 nkp = kpoint_work%nkp
3439 ALLOCATE (xkp(3, nkp))
3440 ALLOCATE (special_pnts(nkp))
3441 special_pnts(:) = ""
3442 xkp(1:3, :) = kpoint_work%xkp(1:3, :)
3443 END IF
3444 CALL timestop(handle)
3445
3446 END SUBROUTINE get_xkp_for_dipole_calc
3447! **************************************************************************************************
3448!> \brief Calculate local moment matrix for a periodic system for all image cells
3449!> \param qs_env ...
3450!> \param moments_rs_img ...
3451!> \param rcc ...
3452!> \author Shridhar Shanbhag
3453! **************************************************************************************************
3454 SUBROUTINE build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
3455
3456 TYPE(qs_environment_type), POINTER :: qs_env
3457 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img
3458 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3459
3460 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix_rs_img'
3461
3462 INTEGER :: handle, i_dir, iatom, ic, ikind, iset, j, jatom, jkind, jset, &
3463 ldsa, ldsb, ldwork, ncoa, ncob, nimg, nkind, nseta, nsetb, nsize, sgfa, sgfb
3464 INTEGER, DIMENSION(3) :: icell
3465 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell
3466 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
3467 LOGICAL :: found
3468 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
3469 REAL(dp), DIMENSION(:, :), POINTER :: dblock, work
3470 REAL(dp), DIMENSION(:, :, :), POINTER :: dipab
3471 REAL(kind=dp) :: dab
3472 TYPE(cell_type), POINTER :: cell
3473 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
3474 TYPE(dft_control_type), POINTER :: dft_control
3475 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3476 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
3477 TYPE(kpoint_type), POINTER :: kpoints_all
3478 TYPE(mp_para_env_type), POINTER :: para_env
3479 TYPE(neighbor_list_iterator_p_type), &
3480 DIMENSION(:), POINTER :: nl_iterator
3481 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3482 POINTER :: sab_all
3483 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3484 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3485 TYPE(qs_kind_type), POINTER :: qs_kind
3486
3487 CALL timeset(routinen, handle)
3488
3489 CALL get_qs_env(qs_env=qs_env, &
3490 dft_control=dft_control, &
3491 qs_kind_set=qs_kind_set, &
3492 matrix_ks_kp=matrix_ks_kp, &
3493 particle_set=particle_set, &
3494 cell=cell, &
3495 para_env=para_env, &
3496 sab_all=sab_all)
3497
3498 NULLIFY (kpoints_all)
3499 CALL kpoint_create(kpoints_all)
3500 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, nimg)
3501
3502 nkind = SIZE(qs_kind_set)
3503 ALLOCATE (basis_set_list(nkind))
3504 DO ikind = 1, nkind
3505 qs_kind => qs_kind_set(ikind)
3506 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set)
3507 IF (ASSOCIATED(basis_set)) THEN
3508 basis_set_list(ikind)%gto_basis_set => basis_set
3509 ELSE
3510 NULLIFY (basis_set_list(ikind)%gto_basis_set)
3511 END IF
3512 END DO
3513
3514 rc(:) = 0._dp
3515 IF (PRESENT(rcc)) rc(:) = rcc(:)
3516
3517 CALL get_particle_set(particle_set, qs_kind_set, basis=basis_set_list)
3518 CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index, index_to_cell=index_to_cell)
3519 nsize = SIZE(index_to_cell, 2)
3520 cpassert(SIZE(moments_rs_img, 2) == nsize)
3521 DO i_dir = 1, 3
3522 DO j = 1, nsize
3523 ALLOCATE (moments_rs_img(i_dir, j)%matrix)
3524 CALL dbcsr_create(matrix=moments_rs_img(i_dir, j)%matrix, &
3525 template=matrix_ks_kp(1, 1)%matrix, &
3526 matrix_type=dbcsr_type_no_symmetry, &
3527 name="DIPMAT")
3528 CALL cp_dbcsr_alloc_block_from_nbl(moments_rs_img(i_dir, j)%matrix, sab_all)
3529 CALL dbcsr_set(moments_rs_img(i_dir, j)%matrix, 0.0_dp)
3530 END DO
3531 END DO
3532
3533 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
3534 ALLOCATE (dipab(ldwork, ldwork, 3))
3535 ALLOCATE (work(ldwork, ldwork))
3536
3537 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
3538 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3539 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
3540 iatom=iatom, jatom=jatom, r=rab, cell=icell)
3541
3542 basis_set_a => basis_set_list(ikind)%gto_basis_set
3543 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
3544 basis_set_b => basis_set_list(jkind)%gto_basis_set
3545 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
3546 associate( &
3547 ! basis ikind
3548 first_sgfa => basis_set_a%first_sgf, &
3549 la_max => basis_set_a%lmax, &
3550 la_min => basis_set_a%lmin, &
3551 npgfa => basis_set_a%npgf, &
3552 nsgfa => basis_set_a%nsgf_set, &
3553 rpgfa => basis_set_a%pgf_radius, &
3554 set_radius_a => basis_set_a%set_radius, &
3555 sphi_a => basis_set_a%sphi, &
3556 zeta => basis_set_a%zet, &
3557 ! basis jkind, &
3558 first_sgfb => basis_set_b%first_sgf, &
3559 lb_max => basis_set_b%lmax, &
3560 lb_min => basis_set_b%lmin, &
3561 npgfb => basis_set_b%npgf, &
3562 nsgfb => basis_set_b%nsgf_set, &
3563 rpgfb => basis_set_b%pgf_radius, &
3564 set_radius_b => basis_set_b%set_radius, &
3565 sphi_b => basis_set_b%sphi, &
3566 zetb => basis_set_b%zet)
3567
3568 nseta = basis_set_a%nset
3569 nsetb = basis_set_b%nset
3570
3571 ldsa = SIZE(sphi_a, 1)
3572 ldsb = SIZE(sphi_b, 1)
3573
3574 NULLIFY (dblock)
3575
3576 ra = pbc(particle_set(iatom)%r(:), cell)
3577 rb(:) = ra(:) + rab(:)
3578 rac = ra - rc
3579 rbc = rb - rc
3580 dab = norm2(rab)
3581
3582 ic = cell_to_index(icell(1), icell(2), icell(3))
3583
3584 DO iset = 1, nseta
3585
3586 ncoa = npgfa(iset)*ncoset(la_max(iset))
3587 sgfa = first_sgfa(1, iset)
3588
3589 DO jset = 1, nsetb
3590
3591 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
3592
3593 ncob = npgfb(jset)*ncoset(lb_max(jset))
3594 sgfb = first_sgfb(1, jset)
3595
3596 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
3597 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), 1, &
3598 rac, rbc, dipab)
3599 DO i_dir = 1, 3
3600 CALL dbcsr_get_block_p(matrix=moments_rs_img(i_dir, ic)%matrix, &
3601 row=iatom, col=jatom, block=dblock, found=found)
3602 cpassert(found)
3603 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
3604 1.0_dp, dipab(1, 1, i_dir), ldwork, &
3605 sphi_b(1, sgfb), ldsb, 0.0_dp, work(1, 1), ldwork)
3606
3607 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
3608 1.0_dp, sphi_a(1, sgfa), ldsa, &
3609 work(1, 1), ldwork, 1.0_dp, dblock(1, 1), SIZE(dblock, 1))
3610 END DO
3611 END DO
3612 END DO
3613 END associate
3614 END DO
3615 CALL neighbor_list_iterator_release(nl_iterator)
3616 CALL kpoint_release(kpoints_all)
3617 DEALLOCATE (dipab, work, basis_set_list)
3618 CALL timestop(handle)
3619
3620 END SUBROUTINE build_local_moment_matrix_rs_img
3621
3622! **************************************************************************************************
3623!> \brief Calculates the dipole moments and berry curvature for periodic systems for kpoints
3624!> \param qs_env ...
3625!> \param xkp list of kpoints
3626!> \param dipole ...
3627!> \param rcc coordinates about which to calculate the dipole
3628!> \param berry_c berry curvature calculated using Ω^γ_n = Σ_m 2*Im[d^α_nm (d^β_mn)*]
3629!> \param do_parallel option to distribute the result in dipole across
3630!> different MPI ranks
3631!> \author Shridhar Shanbhag
3632! **************************************************************************************************
3633 SUBROUTINE qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, do_parallel)
3634 TYPE(qs_environment_type), POINTER :: qs_env
3635 LOGICAL, OPTIONAL :: do_parallel
3636 LOGICAL :: my_do_parallel, calc_bc
3637 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints_deep'
3638 COMPLEX(KIND=dp) :: phase, tmp_max
3639 COMPLEX(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: c_k, h_k, s_k, d_k, cdc, c_dh_c, &
3640 c_ds_c, dh_dk_i, ds_dk_i
3641 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dip
3642 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3643 ALLOCATABLE :: dipole
3644 INTEGER :: handle, i_dir, ikp, nkp, &
3645 n_img_scf, n_img_all, nao, &
3646 num_pe, num_copy, mepos, n, m, mu, &
3647 ispin, nspin
3648 INTEGER, DIMENSION(3) :: periodic
3649 INTEGER, DIMENSION(:, :), POINTER :: index_to_cell_all
3650 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index_all
3651 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3652 REAL(kind=dp), DIMENSION(3) :: my_rcc
3653 REAL(kind=dp), DIMENSION(3, 3) :: hmat
3654 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: eigenvals
3655 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc, xkp
3656 REAL(kind=dp), DIMENSION(:, :, :, :), &
3657 ALLOCATABLE, OPTIONAL :: berry_c
3658 REAL(kind=dp), DIMENSION(:, :, :, :), &
3659 ALLOCATABLE :: d_rs, h_rs, s_rs
3660 TYPE(cell_type), POINTER :: cell
3661 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_rs_img, matrix_ks_kp, &
3662 matrix_s_kp
3663 TYPE(dft_control_type), POINTER :: dft_control
3664 TYPE(kpoint_type), POINTER :: kpoints_all, kpoints_scf
3665 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3666 TYPE(mp_para_env_type), POINTER :: para_env
3667 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3668 POINTER :: sab_all
3669
3670 CALL timeset(routinen, handle)
3671 calc_bc = PRESENT(berry_c)
3672 my_do_parallel = .false.
3673 my_rcc = 0.0_dp
3674 IF (PRESENT(do_parallel)) my_do_parallel = do_parallel
3675 IF (PRESENT(rcc)) my_rcc = rcc
3676
3677 CALL get_qs_env(qs_env, &
3678 matrix_ks_kp=matrix_ks_kp, &
3679 matrix_s_kp=matrix_s_kp, &
3680 sab_all=sab_all, &
3681 cell=cell, &
3682 kpoints=kpoints_scf, &
3683 para_env=para_env, &
3684 dft_control=dft_control, &
3685 mos=mos)
3686
3687 CALL get_mo_set(mo_set=mos(1), nao=nao)
3688 CALL get_cell(cell=cell, h=hmat, periodic=periodic)
3689 nspin = SIZE(matrix_ks_kp, 1)
3690 nkp = SIZE(xkp, 2)
3691
3692 ! create kpoint environment kpoints_all which contains all neighbor cells R
3693 ! without considering any lattice symmetry
3694 NULLIFY (kpoints_all)
3695 CALL kpoint_create(kpoints_all)
3696 CALL kpoint_init_cell_index(kpoints_all, sab_all, para_env, n_img_scf)
3697
3698 CALL get_kpoint_info(kpoints_all, cell_to_index=cell_to_index_all, &
3699 index_to_cell=index_to_cell_all)
3700 n_img_all = SIZE(index_to_cell_all, 2)
3701
3702 NULLIFY (moments_rs_img)
3703 CALL dbcsr_allocate_matrix_set(moments_rs_img, 3, n_img_all)
3704 ! D_μ,ν = <φ_μ|r|φ_ν>
3705 CALL build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc=my_rcc)
3706
3707 ALLOCATE (s_rs(1, nao, nao, n_img_all), h_rs(nspin, nao, nao, n_img_all), source=0.0_dp)
3708 ALLOCATE (d_rs(3, nao, nao, n_img_all), source=0.0_dp)
3709
3710 ! Convert real-space dbcsr matrices into arrays
3711 CALL replicate_rs_matrices(matrix_s_kp, kpoints_scf, s_rs, cell_to_index_all)
3712 CALL replicate_rs_matrices(matrix_ks_kp, kpoints_scf, h_rs, cell_to_index_all)
3713 CALL replicate_rs_matrices(moments_rs_img, kpoints_all, d_rs, cell_to_index_all)
3714
3715 mepos = 0
3716 num_pe = 1
3717 num_copy = nkp
3718 IF (my_do_parallel) THEN
3719 mepos = para_env%mepos
3720 num_pe = para_env%num_pe
3721 num_copy = ceiling(real(nkp)/num_pe)
3722 END IF
3723
3724 ALLOCATE (dipole(nspin, num_copy, 3, nao, nao), source=z_zero)
3725 IF (calc_bc) ALLOCATE (berry_c(nspin, num_copy, 3, nao), source=0.0_dp)
3726
3727!$OMP PARALLEL DEFAULT(NONE) PRIVATE(ikp, S_k, H_k, eigenvals, C_k, ispin, n, m, &
3728!$OMP i_dir, dS_dk_i, dH_dk_i, D_k, dip, bc, C_dS_C, C_dH_C, CDC, tmp_max, phase) &
3729!$OMP SHARED(num_pe, mepos, dipole, berry_c, nao, nspin, periodic, &
3730!$OMP nkp, xkp, S_rs, H_rs, D_rs, index_to_cell_all, hmat, calc_bc)
3731 ALLOCATE (ds_dk_i(nao, nao), c_ds_c(nao, nao), dh_dk_i(nao, nao), c_dh_c(nao, nao), source=z_zero)
3732 ALLOCATE (cdc(nao, nao), dip(3, nao, nao), s_k(nao, nao), h_k(nao, nao), source=z_zero)
3733 ALLOCATE (c_k(nao, nao), d_k(nao, nao), source=z_zero)
3734 ALLOCATE (eigenvals(nao), source=0.0_dp)
3735 IF (calc_bc) ALLOCATE (bc(3, nao), source=0.0_dp)
3736!$OMP DO COLLAPSE(2)
3737 DO ispin = 1, nspin
3738 DO ikp = 1, nkp
3739 IF (mod(ikp - 1, num_pe) /= mepos) cycle
3740
3741 ! S^R -> S(k), H^R -> H(k)
3742 s_k = 0
3743 h_k = 0
3744 CALL rs_to_kp(s_rs(1, :, :, :), s_k, index_to_cell_all, xkp(:, ikp))
3745 CALL rs_to_kp(h_rs(ispin, :, :, :), h_k, index_to_cell_all, xkp(:, ikp))
3746
3747 ! Diagonalize H(k)C(k) = S(k)C(k)ε(k)
3748 CALL geeig_right(h_k, s_k, eigenvals, c_k)
3749
3750 ! To have a smooth complex phase of C(k) as function of k, for every n, we force
3751 ! the largest C_μ,n(k) to be real.
3752 ! This is important to have a continuous dipole moment d_nm(k) as a function of k
3753 DO n = 1, nao
3754 tmp_max = c_k(1, n)
3755 DO mu = 1, nao
3756 IF (abs(c_k(mu, n)) < abs(tmp_max)) cycle
3757 tmp_max = c_k(mu, n)
3758 END DO
3759 phase = tmp_max/abs(tmp_max)
3760 c_k(:, n) = c_k(:, n)/phase
3761 END DO
3762
3763 DO i_dir = 1, 3 ! d^x, d^y, d^z
3764
3765 IF (periodic(i_dir) == 0) cycle
3766 ! ∇ S(k) = Σ_R iR S^R e^(ikR), ∇ H(k) = Σ_R iR H^R e^(ikR)
3767 CALL rs_to_kp(s_rs(1, :, :, :), ds_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3768 CALL rs_to_kp(h_rs(ispin, :, :, :), dh_dk_i, index_to_cell_all, xkp(:, ikp), i_dir, hmat)
3769
3770 ! Σ_R D^R e^(ikR) = D(k), D_μ,ν = <φ_μ|r|φ_ν>
3771 CALL rs_to_kp(d_rs(i_dir, :, :, :), d_k(:, :), index_to_cell_all, xkp(:, ikp))
3772
3773 ! Basis transform to Kohn-Sham basis: (C^H) ∇ S C, (C^H) ∇ H C, (C^H) D C
3774 CALL gemm_square(c_k, 'C', ds_dk_i, 'N', c_k, 'N', c_ds_c)
3775 CALL gemm_square(c_k, 'C', dh_dk_i, 'N', c_k, 'N', c_dh_c)
3776 CALL gemm_square(c_k, 'C', d_k, 'N', c_k, 'N', cdc)
3777
3778 ! Compute the dipole
3779 ! d_nm (k) = - i/(ε(n)-ε(m)) [ (C^H)(dH(k)/dk)C ]_nm
3780 ! + i ε(n)/(ε(n)-ε(m)) [ (C^H)(dS(k)/dk)C ]_nm + [ (C^H)D(k)C ]_nm
3781 DO n = 1, nao
3782 DO m = 1, nao
3783 IF (n == m) cycle ! diagonal elements would need to be computed from
3784 ! a numerical k-derivative which is not implemented
3785 dip(i_dir, n, m) = -gaussi*c_dh_c(n, m)/(eigenvals(n) - eigenvals(m)) &
3786 + gaussi*eigenvals(n)*c_ds_c(n, m)/(eigenvals(n) - eigenvals(m)) &
3787 + cdc(n, m)
3788 END DO
3789 END DO
3790 END DO
3791 ! Compute the Berry curvature from the dipoles
3792 ! Ω^γ_n = Σ_m 2*Im[d^α_nm d^β_mn], where, α, β, γ belong to {x, y, z}
3793 IF (calc_bc) THEN
3794 bc = 0.0_dp
3795 DO i_dir = 1, 3
3796 DO n = 1, nao
3797 DO m = 1, nao
3798 IF (n == m) cycle
3799 bc(i_dir, n) = bc(i_dir, n) &
3800 + 2*aimag(dip(1 + mod(i_dir, 3), n, m)*dip(1 + mod(i_dir + 1, 3), m, n))
3801 END DO
3802 END DO
3803 END DO
3804 END IF
3805 ! Store the dipoles and berry curvature for each MPI rank
3806 dipole(ispin, ceiling(real(ikp)/num_pe), :, :, :) = dip(:, :, :)
3807 IF (calc_bc) berry_c(ispin, ceiling(real(ikp)/num_pe), :, :) = bc(:, :)
3808 END DO
3809 END DO
3810!$OMP END DO
3811 DEALLOCATE (ds_dk_i, c_ds_c, dh_dk_i, c_dh_c, cdc, dip, s_k, h_k, c_k, d_k, eigenvals)
3812 IF (calc_bc) DEALLOCATE (bc)
3813!$OMP END PARALLEL
3814 DEALLOCATE (s_rs, h_rs, d_rs)
3815 CALL dbcsr_deallocate_matrix_set(moments_rs_img)
3816 CALL kpoint_release(kpoints_all)
3817 CALL timestop(handle)
3818 END SUBROUTINE qs_moment_kpoints_deep
3819
3820! **************************************************************************************************
3821!> \brief Calculate and print dipole moment elements d_nm(k) for k-point calculations
3822!> \param qs_env ...
3823!> \param nmoments ...
3824!> \param reference ...
3825!> \param ref_point ...
3826!> \param max_nmo ...
3827!> \param unit_number ...
3828!> \author Shridhar Shanbhag
3829! **************************************************************************************************
3830 SUBROUTINE qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number)
3831 TYPE(qs_environment_type), POINTER :: qs_env
3832 INTEGER, INTENT(IN) :: nmoments, reference, max_nmo
3833 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
3834 INTEGER, INTENT(IN) :: unit_number
3835 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_kpoints'
3836 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
3837 COMPLEX(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: dipole_to_print
3838 COMPLEX(KIND=dp), DIMENSION(:, :, :, :, :), &
3839 ALLOCATABLE :: dipole
3840 INTEGER :: handle, ikp, nkp, nao, &
3841 num_pe, mepos, n, m, &
3842 ispin, nspin, nmin, nmax, homo
3843 REAL(kind=dp), DIMENSION(3) :: rcc
3844 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: xkp
3845 REAL(kind=dp), DIMENSION(:, :), ALLOCATABLE :: bc_to_print
3846 REAL(kind=dp), DIMENSION(:, :, :, :), ALLOCATABLE :: berry_c
3847 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
3848 TYPE(mp_para_env_type), POINTER :: para_env
3849 CHARACTER(LEN=default_string_length), &
3850 DIMENSION(:), POINTER :: special_pnts
3851
3852 CALL timeset(routinen, handle)
3853
3854 IF (nmoments > 1) cpabort("KPOINT quadrupole and higher moments not implemented.")
3855 IF (max_nmo < 0) cpabort("Negative maximum number of molecular orbitals max_nmo provided.")
3856
3857 CALL get_qs_env(qs_env, &
3858 para_env=para_env, &
3859 matrix_ks_kp=matrix_ks_kp, &
3860 mos=mos)
3861
3862 CALL get_mo_set(mo_set=mos(1), nao=nao)
3863 CALL get_xkp_for_dipole_calc(qs_env, xkp, special_pnts)
3864 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
3865 nspin = SIZE(matrix_ks_kp, 1)
3866 nkp = SIZE(xkp, 2)
3867
3868 IF (unit_number > 0) WRITE (unit_number, fmt="(/,T2,A)") &
3869 '!-----------------------------------------------------------------------------!'
3870 IF (unit_number > 0) WRITE (unit_number, "(T22,A)") "Periodic Dipole Matrix Elements"
3871
3872 CALL qs_moment_kpoints_deep(qs_env, &
3873 xkp, &
3874 dipole, &
3875 rcc, &
3876 berry_c, &
3877 do_parallel=.true.)
3878
3879 ALLOCATE (dipole_to_print(3, nao, nao), source=z_zero)
3880 ALLOCATE (bc_to_print(3, nao), source=0.0_dp)
3881
3882 mepos = para_env%mepos
3883 num_pe = para_env%num_pe
3884
3885 DO ikp = 1, nkp
3886 DO ispin = 1, nspin
3887 CALL get_mo_set(mo_set=mos(ispin), homo=homo)
3888 nmin = max(1, homo - (max_nmo - 1)/2)
3889 nmax = min(nao, homo + max_nmo/2)
3890 IF (max_nmo == 0) THEN
3891 nmin = 1
3892 nmax = nao
3893 END IF
3894 dipole_to_print = 0.0_dp
3895 bc_to_print = 0.0_dp
3896 IF (mod(ikp - 1, num_pe) == mepos) THEN
3897 dipole_to_print(:, :, :) = dipole(ispin, ceiling(real(ikp)/num_pe), :, :, :)
3898 bc_to_print(:, :) = berry_c(ispin, ceiling(real(ikp)/num_pe), :, :)
3899 END IF
3900 CALL para_env%sum(dipole_to_print)
3901 CALL para_env%sum(bc_to_print)
3902 IF (unit_number > 0) THEN
3903 IF (special_pnts(ikp) /= "") WRITE (unit_number, "(/,2X,A,A)") &
3904 "Special point: ", adjustl(trim(special_pnts(ikp)))
3905 WRITE (unit_number, "(/,1X,A,I3,1X,3(A,1F12.6))") &
3906 "Kpoint:", ikp, ", kx:", xkp(1, ikp), ", ky:", xkp(2, ikp), ", kz:", xkp(3, ikp)
3907 IF (nspin > 1) WRITE (unit_number, "(/,2X,A,I2)") "Open Shell System. Spin:", ispin
3908 WRITE (unit_number, "(2X,A)") " kp n m Re(dx_nm) Im(dx_nm) &
3909 & Re(dy_nm) Im(dy_nm) Re(dz_nm) Im(dz_nm)"
3910 DO n = nmin, nmax
3911 DO m = nmin, nmax
3912 IF (n == m) cycle
3913 WRITE (unit_number, "(2X,I4,2I4,6(G11.3))") ikp, n, m, dipole_to_print(1:3, n, m)
3914 END DO
3915 END DO
3916 WRITE (unit_number, "(/,1X,A)") "Berry Curvature"
3917 WRITE (unit_number, "(2X,A)") " kp n YZ ZX XY"
3918 DO n = nmin, nmax
3919 WRITE (unit_number, "(2X,2I5,3(1X,G11.3))") &
3920 ikp, n, bc_to_print(1, n), bc_to_print(2, n), bc_to_print(3, n)
3921 END DO
3922 END IF
3923 END DO
3924 END DO
3925 DEALLOCATE (dipole_to_print, bc_to_print, berry_c, dipole)
3926 DEALLOCATE (special_pnts, xkp)
3927
3928 CALL timestop(handle)
3929
3930 END SUBROUTINE qs_moment_kpoints
3931
3932END MODULE qs_moments
ncoset
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
Definition grid_common.h:81
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:218
cell_types::pbc
Definition cell_types.F:97
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
cp_result_methods::put_results
Definition cp_result_methods.F:39
mathlib::gemm_square
Definition mathlib.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
ai_angmom
Calculation of the angular momentum integrals over Cartesian Gaussian-type functions.
Definition ai_angmom.F:17
ai_angmom::angmom
subroutine, public angmom(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, rac, rbc, angab)
...
Definition ai_angmom.F:52
ai_moments
Calculation of the moment integrals over Cartesian Gaussian-type functions.
Definition ai_moments.F:17
ai_moments::diff_momop
subroutine, public diff_momop(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the position of the pr...
Definition ai_moments.F:1583
ai_moments::diff_momop_velocity
subroutine, public diff_momop_velocity(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, lambda, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the primitive on the r...
Definition ai_moments.F:79
ai_moments::contract_cossin
subroutine, public contract_cossin(cos_block, sin_block, iatom, ncoa, nsgfa, sgfa, sphi_a, ldsa, jatom, ncob, nsgfb, sgfb, sphi_b, ldsb, cosab, sinab, ldab, work, ldwork)
...
Definition ai_moments.F:261
ai_moments::cossin
subroutine, public cossin(la_max_set, npgfa, zeta, rpgfa, la_min_set, lb_max, npgfb, zetb, rpgfb, lb_min, rac, rbc, kvec, cosab, sinab, dcosab, dsinab)
...
Definition ai_moments.F:339
ai_moments::moment
subroutine, public moment(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lc_max, rac, rbc, mab)
...
Definition ai_moments.F:986
ai_moments::diff_momop2
subroutine, public diff_momop2(la_max, npgfa, zeta, rpgfa, la_min, lb_max, npgfb, zetb, rpgfb, lb_min, order, rac, rbc, difmab, mab_ext, deltar, iatom, jatom)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the position of the pr...
Definition ai_moments.F:164
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
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.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::mattiat2019
integer, save, public mattiat2019
Definition bibliography.F:36
block_p_types
collect pointers to a block of reals
Definition block_p_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:200
commutator_rpnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition commutator_rpnl.F:14
commutator_rpnl::build_com_mom_nl
subroutine, public build_com_mom_nl(qs_kind_set, sab_all, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr, matrix_r_rxvr, matrix_rxvr_r, matrix_r_doublecom, pseudoatom, ref_point)
Calculate [r,Vnl] (matrix_rv), r x [r,Vnl] (matrix_rxrv) or [rr,Vnl] (matrix_rrv) in AO basis....
Definition commutator_rpnl.F:439
cp_cfm_basic_linalg
Basic linear algebra operations for complex full matrices.
Definition cp_cfm_basic_linalg.F:16
cp_cfm_basic_linalg::cp_cfm_det
subroutine, public cp_cfm_det(matrix_a, det_a)
Computes the determinant (with a correct sign even in parallel environment!) of a complex square matr...
Definition cp_cfm_basic_linalg.F:75
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:630
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:692
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1086
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1194
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_trace
subroutine, public dbcsr_trace(matrix, trace)
Computes the trace of the given matrix, also known as the sum of its diagonal elements.
Definition cp_dbcsr_contrib.F:412
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:552
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_double
subroutine, public cp_fm_struct_double(fmstruct, struct, context, col, row)
creates a struct with twice the number of blocks on each core. If matrix A has to be multiplied with ...
Definition cp_fm_struct.F:523
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1087
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:169
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:556
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::read_float_object
elemental subroutine, public read_float_object(string, object, error_message)
Returns a floating point number read from a string including fraction like z1/z2.
Definition cp_parser_methods.F:1258
cp_result_methods
set of type/routines to handle the storage of results in force_envs
Definition cp_result_methods.F:17
cp_result_methods::cp_results_erase
subroutine, public cp_results_erase(results, description, nval)
erase a part of result_list
Definition cp_result_methods.F:362
cp_result_types
set of type/routines to handle the storage of results in force_envs
Definition cp_result_types.F:18
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:731
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:704
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1047
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::max_line_length
integer, parameter, public max_line_length
Definition kinds.F:59
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_k_r_trafo_simple
Implements transformations from k-space to R-space for Fortran array matrices.
Definition kpoint_k_r_trafo_simple.F:16
kpoint_k_r_trafo_simple::rs_to_kp
subroutine, public rs_to_kp(rs_real, ks_complex, index_to_cell, xkp, deriv_direction, hmat)
Integrate RS matrices (stored as Fortran array) into a kpoint matrix at given kp.
Definition kpoint_k_r_trafo_simple.F:230
kpoint_k_r_trafo_simple::replicate_rs_matrices
subroutine, public replicate_rs_matrices(rs_dbcsr_in, kpoint_in, rs_array_out, cell_to_index_out)
Convert dbcsr matrices representing operators in real-space image cells to arrays.
Definition kpoint_k_r_trafo_simple.F:62
kpoint_methods
Routines needed for kpoint calculation.
Definition kpoint_methods.F:15
kpoint_methods::kpoint_init_cell_index
subroutine, public kpoint_init_cell_index(kpoint, sab_nl, para_env, nimages)
Generates the mapping of cell indices and linear RS index CELL (0,0,0) is always mapped to index 1.
Definition kpoint_methods.F:847
kpoint_methods::kpoint_initialize
subroutine, public kpoint_initialize(kpoint, particle_set, cell)
Generate the kpoints and initialize the kpoint environment.
Definition kpoint_methods.F:119
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::read_kpoint_section
subroutine, public read_kpoint_section(kpoint, kpoint_section, a_vec)
Read the kpoint input section.
Definition kpoint_types.F:724
kpoint_types::get_kpoint_info
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.
Definition kpoint_types.F:522
kpoint_types::kpoint_release
subroutine, public kpoint_release(kpoint)
Release a kpoint environment, deallocate all data.
Definition kpoint_types.F:281
kpoint_types::kpoint_create
subroutine, public kpoint_create(kpoint)
Create a kpoint environment.
Definition kpoint_types.F:230
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::gaussi
complex(kind=dp), parameter, public gaussi
Definition mathconstants.F:142
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::z_zero
complex(kind=dp), parameter, public z_zero
Definition mathconstants.F:143
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::geeig_right
subroutine, public geeig_right(a_in, b_in, eigenvalues, eigenvectors)
Solve the generalized eigenvalue equation for complex matrices A*v = B*v*λ
Definition mathlib.F:2037
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
moments_utils
Calculates the moment integrals <a|r^m|b>
Definition moments_utils.F:15
moments_utils::get_reference_point
subroutine, public get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, ifirst, ilast)
...
Definition moments_utils.F:61
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::current_maxl
integer, save, public current_maxl
Definition orbital_pointers.F:40
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition particle_methods.F:120
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
physcon::debye
real(kind=dp), parameter, public debye
Definition physcon.F:201
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition qs_environment_types.F:530
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
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.
Definition qs_kind_types.F:468
qs_kind_types::get_qs_kind_set
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.
Definition qs_kind_types.F:913
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:341
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:399
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_magmom_matrix
subroutine, public build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:1150
qs_moments::dipole_velocity_deriv
subroutine, public dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
This returns the derivative of the moment integrals [a|\mu|b], with respect to the basis function on ...
Definition qs_moments.F:153
qs_moments::print_moments_nl
subroutine, public print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
...
Definition qs_moments.F:2798
qs_moments::set_label
subroutine, public set_label(label, ix, iy, iz)
...
Definition qs_moments.F:2659
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition qs_moments.F:585
qs_moments::qs_moment_locop
subroutine, public qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
...
Definition qs_moments.F:2258
qs_moments::print_moments
subroutine, public print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
...
Definition qs_moments.F:2690
qs_moments::build_berry_moment_matrix
subroutine, public build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
...
Definition qs_moments.F:1363
qs_moments::op_orbbas
subroutine, public op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
...
Definition qs_moments.F:2136
qs_moments::build_local_moment_matrix_rs_img
subroutine, public build_local_moment_matrix_rs_img(qs_env, moments_rs_img, rcc)
Calculate local moment matrix for a periodic system for all image cells.
Definition qs_moments.F:3455
qs_moments::build_local_moments_der_matrix
subroutine, public build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, ref_point, moments)
Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b > Optionally sto...
Definition qs_moments.F:817
qs_moments::calculate_commutator_nl_terms
subroutine, public calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, ref_point)
Calculate the expectation value of operators related to non-local potential: [r, Vnl],...
Definition qs_moments.F:2866
qs_moments::qs_moment_kpoints_deep
subroutine, public qs_moment_kpoints_deep(qs_env, xkp, dipole, rcc, berry_c, do_parallel)
Calculates the dipole moments and berry curvature for periodic systems for kpoints.
Definition qs_moments.F:3634
qs_moments::qs_moment_kpoints
subroutine, public qs_moment_kpoints(qs_env, nmoments, reference, ref_point, max_nmo, unit_number)
Calculate and print dipole moment elements d_nm(k) for k-point calculations.
Definition qs_moments.F:3831
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
...
Definition qs_moments.F:1530
qs_moments::op_orbbas_rtp
subroutine, public op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
...
Definition qs_moments.F:2166
qs_moments::build_dsdv_moments
subroutine, public build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
Builds the moments for the derivative of the overlap with respect to nuclear velocities.
Definition qs_moments.F:356
qs_moments::qs_moment_berry_phase
subroutine, public qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
...
Definition qs_moments.F:1737
qs_moments::dipole_deriv_ao
subroutine, public dipole_deriv_ao(qs_env, difdip, deltar, order, rcc)
...
Definition qs_moments.F:3121
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition qs_operators_ao.F:74
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
rt_propagation_types
Types and set_get for real time propagation depending on runtype and diagonalization method different...
Definition rt_propagation_types.F:21
rt_propagation_types::get_rtp
subroutine, public get_rtp(rtp, exp_h_old, exp_h_new, h_last_iter, rho_old, rho_next, rho_new, mos, mos_new, mos_old, mos_next, s_inv, s_half, s_minus_half, b_mat, c_mat, propagator_matrix, mixing, mixing_factor, s_der, dt, nsteps, sinvh, sinvh_imag, sinvb, admm_mos)
...
Definition rt_propagation_types.F:416
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
block_p_types::block_p_type
Definition block_p_types.F:23
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::dft_control_type
Definition cp_control_types.F:726
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
cp_result_types::cp_result_type
contains arbitrary information which need to be stored
Definition cp_result_types.F:73
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:168
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:743
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:220
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
rt_propagation_types::rt_prop_type
Definition rt_propagation_types.F:76