(git:374b731)
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-2024 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
39 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
40 dbcsr_allocate_matrix_set,&
41 dbcsr_deallocate_matrix_set
42 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
43 cp_fm_struct_double,&
44 cp_fm_struct_release,&
45 cp_fm_struct_type
46 USE cp_fm_types, ONLY: cp_fm_create,&
47 cp_fm_get_info,&
48 cp_fm_release,&
49 cp_fm_set_all,&
50 cp_fm_type
51 USE cp_result_methods, ONLY: cp_results_erase,&
52 put_results
53 USE cp_result_types, ONLY: cp_result_type
54 USE dbcsr_api, ONLY: &
55 dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_dot, &
56 dbcsr_get_block_p, dbcsr_init_p, dbcsr_multiply, dbcsr_p_type, dbcsr_set, dbcsr_trace, &
57 dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
58 USE distribution_1d_types, ONLY: distribution_1d_type
59 USE kinds, ONLY: default_string_length,&
60 dp
61 USE kpoint_types, ONLY: get_kpoint_info,&
62 kpoint_type
63 USE mathconstants, ONLY: pi,&
64 twopi
65 USE message_passing, ONLY: mp_para_env_type
66 USE moments_utils, ONLY: get_reference_point
67 USE orbital_pointers, ONLY: current_maxl,&
68 indco,&
69 ncoset
70 USE parallel_gemm_api, ONLY: parallel_gemm
71 USE particle_methods, ONLY: get_particle_set
72 USE particle_types, ONLY: particle_type
73 USE physcon, ONLY: bohr,&
74 debye
75 USE qs_environment_types, ONLY: get_qs_env,&
76 qs_environment_type
77 USE qs_kind_types, ONLY: get_qs_kind,&
78 get_qs_kind_set,&
79 qs_kind_type
80 USE qs_ks_types, ONLY: get_ks_env,&
81 qs_ks_env_type
82 USE qs_mo_types, ONLY: get_mo_set,&
83 mo_set_type
84 USE qs_neighbor_list_types, ONLY: get_iterator_info,&
85 neighbor_list_iterate,&
86 neighbor_list_iterator_create,&
87 neighbor_list_iterator_p_type,&
88 neighbor_list_iterator_release,&
89 neighbor_list_set_p_type
90 USE qs_operators_ao, ONLY: build_lin_mom_matrix
91 USE qs_rho_types, ONLY: qs_rho_get,&
92 qs_rho_type
93 USE rt_propagation_types, ONLY: get_rtp,&
94 rt_prop_type
95#include "./base/base_uses.f90"
96
97 IMPLICIT NONE
98
99 PRIVATE
100
101 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_moments'
102
103 ! Public subroutines
104 PUBLIC :: build_berry_moment_matrix, build_local_moment_matrix
105 PUBLIC :: build_berry_kpoint_matrix, build_local_magmom_matrix
106 PUBLIC :: qs_moment_berry_phase, qs_moment_locop
107 PUBLIC :: dipole_deriv_ao
108 PUBLIC :: build_local_moments_der_matrix
109 PUBLIC :: build_dsdv_moments
110 PUBLIC :: dipole_velocity_deriv
111
112CONTAINS
113
114! *****************************************************************************
115!> \brief This returns the derivative of the moment integrals [a|\mu|b], with respect
116!> to the basis function on the right
117!> difdip(beta, alpha) = < mu | r_beta | ∂_alpha nu > * (mu - nu)
118!> \param qs_env ...
119!> \param difdip ...
120!> \param order The order of the derivative (1 for dipole moment)
121!> \param lambda The atom on which we take the derivative
122!> \param rc ...
123!> \author Edward Ditler
124! **************************************************************************************************
125 SUBROUTINE dipole_velocity_deriv(qs_env, difdip, order, lambda, rc)
126 TYPE(qs_environment_type), POINTER :: qs_env
127 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(INOUT) :: difdip
128 INTEGER, INTENT(IN) :: order, lambda
129 REAL(kind=dp), DIMENSION(3) :: rc
130
131 CHARACTER(LEN=*), PARAMETER :: routinen = 'dipole_velocity_deriv'
132
133 INTEGER :: handle, i, iatom, icol, idir, ikind, inode, irow, iset, j, jatom, jkind, jset, &
134 last_jatom, lda, ldab, ldb, m_dim, maxsgf, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, &
135 sgfb
136 LOGICAL :: found
137 REAL(dp) :: dab
138 REAL(dp), DIMENSION(3) :: ra, rab, rac, rb, rbc
139 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
140 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
141 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab, difmab2
142 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mint, mint2
143 TYPE(cell_type), POINTER :: cell
144 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
145 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
146 TYPE(neighbor_list_iterator_p_type), &
147 DIMENSION(:), POINTER :: nl_iterator
148 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
149 POINTER :: sab_all
150 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
151 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
152 TYPE(qs_kind_type), POINTER :: qs_kind
153
154 CALL timeset(routinen, handle)
155
156 NULLIFY (cell, particle_set, qs_kind_set, sab_all)
157 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
158 qs_kind_set=qs_kind_set, sab_all=sab_all)
159 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
160 maxco=ldab, maxsgf=maxsgf)
161
162 nkind = SIZE(qs_kind_set)
163 natom = SIZE(particle_set)
164
165 m_dim = ncoset(order) - 1
166
167 ALLOCATE (basis_set_list(nkind))
168
169 ALLOCATE (mab(ldab, ldab, m_dim))
170 ALLOCATE (difmab2(ldab, ldab, m_dim, 3))
171 ALLOCATE (work(ldab, maxsgf))
172 ALLOCATE (mint(3, 3))
173 ALLOCATE (mint2(3, 3))
174
175 mab(1:ldab, 1:ldab, 1:m_dim) = 0.0_dp
176 difmab2(1:ldab, 1:ldab, 1:m_dim, 1:3) = 0.0_dp
177 work(1:ldab, 1:maxsgf) = 0.0_dp
178
179 DO i = 1, 3
180 DO j = 1, 3
181 NULLIFY (mint(i, j)%block)
182 NULLIFY (mint2(i, j)%block)
183 END DO
184 END DO
185
186 ! Set the basis_set_list(nkind) to point to the corresponding basis sets
187 DO ikind = 1, nkind
188 qs_kind => qs_kind_set(ikind)
189 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
190 IF (ASSOCIATED(basis_set_a)) THEN
191 basis_set_list(ikind)%gto_basis_set => basis_set_a
192 ELSE
193 NULLIFY (basis_set_list(ikind)%gto_basis_set)
194 END IF
195 END DO
196
197 CALL neighbor_list_iterator_create(nl_iterator, sab_all)
198 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
199 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
200 iatom=iatom, jatom=jatom, r=rab)
201
202 basis_set_a => basis_set_list(ikind)%gto_basis_set
203 basis_set_b => basis_set_list(jkind)%gto_basis_set
204 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
205 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
206
207 associate( &
208 ! basis ikind
209 first_sgfa => basis_set_a%first_sgf, &
210 la_max => basis_set_a%lmax, &
211 la_min => basis_set_a%lmin, &
212 npgfa => basis_set_a%npgf, &
213 nsgfa => basis_set_a%nsgf_set, &
214 rpgfa => basis_set_a%pgf_radius, &
215 set_radius_a => basis_set_a%set_radius, &
216 sphi_a => basis_set_a%sphi, &
217 zeta => basis_set_a%zet, &
218 ! basis jkind, &
219 first_sgfb => basis_set_b%first_sgf, &
220 lb_max => basis_set_b%lmax, &
221 lb_min => basis_set_b%lmin, &
222 npgfb => basis_set_b%npgf, &
223 nsgfb => basis_set_b%nsgf_set, &
224 rpgfb => basis_set_b%pgf_radius, &
225 set_radius_b => basis_set_b%set_radius, &
226 sphi_b => basis_set_b%sphi, &
227 zetb => basis_set_b%zet)
228
229 nseta = basis_set_a%nset
230 nsetb = basis_set_b%nset
231
232 IF (inode == 1) last_jatom = 0
233
234 ! this guarantees minimum image convention
235 ! anything else would not make sense
236 IF (jatom == last_jatom) THEN
237 cycle
238 END IF
239
240 last_jatom = jatom
241
242 irow = iatom
243 icol = jatom
244
245 DO i = 1, 3
246 DO j = 1, 3
247 NULLIFY (mint(i, j)%block)
248 CALL dbcsr_get_block_p(matrix=difdip(i, j)%matrix, &
249 row=irow, col=icol, block=mint(i, j)%block, &
250 found=found)
251 cpassert(found)
252 mint(i, j)%block = 0._dp
253 END DO
254 END DO
255
256 ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
257 ra = pbc(particle_set(iatom)%r(:), cell)
258 rb(:) = ra(:) + rab(:)
259 rac = pbc(rc, ra, cell)
260 rbc = rac + rab
261 dab = norm2(rab)
262
263 DO iset = 1, nseta
264 ncoa = npgfa(iset)*ncoset(la_max(iset))
265 sgfa = first_sgfa(1, iset)
266 DO jset = 1, nsetb
267 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
268 ncob = npgfb(jset)*ncoset(lb_max(jset))
269 sgfb = first_sgfb(1, jset)
270 ldab = max(ncoa, ncob)
271 lda = ncoset(la_max(iset))*npgfa(iset)
272 ldb = ncoset(lb_max(jset))*npgfb(jset)
273 ALLOCATE (difmab(lda, ldb, m_dim, 3))
274
275 ! Calculate integral difmab(beta, alpha) = (a|r_beta|db_alpha)
276 ! difmab(beta, alpha) = < a | r_beta | ∂_alpha b >
277 ! difmab(j, idir) = < a | r_j | ∂_idir b >
278 CALL diff_momop_velocity(la_max(iset), npgfa(iset), zeta(:, iset), &
279 rpgfa(:, iset), la_min(iset), lb_max(jset), npgfb(jset), &
280 zetb(:, jset), rpgfb(:, jset), lb_min(jset), order, rac, rbc, &
281 difmab, lambda=lambda, iatom=iatom, jatom=jatom)
282
283 ! *** Contraction step ***
284
285 DO idir = 1, 3 ! derivative of AO function
286 DO j = 1, 3 ! position operator r_j
287 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
288 1.0_dp, difmab(1, 1, j, idir), SIZE(difmab, 1), &
289 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
290 0.0_dp, work(1, 1), SIZE(work, 1))
291
292 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
293 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
294 work(1, 1), SIZE(work, 1), &
295 1.0_dp, mint(j, idir)%block(sgfa, sgfb), &
296 SIZE(mint(j, idir)%block, 1))
297 END DO !j
298 END DO !idir
299 DEALLOCATE (difmab)
300 END DO !jset
301 END DO !iset
302 END associate
303 END do!iterator
304
305 CALL neighbor_list_iterator_release(nl_iterator)
306
307 DO i = 1, 3
308 DO j = 1, 3
309 NULLIFY (mint(i, j)%block)
310 END DO
311 END DO
312
313 DEALLOCATE (mab, difmab2, basis_set_list, work, mint, mint2)
314
315 CALL timestop(handle)
316 END SUBROUTINE dipole_velocity_deriv
317
318! **************************************************************************************************
319!> \brief Builds the moments for the derivative of the overlap with respect to nuclear velocities
320!> \param qs_env ...
321!> \param moments ...
322!> \param nmoments ...
323!> \param ref_point ...
324!> \param ref_points ...
325!> \param basis_type ...
326!> \author Edward Ditler
327! **************************************************************************************************
328 SUBROUTINE build_dsdv_moments(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
329
330 TYPE(qs_environment_type), POINTER :: qs_env
331 TYPE(dbcsr_p_type), DIMENSION(:) :: moments
332 INTEGER, INTENT(IN) :: nmoments
333 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
334 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
335 OPTIONAL :: ref_points
336 CHARACTER(len=*), OPTIONAL :: basis_type
337
338 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_dsdv_moments'
339
340 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
341 maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
342 INTEGER, DIMENSION(3) :: image_cell
343 LOGICAL :: found
344 REAL(kind=dp) :: dab
345 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
346 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
347 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
348 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
349 TYPE(cell_type), POINTER :: cell
350 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
351 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
352 TYPE(neighbor_list_iterator_p_type), &
353 DIMENSION(:), POINTER :: nl_iterator
354 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
355 POINTER :: sab_orb
356 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
357 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
358 TYPE(qs_kind_type), POINTER :: qs_kind
359
360 IF (nmoments < 1) RETURN
361
362 CALL timeset(routinen, handle)
363
364 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
365
366 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
367 cpassert(SIZE(moments) >= nm)
368
369 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
370 CALL get_qs_env(qs_env=qs_env, &
371 qs_kind_set=qs_kind_set, &
372 particle_set=particle_set, cell=cell, &
373 sab_orb=sab_orb)
374
375 nkind = SIZE(qs_kind_set)
376 natom = SIZE(particle_set)
377
378 ! Allocate work storage
379 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
380 maxco=maxco, maxsgf=maxsgf, &
381 basis_type=basis_type)
382
383 ALLOCATE (mab(maxco, maxco, nm))
384 mab(:, :, :) = 0.0_dp
385
386 ALLOCATE (work(maxco, maxsgf))
387 work(:, :) = 0.0_dp
388
389 ALLOCATE (mint(nm))
390 DO i = 1, nm
391 NULLIFY (mint(i)%block)
392 END DO
393
394 ALLOCATE (basis_set_list(nkind))
395 DO ikind = 1, nkind
396 qs_kind => qs_kind_set(ikind)
397 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
398 IF (ASSOCIATED(basis_set_a)) THEN
399 basis_set_list(ikind)%gto_basis_set => basis_set_a
400 ELSE
401 NULLIFY (basis_set_list(ikind)%gto_basis_set)
402 END IF
403 END DO
404 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
405 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
406 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
407 iatom=iatom, jatom=jatom, r=rab, cell=image_cell)
408 basis_set_a => basis_set_list(ikind)%gto_basis_set
409 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
410 basis_set_b => basis_set_list(jkind)%gto_basis_set
411 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
412 ! basis ikind
413 associate( &
414 first_sgfa => basis_set_a%first_sgf, &
415 la_max => basis_set_a%lmax, &
416 la_min => basis_set_a%lmin, &
417 npgfa => basis_set_a%npgf, &
418 nsgfa => basis_set_a%nsgf_set, &
419 rpgfa => basis_set_a%pgf_radius, &
420 set_radius_a => basis_set_a%set_radius, &
421 sphi_a => basis_set_a%sphi, &
422 zeta => basis_set_a%zet, &
423 ! basis jkind, &
424 first_sgfb => basis_set_b%first_sgf, &
425 lb_max => basis_set_b%lmax, &
426 lb_min => basis_set_b%lmin, &
427 npgfb => basis_set_b%npgf, &
428 nsgfb => basis_set_b%nsgf_set, &
429 rpgfb => basis_set_b%pgf_radius, &
430 set_radius_b => basis_set_b%set_radius, &
431 sphi_b => basis_set_b%sphi, &
432 zetb => basis_set_b%zet)
433
434 nseta = basis_set_a%nset
435 nsetb = basis_set_b%nset
436
437 IF (inode == 1) last_jatom = 0
438
439 ! this guarantees minimum image convention
440 ! anything else would not make sense
441 IF (jatom == last_jatom) THEN
442 cycle
443 END IF
444
445 last_jatom = jatom
446
447 IF (iatom <= jatom) THEN
448 irow = iatom
449 icol = jatom
450 ELSE
451 irow = jatom
452 icol = iatom
453 END IF
454
455 DO i = 1, nm
456 NULLIFY (mint(i)%block)
457 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
458 row=irow, col=icol, block=mint(i)%block, found=found)
459 mint(i)%block = 0._dp
460 END DO
461
462 ! fold atomic position back into unit cell
463 IF (PRESENT(ref_points)) THEN
464 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
465 ELSE IF (PRESENT(ref_point)) THEN
466 rc(:) = ref_point(:)
467 ELSE
468 rc(:) = 0._dp
469 END IF
470 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
471 ! by folding around the center, such screwing can be avoided for a proper choice of center.
472 ! we dont use PBC at this point
473
474 ! With and without MIC, rab(:) is the vector giving us the coordinates of rb.
475 ra(:) = particle_set(iatom)%r(:)
476 rb(:) = ra(:) + rab(:)
477 rac = pbc(rc, ra, cell)
478 rbc = rac + rab
479
480 dab = norm2(rab)
481
482 DO iset = 1, nseta
483
484 ncoa = npgfa(iset)*ncoset(la_max(iset))
485 sgfa = first_sgfa(1, iset)
486
487 DO jset = 1, nsetb
488
489 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
490
491 ncob = npgfb(jset)*ncoset(lb_max(jset))
492 sgfb = first_sgfb(1, jset)
493
494 ! Calculate the primitive integrals
495 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
496 rpgfa(:, iset), la_min(iset), &
497 lb_max(jset), npgfb(jset), zetb(:, jset), &
498 rpgfb(:, jset), nmoments, rac, rbc, mab)
499
500 ! Contraction step
501 DO i = 1, nm
502
503 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
504 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
505 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
506 0.0_dp, work(1, 1), SIZE(work, 1))
507
508 IF (iatom <= jatom) THEN
509
510 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
511 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
512 work(1, 1), SIZE(work, 1), &
513 1.0_dp, mint(i)%block(sgfa, sgfb), &
514 SIZE(mint(i)%block, 1))
515
516 ELSE
517
518 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
519 1.0_dp, work(1, 1), SIZE(work, 1), &
520 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
521 1.0_dp, mint(i)%block(sgfb, sgfa), &
522 SIZE(mint(i)%block, 1))
523
524 END IF
525
526 END DO
527
528 END DO
529 END DO
530 END associate
531
532 END DO ! iterator
533
534 CALL neighbor_list_iterator_release(nl_iterator)
535
536 ! Release work storage
537 DEALLOCATE (mab, basis_set_list)
538 DEALLOCATE (work)
539 DO i = 1, nm
540 NULLIFY (mint(i)%block)
541 END DO
542 DEALLOCATE (mint)
543
544 CALL timestop(handle)
545
546 END SUBROUTINE build_dsdv_moments
547
548! **************************************************************************************************
549!> \brief ...
550!> \param qs_env ...
551!> \param moments ...
552!> \param nmoments ...
553!> \param ref_point ...
554!> \param ref_points ...
555!> \param basis_type ...
556! **************************************************************************************************
557 SUBROUTINE build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
558
559 TYPE(qs_environment_type), POINTER :: qs_env
560 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: moments
561 INTEGER, INTENT(IN) :: nmoments
562 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
563 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
564 OPTIONAL :: ref_points
565 CHARACTER(len=*), OPTIONAL :: basis_type
566
567 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moment_matrix'
568
569 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, last_jatom, &
570 maxco, maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
571 LOGICAL :: found
572 REAL(kind=dp) :: dab
573 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
574 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
575 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
576 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
577 TYPE(cell_type), POINTER :: cell
578 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
579 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
580 TYPE(neighbor_list_iterator_p_type), &
581 DIMENSION(:), POINTER :: nl_iterator
582 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
583 POINTER :: sab_orb
584 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
585 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
586 TYPE(qs_kind_type), POINTER :: qs_kind
587
588 IF (nmoments < 1) RETURN
589
590 CALL timeset(routinen, handle)
591
592 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
593
594 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
595 cpassert(SIZE(moments) >= nm)
596
597 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
598 CALL get_qs_env(qs_env=qs_env, &
599 qs_kind_set=qs_kind_set, &
600 particle_set=particle_set, cell=cell, &
601 sab_orb=sab_orb)
602
603 nkind = SIZE(qs_kind_set)
604 natom = SIZE(particle_set)
605
606 ! Allocate work storage
607 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
608 maxco=maxco, maxsgf=maxsgf, &
609 basis_type=basis_type)
610
611 ALLOCATE (mab(maxco, maxco, nm))
612 mab(:, :, :) = 0.0_dp
613
614 ALLOCATE (work(maxco, maxsgf))
615 work(:, :) = 0.0_dp
616
617 ALLOCATE (mint(nm))
618 DO i = 1, nm
619 NULLIFY (mint(i)%block)
620 END DO
621
622 ALLOCATE (basis_set_list(nkind))
623 DO ikind = 1, nkind
624 qs_kind => qs_kind_set(ikind)
625 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
626 IF (ASSOCIATED(basis_set_a)) THEN
627 basis_set_list(ikind)%gto_basis_set => basis_set_a
628 ELSE
629 NULLIFY (basis_set_list(ikind)%gto_basis_set)
630 END IF
631 END DO
632 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
633 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
634 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
635 iatom=iatom, jatom=jatom, r=rab)
636 basis_set_a => basis_set_list(ikind)%gto_basis_set
637 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
638 basis_set_b => basis_set_list(jkind)%gto_basis_set
639 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
640 associate( &
641 ! basis ikind
642 first_sgfa => basis_set_a%first_sgf, &
643 la_max => basis_set_a%lmax, &
644 la_min => basis_set_a%lmin, &
645 npgfa => basis_set_a%npgf, &
646 nsgfa => basis_set_a%nsgf_set, &
647 rpgfa => basis_set_a%pgf_radius, &
648 set_radius_a => basis_set_a%set_radius, &
649 sphi_a => basis_set_a%sphi, &
650 zeta => basis_set_a%zet, &
651 ! basis jkind, &
652 first_sgfb => basis_set_b%first_sgf, &
653 lb_max => basis_set_b%lmax, &
654 lb_min => basis_set_b%lmin, &
655 npgfb => basis_set_b%npgf, &
656 nsgfb => basis_set_b%nsgf_set, &
657 rpgfb => basis_set_b%pgf_radius, &
658 set_radius_b => basis_set_b%set_radius, &
659 sphi_b => basis_set_b%sphi, &
660 zetb => basis_set_b%zet)
661
662 nseta = basis_set_a%nset
663 nsetb = basis_set_b%nset
664
665 IF (inode == 1) last_jatom = 0
666
667 ! this guarantees minimum image convention
668 ! anything else would not make sense
669 IF (jatom == last_jatom) THEN
670 cycle
671 END IF
672
673 last_jatom = jatom
674
675 IF (iatom <= jatom) THEN
676 irow = iatom
677 icol = jatom
678 ELSE
679 irow = jatom
680 icol = iatom
681 END IF
682
683 DO i = 1, nm
684 NULLIFY (mint(i)%block)
685 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
686 row=irow, col=icol, block=mint(i)%block, found=found)
687 mint(i)%block = 0._dp
688 END DO
689
690 ! fold atomic position back into unit cell
691 IF (PRESENT(ref_points)) THEN
692 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
693 ELSE IF (PRESENT(ref_point)) THEN
694 rc(:) = ref_point(:)
695 ELSE
696 rc(:) = 0._dp
697 END IF
698 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
699 ! by folding around the center, such screwing can be avoided for a proper choice of center.
700 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
701 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
702 ! we dont use PBC at this point
703 rab(:) = ra(:) - rb(:)
704 rac(:) = ra(:) - rc(:)
705 rbc(:) = rb(:) - rc(:)
706 dab = norm2(rab)
707
708 DO iset = 1, nseta
709
710 ncoa = npgfa(iset)*ncoset(la_max(iset))
711 sgfa = first_sgfa(1, iset)
712
713 DO jset = 1, nsetb
714
715 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
716
717 ncob = npgfb(jset)*ncoset(lb_max(jset))
718 sgfb = first_sgfb(1, jset)
719
720 ! Calculate the primitive integrals
721 CALL moment(la_max(iset), npgfa(iset), zeta(:, iset), &
722 rpgfa(:, iset), la_min(iset), &
723 lb_max(jset), npgfb(jset), zetb(:, jset), &
724 rpgfb(:, jset), nmoments, rac, rbc, mab)
725
726 ! Contraction step
727 DO i = 1, nm
728
729 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
730 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
731 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
732 0.0_dp, work(1, 1), SIZE(work, 1))
733
734 IF (iatom <= jatom) THEN
735
736 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
737 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
738 work(1, 1), SIZE(work, 1), &
739 1.0_dp, mint(i)%block(sgfa, sgfb), &
740 SIZE(mint(i)%block, 1))
741
742 ELSE
743
744 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
745 1.0_dp, work(1, 1), SIZE(work, 1), &
746 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
747 1.0_dp, mint(i)%block(sgfb, sgfa), &
748 SIZE(mint(i)%block, 1))
749
750 END IF
751
752 END DO
753
754 END DO
755 END DO
756 END associate
757
758 END DO
759 CALL neighbor_list_iterator_release(nl_iterator)
760
761 ! Release work storage
762 DEALLOCATE (mab, basis_set_list)
763 DEALLOCATE (work)
764 DO i = 1, nm
765 NULLIFY (mint(i)%block)
766 END DO
767 DEALLOCATE (mint)
768
769 CALL timestop(handle)
770
771 END SUBROUTINE build_local_moment_matrix
772
773! **************************************************************************************************
774!> \brief Calculate right-hand sided derivatives of multipole moments, e. g. < a | xy d/dz | b >
775!> Optionally stores the multipole moments themselves for free.
776!> Note that the multipole moments are symmetric while their derivatives are anti-symmetric
777!> Only first derivatives are performed, e. g. x d/dy
778!> \param qs_env ...
779!> \param moments_der will contain the derivatives of the multipole moments
780!> \param nmoments_der order of the moments with derivatives
781!> \param nmoments order of the multipole moments (no derivatives, same output as
782!> build_local_moment_matrix, needs moments as arguments to store results)
783!> \param ref_point ...
784!> \param moments contains the multipole moments, optionally for free, up to order nmoments
785!> \note
786!> Adapted from rRc_xyz_der_ao in qs_operators_ao
787! **************************************************************************************************
788 SUBROUTINE build_local_moments_der_matrix(qs_env, moments_der, nmoments_der, nmoments, &
789 ref_point, moments)
790 TYPE(qs_environment_type), POINTER :: qs_env
791 TYPE(dbcsr_p_type), DIMENSION(:, :), &
792 INTENT(INOUT), POINTER :: moments_der
793 INTEGER, INTENT(IN) :: nmoments_der, nmoments
794 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
795 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
796 OPTIONAL, POINTER :: moments
797
798 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_moments_der_matrix'
799
800 INTEGER :: dimders, handle, i, iatom, icol, ider, ii, ikind, inode, ipgf, irow, iset, j, &
801 jatom, jkind, jpgf, jset, last_jatom, m_dim, maxco, maxsgf, na, nb, ncoa, ncob, nda, ndb, &
802 nders, nkind, nm, nmom_build, nseta, nsetb, sgfa, sgfb
803 LOGICAL :: found
804 REAL(kind=dp) :: dab
805 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
806 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
807 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: difmab
808 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
809 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: mab_tmp
810 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mom_block
811 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:, :) :: mom_block_der
812 TYPE(cell_type), POINTER :: cell
813 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
814 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
815 TYPE(neighbor_list_iterator_p_type), &
816 DIMENSION(:), POINTER :: nl_iterator
817 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
818 POINTER :: sab_orb
819 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
820 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
821 TYPE(qs_kind_type), POINTER :: qs_kind
822
823 nmom_build = max(nmoments, nmoments_der) ! build moments up to order nmom_buiod
824 IF (nmom_build < 1) RETURN
825
826 CALL timeset(routinen, handle)
827
828 nders = 1 ! only first order derivatives
829 dimders = ncoset(nders) - 1
830
831 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
832 CALL get_qs_env(qs_env=qs_env, &
833 qs_kind_set=qs_kind_set, &
834 particle_set=particle_set, &
835 cell=cell, &
836 sab_orb=sab_orb)
837
838 nkind = SIZE(qs_kind_set)
839
840 ! Work storage
841 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
842 maxco=maxco, maxsgf=maxsgf)
843
844 IF (nmoments > 0) THEN
845 cpassert(PRESENT(moments))
846 nm = (6 + 11*nmoments + 6*nmoments**2 + nmoments**3)/6 - 1
847 cpassert(SIZE(moments) == nm)
848 ! storage for integrals
849 ALLOCATE (mab(maxco, maxco, nm))
850 ! blocks
851 mab(:, :, :) = 0.0_dp
852 ALLOCATE (mom_block(nm))
853 DO i = 1, nm
854 NULLIFY (mom_block(i)%block)
855 END DO
856 END IF
857
858 IF (nmoments_der > 0) THEN
859 m_dim = ncoset(nmoments_der) - 1
860 cpassert(SIZE(moments_der, dim=1) == m_dim)
861 cpassert(SIZE(moments_der, dim=2) == dimders)
862 ! storage for integrals
863 ALLOCATE (difmab(maxco, maxco, m_dim, dimders))
864 difmab(:, :, :, :) = 0.0_dp
865 ! blocks
866 ALLOCATE (mom_block_der(m_dim, dimders))
867 DO i = 1, m_dim
868 DO ider = 1, dimders
869 NULLIFY (mom_block_der(i, ider)%block)
870 END DO
871 END DO
872 END IF
873
874 ALLOCATE (work(maxco, maxsgf))
875 work(:, :) = 0.0_dp
876
877 NULLIFY (basis_set_a, basis_set_b, basis_set_list)
878 NULLIFY (qs_kind)
879 ALLOCATE (basis_set_list(nkind))
880 DO ikind = 1, nkind
881 qs_kind => qs_kind_set(ikind)
882 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a)
883 IF (ASSOCIATED(basis_set_a)) THEN
884 basis_set_list(ikind)%gto_basis_set => basis_set_a
885 ELSE
886 NULLIFY (basis_set_list(ikind)%gto_basis_set)
887 END IF
888 END DO
889
890 ! Calculate derivatives looping over neighbour list
891 NULLIFY (nl_iterator)
892 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
893 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
894 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
895 iatom=iatom, jatom=jatom, r=rab)
896 basis_set_a => basis_set_list(ikind)%gto_basis_set
897 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
898 basis_set_b => basis_set_list(jkind)%gto_basis_set
899 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
900 associate( &
901 ! basis ikind
902 first_sgfa => basis_set_a%first_sgf, &
903 la_max => basis_set_a%lmax, &
904 la_min => basis_set_a%lmin, &
905 npgfa => basis_set_a%npgf, &
906 nsgfa => basis_set_a%nsgf_set, &
907 rpgfa => basis_set_a%pgf_radius, &
908 set_radius_a => basis_set_a%set_radius, &
909 sphi_a => basis_set_a%sphi, &
910 zeta => basis_set_a%zet, &
911 ! basis jkind, &
912 first_sgfb => basis_set_b%first_sgf, &
913 lb_max => basis_set_b%lmax, &
914 lb_min => basis_set_b%lmin, &
915 npgfb => basis_set_b%npgf, &
916 nsgfb => basis_set_b%nsgf_set, &
917 rpgfb => basis_set_b%pgf_radius, &
918 set_radius_b => basis_set_b%set_radius, &
919 sphi_b => basis_set_b%sphi, &
920 zetb => basis_set_b%zet)
921
922 nseta = basis_set_a%nset
923 nsetb = basis_set_b%nset
924
925 ! reference point
926 IF (PRESENT(ref_point)) THEN
927 rc(:) = ref_point(:)
928 ELSE
929 rc(:) = 0._dp
930 END IF
931 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
932 ! by folding around the center, such screwing can be avoided for a proper choice of center.
933 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
934 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
935 ! we dont use PBC at this point
936 rab(:) = ra(:) - rb(:)
937 rac(:) = ra(:) - rc(:)
938 rbc(:) = rb(:) - rc(:)
939 dab = norm2(rab)
940
941 ! get blocks
942 IF (inode == 1) last_jatom = 0
943
944 IF (jatom == last_jatom) THEN
945 cycle
946 END IF
947
948 last_jatom = jatom
949
950 IF (iatom <= jatom) THEN
951 irow = iatom
952 icol = jatom
953 ELSE
954 irow = jatom
955 icol = iatom
956 END IF
957
958 IF (nmoments > 0) THEN
959 DO i = 1, nm
960 NULLIFY (mom_block(i)%block)
961 ! get block from pre calculated overlap matrix
962 CALL dbcsr_get_block_p(matrix=moments(i)%matrix, &
963 row=irow, col=icol, block=mom_block(i)%block, found=found)
964 cpassert(found .AND. ASSOCIATED(mom_block(i)%block))
965 mom_block(i)%block = 0._dp
966 END DO
967 END IF
968 IF (nmoments_der > 0) THEN
969 DO i = 1, m_dim
970 DO ider = 1, dimders
971 NULLIFY (mom_block_der(i, ider)%block)
972 CALL dbcsr_get_block_p(matrix=moments_der(i, ider)%matrix, &
973 row=irow, col=icol, &
974 block=mom_block_der(i, ider)%block, &
975 found=found)
976 cpassert(found .AND. ASSOCIATED(mom_block_der(i, ider)%block))
977 mom_block_der(i, ider)%block = 0._dp
978 END DO
979 END DO
980 END IF
981
982 DO iset = 1, nseta
983
984 ncoa = npgfa(iset)*ncoset(la_max(iset))
985 sgfa = first_sgfa(1, iset)
986
987 DO jset = 1, nsetb
988
989 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
990
991 ncob = npgfb(jset)*ncoset(lb_max(jset))
992 sgfb = first_sgfb(1, jset)
993
994 NULLIFY (mab_tmp)
995 ALLOCATE (mab_tmp(npgfa(iset)*ncoset(la_max(iset) + 1), &
996 npgfb(jset)*ncoset(lb_max(jset) + 1), ncoset(nmom_build) - 1))
997
998 ! Calculate the primitive integrals (need l+1 for derivatives)
999 CALL moment(la_max(iset) + 1, npgfa(iset), zeta(:, iset), &
1000 rpgfa(:, iset), la_min(iset), &
1001 lb_max(jset) + 1, npgfb(jset), zetb(:, jset), &
1002 rpgfb(:, jset), nmom_build, rac, rbc, mab_tmp)
1003
1004 IF (nmoments_der > 0) THEN
1005 CALL diff_momop(la_max(iset), npgfa(iset), zeta(:, iset), &
1006 rpgfa(:, iset), la_min(iset), &
1007 lb_max(jset), npgfb(jset), zetb(:, jset), &
1008 rpgfb(:, jset), lb_min(jset), &
1009 nmoments_der, rac, rbc, difmab, mab_ext=mab_tmp)
1010 END IF
1011
1012 IF (nmoments > 0) THEN
1013 ! copy subset of mab_tmp (l+1) to mab (l)
1014 mab = 0.0_dp
1015 DO ii = 1, nm
1016 na = 0
1017 nda = 0
1018 DO ipgf = 1, npgfa(iset)
1019 nb = 0
1020 ndb = 0
1021 DO jpgf = 1, npgfb(jset)
1022 DO j = 1, ncoset(lb_max(jset))
1023 DO i = 1, ncoset(la_max(iset))
1024 mab(i + na, j + nb, ii) = mab_tmp(i + nda, j + ndb, ii)
1025 END DO ! i
1026 END DO ! j
1027 nb = nb + ncoset(lb_max(jset))
1028 ndb = ndb + ncoset(lb_max(jset) + 1)
1029 END DO ! jpgf
1030 na = na + ncoset(la_max(iset))
1031 nda = nda + ncoset(la_max(iset) + 1)
1032 END DO ! ipgf
1033 END DO
1034 ! Contraction step
1035 DO i = 1, nm
1036
1037 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1038 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1039 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1040 0.0_dp, work(1, 1), SIZE(work, 1))
1041
1042 IF (iatom <= jatom) THEN
1043 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1044 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1045 work(1, 1), SIZE(work, 1), &
1046 1.0_dp, mom_block(i)%block(sgfa, sgfb), &
1047 SIZE(mom_block(i)%block, 1))
1048 ELSE
1049 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1050 1.0_dp, work(1, 1), SIZE(work, 1), &
1051 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1052 1.0_dp, mom_block(i)%block(sgfb, sgfa), &
1053 SIZE(mom_block(i)%block, 1))
1054 END IF
1055 END DO
1056 END IF
1057
1058 IF (nmoments_der > 0) THEN
1059 DO i = 1, m_dim
1060 DO ider = 1, dimders
1061 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1062 1.0_dp, difmab(1, 1, i, ider), SIZE(difmab, 1), &
1063 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1064 0._dp, work(1, 1), SIZE(work, 1))
1065
1066 IF (iatom <= jatom) THEN
1067 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1068 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1069 work(1, 1), SIZE(work, 1), &
1070 1._dp, mom_block_der(i, ider)%block(sgfa, sgfb), &
1071 SIZE(mom_block_der(i, ider)%block, 1))
1072 ELSE
1073 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1074 -1.0_dp, work(1, 1), SIZE(work, 1), &
1075 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1076 1.0_dp, mom_block_der(i, ider)%block(sgfb, sgfa), &
1077 SIZE(mom_block_der(i, ider)%block, 1))
1078 END IF
1079 END DO
1080 END DO
1081 END IF
1082 DEALLOCATE (mab_tmp)
1083 END DO
1084 END DO
1085 END associate
1086 END DO
1087 CALL neighbor_list_iterator_release(nl_iterator)
1088
1089 ! deallocations
1090 DEALLOCATE (basis_set_list)
1091 DEALLOCATE (work)
1092 IF (nmoments > 0) THEN
1093 DEALLOCATE (mab)
1094 DO i = 1, nm
1095 NULLIFY (mom_block(i)%block)
1096 END DO
1097 DEALLOCATE (mom_block)
1098 END IF
1099 IF (nmoments_der > 0) THEN
1100 DEALLOCATE (difmab)
1101 DO i = 1, m_dim
1102 DO ider = 1, dimders
1103 NULLIFY (mom_block_der(i, ider)%block)
1104 END DO
1105 END DO
1106 DEALLOCATE (mom_block_der)
1107 END IF
1108
1109 CALL timestop(handle)
1110
1111 END SUBROUTINE build_local_moments_der_matrix
1112
1113! **************************************************************************************************
1114!> \brief ...
1115!> \param qs_env ...
1116!> \param magmom ...
1117!> \param nmoments ...
1118!> \param ref_point ...
1119!> \param ref_points ...
1120!> \param basis_type ...
1121! **************************************************************************************************
1122 SUBROUTINE build_local_magmom_matrix(qs_env, magmom, nmoments, ref_point, ref_points, basis_type)
1123
1124 TYPE(qs_environment_type), POINTER :: qs_env
1125 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom
1126 INTEGER, INTENT(IN) :: nmoments
1127 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: ref_point
1128 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
1129 OPTIONAL :: ref_points
1130 CHARACTER(len=*), OPTIONAL :: basis_type
1131
1132 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_local_magmom_matrix'
1133
1134 INTEGER :: handle, i, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, maxco, &
1135 maxsgf, natom, ncoa, ncob, nkind, nm, nseta, nsetb, sgfa, sgfb
1136 LOGICAL :: found
1137 REAL(kind=dp) :: dab
1138 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
1139 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: mab
1140 REAL(kind=dp), DIMENSION(3) :: ra, rab, rac, rb, rbc, rc
1141 TYPE(block_p_type), ALLOCATABLE, DIMENSION(:) :: mint
1142 TYPE(cell_type), POINTER :: cell
1143 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1144 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1145 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1146 TYPE(neighbor_list_iterator_p_type), &
1147 DIMENSION(:), POINTER :: nl_iterator
1148 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1149 POINTER :: sab_orb
1150 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1151 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1152 TYPE(qs_kind_type), POINTER :: qs_kind
1153
1154 IF (nmoments < 1) RETURN
1155
1156 CALL timeset(routinen, handle)
1157
1158 NULLIFY (qs_kind_set, cell, particle_set, sab_orb)
1159 NULLIFY (matrix_s)
1160
1161 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
1162
1163 ! magnetic dipoles/angular moments only
1164 nm = 3
1165
1166 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1167 CALL get_qs_env(qs_env=qs_env, &
1168 qs_kind_set=qs_kind_set, &
1169 particle_set=particle_set, cell=cell, &
1170 sab_orb=sab_orb)
1171
1172 nkind = SIZE(qs_kind_set)
1173 natom = SIZE(particle_set)
1174
1175 ! Allocate work storage
1176 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
1177 maxco=maxco, maxsgf=maxsgf)
1178
1179 ALLOCATE (mab(maxco, maxco, nm))
1180 mab(:, :, :) = 0.0_dp
1181
1182 ALLOCATE (work(maxco, maxsgf))
1183 work(:, :) = 0.0_dp
1184
1185 ALLOCATE (mint(nm))
1186 DO i = 1, nm
1187 NULLIFY (mint(i)%block)
1188 END DO
1189
1190 ALLOCATE (basis_set_list(nkind))
1191 DO ikind = 1, nkind
1192 qs_kind => qs_kind_set(ikind)
1193 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1194 IF (ASSOCIATED(basis_set_a)) THEN
1195 basis_set_list(ikind)%gto_basis_set => basis_set_a
1196 ELSE
1197 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1198 END IF
1199 END DO
1200 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1201 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1202 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1203 iatom=iatom, jatom=jatom, r=rab)
1204 basis_set_a => basis_set_list(ikind)%gto_basis_set
1205 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1206 basis_set_b => basis_set_list(jkind)%gto_basis_set
1207 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1208 associate( &
1209 ! basis ikind
1210 first_sgfa => basis_set_a%first_sgf, &
1211 la_max => basis_set_a%lmax, &
1212 la_min => basis_set_a%lmin, &
1213 npgfa => basis_set_a%npgf, &
1214 nsgfa => basis_set_a%nsgf_set, &
1215 rpgfa => basis_set_a%pgf_radius, &
1216 set_radius_a => basis_set_a%set_radius, &
1217 sphi_a => basis_set_a%sphi, &
1218 zeta => basis_set_a%zet, &
1219 ! basis jkind, &
1220 first_sgfb => basis_set_b%first_sgf, &
1221 lb_max => basis_set_b%lmax, &
1222 lb_min => basis_set_b%lmin, &
1223 npgfb => basis_set_b%npgf, &
1224 nsgfb => basis_set_b%nsgf_set, &
1225 rpgfb => basis_set_b%pgf_radius, &
1226 set_radius_b => basis_set_b%set_radius, &
1227 sphi_b => basis_set_b%sphi, &
1228 zetb => basis_set_b%zet)
1229
1230 nseta = basis_set_a%nset
1231 nsetb = basis_set_b%nset
1232
1233 IF (iatom <= jatom) THEN
1234 irow = iatom
1235 icol = jatom
1236 ELSE
1237 irow = jatom
1238 icol = iatom
1239 END IF
1240
1241 DO i = 1, nm
1242 NULLIFY (mint(i)%block)
1243 CALL dbcsr_get_block_p(matrix=magmom(i)%matrix, &
1244 row=irow, col=icol, block=mint(i)%block, found=found)
1245 mint(i)%block = 0._dp
1246 cpassert(ASSOCIATED(mint(i)%block))
1247 END DO
1248
1249 ! fold atomic position back into unit cell
1250 IF (PRESENT(ref_points)) THEN
1251 rc(:) = 0.5_dp*(ref_points(:, iatom) + ref_points(:, jatom))
1252 ELSE IF (PRESENT(ref_point)) THEN
1253 rc(:) = ref_point(:)
1254 ELSE
1255 rc(:) = 0._dp
1256 END IF
1257 ! using PBC here might screw a molecule that fits the box (but e.g. hasn't been shifted by center_molecule)
1258 ! by folding around the center, such screwing can be avoided for a proper choice of center.
1259 ra(:) = pbc(particle_set(iatom)%r(:) - rc, cell) + rc
1260 rb(:) = pbc(particle_set(jatom)%r(:) - rc, cell) + rc
1261 ! we dont use PBC at this point
1262 rab(:) = ra(:) - rb(:)
1263 rac(:) = ra(:) - rc(:)
1264 rbc(:) = rb(:) - rc(:)
1265 dab = norm2(rab)
1266
1267 DO iset = 1, nseta
1268
1269 ncoa = npgfa(iset)*ncoset(la_max(iset))
1270 sgfa = first_sgfa(1, iset)
1271
1272 DO jset = 1, nsetb
1273
1274 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1275
1276 ncob = npgfb(jset)*ncoset(lb_max(jset))
1277 sgfb = first_sgfb(1, jset)
1278
1279 ! Calculate the primitive integrals
1280 CALL angmom(la_max(iset), npgfa(iset), zeta(:, iset), &
1281 rpgfa(:, iset), la_min(iset), &
1282 lb_max(jset), npgfb(jset), zetb(:, jset), &
1283 rpgfb(:, jset), rac, rbc, mab)
1284
1285 ! Contraction step
1286 DO i = 1, nm
1287 CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
1288 1.0_dp, mab(1, 1, i), SIZE(mab, 1), &
1289 sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1290 0.0_dp, work(1, 1), SIZE(work, 1))
1291
1292 IF (iatom <= jatom) THEN
1293 CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
1294 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1295 work(1, 1), SIZE(work, 1), &
1296 1.0_dp, mint(i)%block(sgfa, sgfb), &
1297 SIZE(mint(i)%block, 1))
1298 ELSE
1299 CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
1300 -1.0_dp, work(1, 1), SIZE(work, 1), &
1301 sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1302 1.0_dp, mint(i)%block(sgfb, sgfa), &
1303 SIZE(mint(i)%block, 1))
1304 END IF
1305
1306 END DO
1307
1308 END DO
1309 END DO
1310 END associate
1311 END DO
1312 CALL neighbor_list_iterator_release(nl_iterator)
1313
1314 ! Release work storage
1315 DEALLOCATE (mab, basis_set_list)
1316 DEALLOCATE (work)
1317 DO i = 1, nm
1318 NULLIFY (mint(i)%block)
1319 END DO
1320 DEALLOCATE (mint)
1321
1322 CALL timestop(handle)
1323
1324 END SUBROUTINE build_local_magmom_matrix
1325
1326! **************************************************************************************************
1327!> \brief ...
1328!> \param qs_env ...
1329!> \param cosmat ...
1330!> \param sinmat ...
1331!> \param kvec ...
1332!> \param sab_orb_external ...
1333!> \param basis_type ...
1334! **************************************************************************************************
1335 SUBROUTINE build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec, sab_orb_external, basis_type)
1336
1337 TYPE(qs_environment_type), POINTER :: qs_env
1338 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
1339 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1340 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1341 OPTIONAL, POINTER :: sab_orb_external
1342 CHARACTER(len=*), OPTIONAL :: basis_type
1343
1344 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_moment_matrix'
1345
1346 INTEGER :: handle, iatom, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, ldsa, &
1347 ldsb, ldwork, natom, ncoa, ncob, nkind, nseta, nsetb, sgfa, sgfb
1348 LOGICAL :: found
1349 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1350 REAL(kind=dp) :: dab
1351 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1352 TYPE(cell_type), POINTER :: cell
1353 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1354 TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
1355 TYPE(neighbor_list_iterator_p_type), &
1356 DIMENSION(:), POINTER :: nl_iterator
1357 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1358 POINTER :: sab_orb
1359 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1360 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1361 TYPE(qs_kind_type), POINTER :: qs_kind
1362
1363 CALL timeset(routinen, handle)
1364
1365 NULLIFY (qs_kind_set, particle_set, sab_orb, cell)
1366 CALL get_qs_env(qs_env=qs_env, &
1367 qs_kind_set=qs_kind_set, &
1368 particle_set=particle_set, cell=cell, &
1369 sab_orb=sab_orb)
1370
1371 IF (PRESENT(sab_orb_external)) sab_orb => sab_orb_external
1372
1373 CALL dbcsr_set(sinmat, 0.0_dp)
1374 CALL dbcsr_set(cosmat, 0.0_dp)
1375
1376 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, maxco=ldwork, basis_type=basis_type)
1377 ldab = ldwork
1378 ALLOCATE (cosab(ldab, ldab))
1379 ALLOCATE (sinab(ldab, ldab))
1380 ALLOCATE (work(ldwork, ldwork))
1381
1382 nkind = SIZE(qs_kind_set)
1383 natom = SIZE(particle_set)
1384
1385 ALLOCATE (basis_set_list(nkind))
1386 DO ikind = 1, nkind
1387 qs_kind => qs_kind_set(ikind)
1388 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
1389 IF (ASSOCIATED(basis_set_a)) THEN
1390 basis_set_list(ikind)%gto_basis_set => basis_set_a
1391 ELSE
1392 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1393 END IF
1394 END DO
1395 CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1396 DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1397 CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, inode=inode, &
1398 iatom=iatom, jatom=jatom, r=rab)
1399 basis_set_a => basis_set_list(ikind)%gto_basis_set
1400 IF (.NOT. ASSOCIATED(basis_set_a)) cycle
1401 basis_set_b => basis_set_list(jkind)%gto_basis_set
1402 IF (.NOT. ASSOCIATED(basis_set_b)) cycle
1403 associate( &
1404 ! basis ikind
1405 first_sgfa => basis_set_a%first_sgf, &
1406 la_max => basis_set_a%lmax, &
1407 la_min => basis_set_a%lmin, &
1408 npgfa => basis_set_a%npgf, &
1409 nsgfa => basis_set_a%nsgf_set, &
1410 rpgfa => basis_set_a%pgf_radius, &
1411 set_radius_a => basis_set_a%set_radius, &
1412 sphi_a => basis_set_a%sphi, &
1413 zeta => basis_set_a%zet, &
1414 ! basis jkind, &
1415 first_sgfb => basis_set_b%first_sgf, &
1416 lb_max => basis_set_b%lmax, &
1417 lb_min => basis_set_b%lmin, &
1418 npgfb => basis_set_b%npgf, &
1419 nsgfb => basis_set_b%nsgf_set, &
1420 rpgfb => basis_set_b%pgf_radius, &
1421 set_radius_b => basis_set_b%set_radius, &
1422 sphi_b => basis_set_b%sphi, &
1423 zetb => basis_set_b%zet)
1424
1425 nseta = basis_set_a%nset
1426 nsetb = basis_set_b%nset
1427
1428 ldsa = SIZE(sphi_a, 1)
1429 ldsb = SIZE(sphi_b, 1)
1430
1431 IF (iatom <= jatom) THEN
1432 irow = iatom
1433 icol = jatom
1434 ELSE
1435 irow = jatom
1436 icol = iatom
1437 END IF
1438
1439 NULLIFY (cblock)
1440 CALL dbcsr_get_block_p(matrix=cosmat, &
1441 row=irow, col=icol, block=cblock, found=found)
1442 NULLIFY (sblock)
1443 CALL dbcsr_get_block_p(matrix=sinmat, &
1444 row=irow, col=icol, block=sblock, found=found)
1445 IF (ASSOCIATED(cblock) .AND. .NOT. ASSOCIATED(sblock) .OR. &
1446 .NOT. ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1447 cpabort("")
1448 END IF
1449
1450 IF (ASSOCIATED(cblock) .AND. ASSOCIATED(sblock)) THEN
1451
1452 ra(:) = pbc(particle_set(iatom)%r(:), cell)
1453 rb(:) = ra + rab
1454 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1455
1456 DO iset = 1, nseta
1457
1458 ncoa = npgfa(iset)*ncoset(la_max(iset))
1459 sgfa = first_sgfa(1, iset)
1460
1461 DO jset = 1, nsetb
1462
1463 IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
1464
1465 ncob = npgfb(jset)*ncoset(lb_max(jset))
1466 sgfb = first_sgfb(1, jset)
1467
1468 ! Calculate the primitive integrals
1469 CALL cossin(la_max(iset), npgfa(iset), zeta(:, iset), rpgfa(:, iset), la_min(iset), &
1470 lb_max(jset), npgfb(jset), zetb(:, jset), rpgfb(:, jset), lb_min(jset), &
1471 ra, rb, kvec, cosab, sinab)
1472 CALL contract_cossin(cblock, sblock, &
1473 iatom, ncoa, nsgfa(iset), sgfa, sphi_a, ldsa, &
1474 jatom, ncob, nsgfb(jset), sgfb, sphi_b, ldsb, &
1475 cosab, sinab, ldab, work, ldwork)
1476
1477 END DO
1478 END DO
1479
1480 END IF
1481 END associate
1482 END DO
1483 CALL neighbor_list_iterator_release(nl_iterator)
1484
1485 DEALLOCATE (cosab)
1486 DEALLOCATE (sinab)
1487 DEALLOCATE (work)
1488 DEALLOCATE (basis_set_list)
1489
1490 CALL timestop(handle)
1491
1492 END SUBROUTINE build_berry_moment_matrix
1493
1494! **************************************************************************************************
1495!> \brief ...
1496!> \param qs_env ...
1497!> \param cosmat ...
1498!> \param sinmat ...
1499!> \param kvec ...
1500!> \param basis_type ...
1501! **************************************************************************************************
1502 SUBROUTINE build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
1503
1504 TYPE(qs_environment_type), POINTER :: qs_env
1505 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: cosmat, sinmat
1506 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: kvec
1507 CHARACTER(len=*), OPTIONAL :: basis_type
1508
1509 CHARACTER(LEN=*), PARAMETER :: routinen = 'build_berry_kpoint_matrix'
1510
1511 INTEGER :: handle, i, iatom, ic, icol, ikind, inode, irow, iset, jatom, jkind, jset, ldab, &
1512 ldsa, ldsb, ldwork, natom, ncoa, ncob, nimg, nkind, nseta, nsetb, sgfa, sgfb
1513 INTEGER, DIMENSION(3) :: icell
1514 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
1515 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1516 LOGICAL :: found, use_cell_mapping
1517 REAL(dp), DIMENSION(:, :), POINTER :: cblock, cosab, sblock, sinab, work
1518 REAL(kind=dp) :: dab
1519 REAL(kind=dp), DIMENSION(3) :: ra, rab, rb
1520 TYPE(cell_type), POINTER :: cell
1521 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
1522 TYPE(dft_control_type), POINTER :: dft_control
1523 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1524 TYPE(gto_basis_set_type), POINTER :: basis_set, basis_set_a, basis_set_b
1525 TYPE(kpoint_type), POINTER :: kpoints
1526 TYPE(neighbor_list_iterator_p_type), &
1527 DIMENSION(:), POINTER :: nl_iterator
1528 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1529 POINTER :: sab_orb
1530 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1531 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1532 TYPE(qs_kind_type), POINTER :: qs_kind
1533 TYPE(qs_ks_env_type), POINTER :: ks_env
1534
1535 CALL timeset(routinen, handle)
1536
1537 CALL get_qs_env(qs_env, &
1538 ks_env=ks_env, &
1539 dft_control=dft_control)
1540 nimg = dft_control%nimages
1541 IF (nimg > 1) THEN
1542 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1543 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1544 use_cell_mapping = .true.
1545 ELSE
1546 use_cell_mapping = .false.
1547 END IF
1548
1549 CALL get_qs_env(qs_env=qs_env, &
1550 qs_kind_set=qs_kind_set, &
1551 particle_set=particle_set, cell=cell, &
1552 sab_orb=sab_orb)
1553
1554 nkind = SIZE(qs_kind_set)
1555 natom = SIZE(particle_set)
1556 ALLOCATE (basis_set_list(nkind))
1557 DO ikind = 1, nkind
1558 qs_kind => qs_kind_set(ikind)
1559 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, basis_type=basis_type)
1560 IF (ASSOCIATED(basis_set)) THEN
1561 basis_set_list(ikind)%gto_basis_set => basis_set
1562 ELSE
1563 NULLIFY (basis_set_list(ikind)%gto_basis_set)
1564 END IF
1565 END DO
1566
1567 ALLOCATE (row_blk_sizes(natom))
1568 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, &
1569 basis=basis_set_list)
1570 CALL get_ks_env(ks_env, dbcsr_dist=dbcsr_dist)
1571 ! (re)allocate matrix sets
1572 CALL dbcsr_allocate_matrix_set(sinmat, 1, nimg)
1573 CALL dbcsr_allocate_matrix_set(cosmat, 1, nimg)
1574 DO i = 1, nimg
1575 ! sin
1576 ALLOCATE (sinmat(1, i)%matrix)
1577 CALL dbcsr_create(matrix=sinmat(1, i)%matrix, &
1578 name="SINMAT", &
1579 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1580 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1581 nze=0)
1582 CALL cp_dbcsr_alloc_block_from_nbl(sinmat(1, i)%matrix, sab_orb)
1583 CALL dbcsr_set(sinmat(1, i)%matrix, 0.0_dp)
1584 ! cos
1585 ALLOCATE (cosmat(1, i)%matrix)
1586 CALL dbcsr_create(matrix=cosmat(1, i)%matrix, &
1587 name="COSMAT", &
1588 dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
1589 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
1590 nze=0)
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(IN) :: 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:73
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:195
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:89
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
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:28
bibliography::mattiat2019
integer, save, public mattiat2019
Definition bibliography.F:43
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:76
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:159
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:607
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_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:744
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:125
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:536
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
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:1016
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:535
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:167
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:333
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_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, 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, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
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, 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_r3d_rs_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_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:451
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)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:864
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_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:330
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:1123
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:126
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:558
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:1336
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:790
qs_moments::build_berry_kpoint_matrix
subroutine, public build_berry_kpoint_matrix(qs_env, cosmat, sinmat, kvec, basis_type)
...
Definition qs_moments.F:1503
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:329
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:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
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:559
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:116
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:138
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:215
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:135
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