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