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