(git:ed6f26b)
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-2025 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 USE commutator_rpnl, ONLY: build_com_mom_nl
32 USE cp_cfm_basic_linalg, ONLY: cp_cfm_det
33 USE cp_cfm_types, ONLY: cp_cfm_create,&
34 cp_cfm_get_info,&
35 cp_cfm_release,&
36 cp_cfm_type
37 USE cp_control_types, ONLY: dft_control_type
38 USE cp_dbcsr_api, ONLY: &
39 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, &
40 dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_type, &
41 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
42 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
43 dbcsr_trace
44 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
45 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
46 dbcsr_allocate_matrix_set,&
47 dbcsr_deallocate_matrix_set
48 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
49 cp_fm_struct_double,&
50 cp_fm_struct_release,&
51 cp_fm_struct_type
52 USE cp_fm_types, ONLY: cp_fm_create,&
53 cp_fm_get_info,&
54 cp_fm_release,&
55 cp_fm_set_all,&
56 cp_fm_type
57 USE cp_result_methods, ONLY: cp_results_erase,&
58 put_results
59 USE cp_result_types, ONLY: cp_result_type
60 USE distribution_1d_types, ONLY: distribution_1d_type
61 USE kinds, ONLY: default_string_length,&
62 dp
63 USE kpoint_types, ONLY: get_kpoint_info,&
64 kpoint_type
65 USE mathconstants, ONLY: pi,&
66 twopi
67 USE message_passing, ONLY: mp_para_env_type
68 USE moments_utils, ONLY: get_reference_point
69 USE orbital_pointers, ONLY: current_maxl,&
70 indco,&
71 ncoset
72 USE parallel_gemm_api, ONLY: parallel_gemm
73 USE particle_methods, ONLY: get_particle_set
74 USE particle_types, ONLY: particle_type
75 USE physcon, ONLY: bohr,&
76 debye
77 USE qs_environment_types, ONLY: get_qs_env,&
78 qs_environment_type
79 USE qs_kind_types, ONLY: get_qs_kind,&
80 get_qs_kind_set,&
81 qs_kind_type
82 USE qs_ks_types, ONLY: get_ks_env,&
83 qs_ks_env_type
84 USE qs_mo_types, ONLY: get_mo_set,&
85 mo_set_type
86 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
87 neighbor_list_iterate,&
88 neighbor_list_iterator_create,&
89 neighbor_list_iterator_p_type,&
90 neighbor_list_iterator_release,&
91 neighbor_list_set_p_type
92 USE qs_operators_ao, ONLY: build_lin_mom_matrix
93 USE qs_rho_types, ONLY: qs_rho_get,&
94 qs_rho_type
95 USE rt_propagation_types, ONLY: get_rtp,&
96 rt_prop_type
97#include "./base/base_uses.f90"
98
99 IMPLICIT NONE
100
101 PRIVATE
102
103 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
104
105 ! Public subroutines
106 PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
107 PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
108 PUBLIC :: qs_moment_berry_phase, qs_moment_locop
109 PUBLIC :: dipole_deriv_ao
110 PUBLIC :: build_local_moments_der_matrix
111 PUBLIC :: build_dsdv_moments
112 PUBLIC :: dipole_velocity_deriv
113
114CONTAINS
115
116! *****************************************************************************
117!> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
118!> to the basis function on the right
119!> difdip(beta, alpha) = < mu | r_beta | ∂_alpha nu > * (mu - nu)
120!> \param qs_env ...
121!> \param difdip ...
122!> \param order The order of the derivative (1 for dipole moment)
123!> \param lambda The atom on which we take the derivative
124!> \param rc ...
125!> \author Edward Ditler
126! **************************************************************************************************
127 SUBROUTINE dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
128 TYPE(qs_environment_type), POINTER :: qs_env
129 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: difdip
130 INTEGER, INTENT(IN) :: order, lambda
131 REAL(kind=dp), DIMENSION(3) :: rc
132
133 CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_velocity_deriv'
134
135 INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
136 last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
137 sgfb
138 LOGICAL :: found
139 REAL(dp) :: dab
140 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
141 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
142 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
143 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab, difmab2
144 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
145 TYPE(cell_type), POINTER :: cell
146 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
147 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
148 TYPE(neighbor_list_iterator_p_type), &
149 DIMENSION(:), POINTER :: nl_iterator
150 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
151 POINTER :: sab_all
152 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
153 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
154 TYPE(qs_kind_type), POINTER :: qs_kind
155
156 CALL timeset(routinen, handle)
157
158 NULLIFY (cell, particle_set, qs_kind_set, sab_all)
159 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
160 qs_kind_set=qs_kind_set, sab_all=sab_all)
161 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
162 maxco=ldab, maxsgf=maxsgf)
163
164 nkind = SIZE(qs_kind_set)
165 natom = SIZE(particle_set)
166
167 m_dim = ncoset(order) - 1
168
169 ALLOCATE (basis_set_list(nkind))
170
171 ALLOCATE (mab(ldab, ldab, m_dim))
172 ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
173 ALLOCATE (work(ldab, maxsgf))
174 ALLOCATE (mint(3, 3))
175 ALLOCATE (mint2(3, 3))
176
177 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
178 difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
179 work(1:ldab, 1:maxsgf) = 0.0_dp
180
181 DO i = 1, 3
182 DO j = 1, 3
183 NULLIFY (mint(i, j)%block)
184 NULLIFY (mint2(i, j)%block)
185 END DO
186 END DO
187
188 ! Set the basis_set_list(nkind) to point to the corresponding basis sets
189 DO ikind = 1, nkind
190 qs_kind => qs_kind_set(ikind)
191 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
192 IF (ASSOCIATED(basis_set_a)) THEN
193 basis_set_list(ikind)%gto_basis_set => basis_set_a
194 ELSE
195 NULLIFY (basis_set_list(ikind)%gto_basis_set)
196 END IF
197 END DO
198
199 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
200 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
201 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
202 iatom=iatom, jatom=jatom, r=rab)
203
204 basis_set_a => basis_set_list(ikind)%gto_basis_set
205 basis_set_b => basis_set_list(jkind)%gto_basis_set
206 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
207 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
208
209 associate( &
210 ! basis ikind
211 first_sgfa => basis_set_a%first_sgf, &
212 la_max => basis_set_a%lmax, &
213 la_min => basis_set_a%lmin, &
214 npgfa => basis_set_a%npgf, &
215 nsgfa => basis_set_a%nsgf_set, &
216 rpgfa => basis_set_a%pgf_radius, &
217 set_radius_a => basis_set_a%set_radius, &
218 sphi_a => basis_set_a%sphi, &
219 zeta => basis_set_a%zet, &
220 ! basis jkind, &
221 first_sgfb => basis_set_b%first_sgf, &
222 lb_max => basis_set_b%lmax, &
223 lb_min => basis_set_b%lmin, &
224 npgfb => basis_set_b%npgf, &
225 nsgfb => basis_set_b%nsgf_set, &
226 rpgfb => basis_set_b%pgf_radius, &
227 set_radius_b => basis_set_b%set_radius, &
228 sphi_b => basis_set_b%sphi, &
229 zetb => basis_set_b%zet)
230
231 nseta = basis_set_a%nset
232 nsetb = basis_set_b%nset
233
234 IF (inode == 1) last_jatom = 0
235
236 ! this guarantees minimum image convention
237 ! anything else would not make sense
238 IF (jatom == last_jatom) THEN
239 cycle
240 END IF
241
242 last_jatom = jatom
243
244 irow = iatom
245 icol = jatom
246
247 DO i = 1, 3
248 DO j = 1, 3
249 NULLIFY (mint(i, j)%block)
250 CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
251 row=irow, col=icol, block=mint(i, j)%block, &
252 found=found)
253 cpassert(found)
254 mint(i, j)%block = 0._dp
255 END DO
256 END DO
257
258 ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
259 ra = pbc(particle_set(iatom)%r(:), cell)
260 rb(:) = ra(:) + rab(:)
261 rac = pbc(rc, ra, cell)
262 rbc = rac + rab
263 dab = norm2(rab)
264
265 DO iset = 1, nseta
266 ncoa = npgfa(iset)*ncoset(la_max(iset))
267 sgfa = first_sgfa(1, iset)
268 DO jset = 1, nsetb
269 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
270 ncob = npgfb(jset)*ncoset(lb_max(jset))
271 sgfb = first_sgfb(1, jset)
272 ldab = max(ncoa, ncob)
273 lda = ncoset(la_max(iset))*npgfa(iset)
274 ldb = ncoset(lb_max(jset))*npgfb(jset)
275 ALLOCATE (difmab(lda, ldb, m_dim, 3))
276
277 ! Calculate integral difmab(beta, alpha) = (a|r_beta|db_alpha)
278 ! difmab(beta, alpha) = < a | r_beta | ∂_alpha b >
279 ! difmab(j, idir) = < a | r_j | ∂_idir b >
280 CALL diff_momop_velocity(la_max(iset), npgfa(iset), zeta(:, iset), &
281 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
282 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
283 difmab, lambda=lambda, iatom=iatom, jatom=jatom)
284
285 ! *** Contraction step ***
286
287 DO idir = 1, 3 ! derivative of AO function
288 DO j = 1, 3 ! position operator r_j
289 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
290 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
291 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
292 0.0_dp, work(1, 1), SIZE(work, 1))
293
294 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
295 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
296 work(1, 1), SIZE(work, 1), &
297 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
298 SIZE(mint(j, idir)%block, 1))
299 END DO !j
300 END DO !idir
301 DEALLOCATE (difmab)
302 END DO !jset
303 END DO !iset
304 END associate
305 END do!iterator
306
307 CALL neighbor_list_iterator_release(nl_iterator)
308
309 DO i = 1, 3
310 DO j = 1, 3
311 NULLIFY (mint(i, j)%block)
312 END DO
313 END DO
314
315 DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
316
317 CALL timestop(handle)
318 END SUBROUTINE dipole_velocity_deriv
319
320! **************************************************************************************************
321!> \brief Builds the moments for the derivative of the overlap with respect to nuclear velocities
322!> \param qs_env ...
323!> \param moments ...
324!> \param nmoments ...
325!> \param ref_point ...
326!> \param ref_points ...
327!> \param basis_type ...
328!> \author Edward Ditler
329! **************************************************************************************************
330 SUBROUTINE build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
331
332 TYPE(qs_environment_type), POINTER :: qs_env
333 TYPE(dbcsr_p_type), DIMENSION(:) :: moments
334 INTEGER, INTENT(IN) :: nmoments
335 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
336 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
337 OPTIONAL :: ref_points
338 CHARACTER(len=*), OPTIONAL :: basis_type
339
340 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_dsdv_moments'
341
342 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
343 maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
344 INTEGER, DIMENSION(3) :: image_cell
345 LOGICAL :: found
346 REAL(kind=dp) :: dab
347 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
348 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
349 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
350 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
351 TYPE(cell_type), POINTER :: cell
352 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
353 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
354 TYPE(neighbor_list_iterator_p_type), &
355 DIMENSION(:), POINTER :: nl_iterator
356 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
357 POINTER :: sab_orb
358 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
359 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
360 TYPE(qs_kind_type), POINTER :: qs_kind
361
362 IF (nmoments < 1) RETURN
363
364 CALL timeset(routinen, handle)
365
366 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
367
368 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
369 cpassert(SIZE(moments) >= nm)
370
371 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
372 CALL get_qs_env(qs_env=qs_env, &
373 qs_kind_set=qs_kind_set, &
374 particle_set=particle_set, cell=cell, &
375 sab_orb=sab_orb)
376
377 nkind = SIZE(qs_kind_set)
378 natom = SIZE(particle_set)
379
380 ! Allocate work storage
381 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
382 maxco=maxco, maxsgf=maxsgf, &
383 basis_type=basis_type)
384
385 ALLOCATE (mab(maxco, maxco, nm))
386 mab(:, :, :) = 0.0_dp
387
388 ALLOCATE (work(maxco, maxsgf))
389 work(:, :) = 0.0_dp
390
391 ALLOCATE (mint(nm))
392 DO i = 1, nm
393 NULLIFY (mint(i)%block)
394 END DO
395
396 ALLOCATE (basis_set_list(nkind))
397 DO ikind = 1, nkind
398 qs_kind => qs_kind_set(ikind)
399 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
400 IF (ASSOCIATED(basis_set_a)) THEN
401 basis_set_list(ikind)%gto_basis_set => basis_set_a
402 ELSE
403 NULLIFY (basis_set_list(ikind)%gto_basis_set)
404 END IF
405 END DO
406 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
407 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
408 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
409 iatom=iatom, jatom=jatom, r=rab, cell=image_cell)
410 basis_set_a => basis_set_list(ikind)%gto_basis_set
411 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
412 basis_set_b => basis_set_list(jkind)%gto_basis_set
413 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
414 ! basis ikind
415 associate( &
416 first_sgfa => basis_set_a%first_sgf, &
417 la_max => basis_set_a%lmax, &
418 la_min => basis_set_a%lmin, &
419 npgfa => basis_set_a%npgf, &
420 nsgfa => basis_set_a%nsgf_set, &
421 rpgfa => basis_set_a%pgf_radius, &
422 set_radius_a => basis_set_a%set_radius, &
423 sphi_a => basis_set_a%sphi, &
424 zeta => basis_set_a%zet, &
425 ! basis jkind, &
426 first_sgfb => basis_set_b%first_sgf, &
427 lb_max => basis_set_b%lmax, &
428 lb_min => basis_set_b%lmin, &
429 npgfb => basis_set_b%npgf, &
430 nsgfb => basis_set_b%nsgf_set, &
431 rpgfb => basis_set_b%pgf_radius, &
432 set_radius_b => basis_set_b%set_radius, &
433 sphi_b => basis_set_b%sphi, &
434 zetb => basis_set_b%zet)
435
436 nseta = basis_set_a%nset
437 nsetb = basis_set_b%nset
438
439 IF (inode == 1) last_jatom = 0
440
441 ! this guarantees minimum image convention
442 ! anything else would not make sense
443 IF (jatom == last_jatom) THEN
444 cycle
445 END IF
446
447 last_jatom = jatom
448
449 IF (iatom <= jatom) THEN
450 irow = iatom
451 icol = jatom
452 ELSE
453 irow = jatom
454 icol = iatom
455 END IF
456
457 DO i = 1, nm
458 NULLIFY (mint(i)%block)
459 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
460 row=irow, col=icol, block=mint(i)%block, found=found)
461 mint(i)%block = 0._dp
462 END DO
463
464 ! fold atomic position back into unit cell
465 IF (PRESENT(ref_points)) THEN
466 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
467 ELSE IF (PRESENT(ref_point)) THEN
468 rc(:) = ref_point(:)
469 ELSE
470 rc(:) = 0._dp
471 END IF
472 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
473 ! by folding around the center, such screwing can be avoided for a proper choice of center.
474 ! we dont use PBC at this point
475
476 ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
477 ra(:) = particle_set(iatom)%r(:)
478 rb(:) = ra(:) + rab(:)
479 rac = pbc(rc, ra, cell)
480 rbc = rac + rab
481
482 dab = norm2(rab)
483
484 DO iset = 1, nseta
485
486 ncoa = npgfa(iset)*ncoset(la_max(iset))
487 sgfa = first_sgfa(1, iset)
488
489 DO jset = 1, nsetb
490
491 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
492
493 ncob = npgfb(jset)*ncoset(lb_max(jset))
494 sgfb = first_sgfb(1, jset)
495
496 ! Calculate the primitive integrals
497 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
498 rpgfa(:, iset), la_min(iset), &
499 lb_max(jset), npgfb(jset), zetb(:, jset), &
500 rpgfb(:, jset), nmoments, rac, rbc, mab)
501
502 ! Contraction step
503 DO i = 1, nm
504
505 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
506 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
507 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
508 0.0_dp, work(1, 1), SIZE(work, 1))
509
510 IF (iatom <= jatom) THEN
511
512 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
513 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
514 work(1, 1), SIZE(work, 1), &
515 1.0_dp, mint(i)%block(sgfa, sgfb), &
516 SIZE(mint(i)%block, 1))
517
518 ELSE
519
520 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
521 1.0_dp, work(1, 1), SIZE(work, 1), &
522 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
523 1.0_dp, mint(i)%block(sgfb, sgfa), &
524 SIZE(mint(i)%block, 1))
525
526 END IF
527
528 END DO
529
530 END DO
531 END DO
532 END associate
533
534 END DO ! iterator
535
536 CALL neighbor_list_iterator_release(nl_iterator)
537
538 ! Release work storage
539 DEALLOCATE (mab, basis_set_list)
540 DEALLOCATE (work)
541 DO i = 1, nm
542 NULLIFY (mint(i)%block)
543 END DO
544 DEALLOCATE (mint)
545
546 CALL timestop(handle)
547
548 END SUBROUTINE build_dsdv_moments
549
550! **************************************************************************************************
551!> \brief ...
552!> \param qs_env ...
553!> \param moments ...
554!> \param nmoments ...
555!> \param ref_point ...
556!> \param ref_points ...
557!> \param basis_type ...
558! **************************************************************************************************
559 SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
560
561 TYPE(qs_environment_type), POINTER :: qs_env
562 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
563 INTEGER, INTENT(IN) :: nmoments
564 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
565 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
566 OPTIONAL :: ref_points
567 CHARACTER(len=*), OPTIONAL :: basis_type
568
569 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix'
570
571 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
572 maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
573 LOGICAL :: found
574 REAL(kind=dp) :: dab
575 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
576 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
577 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
578 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
579 TYPE(cell_type), POINTER :: cell
580 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
581 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
582 TYPE(neighbor_list_iterator_p_type), &
583 DIMENSION(:), POINTER :: nl_iterator
584 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
585 POINTER :: sab_orb
586 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
587 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
588 TYPE(qs_kind_type), POINTER :: qs_kind
589
590 IF (nmoments < 1) RETURN
591
592 CALL timeset(routinen, handle)
593
594 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
595
596 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
597 cpassert(SIZE(moments) >= nm)
598
599 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
600 CALL get_qs_env(qs_env=qs_env, &
601 qs_kind_set=qs_kind_set, &
602 particle_set=particle_set, cell=cell, &
603 sab_orb=sab_orb)
604
605 nkind = SIZE(qs_kind_set)
606 natom = SIZE(particle_set)
607
608 ! Allocate work storage
609 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
610 maxco=maxco, maxsgf=maxsgf, &
611 basis_type=basis_type)
612
613 ALLOCATE (mab(maxco, maxco, nm))
614 mab(:, :, :) = 0.0_dp
615
616 ALLOCATE (work(maxco, maxsgf))
617 work(:, :) = 0.0_dp
618
619 ALLOCATE (mint(nm))
620 DO i = 1, nm
621 NULLIFY (mint(i)%block)
622 END DO
623
624 ALLOCATE (basis_set_list(nkind))
625 DO ikind = 1, nkind
626 qs_kind => qs_kind_set(ikind)
627 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
628 IF (ASSOCIATED(basis_set_a)) THEN
629 basis_set_list(ikind)%gto_basis_set => basis_set_a
630 ELSE
631 NULLIFY (basis_set_list(ikind)%gto_basis_set)
632 END IF
633 END DO
634 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
635 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
636 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
637 iatom=iatom, jatom=jatom, r=rab)
638 basis_set_a => basis_set_list(ikind)%gto_basis_set
639 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
640 basis_set_b => basis_set_list(jkind)%gto_basis_set
641 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
642 associate( &
643 ! basis ikind
644 first_sgfa => basis_set_a%first_sgf, &
645 la_max => basis_set_a%lmax, &
646 la_min => basis_set_a%lmin, &
647 npgfa => basis_set_a%npgf, &
648 nsgfa => basis_set_a%nsgf_set, &
649 rpgfa => basis_set_a%pgf_radius, &
650 set_radius_a => basis_set_a%set_radius, &
651 sphi_a => basis_set_a%sphi, &
652 zeta => basis_set_a%zet, &
653 ! basis jkind, &
654 first_sgfb => basis_set_b%first_sgf, &
655 lb_max => basis_set_b%lmax, &
656 lb_min => basis_set_b%lmin, &
657 npgfb => basis_set_b%npgf, &
658 nsgfb => basis_set_b%nsgf_set, &
659 rpgfb => basis_set_b%pgf_radius, &
660 set_radius_b => basis_set_b%set_radius, &
661 sphi_b => basis_set_b%sphi, &
662 zetb => basis_set_b%zet)
663
664 nseta = basis_set_a%nset
665 nsetb = basis_set_b%nset
666
667 IF (inode == 1) last_jatom = 0
668
669 ! this guarantees minimum image convention
670 ! anything else would not make sense
671 IF (jatom == last_jatom) THEN
672 cycle
673 END IF
674
675 last_jatom = jatom
676
677 IF (iatom <= jatom) THEN
678 irow = iatom
679 icol = jatom
680 ELSE
681 irow = jatom
682 icol = iatom
683 END IF
684
685 DO i = 1, nm
686 NULLIFY (mint(i)%block)
687 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
688 row=irow, col=icol, block=mint(i)%block, found=found)
689 mint(i)%block = 0._dp
690 END DO
691
692 ! fold atomic position back into unit cell
693 IF (PRESENT(ref_points)) THEN
694 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
695 ELSE IF (PRESENT(ref_point)) THEN
696 rc(:) = ref_point(:)
697 ELSE
698 rc(:) = 0._dp
699 END IF
700 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
701 ! by folding around the center, such screwing can be avoided for a proper choice of center.
702 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
703 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
704 ! we dont use PBC at this point
705 rab(:) = ra(:) - rb(:)
706 rac(:) = ra(:) - rc(:)
707 rbc(:) = rb(:) - rc(:)
708 dab = norm2(rab)
709
710 DO iset = 1, nseta
711
712 ncoa = npgfa(iset)*ncoset(la_max(iset))
713 sgfa = first_sgfa(1, iset)
714
715 DO jset = 1, nsetb
716
717 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
718
719 ncob = npgfb(jset)*ncoset(lb_max(jset))
720 sgfb = first_sgfb(1, jset)
721
722 ! Calculate the primitive integrals
723 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
724 rpgfa(:, iset), la_min(iset), &
725 lb_max(jset), npgfb(jset), zetb(:, jset), &
726 rpgfb(:, jset), nmoments, rac, rbc, mab)
727
728 ! Contraction step
729 DO i = 1, nm
730
731 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
732 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
733 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
734 0.0_dp, work(1, 1), SIZE(work, 1))
735
736 IF (iatom <= jatom) THEN
737
738 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
739 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
740 work(1, 1), SIZE(work, 1), &
741 1.0_dp, mint(i)%block(sgfa, sgfb), &
742 SIZE(mint(i)%block, 1))
743
744 ELSE
745
746 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
747 1.0_dp, work(1, 1), SIZE(work, 1), &
748 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
749 1.0_dp, mint(i)%block(sgfb, sgfa), &
750 SIZE(mint(i)%block, 1))
751
752 END IF
753
754 END DO
755
756 END DO
757 END DO
758 END associate
759
760 END DO
761 CALL neighbor_list_iterator_release(nl_iterator)
762
763 ! Release work storage
764 DEALLOCATE (mab, basis_set_list)
765 DEALLOCATE (work)
766 DO i = 1, nm
767 NULLIFY (mint(i)%block)
768 END DO
769 DEALLOCATE (mint)
770
771 CALL timestop(handle)
772
773 END SUBROUTINE build_local_moment_matrix
774
775! **************************************************************************************************
776!> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
777!> Optionally stores the multipole moments themselves for free.
778!> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
779!> Only first derivatives are performed, e. g. x d/dy
780!> \param qs_env ...
781!> \param moments_der will contain the derivatives of the multipole moments
782!> \param nmoments_der order of the moments with derivatives
783!> \param nmoments order of the multipole moments (no derivatives, same output as
784!> build_local_moment_matrix, needs moments as arguments to store results)
785!> \param ref_point ...
786!> \param moments contains the multipole moments, optionally for free, up to order nmoments
787!> \note
788!> Adapted from rRc_xyz_der_ao in qs_operators_ao
789! **************************************************************************************************
790 SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
791 ref_point, moments)
792 TYPE(qs_environment_type), POINTER :: qs_env
793 TYPE(dbcsr_p_type), DIMENSION(:, :), &
794 INTENT(INOUT), POINTER :: moments_der
795 INTEGER, INTENT(IN) :: nmoments_der, nmoments
796 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
797 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
798 OPTIONAL, POINTER :: moments
799
800 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moments_der_matrix'
801
802 INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
803 jatom, jkind, jpgf, jset, last_jatom, m_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
804 nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
805 LOGICAL :: found
806 REAL(kind=dp) :: dab
807 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
808 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
809 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
810 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
811 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
812 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
813 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
814 TYPE(cell_type), POINTER :: cell
815 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
816 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
817 TYPE(neighbor_list_iterator_p_type), &
818 DIMENSION(:), POINTER :: nl_iterator
819 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
820 POINTER :: sab_orb
821 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
822 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
823 TYPE(qs_kind_type), POINTER :: qs_kind
824
825 nmom_build = max(nmoments, nmoments_der) ! build moments up to order nmom_buiod
826 IF (nmom_build < 1) RETURN
827
828 CALL timeset(routinen, handle)
829
830 nders = 1 ! only first order derivatives
831 dimders = ncoset(nders) - 1
832
833 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
834 CALL get_qs_env(qs_env=qs_env, &
835 qs_kind_set=qs_kind_set, &
836 particle_set=particle_set, &
837 cell=cell, &
838 sab_orb=sab_orb)
839
840 nkind = SIZE(qs_kind_set)
841
842 ! Work storage
843 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
844 maxco=maxco, maxsgf=maxsgf)
845
846 IF (nmoments > 0) THEN
847 cpassert(PRESENT(moments))
848 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
849 cpassert(SIZE(moments) == nm)
850 ! storage for integrals
851 ALLOCATE (mab(maxco, maxco, nm))
852 ! blocks
853 mab(:, :, :) = 0.0_dp
854 ALLOCATE (mom_block(nm))
855 DO i = 1, nm
856 NULLIFY (mom_block(i)%block)
857 END DO
858 END IF
859
860 IF (nmoments_der > 0) THEN
861 m_dim = ncoset(nmoments_der) - 1
862 cpassert(SIZE(moments_der, dim=1) == m_dim)
863 cpassert(SIZE(moments_der, dim=2) == dimders)
864 ! storage for integrals
865 ALLOCATE (difmab(maxco, maxco, m_dim, dimders))
866 difmab(:, :, :, :) = 0.0_dp
867 ! blocks
868 ALLOCATE (mom_block_der(m_dim, dimders))
869 DO i = 1, m_dim
870 DO ider = 1, dimders
871 NULLIFY (mom_block_der(i, ider)%block)
872 END DO
873 END DO
874 END IF
875
876 ALLOCATE (work(maxco, maxsgf))
877 work(:, :) = 0.0_dp
878
879 NULLIFY (basis_set_a, basis_set_b, basis_set_list)
880 NULLIFY (qs_kind)
881 ALLOCATE (basis_set_list(nkind))
882 DO ikind = 1, nkind
883 qs_kind => qs_kind_set(ikind)
884 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
885 IF (ASSOCIATED(basis_set_a)) THEN
886 basis_set_list(ikind)%gto_basis_set => basis_set_a
887 ELSE
888 NULLIFY (basis_set_list(ikind)%gto_basis_set)
889 END IF
890 END DO
891
892 ! Calculate derivatives looping over neighbour list
893 NULLIFY (nl_iterator)
894 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
895 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
896 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
897 iatom=iatom, jatom=jatom, r=rab)
898 basis_set_a => basis_set_list(ikind)%gto_basis_set
899 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
900 basis_set_b => basis_set_list(jkind)%gto_basis_set
901 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
902 associate( &
903 ! basis ikind
904 first_sgfa => basis_set_a%first_sgf, &
905 la_max => basis_set_a%lmax, &
906 la_min => basis_set_a%lmin, &
907 npgfa => basis_set_a%npgf, &
908 nsgfa => basis_set_a%nsgf_set, &
909 rpgfa => basis_set_a%pgf_radius, &
910 set_radius_a => basis_set_a%set_radius, &
911 sphi_a => basis_set_a%sphi, &
912 zeta => basis_set_a%zet, &
913 ! basis jkind, &
914 first_sgfb => basis_set_b%first_sgf, &
915 lb_max => basis_set_b%lmax, &
916 lb_min => basis_set_b%lmin, &
917 npgfb => basis_set_b%npgf, &
918 nsgfb => basis_set_b%nsgf_set, &
919 rpgfb => basis_set_b%pgf_radius, &
920 set_radius_b => basis_set_b%set_radius, &
921 sphi_b => basis_set_b%sphi, &
922 zetb => basis_set_b%zet)
923
924 nseta = basis_set_a%nset
925 nsetb = basis_set_b%nset
926
927 ! reference point
928 IF (PRESENT(ref_point)) THEN
929 rc(:) = ref_point(:)
930 ELSE
931 rc(:) = 0._dp
932 END IF
933 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
934 ! by folding around the center, such screwing can be avoided for a proper choice of center.
935 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
936 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
937 ! we dont use PBC at this point
938 rab(:) = ra(:) - rb(:)
939 rac(:) = ra(:) - rc(:)
940 rbc(:) = rb(:) - rc(:)
941 dab = norm2(rab)
942
943 ! get blocks
944 IF (inode == 1) last_jatom = 0
945
946 IF (jatom == last_jatom) THEN
947 cycle
948 END IF
949
950 last_jatom = jatom
951
952 IF (iatom <= jatom) THEN
953 irow = iatom
954 icol = jatom
955 ELSE
956 irow = jatom
957 icol = iatom
958 END IF
959
960 IF (nmoments > 0) THEN
961 DO i = 1, nm
962 NULLIFY (mom_block(i)%block)
963 ! get block from pre calculated overlap matrix
964 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
965 row=irow, col=icol, block=mom_block(i)%block, found=found)
966 cpassert(found .AND. ASSOCIATED(mom_block(i)%block))
967 mom_block(i)%block = 0._dp
968 END DO
969 END IF
970 IF (nmoments_der > 0) THEN
971 DO i = 1, m_dim
972 DO ider = 1, dimders
973 NULLIFY (mom_block_der(i, ider)%block)
974 CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
975 row=irow, col=icol, &
976 block=mom_block_der(i, ider)%block, &
977 found=found)
978 cpassert(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
979 mom_block_der(i, ider)%block = 0._dp
980 END DO
981 END DO
982 END IF
983
984 DO iset = 1, nseta
985
986 ncoa = npgfa(iset)*ncoset(la_max(iset))
987 sgfa = first_sgfa(1, iset)
988
989 DO jset = 1, nsetb
990
991 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
992
993 ncob = npgfb(jset)*ncoset(lb_max(jset))
994 sgfb = first_sgfb(1, jset)
995
996 NULLIFY (mab_tmp)
997 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
998 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
999
1000 ! Calculate the primitive integrals (need l+1 for derivatives)
1001 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1002 rpgfa(:, iset), la_min(iset), &
1003 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
1004 rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
1005
1006 IF (nmoments_der > 0) THEN
1007 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1008 rpgfa(:, iset), la_min(iset), &
1009 lb_max(jset), npgfb(jset), zetb(:, jset), &
1010 rpgfb(:, jset), lb_min(jset), &
1011 nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp)
1012 END IF
1013
1014 IF (nmoments > 0) THEN
1015 ! copy subset of mab_tmp (l+1) to mab (l)
1016 mab = 0.0_dp
1017 DO ii = 1, nm
1018 na = 0
1019 nda = 0
1020 DO ipgf = 1, npgfa(iset)
1021 nb = 0
1022 ndb = 0
1023 DO jpgf = 1, npgfb(jset)
1024 DO j = 1, ncoset(lb_max(jset))
1025 DO i = 1, ncoset(la_max(iset))
1026 mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
1027 END DO ! i
1028 END DO ! j
1029 nb = nb + ncoset(lb_max(jset))
1030 ndb = ndb + ncoset(lb_max(jset) + 1)
1031 END DO ! jpgf
1032 na = na + ncoset(la_max(iset))
1033 nda = nda + ncoset(la_max(iset) + 1)
1034 END DO ! ipgf
1035 END DO
1036 ! Contraction step
1037 DO i = 1, nm
1038
1039 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1040 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1041 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1042 0.0_dp, work(1, 1), SIZE(work, 1))
1043
1044 IF (iatom <= jatom) THEN
1045 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1046 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1047 work(1, 1), SIZE(work, 1), &
1048 1.0_dp, mom_block(i)%block(sgfa, sgfb), &
1049 SIZE(mom_block(i)%block, 1))
1050 ELSE
1051 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1052 1.0_dp, work(1, 1), SIZE(work, 1), &
1053 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1054 1.0_dp, mom_block(i)%block(sgfb, sgfa), &
1055 SIZE(mom_block(i)%block, 1))
1056 END IF
1057 END DO
1058 END IF
1059
1060 IF (nmoments_der > 0) THEN
1061 DO i = 1, m_dim
1062 DO ider = 1, dimders
1063 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1064 1.0_dp, difmab(1, 1, i, ider), SIZE(difmab, 1), &
1065 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1066 0._dp, work(1, 1), SIZE(work, 1))
1067
1068 IF (iatom <= jatom) THEN
1069 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1070 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1071 work(1, 1), SIZE(work, 1), &
1072 1._dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
1073 SIZE(mom_block_der(i, ider)%block, 1))
1074 ELSE
1075 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1076 -1.0_dp, work(1, 1), SIZE(work, 1), &
1077 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1078 1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
1079 SIZE(mom_block_der(i, ider)%block, 1))
1080 END IF
1081 END DO
1082 END DO
1083 END IF
1084 DEALLOCATE (mab_tmp)
1085 END DO
1086 END DO
1087 END associate
1088 END DO
1089 CALL neighbor_list_iterator_release(nl_iterator)
1090
1091 ! deallocations
1092 DEALLOCATE (basis_set_list)
1093 DEALLOCATE (work)
1094 IF (nmoments > 0) THEN
1095 DEALLOCATE (mab)
1096 DO i = 1, nm
1097 NULLIFY (mom_block(i)%block)
1098 END DO
1099 DEALLOCATE (mom_block)
1100 END IF
1101 IF (nmoments_der > 0) THEN
1102 DEALLOCATE (difmab)
1103 DO i = 1, m_dim
1104 DO ider = 1, dimders
1105 NULLIFY (mom_block_der(i, ider)%block)
1106 END DO
1107 END DO
1108 DEALLOCATE (mom_block_der)
1109 END IF
1110
1111 CALL timestop(handle)
1112
1113 END SUBROUTINE build_local_moments_der_matrix
1114
1115! **************************************************************************************************
1116!> \brief ...
1117!> \param qs_env ...
1118!> \param magmom ...
1119!> \param nmoments ...
1120!> \param ref_point ...
1121!> \param ref_points ...
1122!> \param basis_type ...
1123! **************************************************************************************************
1124 SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
1125
1126 TYPE(qs_environment_type), POINTER :: qs_env
1127 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
1128 INTEGER, INTENT(IN) :: nmoments
1129 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
1130 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
1131 OPTIONAL :: ref_points
1132 CHARACTER(len=*), OPTIONAL :: basis_type
1133
1134 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_magmom_matrix'
1135
1136 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
1137 maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
1138 LOGICAL :: found
1139 REAL(kind=dp) :: dab
1140 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1141 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
1142 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
1143 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
1144 TYPE(cell_type), POINTER :: cell
1145 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1146 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1147 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1148 TYPE(neighbor_list_iterator_p_type), &
1149 DIMENSION(:), POINTER :: nl_iterator
1150 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1151 POINTER :: sab_orb
1152 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1153 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1154 TYPE(qs_kind_type), POINTER :: qs_kind
1155
1156 IF (nmoments < 1) RETURN
1157
1158 CALL timeset(routinen, handle)
1159
1160 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
1161 NULLIFY (matrix_s)
1162
1163 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1164
1165 ! magnetic dipoles/angular moments only
1166 nm = 3
1167
1168 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1169 CALL get_qs_env(qs_env=qs_env, &
1170 qs_kind_set=qs_kind_set, &
1171 particle_set=particle_set, cell=cell, &
1172 sab_orb=sab_orb)
1173
1174 nkind = SIZE(qs_kind_set)
1175 natom = SIZE(particle_set)
1176
1177 ! Allocate work storage
1178 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1179 maxco=maxco, maxsgf=maxsgf)
1180
1181 ALLOCATE (mab(maxco, maxco, nm))
1182 mab(:, :, :) = 0.0_dp
1183
1184 ALLOCATE (work(maxco, maxsgf))
1185 work(:, :) = 0.0_dp
1186
1187 ALLOCATE (mint(nm))
1188 DO i = 1, nm
1189 NULLIFY (mint(i)%block)
1190 END DO
1191
1192 ALLOCATE (basis_set_list(nkind))
1193 DO ikind = 1, nkind
1194 qs_kind => qs_kind_set(ikind)
1195 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1196 IF (ASSOCIATED(basis_set_a)) THEN
1197 basis_set_list(ikind)%gto_basis_set => basis_set_a
1198 ELSE
1199 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1200 END IF
1201 END DO
1202 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1203 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1204 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1205 iatom=iatom, jatom=jatom, r=rab)
1206 basis_set_a => basis_set_list(ikind)%gto_basis_set
1207 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1208 basis_set_b => basis_set_list(jkind)%gto_basis_set
1209 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1210 associate( &
1211 ! basis ikind
1212 first_sgfa => basis_set_a%first_sgf, &
1213 la_max => basis_set_a%lmax, &
1214 la_min => basis_set_a%lmin, &
1215 npgfa => basis_set_a%npgf, &
1216 nsgfa => basis_set_a%nsgf_set, &
1217 rpgfa => basis_set_a%pgf_radius, &
1218 set_radius_a => basis_set_a%set_radius, &
1219 sphi_a => basis_set_a%sphi, &
1220 zeta => basis_set_a%zet, &
1221 ! basis jkind, &
1222 first_sgfb => basis_set_b%first_sgf, &
1223 lb_max => basis_set_b%lmax, &
1224 lb_min => basis_set_b%lmin, &
1225 npgfb => basis_set_b%npgf, &
1226 nsgfb => basis_set_b%nsgf_set, &
1227 rpgfb => basis_set_b%pgf_radius, &
1228 set_radius_b => basis_set_b%set_radius, &
1229 sphi_b => basis_set_b%sphi, &
1230 zetb => basis_set_b%zet)
1231
1232 nseta = basis_set_a%nset
1233 nsetb = basis_set_b%nset
1234
1235 IF (iatom <= jatom) THEN
1236 irow = iatom
1237 icol = jatom
1238 ELSE
1239 irow = jatom
1240 icol = iatom
1241 END IF
1242
1243 DO i = 1, nm
1244 NULLIFY (mint(i)%block)
1245 CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
1246 row=irow, col=icol, block=mint(i)%block, found=found)
1247 mint(i)%block = 0._dp
1248 cpassert(ASSOCIATED(mint(i)%block))
1249 END DO
1250
1251 ! fold atomic position back into unit cell
1252 IF (PRESENT(ref_points)) THEN
1253 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
1254 ELSE IF (PRESENT(ref_point)) THEN
1255 rc(:) = ref_point(:)
1256 ELSE
1257 rc(:) = 0._dp
1258 END IF
1259 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
1260 ! by folding around the center, such screwing can be avoided for a proper choice of center.
1261 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
1262 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
1263 ! we dont use PBC at this point
1264 rab(:) = ra(:) - rb(:)
1265 rac(:) = ra(:) - rc(:)
1266 rbc(:) = rb(:) - rc(:)
1267 dab = norm2(rab)
1268
1269 DO iset = 1, nseta
1270
1271 ncoa = npgfa(iset)*ncoset(la_max(iset))
1272 sgfa = first_sgfa(1, iset)
1273
1274 DO jset = 1, nsetb
1275
1276 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1277
1278 ncob = npgfb(jset)*ncoset(lb_max(jset))
1279 sgfb = first_sgfb(1, jset)
1280
1281 ! Calculate the primitive integrals
1282 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
1283 rpgfa(:, iset), la_min(iset), &
1284 lb_max(jset), npgfb(jset), zetb(:, jset), &
1285 rpgfb(:, jset), rac, rbc, mab)
1286
1287 ! Contraction step
1288 DO i = 1, nm
1289 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1290 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1291 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1292 0.0_dp, work(1, 1), SIZE(work, 1))
1293
1294 IF (iatom <= jatom) THEN
1295 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1296 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1297 work(1, 1), SIZE(work, 1), &
1298 1.0_dp, mint(i)%block(sgfa, sgfb), &
1299 SIZE(mint(i)%block, 1))
1300 ELSE
1301 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1302 -1.0_dp, work(1, 1), SIZE(work, 1), &
1303 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1304 1.0_dp, mint(i)%block(sgfb, sgfa), &
1305 SIZE(mint(i)%block, 1))
1306 END IF
1307
1308 END DO
1309
1310 END DO
1311 END DO
1312 END associate
1313 END DO
1314 CALL neighbor_list_iterator_release(nl_iterator)
1315
1316 ! Release work storage
1317 DEALLOCATE (mab, basis_set_list)
1318 DEALLOCATE (work)
1319 DO i = 1, nm
1320 NULLIFY (mint(i)%block)
1321 END DO
1322 DEALLOCATE (mint)
1323
1324 CALL timestop(handle)
1325
1326 END SUBROUTINE build_local_magmom_matrix
1327
1328! **************************************************************************************************
1329!> \brief ...
1330!> \param qs_env ...
1331!> \param cosmat ...
1332!> \param sinmat ...
1333!> \param kvec ...
1334!> \param sab_orb_external ...
1335!> \param basis_type ...
1336! **************************************************************************************************
1337 SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1338
1339 TYPE(qs_environment_type), POINTER :: qs_env
1340 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1341 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1342 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1343 OPTIONAL, POINTER :: sab_orb_external
1344 CHARACTER(len=*), OPTIONAL :: basis_type
1345
1346 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_moment_matrix'
1347
1348 INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, ldsa, &
1349 ldsb, ldwork, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
1350 LOGICAL :: found
1351 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1352 REAL(kind=dp) :: dab
1353 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1354 TYPE(cell_type), POINTER :: cell
1355 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1356 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1357 TYPE(neighbor_list_iterator_p_type), &
1358 DIMENSION(:), POINTER :: nl_iterator
1359 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1360 POINTER :: sab_orb
1361 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1362 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1363 TYPE(qs_kind_type), POINTER :: qs_kind
1364
1365 CALL timeset(routinen, handle)
1366
1367 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1368 CALL get_qs_env(qs_env=qs_env, &
1369 qs_kind_set=qs_kind_set, &
1370 particle_set=particle_set, cell=cell, &
1371 sab_orb=sab_orb)
1372
1373 IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
1374
1375 CALL dbcsr_set(sinmat, 0.0_dp)
1376 CALL dbcsr_set(cosmat, 0.0_dp)
1377
1378 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
1379 ldab = ldwork
1380 ALLOCATE (cosab(ldab, ldab))
1381 ALLOCATE (sinab(ldab, ldab))
1382 ALLOCATE (work(ldwork, ldwork))
1383
1384 nkind = SIZE(qs_kind_set)
1385 natom = SIZE(particle_set)
1386
1387 ALLOCATE (basis_set_list(nkind))
1388 DO ikind = 1, nkind
1389 qs_kind => qs_kind_set(ikind)
1390 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1391 IF (ASSOCIATED(basis_set_a)) THEN
1392 basis_set_list(ikind)%gto_basis_set => basis_set_a
1393 ELSE
1394 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1395 END IF
1396 END DO
1397 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1398 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1399 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1400 iatom=iatom, jatom=jatom, r=rab)
1401 basis_set_a => basis_set_list(ikind)%gto_basis_set
1402 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1403 basis_set_b => basis_set_list(jkind)%gto_basis_set
1404 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1405 associate( &
1406 ! basis ikind
1407 first_sgfa => basis_set_a%first_sgf, &
1408 la_max => basis_set_a%lmax, &
1409 la_min => basis_set_a%lmin, &
1410 npgfa => basis_set_a%npgf, &
1411 nsgfa => basis_set_a%nsgf_set, &
1412 rpgfa => basis_set_a%pgf_radius, &
1413 set_radius_a => basis_set_a%set_radius, &
1414 sphi_a => basis_set_a%sphi, &
1415 zeta => basis_set_a%zet, &
1416 ! basis jkind, &
1417 first_sgfb => basis_set_b%first_sgf, &
1418 lb_max => basis_set_b%lmax, &
1419 lb_min => basis_set_b%lmin, &
1420 npgfb => basis_set_b%npgf, &
1421 nsgfb => basis_set_b%nsgf_set, &
1422 rpgfb => basis_set_b%pgf_radius, &
1423 set_radius_b => basis_set_b%set_radius, &
1424 sphi_b => basis_set_b%sphi, &
1425 zetb => basis_set_b%zet)
1426
1427 nseta = basis_set_a%nset
1428 nsetb = basis_set_b%nset
1429
1430 ldsa = SIZE(sphi_a, 1)
1431 ldsb = SIZE(sphi_b, 1)
1432
1433 IF (iatom <= jatom) THEN
1434 irow = iatom
1435 icol = jatom
1436 ELSE
1437 irow = jatom
1438 icol = iatom
1439 END IF
1440
1441 NULLIFY (cblock)
1442 CALL dbcsr_get_block_p(matrix=cosmat, &
1443 row=irow, col=icol, block=cblock, found=found)
1444 NULLIFY (sblock)
1445 CALL dbcsr_get_block_p(matrix=sinmat, &
1446 row=irow, col=icol, block=sblock, found=found)
1447 IF (ASSOCIATED(cblock) .AND. .NOT. ASSOCIATED(sblock) .OR. &
1448 .NOT. ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1449 cpabort("")
1450 END IF
1451
1452 IF (ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1453
1454 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1455 rb(:) = ra + rab
1456 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1457
1458 DO iset = 1, nseta
1459
1460 ncoa = npgfa(iset)*ncoset(la_max(iset))
1461 sgfa = first_sgfa(1, iset)
1462
1463 DO jset = 1, nsetb
1464
1465 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1466
1467 ncob = npgfb(jset)*ncoset(lb_max(jset))
1468 sgfb = first_sgfb(1, jset)
1469
1470 ! Calculate the primitive integrals
1471 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1472 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1473 ra, rb, kvec, cosab, sinab)
1474 CALL contract_cossin(cblock, sblock, &
1475 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1476 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1477 cosab, sinab, ldab, work, ldwork)
1478
1479 END DO
1480 END DO
1481
1482 END IF
1483 END associate
1484 END DO
1485 CALL neighbor_list_iterator_release(nl_iterator)
1486
1487 DEALLOCATE (cosab)
1488 DEALLOCATE (sinab)
1489 DEALLOCATE (work)
1490 DEALLOCATE (basis_set_list)
1491
1492 CALL timestop(handle)
1493
1494 END SUBROUTINE build_berry_moment_matrix
1495
1496! **************************************************************************************************
1497!> \brief ...
1498!> \param qs_env ...
1499!> \param cosmat ...
1500!> \param sinmat ...
1501!> \param kvec ...
1502!> \param basis_type ...
1503! **************************************************************************************************
1504 SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
1505
1506 TYPE(qs_environment_type), POINTER :: qs_env
1507 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1508 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1509 CHARACTER(len=*), OPTIONAL :: basis_type
1510
1511 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_kpoint_matrix'
1512
1513 INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1514 ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1515 INTEGER, DIMENSION(3) :: icell
1516 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1517 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1518 LOGICAL :: found, use_cell_mapping
1519 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1520 REAL(kind=dp) :: dab
1521 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1522 TYPE(cell_type), POINTER :: cell
1523 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1524 TYPE(dft_control_type), POINTER :: dft_control
1525 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1526 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1527 TYPE(kpoint_type), POINTER :: kpoints
1528 TYPE(neighbor_list_iterator_p_type), &
1529 DIMENSION(:), POINTER :: nl_iterator
1530 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1531 POINTER :: sab_orb
1532 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1533 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1534 TYPE(qs_kind_type), POINTER :: qs_kind
1535 TYPE(qs_ks_env_type), POINTER :: ks_env
1536
1537 CALL timeset(routinen, handle)
1538
1539 CALL get_qs_env(qs_env, &
1540 ks_env=ks_env, &
1541 dft_control=dft_control)
1542 nimg = dft_control%nimages
1543 IF (nimg > 1) THEN
1544 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1545 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1546 use_cell_mapping = .true.
1547 ELSE
1548 use_cell_mapping = .false.
1549 END IF
1550
1551 CALL get_qs_env(qs_env=qs_env, &
1552 qs_kind_set=qs_kind_set, &
1553 particle_set=particle_set, cell=cell, &
1554 sab_orb=sab_orb)
1555
1556 nkind = SIZE(qs_kind_set)
1557 natom = SIZE(particle_set)
1558 ALLOCATE (basis_set_list(nkind))
1559 DO ikind = 1, nkind
1560 qs_kind => qs_kind_set(ikind)
1561 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1562 IF (ASSOCIATED(basis_set)) THEN
1563 basis_set_list(ikind)%gto_basis_set => basis_set
1564 ELSE
1565 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1566 END IF
1567 END DO
1568
1569 ALLOCATE (row_blk_sizes(natom))
1570 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1571 basis=basis_set_list)
1572 CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1573 ! (re)allocate matrix sets
1574 CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1575 CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1576 DO i = 1, nimg
1577 ! sin
1578 ALLOCATE (sinmat(1, i)%matrix)
1579 CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1580 name="SINMAT", &
1581 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1582 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1583 CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1584 CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1585 ! cos
1586 ALLOCATE (cosmat(1, i)%matrix)
1587 CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1588 name="COSMAT", &
1589 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1590 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes)
1591 CALL cp_dbcsr_alloc_block_from_nbl(cosmat(1, i)%matrix, sab_orb)
1592 CALL dbcsr_set(cosmat(1, i)%matrix, 0.0_dp)
1593 END DO
1594
1595 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork)
1596 ldab = ldwork
1597 ALLOCATE (cosab(ldab, ldab))
1598 ALLOCATE (sinab(ldab, ldab))
1599 ALLOCATE (work(ldwork, ldwork))
1600
1601 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1602 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1603 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1604 iatom=iatom, jatom=jatom, r=rab, cell=icell)
1605 basis_set_a => basis_set_list(ikind)%gto_basis_set
1606 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1607 basis_set_b => basis_set_list(jkind)%gto_basis_set
1608 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1609 associate( &
1610 ! basis ikind
1611 first_sgfa => basis_set_a%first_sgf, &
1612 la_max => basis_set_a%lmax, &
1613 la_min => basis_set_a%lmin, &
1614 npgfa => basis_set_a%npgf, &
1615 nsgfa => basis_set_a%nsgf_set, &
1616 rpgfa => basis_set_a%pgf_radius, &
1617 set_radius_a => basis_set_a%set_radius, &
1618 sphi_a => basis_set_a%sphi, &
1619 zeta => basis_set_a%zet, &
1620 ! basis jkind, &
1621 first_sgfb => basis_set_b%first_sgf, &
1622 lb_max => basis_set_b%lmax, &
1623 lb_min => basis_set_b%lmin, &
1624 npgfb => basis_set_b%npgf, &
1625 nsgfb => basis_set_b%nsgf_set, &
1626 rpgfb => basis_set_b%pgf_radius, &
1627 set_radius_b => basis_set_b%set_radius, &
1628 sphi_b => basis_set_b%sphi, &
1629 zetb => basis_set_b%zet)
1630
1631 nseta = basis_set_a%nset
1632 nsetb = basis_set_b%nset
1633
1634 ldsa = SIZE(sphi_a, 1)
1635 ldsb = SIZE(sphi_b, 1)
1636
1637 IF (iatom <= jatom) THEN
1638 irow = iatom
1639 icol = jatom
1640 ELSE
1641 irow = jatom
1642 icol = iatom
1643 END IF
1644
1645 IF (use_cell_mapping) THEN
1646 ic = cell_to_index(icell(1), icell(2), icell(3))
1647 cpassert(ic > 0)
1648 ELSE
1649 ic = 1
1650 END IF
1651
1652 NULLIFY (sblock)
1653 CALL dbcsr_get_block_p(matrix=sinmat(1, ic)%matrix, &
1654 row=irow, col=icol, block=sblock, found=found)
1655 cpassert(found)
1656 NULLIFY (cblock)
1657 CALL dbcsr_get_block_p(matrix=cosmat(1, ic)%matrix, &
1658 row=irow, col=icol, block=cblock, found=found)
1659 cpassert(found)
1660
1661 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1662 rb(:) = ra + rab
1663 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1664
1665 DO iset = 1, nseta
1666
1667 ncoa = npgfa(iset)*ncoset(la_max(iset))
1668 sgfa = first_sgfa(1, iset)
1669
1670 DO jset = 1, nsetb
1671
1672 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1673
1674 ncob = npgfb(jset)*ncoset(lb_max(jset))
1675 sgfb = first_sgfb(1, jset)
1676
1677 ! Calculate the primitive integrals
1678 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1679 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1680 ra, rb, kvec, cosab, sinab)
1681 CALL contract_cossin(cblock, sblock, &
1682 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1683 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1684 cosab, sinab, ldab, work, ldwork)
1685
1686 END DO
1687 END DO
1688 END associate
1689 END DO
1690 CALL neighbor_list_iterator_release(nl_iterator)
1691
1692 DEALLOCATE (cosab)
1693 DEALLOCATE (sinab)
1694 DEALLOCATE (work)
1695 DEALLOCATE (basis_set_list)
1696 DEALLOCATE (row_blk_sizes)
1697
1698 CALL timestop(handle)
1699
1700 END SUBROUTINE build_berry_kpoint_matrix
1701
1702! **************************************************************************************************
1703!> \brief ...
1704!> \param qs_env ...
1705!> \param magnetic ...
1706!> \param nmoments ...
1707!> \param reference ...
1708!> \param ref_point ...
1709!> \param unit_number ...
1710! **************************************************************************************************
1711 SUBROUTINE qs_moment_berry_phase(qs_env, magnetic, nmoments, reference, ref_point, unit_number)
1712
1713 TYPE(qs_environment_type), POINTER :: qs_env
1714 LOGICAL, INTENT(IN) :: magnetic
1715 INTEGER, INTENT(IN) :: nmoments, reference
1716 REAL(dp), DIMENSION(:), POINTER :: ref_point
1717 INTEGER, INTENT(IN) :: unit_number
1718
1719 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_berry_phase'
1720
1721 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
1722 CHARACTER(LEN=default_string_length) :: description
1723 COMPLEX(dp) :: xphase(3), zdet, zdeta, zi(3), &
1724 zij(3, 3), zijk(3, 3, 3), &
1725 zijkl(3, 3, 3, 3), zphase(3), zz
1726 INTEGER :: handle, i, ia, idim, ikind, ispin, ix, &
1727 iy, iz, j, k, l, nao, nm, nmo, nmom, &
1728 nmotot, tmp_dim
1729 LOGICAL :: floating, ghost, uniform
1730 REAL(dp) :: charge, ci(3), cij(3, 3), dd, occ, trace
1731 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom
1732 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
1733 REAL(dp), DIMENSION(3) :: kvec, qq, rcc, ria
1734 TYPE(atomic_kind_type), POINTER :: atomic_kind
1735 TYPE(cell_type), POINTER :: cell
1736 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
1737 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
1738 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: opvec
1739 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set
1740 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
1741 TYPE(cp_fm_type), POINTER :: mo_coeff
1742 TYPE(cp_result_type), POINTER :: results
1743 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ao
1744 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1745 TYPE(dft_control_type), POINTER :: dft_control
1746 TYPE(distribution_1d_type), POINTER :: local_particles
1747 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1748 TYPE(mp_para_env_type), POINTER :: para_env
1749 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1750 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1751 TYPE(qs_rho_type), POINTER :: rho
1752 TYPE(rt_prop_type), POINTER :: rtp
1753
1754 cpassert(ASSOCIATED(qs_env))
1755
1756 IF (ASSOCIATED(qs_env%ls_scf_env)) THEN
1757 IF (unit_number > 0) WRITE (unit_number, *) "Periodic moment calculation not implemented in linear scaling code"
1758 RETURN
1759 END IF
1760
1761 CALL timeset(routinen, handle)
1762
1763 ! restrict maximum moment available
1764 nmom = min(nmoments, 2)
1765
1766 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
1767 ! rmom(:,1)=electronic
1768 ! rmom(:,2)=nuclear
1769 ! rmom(:,1)=total
1770 ALLOCATE (rmom(nm + 1, 3))
1771 ALLOCATE (rlab(nm + 1))
1772 rmom = 0.0_dp
1773 rlab = ""
1774 IF (magnetic) THEN
1775 nm = 3
1776 ALLOCATE (mmom(nm))
1777 mmom = 0._dp
1778 END IF
1779
1780 NULLIFY (dft_control, rho, cell, particle_set, results, para_env, &
1781 local_particles, matrix_s, mos, rho_ao)
1782
1783 CALL get_qs_env(qs_env, &
1784 dft_control=dft_control, &
1785 rho=rho, &
1786 cell=cell, &
1787 results=results, &
1788 particle_set=particle_set, &
1789 qs_kind_set=qs_kind_set, &
1790 para_env=para_env, &
1791 local_particles=local_particles, &
1792 matrix_s=matrix_s, &
1793 mos=mos)
1794
1795 CALL qs_rho_get(rho, rho_ao=rho_ao)
1796
1797 NULLIFY (cosmat, sinmat)
1798 ALLOCATE (cosmat, sinmat)
1799 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
1800 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
1801 CALL dbcsr_set(cosmat, 0.0_dp)
1802 CALL dbcsr_set(sinmat, 0.0_dp)
1803
1804 ALLOCATE (op_fm_set(2, dft_control%nspins))
1805 ALLOCATE (opvec(dft_control%nspins))
1806 ALLOCATE (eigrmat(dft_control%nspins))
1807 nmotot = 0
1808 DO ispin = 1, dft_control%nspins
1809 NULLIFY (tmp_fm_struct, mo_coeff)
1810 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
1811 nmotot = nmotot + nmo
1812 CALL cp_fm_create(opvec(ispin), mo_coeff%matrix_struct)
1813 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
1814 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
1815 DO i = 1, SIZE(op_fm_set, 1)
1816 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
1817 END DO
1818 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
1819 CALL cp_fm_struct_release(tmp_fm_struct)
1820 END DO
1821
1822 ! occupation
1823 DO ispin = 1, dft_control%nspins
1824 CALL get_mo_set(mo_set=mos(ispin), maxocc=occ, uniform_occupation=uniform)
1825 IF (.NOT. uniform) THEN
1826 cpwarn("Berry phase moments for non uniform MOs' occupation numbers not implemented")
1827 END IF
1828 END DO
1829
1830 ! reference point
1831 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
1832 rcc = pbc(rcc, cell)
1833
1834 ! label
1835 DO l = 1, nm
1836 ix = indco(1, l + 1)
1837 iy = indco(2, l + 1)
1838 iz = indco(3, l + 1)
1839 CALL set_label(rlab(l + 1), ix, iy, iz)
1840 END DO
1841
1842 ! nuclear contribution
1843 DO ia = 1, SIZE(particle_set)
1844 atomic_kind => particle_set(ia)%atomic_kind
1845 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1846 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1847 IF (.NOT. ghost .AND. .NOT. floating) THEN
1848 rmom(1, 2) = rmom(1, 2) - charge
1849 END IF
1850 END DO
1851 ria = twopi*matmul(cell%h_inv, rcc)
1852 zphase = cmplx(cos(ria), sin(ria), dp)**rmom(1, 2)
1853
1854 zi = 0._dp
1855 zij = 0._dp
1856 zijk = 0._dp
1857 zijkl = 0._dp
1858
1859 DO l = 1, nmom
1860 SELECT CASE (l)
1861 CASE (1)
1862 ! Dipole
1863 zi(:) = cmplx(1._dp, 0._dp, dp)
1864 DO ia = 1, SIZE(particle_set)
1865 atomic_kind => particle_set(ia)%atomic_kind
1866 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
1867 CALL get_qs_kind(qs_kind_set(ikind), core_charge=charge, ghost=ghost, floating=floating)
1868 IF (.NOT. ghost .AND. .NOT. floating) THEN
1869 ria = particle_set(ia)%r
1870 ria = pbc(ria, cell)
1871 DO i = 1, 3
1872 kvec(:) = twopi*cell%h_inv(i, :)
1873 dd = sum(kvec(:)*ria(:))
1874 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1875 zi(i) = zi(i)*zdeta
1876 END DO
1877 END IF
1878 END DO
1879 zi = zi*zphase
1880 ci = aimag(log(zi))/twopi
1881 qq = aimag(log(zi))
1882 rmom(2:4, 2) = matmul(cell%hmat, ci)
1883 CASE (2)
1884 ! Quadrupole
1885 cpabort("Berry phase moments bigger than 1 not implemented")
1886 zij(:, :) = 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)
1891 ria = particle_set(ia)%r
1892 ria = pbc(ria, cell)
1893 DO i = 1, 3
1894 DO j = i, 3
1895 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1896 dd = sum(kvec(:)*ria(:))
1897 zdeta = cmplx(cos(dd), sin(dd), kind=dp)**charge
1898 zij(i, j) = zij(i, j)*zdeta
1899 zij(j, i) = zij(i, j)
1900 END DO
1901 END DO
1902 END DO
1903 DO i = 1, 3
1904 DO j = 1, 3
1905 zij(i, j) = zij(i, j)*zphase(i)*zphase(j)
1906 zz = zij(i, j)/zi(i)/zi(j)
1907 cij(i, j) = aimag(log(zz))/twopi
1908 END DO
1909 END DO
1910 cij = 0.5_dp*cij/twopi/twopi
1911 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
1912 DO k = 4, 9
1913 ix = indco(1, k + 1)
1914 iy = indco(2, k + 1)
1915 iz = indco(3, k + 1)
1916 IF (ix == 0) THEN
1917 rmom(k + 1, 2) = cij(iy, iz)
1918 ELSE IF (iy == 0) THEN
1919 rmom(k + 1, 2) = cij(ix, iz)
1920 ELSE IF (iz == 0) THEN
1921 rmom(k + 1, 2) = cij(ix, iy)
1922 END IF
1923 END DO
1924 CASE (3)
1925 ! Octapole
1926 cpabort("Berry phase moments bigger than 2 not implemented")
1927 CASE (4)
1928 ! Hexadecapole
1929 cpabort("Berry phase moments bigger than 3 not implemented")
1930 CASE DEFAULT
1931 cpabort("Berry phase moments bigger than 4 not implemented")
1932 END SELECT
1933 END DO
1934
1935 ! electronic contribution
1936
1937 ria = twopi*real(nmotot, dp)*occ*matmul(cell%h_inv, rcc)
1938 xphase = cmplx(cos(ria), sin(ria), dp)
1939
1940 ! charge
1941 trace = 0.0_dp
1942 DO ispin = 1, dft_control%nspins
1943 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
1944 rmom(1, 1) = rmom(1, 1) + trace
1945 END DO
1946
1947 zi = 0._dp
1948 zij = 0._dp
1949 zijk = 0._dp
1950 zijkl = 0._dp
1951
1952 DO l = 1, nmom
1953 SELECT CASE (l)
1954 CASE (1)
1955 ! Dipole
1956 DO i = 1, 3
1957 kvec(:) = twopi*cell%h_inv(i, :)
1958 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1959 IF (qs_env%run_rtp) THEN
1960 CALL get_qs_env(qs_env, rtp=rtp)
1961 CALL get_rtp(rtp, mos_new=mos_new)
1962 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1963 ELSE
1964 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1965 END IF
1966 zdet = cmplx(1._dp, 0._dp, dp)
1967 DO ispin = 1, dft_control%nspins
1968 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
1969 DO idim = 1, tmp_dim
1970 eigrmat(ispin)%local_data(:, idim) = &
1971 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
1972 -op_fm_set(2, ispin)%local_data(:, idim), dp)
1973 END DO
1974 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
1975 CALL cp_cfm_det(eigrmat(ispin), zdeta)
1976 zdet = zdet*zdeta
1977 IF (dft_control%nspins == 1) THEN
1978 zdet = zdet*zdeta
1979 END IF
1980 END DO
1981 zi(i) = zdet
1982 END DO
1983 zi = zi*xphase
1984 CASE (2)
1985 ! Quadrupole
1986 cpabort("Berry phase moments bigger than 1 not implemented")
1987 DO i = 1, 3
1988 DO j = i, 3
1989 kvec(:) = twopi*(cell%h_inv(i, :) + cell%h_inv(j, :))
1990 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
1991 IF (qs_env%run_rtp) THEN
1992 CALL get_qs_env(qs_env, rtp=rtp)
1993 CALL get_rtp(rtp, mos_new=mos_new)
1994 CALL op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
1995 ELSE
1996 CALL op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
1997 END IF
1998 zdet = cmplx(1._dp, 0._dp, dp)
1999 DO ispin = 1, dft_control%nspins
2000 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
2001 DO idim = 1, tmp_dim
2002 eigrmat(ispin)%local_data(:, idim) = &
2003 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
2004 -op_fm_set(2, ispin)%local_data(:, idim), dp)
2005 END DO
2006 ! CALL cp_cfm_lu_decompose(eigrmat(ispin), zdeta)
2007 CALL cp_cfm_det(eigrmat(ispin), zdeta)
2008 zdet = zdet*zdeta
2009 IF (dft_control%nspins == 1) THEN
2010 zdet = zdet*zdeta
2011 END IF
2012 END DO
2013 zij(i, j) = zdet*xphase(i)*xphase(j)
2014 zij(j, i) = zdet*xphase(i)*xphase(j)
2015 END DO
2016 END DO
2017 CASE (3)
2018 ! Octapole
2019 cpabort("Berry phase moments bigger than 2 not implemented")
2020 CASE (4)
2021 ! Hexadecapole
2022 cpabort("Berry phase moments bigger than 3 not implemented")
2023 CASE DEFAULT
2024 cpabort("Berry phase moments bigger than 4 not implemented")
2025 END SELECT
2026 END DO
2027 DO l = 1, nmom
2028 SELECT CASE (l)
2029 CASE (1)
2030 ! Dipole (apply periodic (2 Pi) boundary conditions)
2031 ci = aimag(log(zi))
2032 DO i = 1, 3
2033 IF (qq(i) + ci(i) > pi) ci(i) = ci(i) - twopi
2034 IF (qq(i) + ci(i) < -pi) ci(i) = ci(i) + twopi
2035 END DO
2036 rmom(2:4, 1) = matmul(cell%hmat, ci)/twopi
2037 CASE (2)
2038 ! Quadrupole
2039 cpabort("Berry phase moments bigger than 1 not implemented")
2040 DO i = 1, 3
2041 DO j = 1, 3
2042 zz = zij(i, j)/zi(i)/zi(j)
2043 cij(i, j) = aimag(log(zz))/twopi
2044 END DO
2045 END DO
2046 cij = 0.5_dp*cij/twopi/twopi
2047 cij = matmul(matmul(cell%hmat, cij), transpose(cell%hmat))
2048 DO k = 4, 9
2049 ix = indco(1, k + 1)
2050 iy = indco(2, k + 1)
2051 iz = indco(3, k + 1)
2052 IF (ix == 0) THEN
2053 rmom(k + 1, 1) = cij(iy, iz)
2054 ELSE IF (iy == 0) THEN
2055 rmom(k + 1, 1) = cij(ix, iz)
2056 ELSE IF (iz == 0) THEN
2057 rmom(k + 1, 1) = cij(ix, iy)
2058 END IF
2059 END DO
2060 CASE (3)
2061 ! Octapole
2062 cpabort("Berry phase moments bigger than 2 not implemented")
2063 CASE (4)
2064 ! Hexadecapole
2065 cpabort("Berry phase moments bigger than 3 not implemented")
2066 CASE DEFAULT
2067 cpabort("Berry phase moments bigger than 4 not implemented")
2068 END SELECT
2069 END DO
2070
2071 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2072 description = "[DIPOLE]"
2073 CALL cp_results_erase(results=results, description=description)
2074 CALL put_results(results=results, description=description, &
2075 values=rmom(2:4, 3))
2076 IF (magnetic) THEN
2077 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true., mmom=mmom)
2078 ELSE
2079 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.true.)
2080 END IF
2081
2082 DEALLOCATE (rmom)
2083 DEALLOCATE (rlab)
2084 IF (magnetic) THEN
2085 DEALLOCATE (mmom)
2086 END IF
2087
2088 CALL dbcsr_deallocate_matrix(cosmat)
2089 CALL dbcsr_deallocate_matrix(sinmat)
2090
2091 CALL cp_fm_release(opvec)
2092 CALL cp_fm_release(op_fm_set)
2093 DO ispin = 1, dft_control%nspins
2094 CALL cp_cfm_release(eigrmat(ispin))
2095 END DO
2096 DEALLOCATE (eigrmat)
2097
2098 CALL timestop(handle)
2099
2100 END SUBROUTINE qs_moment_berry_phase
2101
2102! **************************************************************************************************
2103!> \brief ...
2104!> \param cosmat ...
2105!> \param sinmat ...
2106!> \param mos ...
2107!> \param op_fm_set ...
2108!> \param opvec ...
2109! **************************************************************************************************
2110 SUBROUTINE op_orbbas(cosmat, sinmat, mos, op_fm_set, opvec)
2111
2112 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2113 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2114 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2115 TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT) :: opvec
2116
2117 INTEGER :: i, nao, nmo
2118 TYPE(cp_fm_type), POINTER :: mo_coeff
2119
2120 DO i = 1, SIZE(op_fm_set, 2) ! spin
2121 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2122 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(i), ncol=nmo)
2123 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2124 op_fm_set(1, i))
2125 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(i), ncol=nmo)
2126 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(i), 0.0_dp, &
2127 op_fm_set(2, i))
2128 END DO
2129
2130 END SUBROUTINE op_orbbas
2131
2132! **************************************************************************************************
2133!> \brief ...
2134!> \param cosmat ...
2135!> \param sinmat ...
2136!> \param mos ...
2137!> \param op_fm_set ...
2138!> \param mos_new ...
2139! **************************************************************************************************
2140 SUBROUTINE op_orbbas_rtp(cosmat, sinmat, mos, op_fm_set, mos_new)
2141
2142 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
2143 TYPE(mo_set_type), DIMENSION(:), INTENT(IN) :: mos
2144 TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: op_fm_set
2145 TYPE(cp_fm_type), DIMENSION(:), POINTER :: mos_new
2146
2147 INTEGER :: i, icol, lcol, nao, newdim, nmo
2148 LOGICAL :: double_col, double_row
2149 TYPE(cp_fm_struct_type), POINTER :: newstruct, newstruct1
2150 TYPE(cp_fm_type) :: work, work1, work2
2151 TYPE(cp_fm_type), POINTER :: mo_coeff
2152
2153 DO i = 1, SIZE(op_fm_set, 2) ! spin
2154 CALL get_mo_set(mo_set=mos(i), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
2155 CALL cp_fm_get_info(mos_new(2*i), ncol_local=lcol, ncol_global=nmo)
2156 double_col = .true.
2157 double_row = .false.
2158 CALL cp_fm_struct_double(newstruct, &
2159 mos_new(2*i)%matrix_struct, &
2160 mos_new(2*i)%matrix_struct%context, &
2161 double_col, &
2162 double_row)
2163
2164 CALL cp_fm_create(work, matrix_struct=newstruct)
2165 CALL cp_fm_create(work1, matrix_struct=newstruct)
2166 CALL cp_fm_create(work2, matrix_struct=newstruct)
2167 CALL cp_fm_get_info(work, ncol_global=newdim)
2168
2169 CALL cp_fm_set_all(work, 0.0_dp, 0.0_dp)
2170 DO icol = 1, lcol
2171 work%local_data(:, icol) = mos_new(2*i - 1)%local_data(:, icol)
2172 work%local_data(:, icol + lcol) = mos_new(2*i)%local_data(:, icol)
2173 END DO
2174
2175 CALL cp_dbcsr_sm_fm_multiply(cosmat, work, work1, ncol=newdim)
2176 CALL cp_dbcsr_sm_fm_multiply(sinmat, work, work2, ncol=newdim)
2177
2178 DO icol = 1, lcol
2179 work%local_data(:, icol) = work1%local_data(:, icol) - work2%local_data(:, icol + lcol)
2180 work%local_data(:, icol + lcol) = work1%local_data(:, icol + lcol) + work2%local_data(:, icol)
2181 END DO
2182
2183 CALL cp_fm_release(work1)
2184 CALL cp_fm_release(work2)
2185
2186 CALL cp_fm_struct_double(newstruct1, &
2187 op_fm_set(1, i)%matrix_struct, &
2188 op_fm_set(1, i)%matrix_struct%context, &
2189 double_col, &
2190 double_row)
2191
2192 CALL cp_fm_create(work1, matrix_struct=newstruct1)
2193
2194 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i - 1), &
2195 work, 0.0_dp, work1)
2196
2197 DO icol = 1, lcol
2198 op_fm_set(1, i)%local_data(:, icol) = work1%local_data(:, icol)
2199 op_fm_set(2, i)%local_data(:, icol) = work1%local_data(:, icol + lcol)
2200 END DO
2201
2202 CALL parallel_gemm("T", "N", nmo, newdim, nao, 1.0_dp, mos_new(2*i), &
2203 work, 0.0_dp, work1)
2204
2205 DO icol = 1, lcol
2206 op_fm_set(1, i)%local_data(:, icol) = &
2207 op_fm_set(1, i)%local_data(:, icol) + work1%local_data(:, icol + lcol)
2208 op_fm_set(2, i)%local_data(:, icol) = &
2209 op_fm_set(2, i)%local_data(:, icol) - work1%local_data(:, icol)
2210 END DO
2211
2212 CALL cp_fm_release(work)
2213 CALL cp_fm_release(work1)
2214 CALL cp_fm_struct_release(newstruct)
2215 CALL cp_fm_struct_release(newstruct1)
2216
2217 END DO
2218
2219 END SUBROUTINE op_orbbas_rtp
2220
2221! **************************************************************************************************
2222!> \brief ...
2223!> \param qs_env ...
2224!> \param magnetic ...
2225!> \param nmoments ...
2226!> \param reference ...
2227!> \param ref_point ...
2228!> \param unit_number ...
2229!> \param vel_reprs ...
2230!> \param com_nl ...
2231! **************************************************************************************************
2232 SUBROUTINE qs_moment_locop(qs_env, magnetic, nmoments, reference, ref_point, unit_number, vel_reprs, com_nl)
2233
2234 TYPE(qs_environment_type), POINTER :: qs_env
2235 LOGICAL, INTENT(IN) :: magnetic
2236 INTEGER, INTENT(IN) :: nmoments, reference
2237 REAL(dp), DIMENSION(:), INTENT(IN), POINTER :: ref_point
2238 INTEGER, INTENT(IN) :: unit_number
2239 LOGICAL, INTENT(IN), OPTIONAL :: vel_reprs, com_nl
2240
2241 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_moment_locop'
2242
2243 CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
2244 CHARACTER(LEN=default_string_length) :: description
2245 INTEGER :: akind, handle, i, ia, iatom, idir, &
2246 ikind, ispin, ix, iy, iz, l, nm, nmom, &
2247 order
2248 LOGICAL :: my_com_nl, my_velreprs
2249 REAL(dp) :: charge, dd, strace, trace
2250 REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
2251 nlcom_rv, nlcom_rvr, nlcom_rxrv, &
2252 qupole_der, rmom_vel
2253 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
2254 REAL(dp), DIMENSION(3) :: rcc, ria
2255 TYPE(atomic_kind_type), POINTER :: atomic_kind
2256 TYPE(cell_type), POINTER :: cell
2257 TYPE(cp_result_type), POINTER :: results
2258 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
2259 rho_ao
2260 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
2261 TYPE(dbcsr_type), POINTER :: tmp_ao
2262 TYPE(dft_control_type), POINTER :: dft_control
2263 TYPE(distribution_1d_type), POINTER :: local_particles
2264 TYPE(mp_para_env_type), POINTER :: para_env
2265 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2266 POINTER :: sab_all, sab_orb
2267 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2268 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2269 TYPE(qs_rho_type), POINTER :: rho
2270
2271 cpassert(ASSOCIATED(qs_env))
2272
2273 CALL timeset(routinen, handle)
2274
2275 my_velreprs = .false.
2276 IF (PRESENT(vel_reprs)) my_velreprs = vel_reprs
2277 IF (PRESENT(com_nl)) my_com_nl = com_nl
2278 IF (my_velreprs) CALL cite_reference(mattiat2019)
2279
2280 ! reference point
2281 CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
2282
2283 ! only allow for moments up to maxl set by basis
2284 nmom = min(nmoments, current_maxl)
2285 ! electronic contribution
2286 NULLIFY (dft_control, rho, cell, particle_set, qs_kind_set, results, para_env, matrix_s, rho_ao, sab_all, sab_orb)
2287 CALL get_qs_env(qs_env, &
2288 dft_control=dft_control, &
2289 rho=rho, &
2290 cell=cell, &
2291 results=results, &
2292 particle_set=particle_set, &
2293 qs_kind_set=qs_kind_set, &
2294 para_env=para_env, &
2295 matrix_s=matrix_s, &
2296 sab_all=sab_all, &
2297 sab_orb=sab_orb)
2298
2299 IF (my_com_nl) THEN
2300 IF ((nmom >= 1) .AND. my_velreprs) THEN
2301 ALLOCATE (nlcom_rv(3))
2302 nlcom_rv(:) = 0._dp
2303 END IF
2304 IF ((nmom >= 2) .AND. my_velreprs) THEN
2305 ALLOCATE (nlcom_rrv(6))
2306 nlcom_rrv(:) = 0._dp
2307 ALLOCATE (nlcom_rvr(6))
2308 nlcom_rvr(:) = 0._dp
2309 ALLOCATE (nlcom_rrv_vrr(6))
2310 nlcom_rrv_vrr(:) = 0._dp
2311 END IF
2312 IF (magnetic) THEN
2313 ALLOCATE (nlcom_rxrv(3))
2314 nlcom_rxrv = 0._dp
2315 END IF
2316 ! Calculate non local correction terms
2317 CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr, rcc)
2318 END IF
2319
2320 NULLIFY (moments)
2321 nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
2322 CALL dbcsr_allocate_matrix_set(moments, nm)
2323 DO i = 1, nm
2324 ALLOCATE (moments(i)%matrix)
2325 IF (my_velreprs .AND. (nmom >= 2)) THEN
2326 CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
2327 matrix_type=dbcsr_type_symmetric)
2328 CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
2329 ELSE
2330 CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
2331 END IF
2332 CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
2333 END DO
2334
2335 ! calculate derivatives if quadrupole in vel. reprs. is requested
2336 IF (my_velreprs .AND. (nmom >= 2)) THEN
2337 NULLIFY (moments_der)
2338 CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
2339 DO i = 1, 3 ! x, y, z
2340 DO idir = 1, 3 ! d/dx, d/dy, d/dz
2341 CALL dbcsr_init_p(moments_der(i, idir)%matrix)
2342 CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
2343 matrix_type=dbcsr_type_antisymmetric)
2344 CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
2345 CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
2346 END DO
2347 END DO
2348 CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, moments=moments)
2349 ELSE
2350 CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
2351 END IF
2352
2353 CALL qs_rho_get(rho, rho_ao=rho_ao)
2354
2355 ALLOCATE (rmom(nm + 1, 3))
2356 ALLOCATE (rlab(nm + 1))
2357 rmom = 0.0_dp
2358 rlab = ""
2359
2360 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2361 ! Allocate matrix to store the matrix product to be traced (dbcsr_dot only works for products of
2362 ! symmetric matrices)
2363 NULLIFY (tmp_ao)
2364 CALL dbcsr_init_p(tmp_ao)
2365 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2366 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2367 CALL dbcsr_set(tmp_ao, 0.0_dp)
2368 END IF
2369
2370 trace = 0.0_dp
2371 DO ispin = 1, dft_control%nspins
2372 CALL dbcsr_dot(rho_ao(ispin)%matrix, matrix_s(1)%matrix, trace)
2373 rmom(1, 1) = rmom(1, 1) + trace
2374 END DO
2375
2376 DO i = 1, SIZE(moments)
2377 strace = 0._dp
2378 DO ispin = 1, dft_control%nspins
2379 IF (my_velreprs .AND. nmoments >= 2) THEN
2380 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, moments(i)%matrix, &
2381 0.0_dp, tmp_ao)
2382 CALL dbcsr_trace(tmp_ao, trace)
2383 ELSE
2384 CALL dbcsr_dot(rho_ao(ispin)%matrix, moments(i)%matrix, trace)
2385 END IF
2386 strace = strace + trace
2387 END DO
2388 rmom(i + 1, 1) = strace
2389 END DO
2390
2391 CALL dbcsr_deallocate_matrix_set(moments)
2392
2393 ! nuclear contribution
2394 CALL get_qs_env(qs_env=qs_env, &
2395 local_particles=local_particles)
2396 DO ikind = 1, SIZE(local_particles%n_el)
2397 DO ia = 1, local_particles%n_el(ikind)
2398 iatom = local_particles%list(ikind)%array(ia)
2399 ! fold atomic positions back into unit cell
2400 ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
2401 ria = ria - rcc
2402 atomic_kind => particle_set(iatom)%atomic_kind
2403 CALL get_atomic_kind(atomic_kind, kind_number=akind)
2404 CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
2405 rmom(1, 2) = rmom(1, 2) - charge
2406 DO l = 1, nm
2407 ix = indco(1, l + 1)
2408 iy = indco(2, l + 1)
2409 iz = indco(3, l + 1)
2410 dd = 1._dp
2411 IF (ix > 0) dd = dd*ria(1)**ix
2412 IF (iy > 0) dd = dd*ria(2)**iy
2413 IF (iz > 0) dd = dd*ria(3)**iz
2414 rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
2415 CALL set_label(rlab(l + 1), ix, iy, iz)
2416 END DO
2417 END DO
2418 END DO
2419 CALL para_env%sum(rmom(:, 2))
2420 rmom(:, :) = -rmom(:, :)
2421 rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
2422
2423 ! magnetic moments
2424 IF (magnetic) THEN
2425 NULLIFY (magmom)
2426 CALL dbcsr_allocate_matrix_set(magmom, 3)
2427 DO i = 1, SIZE(magmom)
2428 CALL dbcsr_init_p(magmom(i)%matrix)
2429 CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
2430 matrix_type=dbcsr_type_antisymmetric)
2431 CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
2432 CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
2433 END DO
2434
2435 CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
2436
2437 ALLOCATE (mmom(SIZE(magmom)))
2438 mmom(:) = 0.0_dp
2439 IF (qs_env%run_rtp) THEN
2440 ! get imaginary part of the density in rho_ao (the real part is not needed since the trace of the product
2441 ! of a symmetric (REAL(rho_ao)) and an anti-symmetric (L_AO) matrix is zero)
2442 ! There may be other cases, where the imaginary part of the density is relevant
2443 NULLIFY (rho_ao)
2444 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2445 END IF
2446 ! if the density is purely real this is an expensive way to calculate zero
2447 DO i = 1, SIZE(magmom)
2448 strace = 0._dp
2449 DO ispin = 1, dft_control%nspins
2450 CALL dbcsr_set(tmp_ao, 0.0_dp)
2451 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, magmom(i)%matrix, &
2452 0.0_dp, tmp_ao)
2453 CALL dbcsr_trace(tmp_ao, trace)
2454 strace = strace + trace
2455 END DO
2456 mmom(i) = strace
2457 END DO
2458
2459 CALL dbcsr_deallocate_matrix_set(magmom)
2460 END IF
2461
2462 ! velocity representations
2463 IF (my_velreprs) THEN
2464 ALLOCATE (rmom_vel(nm))
2465 rmom_vel = 0.0_dp
2466
2467 DO order = 1, nmom
2468 SELECT CASE (order)
2469
2470 CASE (1) ! expectation value of momentum
2471 NULLIFY (momentum)
2472 CALL dbcsr_allocate_matrix_set(momentum, 3)
2473 DO i = 1, 3
2474 CALL dbcsr_init_p(momentum(i)%matrix)
2475 CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
2476 matrix_type=dbcsr_type_antisymmetric)
2477 CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
2478 CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
2479 END DO
2480 CALL build_lin_mom_matrix(qs_env, momentum)
2481
2482 ! imaginary part of the density for RTP, real part gives 0 since momentum is antisymmetric
2483 IF (qs_env%run_rtp) THEN
2484 NULLIFY (rho_ao)
2485 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2486 DO idir = 1, SIZE(momentum)
2487 strace = 0._dp
2488 DO ispin = 1, dft_control%nspins
2489 CALL dbcsr_set(tmp_ao, 0.0_dp)
2490 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, momentum(idir)%matrix, &
2491 0.0_dp, tmp_ao)
2492 CALL dbcsr_trace(tmp_ao, trace)
2493 strace = strace + trace
2494 END DO
2495 rmom_vel(idir) = rmom_vel(idir) + strace
2496 END DO
2497 END IF
2498
2499 CALL dbcsr_deallocate_matrix_set(momentum)
2500
2501 CASE (2) ! expectation value of quadrupole moment in vel. reprs.
2502 ALLOCATE (qupole_der(9)) ! will contain the expectation values of r_\alpha * d/d r_\beta
2503 qupole_der = 0._dp
2504
2505 NULLIFY (rho_ao)
2506 CALL qs_rho_get(rho, rho_ao=rho_ao)
2507
2508 ! Calculate expectation value over real part
2509 trace = 0._dp
2510 DO i = 1, 3
2511 DO idir = 1, 3
2512 strace = 0._dp
2513 DO ispin = 1, dft_control%nspins
2514 CALL dbcsr_set(tmp_ao, 0._dp)
2515 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2516 CALL dbcsr_trace(tmp_ao, trace)
2517 strace = strace + trace
2518 END DO
2519 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2520 END DO
2521 END DO
2522
2523 IF (qs_env%run_rtp) THEN
2524 NULLIFY (rho_ao)
2525 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
2526
2527 ! Calculate expectation value over imaginary part
2528 trace = 0._dp
2529 DO i = 1, 3
2530 DO idir = 1, 3
2531 strace = 0._dp
2532 DO ispin = 1, dft_control%nspins
2533 CALL dbcsr_set(tmp_ao, 0._dp)
2534 CALL dbcsr_multiply("T", "N", 1._dp, rho_ao(ispin)%matrix, moments_der(i, idir)%matrix, 0._dp, tmp_ao)
2535 CALL dbcsr_trace(tmp_ao, trace)
2536 strace = strace + trace
2537 END DO
2538 qupole_der((i - 1)*3 + idir) = qupole_der((i - 1)*3 + idir) + strace
2539 END DO
2540 END DO
2541 END IF
2542
2543 ! calculate vel. reprs. of quadrupole moment from derivatives
2544 rmom_vel(4) = -2*qupole_der(1) - rmom(1, 1)
2545 rmom_vel(5) = -qupole_der(2) - qupole_der(4)
2546 rmom_vel(6) = -qupole_der(3) - qupole_der(7)
2547 rmom_vel(7) = -2*qupole_der(5) - rmom(1, 1)
2548 rmom_vel(8) = -qupole_der(6) - qupole_der(8)
2549 rmom_vel(9) = -2*qupole_der(9) - rmom(1, 1)
2550
2551 DEALLOCATE (qupole_der)
2552 CASE DEFAULT
2553 END SELECT
2554 END DO
2555 END IF
2556
2557 IF ((my_velreprs .AND. (nmoments >= 1)) .OR. magnetic) THEN
2558 CALL dbcsr_deallocate_matrix(tmp_ao)
2559 END IF
2560 IF (my_velreprs .AND. (nmoments >= 2)) THEN
2561 CALL dbcsr_deallocate_matrix_set(moments_der)
2562 END IF
2563
2564 description = "[DIPOLE]"
2565 CALL cp_results_erase(results=results, description=description)
2566 CALL put_results(results=results, description=description, &
2567 values=rmom(2:4, 3))
2568
2569 IF (magnetic .AND. my_velreprs) THEN
2570 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom, rmom_vel=rmom_vel)
2571 ELSE IF (magnetic .AND. .NOT. my_velreprs) THEN
2572 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., mmom=mmom)
2573 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2574 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false., rmom_vel=rmom_vel)
2575 ELSE
2576 CALL print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic=.false.)
2577 END IF
2578
2579 IF (my_com_nl) THEN
2580 IF (magnetic) THEN
2581 mmom(:) = nlcom_rxrv(:)
2582 END IF
2583 IF (my_velreprs .AND. (nmom >= 1)) THEN
2584 DEALLOCATE (rmom_vel)
2585 ALLOCATE (rmom_vel(21))
2586 rmom_vel(1:3) = nlcom_rv
2587 END IF
2588 IF (my_velreprs .AND. (nmom >= 2)) THEN
2589 rmom_vel(4:9) = nlcom_rrv
2590 rmom_vel(10:15) = nlcom_rvr
2591 rmom_vel(16:21) = nlcom_rrv_vrr
2592 END IF
2593 IF (magnetic .AND. .NOT. my_velreprs) THEN
2594 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom)
2595 ELSE IF (my_velreprs .AND. .NOT. magnetic) THEN
2596 CALL print_moments_nl(unit_number, nmom, rlab, rmom_vel=rmom_vel)
2597 ELSE IF (my_velreprs .AND. magnetic) THEN
2598 CALL print_moments_nl(unit_number, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
2599 END IF
2600
2601 END IF
2602
2603 IF (my_com_nl) THEN
2604 IF (nmom >= 1 .AND. my_velreprs) DEALLOCATE (nlcom_rv)
2605 IF (nmom >= 2 .AND. my_velreprs) THEN
2606 DEALLOCATE (nlcom_rrv)
2607 DEALLOCATE (nlcom_rvr)
2608 DEALLOCATE (nlcom_rrv_vrr)
2609 END IF
2610 IF (magnetic) DEALLOCATE (nlcom_rxrv)
2611 END IF
2612
2613 DEALLOCATE (rmom)
2614 DEALLOCATE (rlab)
2615 IF (magnetic) THEN
2616 DEALLOCATE (mmom)
2617 END IF
2618 IF (my_velreprs) THEN
2619 DEALLOCATE (rmom_vel)
2620 END IF
2621
2622 CALL timestop(handle)
2623
2624 END SUBROUTINE qs_moment_locop
2625
2626! **************************************************************************************************
2627!> \brief ...
2628!> \param label ...
2629!> \param ix ...
2630!> \param iy ...
2631!> \param iz ...
2632! **************************************************************************************************
2633 SUBROUTINE set_label(label, ix, iy, iz)
2634 CHARACTER(LEN=*), INTENT(OUT) :: label
2635 INTEGER, INTENT(IN) :: ix, iy, iz
2636
2637 INTEGER :: i
2638
2639 label = ""
2640 DO i = 1, ix
2641 WRITE (label(i:), "(A1)") "X"
2642 END DO
2643 DO i = ix + 1, ix + iy
2644 WRITE (label(i:), "(A1)") "Y"
2645 END DO
2646 DO i = ix + iy + 1, ix + iy + iz
2647 WRITE (label(i:), "(A1)") "Z"
2648 END DO
2649
2650 END SUBROUTINE set_label
2651
2652! **************************************************************************************************
2653!> \brief ...
2654!> \param unit_number ...
2655!> \param nmom ...
2656!> \param rmom ...
2657!> \param rlab ...
2658!> \param rcc ...
2659!> \param cell ...
2660!> \param periodic ...
2661!> \param mmom ...
2662!> \param rmom_vel ...
2663! **************************************************************************************************
2664 SUBROUTINE print_moments(unit_number, nmom, rmom, rlab, rcc, cell, periodic, mmom, rmom_vel)
2665 INTEGER, INTENT(IN) :: unit_number, nmom
2666 REAL(dp), DIMENSION(:, :), INTENT(IN) :: rmom
2667 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2668 REAL(dp), DIMENSION(3), INTENT(IN) :: rcc
2669 TYPE(cell_type), POINTER :: cell
2670 LOGICAL :: periodic
2671 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2672
2673 INTEGER :: i, i0, i1, j, l
2674 REAL(dp) :: dd
2675
2676 IF (unit_number > 0) THEN
2677 DO l = 0, nmom
2678 SELECT CASE (l)
2679 CASE (0)
2680 WRITE (unit_number, "(T3,A,T33,3F16.8)") "Reference Point [Bohr]", rcc
2681 WRITE (unit_number, "(T3,A)") "Charges"
2682 WRITE (unit_number, "(T5,A,T18,F14.8,T36,A,T42,F14.8,T60,A,T67,F14.8)") &
2683 "Electronic=", rmom(1, 1), "Core=", rmom(1, 2), "Total=", rmom(1, 3)
2684 CASE (1)
2685 IF (periodic) THEN
2686 WRITE (unit_number, "(T3,A)") "Dipole vectors are based on the periodic (Berry phase) operator."
2687 WRITE (unit_number, "(T3,A)") "They are defined modulo integer multiples of the cell matrix [Debye]."
2688 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[X] [", cell%hmat(1, :)*debye, "] [i]"
2689 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Y]=[", cell%hmat(2, :)*debye, "]*[j]"
2690 WRITE (unit_number, "(T3,A,3(F14.8,1X),A)") "[Z] [", cell%hmat(3, :)*debye, "] [k]"
2691 ELSE
2692 WRITE (unit_number, "(T3,A)") "Dipoles are based on the traditional operator."
2693 END IF
2694 dd = sqrt(sum(rmom(2:4, 3)**2))*debye
2695 WRITE (unit_number, "(T3,A)") "Dipole moment [Debye]"
2696 WRITE (unit_number, "(T5,3(A,A,E15.7,1X),T60,A,T68,F13.7)") &
2697 (trim(rlab(i)), "=", rmom(i, 3)*debye, i=2, 4), "Total=", dd
2698 CASE (2)
2699 WRITE (unit_number, "(T3,A)") "Quadrupole moment [Debye*Angstrom]"
2700 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2701 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=5, 7)
2702 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2703 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr, i=8, 10)
2704 CASE (3)
2705 WRITE (unit_number, "(T3,A)") "Octapole moment [Debye*Angstrom**2]"
2706 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2707 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=11, 14)
2708 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2709 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=15, 18)
2710 WRITE (unit_number, "(T7,4(A,A,F14.8,3X))") &
2711 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr, i=19, 20)
2712 CASE (4)
2713 WRITE (unit_number, "(T3,A)") "Hexadecapole moment [Debye*Angstrom**3]"
2714 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2715 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=21, 24)
2716 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2717 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=25, 28)
2718 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2719 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=29, 32)
2720 WRITE (unit_number, "(T6,4(A,A,F14.8,2X))") &
2721 (trim(rlab(i)), "=", rmom(i, 3)*debye/bohr/bohr/bohr, i=32, 35)
2722 CASE DEFAULT
2723 WRITE (unit_number, "(T3,A,A,I2)") "Higher moment [Debye*Angstrom**(L-1)]", &
2724 " L=", l
2725 i0 = (6 + 11*(l - 1) + 6*(l - 1)**2 + (l - 1)**3)/6
2726 i1 = (6 + 11*l + 6*l**2 + l**3)/6 - 1
2727 dd = debye/(bohr)**(l - 1)
2728 DO i = i0, i1, 3
2729 WRITE (unit_number, "(T18,3(A,A,F14.8,4X))") &
2730 (trim(rlab(j + 1)), "=", rmom(j + 1, 3)*dd, j=i, min(i1, i + 2))
2731 END DO
2732 END SELECT
2733 END DO
2734 IF (PRESENT(mmom)) THEN
2735 IF (nmom >= 1) THEN
2736 dd = sqrt(sum(mmom(1:3)**2))
2737 WRITE (unit_number, "(T3,A)") "Orbital angular momentum [a. u.]"
2738 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2739 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2740 END IF
2741 END IF
2742 IF (PRESENT(rmom_vel)) THEN
2743 DO l = 1, nmom
2744 SELECT CASE (l)
2745 CASE (1)
2746 dd = sqrt(sum(rmom_vel(1:3)**2))
2747 WRITE (unit_number, "(T3,A)") "Expectation value of momentum operator [a. u.]"
2748 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2749 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2750 CASE (2)
2751 WRITE (unit_number, "(T3,A)") "Expectation value of quadrupole operator in vel. repr. [a. u.]"
2752 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2753 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2754 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2755 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2756 CASE DEFAULT
2757 END SELECT
2758 END DO
2759 END IF
2760 END IF
2761
2762 END SUBROUTINE print_moments
2763
2764! **************************************************************************************************
2765!> \brief ...
2766!> \param unit_number ...
2767!> \param nmom ...
2768!> \param rlab ...
2769!> \param mmom ...
2770!> \param rmom_vel ...
2771! **************************************************************************************************
2772 SUBROUTINE print_moments_nl(unit_number, nmom, rlab, mmom, rmom_vel)
2773 INTEGER, INTENT(IN) :: unit_number, nmom
2774 CHARACTER(LEN=8), DIMENSION(:) :: rlab
2775 REAL(dp), DIMENSION(:), INTENT(IN), OPTIONAL :: mmom, rmom_vel
2776
2777 INTEGER :: i, l
2778 REAL(dp) :: dd
2779
2780 IF (unit_number > 0) THEN
2781 IF (PRESENT(mmom)) THEN
2782 IF (nmom >= 1) THEN
2783 dd = sqrt(sum(mmom(1:3)**2))
2784 WRITE (unit_number, "(T3,A)") "Expectation value of rx[r,V_nl] [a. u.]"
2785 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2786 (trim(rlab(i + 1)), "=", mmom(i), i=1, 3), "Total=", dd
2787 END IF
2788 END IF
2789 IF (PRESENT(rmom_vel)) THEN
2790 DO l = 1, nmom
2791 SELECT CASE (l)
2792 CASE (1)
2793 dd = sqrt(sum(rmom_vel(1:3)**2))
2794 WRITE (unit_number, "(T3,A)") "Expectation value of [r,V_nl] [a. u.]"
2795 WRITE (unit_number, "(T5,3(A,A,E16.8,1X),T64,A,T71,F14.8)") &
2796 (trim(rlab(i + 1)), "=", rmom_vel(i), i=1, 3), "Total=", dd
2797 CASE (2)
2798 WRITE (unit_number, "(T3,A)") "Expectation value of [rr,V_nl] [a. u.]"
2799 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2800 (trim(rlab(i + 1)), "=", rmom_vel(i), i=4, 6)
2801 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2802 (trim(rlab(i + 1)), "=", rmom_vel(i), i=7, 9)
2803 WRITE (unit_number, "(T3,A)") "Expectation value of r x V_nl x r [a. u.]"
2804 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2805 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=10, 12)
2806 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2807 (trim(rlab(i + 1 - 6)), "=", rmom_vel(i), i=13, 15)
2808 WRITE (unit_number, "(T3,A)") "Expectation value of r x r x V_nl + V_nl x r x r [a. u.]"
2809 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2810 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=16, 18)
2811 WRITE (unit_number, "(T17,3(A,A,E16.8,9X))") &
2812 (trim(rlab(i + 1 - 12)), "=", rmom_vel(i), i=19, 21)
2813 CASE DEFAULT
2814 END SELECT
2815 END DO
2816 END IF
2817 END IF
2818
2819 END SUBROUTINE print_moments_nl
2820
2821! **************************************************************************************************
2822!> \brief Calculate the expectation value of operators related to non-local potential:
2823!> [r, Vnl], noted rv
2824!> r x [r,Vnl], noted rxrv
2825!> [rr,Vnl], noted rrv
2826!> r x Vnl x r, noted rvr
2827!> r x r x Vnl + Vnl x r x r, noted rrv_vrr
2828!> Note that the 3 first operator are commutator while the 2 last
2829!> are not. For reading clarity the same notation is used for all 5
2830!> operators.
2831!> \param qs_env ...
2832!> \param nlcom_rv ...
2833!> \param nlcom_rxrv ...
2834!> \param nlcom_rrv ...
2835!> \param nlcom_rvr ...
2836!> \param nlcom_rrv_vrr ...
2837!> \param ref_point ...
2838! **************************************************************************************************
2839 SUBROUTINE calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, nlcom_rvr, &
2840 nlcom_rrv_vrr, ref_point)
2841
2842 TYPE(qs_environment_type), POINTER :: qs_env
2843 REAL(dp), ALLOCATABLE, DIMENSION(:), OPTIONAL :: nlcom_rv, nlcom_rxrv, nlcom_rrv, &
2844 nlcom_rvr, nlcom_rrv_vrr
2845 REAL(dp), DIMENSION(3) :: ref_point
2846
2847 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_commutator_nl_terms'
2848
2849 INTEGER :: handle, ind, ispin
2850 LOGICAL :: calc_rrv, calc_rrv_vrr, calc_rv, &
2851 calc_rvr, calc_rxrv
2852 REAL(dp) :: eps_ppnl, strace, trace
2853 TYPE(cell_type), POINTER :: cell
2854 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_rrv, matrix_rrv_vrr, matrix_rv, &
2855 matrix_rvr, matrix_rxrv, matrix_s, &
2856 rho_ao
2857 TYPE(dbcsr_type), POINTER :: tmp_ao
2858 TYPE(dft_control_type), POINTER :: dft_control
2859 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2860 POINTER :: sab_all, sab_orb, sap_ppnl
2861 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2862 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2863 TYPE(qs_rho_type), POINTER :: rho
2864
2865 CALL timeset(routinen, handle)
2866
2867 calc_rv = .false.
2868 calc_rxrv = .false.
2869 calc_rrv = .false.
2870 calc_rvr = .false.
2871 calc_rrv_vrr = .false.
2872
2873 ! rv, rxrv and rrv are commutator matrices: anti-symmetric.
2874 ! The real part of the density matrix rho_ao is symmetric so that
2875 ! the expectation value of real density matrix is zero. Hence, if
2876 ! the density matrix is real, no need to compute these quantities.
2877 ! This is not the case for rvr and rrv_vrr which are symmetric.
2878
2879 IF (ALLOCATED(nlcom_rv)) THEN
2880 nlcom_rv(:) = 0._dp
2881 IF (qs_env%run_rtp) calc_rv = .true.
2882 END IF
2883 IF (ALLOCATED(nlcom_rxrv)) THEN
2884 nlcom_rxrv(:) = 0._dp
2885 IF (qs_env%run_rtp) calc_rxrv = .true.
2886 END IF
2887 IF (ALLOCATED(nlcom_rrv)) THEN
2888 nlcom_rrv(:) = 0._dp
2889 IF (qs_env%run_rtp) calc_rrv = .true.
2890 END IF
2891 IF (ALLOCATED(nlcom_rvr)) THEN
2892 nlcom_rvr(:) = 0._dp
2893 calc_rvr = .true.
2894 END IF
2895 IF (ALLOCATED(nlcom_rrv_vrr)) THEN
2896 nlcom_rrv_vrr(:) = 0._dp
2897 calc_rrv_vrr = .true.
2898 END IF
2899
2900 IF (.NOT. (calc_rv .OR. calc_rrv .OR. calc_rxrv &
2901 .OR. calc_rvr .OR. calc_rrv_vrr)) RETURN
2902
2903 NULLIFY (cell, matrix_s, particle_set, qs_kind_set, rho, sab_all, sab_orb, sap_ppnl)
2904 CALL get_qs_env(qs_env, &
2905 cell=cell, &
2906 dft_control=dft_control, &
2907 matrix_s=matrix_s, &
2908 particle_set=particle_set, &
2909 qs_kind_set=qs_kind_set, &
2910 rho=rho, &
2911 sab_orb=sab_orb, &
2912 sab_all=sab_all, &
2913 sap_ppnl=sap_ppnl)
2914
2915 eps_ppnl = dft_control%qs_control%eps_ppnl
2916
2917 ! Allocate storage
2918 NULLIFY (matrix_rv, matrix_rxrv, matrix_rrv, matrix_rvr, matrix_rrv_vrr)
2919 IF (calc_rv) THEN
2920 CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
2921 DO ind = 1, 3
2922 CALL dbcsr_init_p(matrix_rv(ind)%matrix)
2923 CALL dbcsr_create(matrix_rv(ind)%matrix, template=matrix_s(1)%matrix, &
2924 matrix_type=dbcsr_type_antisymmetric)
2925 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(ind)%matrix, sab_orb)
2926 CALL dbcsr_set(matrix_rv(ind)%matrix, 0._dp)
2927 END DO
2928 END IF
2929
2930 IF (calc_rxrv) THEN
2931 CALL dbcsr_allocate_matrix_set(matrix_rxrv, 3)
2932 DO ind = 1, 3
2933 CALL dbcsr_init_p(matrix_rxrv(ind)%matrix)
2934 CALL dbcsr_create(matrix_rxrv(ind)%matrix, template=matrix_s(1)%matrix, &
2935 matrix_type=dbcsr_type_antisymmetric)
2936 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rxrv(ind)%matrix, sab_orb)
2937 CALL dbcsr_set(matrix_rxrv(ind)%matrix, 0._dp)
2938 END DO
2939 END IF
2940
2941 IF (calc_rrv) THEN
2942 CALL dbcsr_allocate_matrix_set(matrix_rrv, 6)
2943 DO ind = 1, 6
2944 CALL dbcsr_init_p(matrix_rrv(ind)%matrix)
2945 CALL dbcsr_create(matrix_rrv(ind)%matrix, template=matrix_s(1)%matrix, &
2946 matrix_type=dbcsr_type_antisymmetric)
2947 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv(ind)%matrix, sab_orb)
2948 CALL dbcsr_set(matrix_rrv(ind)%matrix, 0._dp)
2949 END DO
2950 END IF
2951
2952 IF (calc_rvr) THEN
2953 CALL dbcsr_allocate_matrix_set(matrix_rvr, 6)
2954 DO ind = 1, 6
2955 CALL dbcsr_init_p(matrix_rvr(ind)%matrix)
2956 CALL dbcsr_create(matrix_rvr(ind)%matrix, template=matrix_s(1)%matrix, &
2957 matrix_type=dbcsr_type_symmetric)
2958 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rvr(ind)%matrix, sab_orb)
2959 CALL dbcsr_set(matrix_rvr(ind)%matrix, 0._dp)
2960 END DO
2961 END IF
2962 IF (calc_rrv_vrr) THEN
2963 CALL dbcsr_allocate_matrix_set(matrix_rrv_vrr, 6)
2964 DO ind = 1, 6
2965 CALL dbcsr_init_p(matrix_rrv_vrr(ind)%matrix)
2966 CALL dbcsr_create(matrix_rrv_vrr(ind)%matrix, template=matrix_s(1)%matrix, &
2967 matrix_type=dbcsr_type_symmetric)
2968 CALL cp_dbcsr_alloc_block_from_nbl(matrix_rrv_vrr(ind)%matrix, sab_orb)
2969 CALL dbcsr_set(matrix_rrv_vrr(ind)%matrix, 0._dp)
2970 END DO
2971 END IF
2972
2973 ! calculate evaluation of operators in AO basis set
2974 CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv, &
2975 matrix_rxrv=matrix_rxrv, matrix_rrv=matrix_rrv, matrix_rvr=matrix_rvr, &
2976 matrix_rrv_vrr=matrix_rrv_vrr, ref_point=ref_point)
2977
2978 ! Calculate expectation values
2979 ! Real part
2980 NULLIFY (tmp_ao)
2981 CALL dbcsr_init_p(tmp_ao)
2982 CALL dbcsr_create(tmp_ao, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry, name="tmp")
2983 CALL cp_dbcsr_alloc_block_from_nbl(tmp_ao, sab_all)
2984 CALL dbcsr_set(tmp_ao, 0.0_dp)
2985
2986 IF (calc_rvr .OR. calc_rrv_vrr) THEN
2987 NULLIFY (rho_ao)
2988 CALL qs_rho_get(rho, rho_ao=rho_ao)
2989
2990 IF (calc_rvr) THEN
2991 trace = 0._dp
2992 DO ind = 1, SIZE(matrix_rvr)
2993 strace = 0._dp
2994 DO ispin = 1, dft_control%nspins
2995 CALL dbcsr_set(tmp_ao, 0.0_dp)
2996 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rvr(ind)%matrix, &
2997 0.0_dp, tmp_ao)
2998 CALL dbcsr_trace(tmp_ao, trace)
2999 strace = strace + trace
3000 END DO
3001 nlcom_rvr(ind) = nlcom_rvr(ind) + strace
3002 END DO
3003 END IF
3004
3005 IF (calc_rrv_vrr) THEN
3006 trace = 0._dp
3007 DO ind = 1, SIZE(matrix_rrv_vrr)
3008 strace = 0._dp
3009 DO ispin = 1, dft_control%nspins
3010 CALL dbcsr_set(tmp_ao, 0.0_dp)
3011 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv_vrr(ind)%matrix, &
3012 0.0_dp, tmp_ao)
3013 CALL dbcsr_trace(tmp_ao, trace)
3014 strace = strace + trace
3015 END DO
3016 nlcom_rrv_vrr(ind) = nlcom_rrv_vrr(ind) + strace
3017 END DO
3018 END IF
3019 END IF
3020
3021 ! imagninary part of the density matrix
3022 NULLIFY (rho_ao)
3023 CALL qs_rho_get(rho, rho_ao_im=rho_ao)
3024
3025 IF (calc_rv) THEN
3026 trace = 0._dp
3027 DO ind = 1, SIZE(matrix_rv)
3028 strace = 0._dp
3029 DO ispin = 1, dft_control%nspins
3030 CALL dbcsr_set(tmp_ao, 0.0_dp)
3031 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rv(ind)%matrix, &
3032 0.0_dp, tmp_ao)
3033 CALL dbcsr_trace(tmp_ao, trace)
3034 strace = strace + trace
3035 END DO
3036 nlcom_rv(ind) = nlcom_rv(ind) + strace
3037 END DO
3038 END IF
3039
3040 IF (calc_rrv) THEN
3041 trace = 0._dp
3042 DO ind = 1, SIZE(matrix_rrv)
3043 strace = 0._dp
3044 DO ispin = 1, dft_control%nspins
3045 CALL dbcsr_set(tmp_ao, 0.0_dp)
3046 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rrv(ind)%matrix, &
3047 0.0_dp, tmp_ao)
3048 CALL dbcsr_trace(tmp_ao, trace)
3049 strace = strace + trace
3050 END DO
3051 nlcom_rrv(ind) = nlcom_rrv(ind) + strace
3052 END DO
3053 END IF
3054
3055 IF (calc_rxrv) THEN
3056 trace = 0._dp
3057 DO ind = 1, SIZE(matrix_rxrv)
3058 strace = 0._dp
3059 DO ispin = 1, dft_control%nspins
3060 CALL dbcsr_set(tmp_ao, 0.0_dp)
3061 CALL dbcsr_multiply("T", "N", 1.0_dp, rho_ao(ispin)%matrix, matrix_rxrv(ind)%matrix, &
3062 0.0_dp, tmp_ao)
3063 CALL dbcsr_trace(tmp_ao, trace)
3064 strace = strace + trace
3065 END DO
3066 nlcom_rxrv(ind) = nlcom_rxrv(ind) + strace
3067 END DO
3068 END IF
3069 CALL dbcsr_deallocate_matrix(tmp_ao)
3070 IF (calc_rv) CALL dbcsr_deallocate_matrix_set(matrix_rv)
3071 IF (calc_rxrv) CALL dbcsr_deallocate_matrix_set(matrix_rxrv)
3072 IF (calc_rrv) CALL dbcsr_deallocate_matrix_set(matrix_rrv)
3073 IF (calc_rvr) CALL dbcsr_deallocate_matrix_set(matrix_rvr)
3074 IF (calc_rrv_vrr) CALL dbcsr_deallocate_matrix_set(matrix_rrv_vrr)
3075
3076 CALL timestop(handle)
3077 END SUBROUTINE calculate_commutator_nl_terms
3078
3079! *****************************************************************************
3080!> \brief ...
3081!> \param qs_env ...
3082!> \param difdip ...
3083!> \param deltaR ...
3084!> \param order ...
3085!> \param rcc ...
3086!> \note calculate matrix elements <a|r_beta |db/dR_alpha > + <da/dR_alpha | r_beta | b >
3087!> be aware: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > only valid
3088!> if alpha .neq.beta
3089!> if alpha=beta: < a | r_beta| db/dR_alpha > = - < da/dR_alpha | r_beta | b > - < a | b >
3090!> modified from qs_efield_mo_derivatives
3091!> SL July 2015
3092! **************************************************************************************************
3093 SUBROUTINE dipole_deriv_ao(qs_env, difdip, deltaR, order, rcc)
3094 TYPE(qs_environment_type), POINTER :: qs_env
3095 TYPE(dbcsr_p_type), DIMENSION(:, :), &
3096 INTENT(INOUT), POINTER :: difdip
3097 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
3098 POINTER :: deltar
3099 INTEGER, INTENT(IN) :: order
3100 REAL(kind=dp), DIMENSION(3), OPTIONAL :: rcc
3101
3102 CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_deriv_ao'
3103
3104 INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
3105 last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
3106 sgfb
3107 INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3108 npgfb, nsgfa, nsgfb
3109 INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
3110 LOGICAL :: found
3111 REAL(dp) :: dab
3112 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
3113 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
3114 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
3115 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
3116 REAL(kind=dp), DIMENSION(:, :), POINTER :: rpgfa, rpgfb, sphi_a, sphi_b, zeta, zetb
3117 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab
3118 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: difmab2
3119 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
3120 TYPE(cell_type), POINTER :: cell
3121 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
3122 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
3123 TYPE(neighbor_list_iterator_p_type), &
3124 DIMENSION(:), POINTER :: nl_iterator
3125 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
3126 POINTER :: sab_all, sab_orb
3127 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
3128 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
3129 TYPE(qs_kind_type), POINTER :: qs_kind
3130
3131 CALL timeset(routinen, handle)
3132
3133 NULLIFY (cell, particle_set, qs_kind_set, sab_orb, sab_all)
3134 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
3135 qs_kind_set=qs_kind_set, sab_orb=sab_orb, sab_all=sab_all)
3136 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
3137 maxco=ldab, maxsgf=maxsgf)
3138
3139 nkind = SIZE(qs_kind_set)
3140 natom = SIZE(particle_set)
3141
3142 m_dim = ncoset(order) - 1
3143
3144 IF (PRESENT(rcc)) THEN
3145 rc = rcc
3146 ELSE
3147 rc = 0._dp
3148 END IF
3149
3150 ALLOCATE (basis_set_list(nkind))
3151
3152 ALLOCATE (mab(ldab, ldab, m_dim))
3153 ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
3154 ALLOCATE (work(ldab, maxsgf))
3155 ALLOCATE (mint(3, 3))
3156 ALLOCATE (mint2(3, 3))
3157
3158 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
3159 difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
3160 work(1:ldab, 1:maxsgf) = 0.0_dp
3161
3162 DO i = 1, 3
3163 DO j = 1, 3
3164 NULLIFY (mint(i, j)%block)
3165 NULLIFY (mint2(i, j)%block)
3166 END DO
3167 END DO
3168
3169 ! Set the basis_set_list(nkind) to point to the corresponding basis sets
3170 DO ikind = 1, nkind
3171 qs_kind => qs_kind_set(ikind)
3172 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
3173 IF (ASSOCIATED(basis_set_a)) THEN
3174 basis_set_list(ikind)%gto_basis_set => basis_set_a
3175 ELSE
3176 NULLIFY (basis_set_list(ikind)%gto_basis_set)
3177 END IF
3178 END DO
3179
3180 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
3181 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
3182 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
3183 iatom=iatom, jatom=jatom, r=rab)
3184
3185 basis_set_a => basis_set_list(ikind)%gto_basis_set
3186 basis_set_b => basis_set_list(jkind)%gto_basis_set
3187 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
3188 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
3189
3190 ! basis ikind
3191 first_sgfa => basis_set_a%first_sgf
3192 la_max => basis_set_a%lmax
3193 la_min => basis_set_a%lmin
3194 npgfa => basis_set_a%npgf
3195 nseta = basis_set_a%nset
3196 nsgfa => basis_set_a%nsgf_set
3197 rpgfa => basis_set_a%pgf_radius
3198 set_radius_a => basis_set_a%set_radius
3199 sphi_a => basis_set_a%sphi
3200 zeta => basis_set_a%zet
3201 ! basis jkind
3202 first_sgfb => basis_set_b%first_sgf
3203 lb_max => basis_set_b%lmax
3204 lb_min => basis_set_b%lmin
3205 npgfb => basis_set_b%npgf
3206 nsetb = basis_set_b%nset
3207 nsgfb => basis_set_b%nsgf_set
3208 rpgfb => basis_set_b%pgf_radius
3209 set_radius_b => basis_set_b%set_radius
3210 sphi_b => basis_set_b%sphi
3211 zetb => basis_set_b%zet
3212
3213 IF (inode == 1) last_jatom = 0
3214
3215 ! this guarentees minimum image convention
3216 ! anything else would not make sense
3217 IF (jatom == last_jatom) THEN
3218 cycle
3219 END IF
3220
3221 last_jatom = jatom
3222
3223 irow = iatom
3224 icol = jatom
3225
3226 DO i = 1, 3
3227 DO j = 1, 3
3228 NULLIFY (mint(i, j)%block)
3229 CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
3230 row=irow, col=icol, block=mint(i, j)%block, &
3231 found=found)
3232 cpassert(found)
3233 END DO
3234 END DO
3235
3236 ra(:) = particle_set(iatom)%r(:)
3237 rb(:) = particle_set(jatom)%r(:)
3238 rab(:) = pbc(rb, ra, cell)
3239 rac(:) = pbc(ra - rc, cell)
3240 rbc(:) = pbc(rb - rc, cell)
3241 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
3242
3243 DO iset = 1, nseta
3244 ncoa = npgfa(iset)*ncoset(la_max(iset))
3245 sgfa = first_sgfa(1, iset)
3246 DO jset = 1, nsetb
3247 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
3248 ncob = npgfb(jset)*ncoset(lb_max(jset))
3249 sgfb = first_sgfb(1, jset)
3250 ldab = max(ncoa, ncob)
3251 lda = ncoset(la_max(iset))*npgfa(iset)
3252 ldb = ncoset(lb_max(jset))*npgfb(jset)
3253 ALLOCATE (difmab(lda, ldb, m_dim, 3))
3254
3255 ! Calculate integral (da|r|b)
3256 CALL diff_momop2(la_max(iset), npgfa(iset), zeta(:, iset), &
3257 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
3258 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
3259 difmab, deltar=deltar, iatom=iatom, jatom=jatom)
3260
3261! *** Contraction step ***
3262
3263 DO idir = 1, 3 ! derivative of AO function
3264 DO j = 1, 3 ! position operator r_j
3265 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
3266 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
3267 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
3268 0.0_dp, work(1, 1), SIZE(work, 1))
3269
3270 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
3271 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
3272 work(1, 1), SIZE(work, 1), &
3273 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
3274 SIZE(mint(j, idir)%block, 1))
3275 END DO !j
3276 END DO !idir
3277 DEALLOCATE (difmab)
3278 END DO !jset
3279 END DO !iset
3280 END do!iterator
3281
3282 CALL neighbor_list_iterator_release(nl_iterator)
3283
3284 DO i = 1, 3
3285 DO j = 1, 3
3286 NULLIFY (mint(i, j)%block)
3287 END DO
3288 END DO
3289
3290 DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
3291
3292 CALL timestop(handle)
3293 END SUBROUTINE dipole_deriv_ao
3294
3295END 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:206
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:192
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:87
cp_result_methods::put_results
Definition cp_result_methods.F:39
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
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
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_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:121
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:155
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:603
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:233
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
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:679
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:1073
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:205
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
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:1009
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:528
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:164
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
kinds
Defines the basic variable types.
Definition kinds.F:23
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_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:365
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
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:98
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::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, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, 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, 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:454
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)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:880
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, 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, 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, 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:333
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:397
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:1125
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:128
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:560
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:2233
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:1338
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:792
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
...
Definition qs_moments.F:1505
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:331
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:1712
qs_moments::dipole_deriv_ao
subroutine, public dipole_deriv_ao(qs_env, difdip, deltar, order, rcc)
...
Definition qs_moments.F:3094
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:408
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:55
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:68
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:180
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:174
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:113
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
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:150
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
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:136
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
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