Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines for the calculation of moments from Wannier functions
10 : !> \author Lukas Schreder
11 : ! **************************************************************************************************
12 : MODULE localized_moments
13 :
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE bibliography, ONLY: Schreder2021,&
17 : Schreder2024_1,&
18 : cite_reference
19 : USE cell_types, ONLY: cell_type,&
20 : pbc
21 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm
22 : USE cp_cfm_types, ONLY: cp_cfm_create,&
23 : cp_cfm_get_element,&
24 : cp_cfm_release,&
25 : cp_cfm_type,&
26 : cp_fm_to_cfm
27 : USE cp_control_types, ONLY: dft_control_type
28 : USE cp_dbcsr_api, ONLY: &
29 : dbcsr_copy, dbcsr_create, dbcsr_get_readonly_block_p, dbcsr_init_p, &
30 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, &
31 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_set, dbcsr_type, &
32 : dbcsr_type_antisymmetric, dbcsr_type_symmetric
33 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
34 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
35 : dbcsr_allocate_matrix_set,&
36 : dbcsr_deallocate_matrix_set
37 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
38 : cp_fm_struct_release,&
39 : cp_fm_struct_type
40 : USE cp_fm_types, ONLY: cp_fm_create,&
41 : cp_fm_get_info,&
42 : cp_fm_release,&
43 : cp_fm_set_all,&
44 : cp_fm_to_fm,&
45 : cp_fm_type
46 : USE cp_log_handling, ONLY: cp_get_default_logger,&
47 : cp_logger_type
48 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
49 : cp_print_key_unit_nr
50 : USE input_constants, ONLY: use_mom_ref_com
51 : USE input_section_types, ONLY: section_vals_type
52 : USE kg_environment_types, ONLY: kg_environment_type
53 : USE kinds, ONLY: default_string_length,&
54 : dp
55 : USE message_passing, ONLY: mp_para_env_type
56 : USE molecule_kind_types, ONLY: get_molecule_kind,&
57 : molecule_kind_type
58 : USE molecule_types, ONLY: molecule_type
59 : USE moments_utils, ONLY: get_reference_point
60 : USE orbital_pointers, ONLY: current_maxl,&
61 : indco
62 : USE particle_types, ONLY: particle_type
63 : USE qs_environment_types, ONLY: get_qs_env,&
64 : qs_environment_type
65 : USE qs_kind_types, ONLY: get_qs_kind,&
66 : qs_kind_type
67 : USE qs_loc_types, ONLY: qs_loc_env_type
68 : USE qs_moments, ONLY: build_local_magmom_matrix,&
69 : build_local_moment_matrix,&
70 : build_local_moments_der_matrix,&
71 : calculate_commutator_nl_terms,&
72 : print_moments,&
73 : print_moments_nl,&
74 : set_label
75 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
76 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
77 : USE qs_rho_types, ONLY: qs_rho_get,&
78 : qs_rho_type
79 : #include "./base/base_uses.f90"
80 :
81 : IMPLICIT NONE
82 :
83 : PRIVATE
84 :
85 : ! *** public subroutines ***
86 :
87 : PUBLIC :: calculate_kg_moments, &
88 : calculate_localized_moments
89 :
90 : ! *** body ***
91 :
92 : CONTAINS
93 :
94 : ! **************************************************************************************************
95 : !> \brief Calculates localized multipole moments per molecule in the Wannier MO basis
96 : !> \param qs_env ...
97 : !> \param qs_loc_env ...
98 : !> \param mo_local ...
99 : !> \param max_moment ...
100 : !> \param magnetic ...
101 : !> \param vel_reprs ...
102 : !> \param com_nl ...
103 : !> \param loc_print_section ...
104 : ! **************************************************************************************************
105 12 : SUBROUTINE calculate_localized_moments(qs_env, qs_loc_env, mo_local, max_moment, magnetic, &
106 : vel_reprs, com_nl, loc_print_section)
107 : TYPE(qs_environment_type), POINTER :: qs_env
108 : TYPE(qs_loc_env_type), INTENT(IN) :: qs_loc_env
109 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_local
110 : INTEGER, INTENT(IN) :: max_moment
111 : LOGICAL, INTENT(IN) :: magnetic, vel_reprs, com_nl
112 : TYPE(section_vals_type), POINTER :: loc_print_section
113 :
114 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_localized_moments'
115 :
116 12 : CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
117 : CHARACTER(LEN=default_string_length) :: mol_name
118 : INTEGER :: akind, first_atom, handle, i, iatom, idir, imo_im, imo_re, imol, iounit, iproc, &
119 : ispin, ix, iy, iz, l, last_atom, nao, natom, nm, nmols, nmom, ns, nspins, order, &
120 : rmom_vel_size
121 12 : INTEGER, ALLOCATABLE, DIMENSION(:) :: states
122 : INTEGER, DIMENSION(2) :: nstates
123 : LOGICAL :: do_rtp, floating, ghost
124 : REAL(dp) :: charge, dd, im_part, re_part, zwfc
125 12 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
126 12 : nlcom_rv, nlcom_rvr, nlcom_rxrv, &
127 12 : qupole_der, rmom_vel
128 12 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
129 : REAL(dp), DIMENSION(3) :: rcc, ria
130 12 : REAL(dp), DIMENSION(:), POINTER :: ref_point
131 : TYPE(atomic_kind_type), POINTER :: atomic_kind
132 : TYPE(cell_type), POINTER :: cell
133 : TYPE(cp_cfm_type) :: ov_cfm, v_cfm, vov_cfm
134 : TYPE(cp_fm_struct_type), POINTER :: o_struct, v_struct, vov_struct
135 : TYPE(cp_fm_type) :: or_fm, vi_fm, vr_fm
136 : TYPE(cp_logger_type), POINTER :: logger
137 12 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum
138 12 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
139 : TYPE(dft_control_type), POINTER :: dft_control
140 : TYPE(molecule_kind_type), POINTER :: molecule_kind
141 12 : TYPE(molecule_type), POINTER :: molecule_set(:)
142 : TYPE(mp_para_env_type), POINTER :: para_env
143 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
144 12 : POINTER :: sab_orb
145 12 : TYPE(particle_type), POINTER :: particle_set(:)
146 12 : TYPE(qs_kind_type), POINTER :: qs_kind_set(:)
147 :
148 12 : CALL timeset(routineN, handle)
149 :
150 12 : CALL cite_reference(Schreder2021)
151 12 : IF (max_moment > 1 .OR. magnetic .OR. vel_reprs .OR. com_nl) THEN
152 12 : CALL cite_reference(Schreder2024_1)
153 : END IF
154 :
155 12 : logger => cp_get_default_logger()
156 :
157 12 : NULLIFY (cell, dft_control, matrix_s, molecule_set, qs_kind_set, sab_orb)
158 : CALL get_qs_env(qs_env, &
159 : cell=cell, &
160 : dft_control=dft_control, &
161 : matrix_s=matrix_s, &
162 : molecule_set=molecule_set, &
163 : qs_kind_set=qs_kind_set, &
164 12 : sab_orb=sab_orb)
165 12 : particle_set => qs_loc_env%particle_set
166 12 : para_env => qs_loc_env%para_env
167 :
168 12 : nspins = dft_control%nspins
169 12 : zwfc = 3.0_dp - REAL(nspins, KIND=dp)
170 : ! mo_local is sized nspins for SCF and 2*nspins for RTP (real, imag interleaved per spin)
171 12 : do_rtp = qs_env%run_rtp .OR. (SIZE(mo_local) == 2*nspins)
172 :
173 12 : nmom = MIN(max_moment, current_maxl)
174 : ! number of independent multipole components for orders 1..nmom
175 12 : nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
176 12 : nmols = SIZE(molecule_set)
177 :
178 : iounit = cp_print_key_unit_nr(logger, loc_print_section, "LOCALIZED_MOMENTS", &
179 12 : extension=".LocMom", middle_name="LOCALIZED_MOMENTS")
180 :
181 72 : ALLOCATE (rmom(nm + 1, 3), rlab(nm + 1))
182 12 : IF (magnetic) ALLOCATE (mmom(3))
183 12 : IF (vel_reprs) THEN
184 12 : IF (com_nl .AND. (nmom >= 2)) THEN
185 : rmom_vel_size = 21
186 : ELSE
187 12 : rmom_vel_size = MAX(nm, 9)
188 : END IF
189 36 : ALLOCATE (rmom_vel(rmom_vel_size))
190 : END IF
191 :
192 48 : DO imol = 1, nmols
193 36 : molecule_kind => molecule_set(imol)%molecule_kind
194 36 : first_atom = molecule_set(imol)%first_atom
195 36 : last_atom = molecule_set(imol)%last_atom
196 36 : CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom, name=mol_name)
197 36 : NULLIFY (ref_point)
198 : CALL get_reference_point(rcc, qs_env=qs_env, reference=use_mom_ref_com, &
199 : ref_point=ref_point, ifirst=first_atom, &
200 36 : ilast=last_atom)
201 :
202 36 : rmom = 0.0_dp
203 180 : rlab = ""
204 :
205 36 : NULLIFY (moments)
206 36 : CALL dbcsr_allocate_matrix_set(moments, nm)
207 144 : DO i = 1, nm
208 108 : ALLOCATE (moments(i)%matrix)
209 108 : IF (vel_reprs .AND. (nmom >= 2)) THEN
210 : CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
211 0 : matrix_type=dbcsr_type_symmetric)
212 0 : CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
213 : ELSE
214 108 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
215 : END IF
216 144 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
217 : END DO
218 36 : IF (vel_reprs .AND. (nmom >= 2)) THEN
219 0 : NULLIFY (moments_der)
220 0 : CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
221 0 : DO i = 1, 3
222 0 : DO idir = 1, 3
223 0 : CALL dbcsr_init_p(moments_der(i, idir)%matrix)
224 : CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
225 0 : matrix_type=dbcsr_type_antisymmetric)
226 0 : CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
227 0 : CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
228 : END DO
229 : END DO
230 : CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, &
231 0 : moments=moments)
232 : ELSE
233 36 : CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
234 : END IF
235 :
236 36 : IF (magnetic) THEN
237 36 : NULLIFY (magmom)
238 36 : CALL dbcsr_allocate_matrix_set(magmom, 3)
239 144 : DO i = 1, 3
240 108 : CALL dbcsr_init_p(magmom(i)%matrix)
241 : CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
242 108 : matrix_type=dbcsr_type_antisymmetric)
243 108 : CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
244 144 : CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
245 : END DO
246 36 : CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
247 36 : mmom = 0.0_dp
248 : END IF
249 :
250 36 : IF (vel_reprs) THEN
251 36 : NULLIFY (momentum)
252 36 : CALL dbcsr_allocate_matrix_set(momentum, 3)
253 144 : DO i = 1, 3
254 108 : CALL dbcsr_init_p(momentum(i)%matrix)
255 : CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
256 108 : matrix_type=dbcsr_type_antisymmetric)
257 108 : CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
258 144 : CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
259 : END DO
260 36 : CALL build_lin_mom_matrix(qs_env, momentum)
261 36 : rmom_vel = 0.0_dp
262 : END IF
263 :
264 36 : CALL cp_fm_get_info(mo_local(1), nrow_global=nao)
265 : CALL cp_fm_struct_create(o_struct, nrow_global=nao, ncol_global=nao, &
266 36 : template_fmstruct=mo_local(1)%matrix_struct)
267 36 : CALL cp_fm_create(or_fm, o_struct)
268 :
269 72 : DO ispin = 1, nspins
270 36 : IF (ASSOCIATED(molecule_set(imol)%lmi)) THEN
271 18 : nstates(1) = molecule_set(imol)%lmi(ispin)%nstates
272 : ELSE
273 18 : nstates(1) = 0
274 : END IF
275 36 : nstates(2) = para_env%mepos
276 36 : CALL para_env%maxloc(nstates)
277 36 : IF (nstates(1) == 0) CYCLE
278 36 : ns = nstates(1)
279 36 : iproc = nstates(2)
280 108 : ALLOCATE (states(ns))
281 36 : IF (iproc == para_env%mepos) THEN
282 222 : states(:) = molecule_set(imol)%lmi(ispin)%states(:)
283 : ELSE
284 18 : states(:) = 0
285 : END IF
286 36 : CALL para_env%bcast(states, iproc)
287 :
288 36 : IF (do_rtp) THEN
289 24 : imo_re = 2*ispin - 1
290 24 : imo_im = 2*ispin
291 : ELSE
292 : imo_re = ispin
293 : imo_im = 0
294 : END IF
295 :
296 : CALL cp_fm_struct_create(v_struct, nrow_global=nao, ncol_global=ns, &
297 36 : template_fmstruct=mo_local(imo_re)%matrix_struct)
298 : CALL cp_fm_struct_create(vov_struct, nrow_global=ns, ncol_global=ns, &
299 36 : template_fmstruct=mo_local(imo_re)%matrix_struct)
300 36 : CALL cp_fm_create(vr_fm, v_struct, name="vr")
301 36 : CALL cp_fm_set_all(vr_fm, 0.0_dp)
302 36 : IF (do_rtp) THEN
303 24 : CALL cp_fm_create(vi_fm, v_struct, name="vi")
304 24 : CALL cp_fm_set_all(vi_fm, 0.0_dp)
305 : END IF
306 36 : CALL cp_cfm_create(v_cfm, v_struct, name="v")
307 36 : CALL cp_cfm_create(ov_cfm, v_struct, name="ov")
308 36 : CALL cp_cfm_create(vov_cfm, vov_struct, name="vov")
309 36 : CALL cp_fm_struct_release(v_struct)
310 36 : CALL cp_fm_struct_release(vov_struct)
311 :
312 444 : DO i = 1, ns
313 408 : CALL cp_fm_to_fm(mo_local(imo_re), vr_fm, 1, states(i), i)
314 444 : IF (do_rtp) THEN
315 272 : CALL cp_fm_to_fm(mo_local(imo_im), vi_fm, 1, states(i), i)
316 : END IF
317 : END DO
318 36 : IF (do_rtp) THEN
319 24 : CALL cp_fm_to_cfm(vr_fm, vi_fm, v_cfm)
320 24 : CALL cp_fm_release(vi_fm)
321 : ELSE
322 12 : CALL cp_fm_to_cfm(vr_fm, mtarget=v_cfm)
323 : END IF
324 36 : CALL cp_fm_release(vr_fm)
325 :
326 144 : DO i = 1, nm
327 : CALL trace_op_in_mo(moments(i)%matrix, or_fm, v_cfm, ov_cfm, vov_cfm, &
328 108 : o_struct, ns, nao, re_part, im_part)
329 144 : rmom(i + 1, 1) = rmom(i + 1, 1) + zwfc*re_part
330 : END DO
331 36 : rmom(1, 1) = rmom(1, 1) + zwfc*REAL(ns, KIND=dp)
332 :
333 36 : IF (magnetic) THEN
334 144 : DO i = 1, 3
335 : CALL trace_op_in_mo(magmom(i)%matrix, or_fm, v_cfm, ov_cfm, vov_cfm, &
336 108 : o_struct, ns, nao, re_part, im_part)
337 144 : mmom(i) = mmom(i) - zwfc*im_part
338 : END DO
339 : END IF
340 :
341 36 : IF (vel_reprs) THEN
342 72 : DO order = 1, nmom
343 36 : SELECT CASE (order)
344 : CASE (1)
345 144 : DO i = 1, 3
346 : CALL trace_op_in_mo(momentum(i)%matrix, or_fm, v_cfm, ov_cfm, vov_cfm, &
347 108 : o_struct, ns, nao, re_part, im_part)
348 144 : rmom_vel(i) = rmom_vel(i) - zwfc*im_part
349 : END DO
350 : CASE (2)
351 0 : ALLOCATE (qupole_der(9))
352 0 : qupole_der = 0.0_dp
353 0 : DO i = 1, 3
354 0 : DO idir = 1, 3
355 : CALL trace_op_in_mo(moments_der(i, idir)%matrix, or_fm, v_cfm, ov_cfm, &
356 0 : vov_cfm, o_struct, ns, nao, re_part, im_part)
357 : qupole_der((i - 1)*3 + idir) = &
358 0 : qupole_der((i - 1)*3 + idir) + zwfc*(re_part - im_part)
359 : END DO
360 : END DO
361 0 : rmom_vel(4) = -2.0_dp*qupole_der(1) - rmom(1, 1)
362 0 : rmom_vel(5) = -qupole_der(2) - qupole_der(4)
363 0 : rmom_vel(6) = -qupole_der(3) - qupole_der(7)
364 0 : rmom_vel(7) = -2.0_dp*qupole_der(5) - rmom(1, 1)
365 0 : rmom_vel(8) = -qupole_der(6) - qupole_der(8)
366 0 : rmom_vel(9) = -2.0_dp*qupole_der(9) - rmom(1, 1)
367 36 : DEALLOCATE (qupole_der)
368 : END SELECT
369 : END DO
370 : END IF
371 :
372 36 : CALL cp_cfm_release(v_cfm)
373 36 : CALL cp_cfm_release(ov_cfm)
374 36 : CALL cp_cfm_release(vov_cfm)
375 144 : DEALLOCATE (states)
376 : END DO
377 :
378 36 : CALL cp_fm_release(or_fm)
379 36 : CALL cp_fm_struct_release(o_struct)
380 :
381 36 : CALL para_env%sum(rmom(:, 1))
382 36 : IF (magnetic) CALL para_env%sum(mmom)
383 36 : IF (vel_reprs) CALL para_env%sum(rmom_vel)
384 :
385 264 : DO iatom = first_atom, first_atom + natom - 1
386 1824 : ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
387 912 : ria = ria - rcc
388 228 : atomic_kind => particle_set(iatom)%atomic_kind
389 228 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
390 : CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating, &
391 228 : core_charge=charge)
392 228 : IF (ghost .OR. floating) CYCLE
393 228 : rmom(1, 2) = rmom(1, 2) - charge
394 1176 : DO l = 1, nm
395 684 : ix = indco(1, l + 1)
396 684 : iy = indco(2, l + 1)
397 684 : iz = indco(3, l + 1)
398 684 : dd = 1.0_dp
399 684 : IF (ix > 0) dd = dd*ria(1)**ix
400 684 : IF (iy > 0) dd = dd*ria(2)**iy
401 684 : IF (iz > 0) dd = dd*ria(3)**iz
402 684 : rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
403 912 : CALL set_label(rlab(l + 1), ix, iy, iz)
404 : END DO
405 : END DO
406 576 : rmom(:, :) = -rmom(:, :)
407 180 : rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
408 :
409 36 : IF (com_nl) THEN
410 36 : IF ((nmom >= 1) .AND. vel_reprs) THEN
411 36 : ALLOCATE (nlcom_rv(3))
412 36 : nlcom_rv(:) = 0.0_dp
413 : END IF
414 36 : IF ((nmom >= 2) .AND. vel_reprs) THEN
415 0 : ALLOCATE (nlcom_rrv(6), nlcom_rvr(6), nlcom_rrv_vrr(6))
416 0 : nlcom_rrv(:) = 0.0_dp
417 0 : nlcom_rvr(:) = 0.0_dp
418 0 : nlcom_rrv_vrr(:) = 0.0_dp
419 : END IF
420 36 : IF (magnetic) THEN
421 36 : ALLOCATE (nlcom_rxrv(3))
422 36 : nlcom_rxrv(:) = 0.0_dp
423 : END IF
424 : CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, &
425 36 : nlcom_rvr, nlcom_rrv_vrr, rcc)
426 : END IF
427 :
428 36 : IF (iounit > 0) THEN
429 : WRITE (UNIT=iounit, FMT='(/,T3,A,I6,2X,A)') &
430 18 : "# Molecule:", imol, TRIM(mol_name)
431 18 : IF (magnetic .AND. vel_reprs) THEN
432 : CALL print_moments(iounit, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., &
433 18 : mmom=mmom, rmom_vel=rmom_vel)
434 0 : ELSE IF (magnetic) THEN
435 : CALL print_moments(iounit, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., &
436 0 : mmom=mmom)
437 0 : ELSE IF (vel_reprs) THEN
438 : CALL print_moments(iounit, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., &
439 0 : rmom_vel=rmom_vel)
440 : ELSE
441 0 : CALL print_moments(iounit, nmom, rmom, rlab, rcc, cell, periodic=.FALSE.)
442 : END IF
443 :
444 18 : IF (com_nl) THEN
445 72 : IF (magnetic) mmom(:) = nlcom_rxrv(:)
446 72 : IF (vel_reprs .AND. (nmom >= 1)) rmom_vel(1:3) = nlcom_rv
447 18 : IF (vel_reprs .AND. (nmom >= 2)) THEN
448 0 : rmom_vel(4:9) = nlcom_rrv
449 0 : rmom_vel(10:15) = nlcom_rvr
450 0 : rmom_vel(16:21) = nlcom_rrv_vrr
451 : END IF
452 18 : IF (magnetic .AND. vel_reprs) THEN
453 18 : CALL print_moments_nl(iounit, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
454 0 : ELSE IF (magnetic) THEN
455 0 : CALL print_moments_nl(iounit, nmom, rlab, mmom=mmom)
456 0 : ELSE IF (vel_reprs) THEN
457 0 : CALL print_moments_nl(iounit, nmom, rlab, rmom_vel=rmom_vel)
458 : END IF
459 : END IF
460 : END IF
461 :
462 36 : IF (com_nl) THEN
463 36 : IF ((nmom >= 1) .AND. vel_reprs) DEALLOCATE (nlcom_rv)
464 36 : IF ((nmom >= 2) .AND. vel_reprs) DEALLOCATE (nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr)
465 36 : IF (magnetic) DEALLOCATE (nlcom_rxrv)
466 : END IF
467 36 : CALL dbcsr_deallocate_matrix_set(moments)
468 36 : IF (vel_reprs .AND. (nmom >= 2)) CALL dbcsr_deallocate_matrix_set(moments_der)
469 36 : IF (magnetic) CALL dbcsr_deallocate_matrix_set(magmom)
470 156 : IF (vel_reprs) CALL dbcsr_deallocate_matrix_set(momentum)
471 : END DO
472 :
473 12 : DEALLOCATE (rmom, rlab)
474 12 : IF (magnetic) DEALLOCATE (mmom)
475 12 : IF (vel_reprs) DEALLOCATE (rmom_vel)
476 :
477 12 : CALL cp_print_key_finished_output(iounit, logger, loc_print_section, "LOCALIZED_MOMENTS")
478 :
479 12 : CALL timestop(handle)
480 :
481 24 : END SUBROUTINE calculate_localized_moments
482 :
483 : ! **************************************************************************************************
484 : !> \brief Trace a real dbcsr operator against a complex MO matrix
485 : !> \param op_dbcsr ...
486 : !> \param or_fm ...
487 : !> \param v_cfm ...
488 : !> \param ov_cfm ...
489 : !> \param vov_cfm ...
490 : !> \param o_struct ...
491 : !> \param ns ...
492 : !> \param nao ...
493 : !> \param re_sum ...
494 : !> \param im_sum ...
495 : ! **************************************************************************************************
496 648 : SUBROUTINE trace_op_in_mo(op_dbcsr, or_fm, v_cfm, ov_cfm, vov_cfm, o_struct, ns, nao, &
497 : re_sum, im_sum)
498 : TYPE(dbcsr_type), POINTER :: op_dbcsr
499 : TYPE(cp_fm_type), INTENT(INOUT) :: or_fm
500 : TYPE(cp_cfm_type), INTENT(IN) :: v_cfm
501 : TYPE(cp_cfm_type), INTENT(INOUT) :: ov_cfm, vov_cfm
502 : TYPE(cp_fm_struct_type), POINTER :: o_struct
503 : INTEGER, INTENT(IN) :: ns, nao
504 : REAL(dp), INTENT(OUT) :: re_sum, im_sum
505 :
506 : COMPLEX(dp) :: di
507 : INTEGER :: j
508 : TYPE(cp_cfm_type) :: o_cfm
509 :
510 324 : CALL copy_dbcsr_to_fm(op_dbcsr, or_fm)
511 324 : CALL cp_cfm_create(o_cfm, o_struct)
512 324 : CALL cp_fm_to_cfm(or_fm, mtarget=o_cfm)
513 : CALL cp_cfm_gemm("N", "N", nao, ns, nao, (1._dp, 0._dp), o_cfm, v_cfm, &
514 324 : (0._dp, 0._dp), ov_cfm)
515 : CALL cp_cfm_gemm("C", "N", ns, ns, nao, (1._dp, 0._dp), v_cfm, ov_cfm, &
516 324 : (0._dp, 0._dp), vov_cfm)
517 324 : re_sum = 0.0_dp
518 324 : im_sum = 0.0_dp
519 3996 : DO j = 1, ns
520 3672 : CALL cp_cfm_get_element(vov_cfm, j, j, di)
521 3672 : re_sum = re_sum + REAL(di, KIND=dp)
522 3996 : im_sum = im_sum + AIMAG(di)
523 : END DO
524 324 : CALL cp_cfm_release(o_cfm)
525 :
526 324 : END SUBROUTINE trace_op_in_mo
527 :
528 : ! **************************************************************************************************
529 : !> \brief Calculates multipole moments per molecule from the Kim-Gordon AO density matrix
530 : !> \param qs_env ...
531 : !> \param unit_nr output unit (the open DFT%PRINT%MOMENTS unit from the caller)
532 : !> \param max_moment ...
533 : !> \param magnetic ...
534 : !> \param vel_reprs ...
535 : !> \param com_nl ...
536 : ! **************************************************************************************************
537 0 : SUBROUTINE calculate_kg_moments(qs_env, unit_nr, max_moment, magnetic, vel_reprs, com_nl)
538 : TYPE(qs_environment_type), POINTER :: qs_env
539 : INTEGER, INTENT(IN) :: unit_nr, max_moment
540 : LOGICAL, INTENT(IN) :: magnetic, vel_reprs, com_nl
541 :
542 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_kg_moments'
543 :
544 0 : CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: rlab
545 : CHARACTER(LEN=default_string_length) :: mol_name
546 : INTEGER :: akind, first_atom, handle, i, iatom, &
547 : idir, imol, ispin, ix, iy, iz, jatom, &
548 : l, last_atom, natom, nm, nmol, nmom, &
549 : nspins, rmom_vel_size
550 0 : INTEGER, DIMENSION(:), POINTER :: atom_to_molecule
551 : LOGICAL :: do_rtp, floating, found, ghost
552 : REAL(dp) :: charge, dd, factor
553 0 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: mmom, nlcom_rrv, nlcom_rrv_vrr, &
554 0 : nlcom_rv, nlcom_rvr, nlcom_rxrv, &
555 0 : qupole_der, rmom_vel
556 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: rmom
557 : REAL(dp), DIMENSION(3) :: rcc, ria
558 0 : REAL(dp), DIMENSION(:), POINTER :: ref_point
559 0 : REAL(dp), DIMENSION(:, :), POINTER :: oblock, pblock, sblock
560 : TYPE(atomic_kind_type), POINTER :: atomic_kind
561 : TYPE(cell_type), POINTER :: cell
562 : TYPE(dbcsr_iterator_type) :: iter
563 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: magmom, matrix_s, moments, momentum, &
564 0 : rho_ao, rho_ao_im
565 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: moments_der
566 : TYPE(dft_control_type), POINTER :: dft_control
567 : TYPE(kg_environment_type), POINTER :: kg_env
568 : TYPE(molecule_kind_type), POINTER :: molecule_kind
569 0 : TYPE(molecule_type), POINTER :: molecule_set(:)
570 : TYPE(mp_para_env_type), POINTER :: para_env
571 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
572 0 : POINTER :: sab_orb
573 0 : TYPE(particle_type), POINTER :: particle_set(:)
574 0 : TYPE(qs_kind_type), POINTER :: qs_kind_set(:)
575 : TYPE(qs_rho_type), POINTER :: rho
576 :
577 0 : CALL timeset(routineN, handle)
578 :
579 0 : CALL cite_reference(Schreder2021)
580 0 : IF (max_moment > 1 .OR. magnetic .OR. vel_reprs .OR. com_nl) THEN
581 0 : CALL cite_reference(Schreder2024_1)
582 : END IF
583 :
584 0 : NULLIFY (cell, dft_control, kg_env, matrix_s, para_env, particle_set, &
585 0 : qs_kind_set, rho, sab_orb)
586 : CALL get_qs_env(qs_env, &
587 : cell=cell, &
588 : dft_control=dft_control, &
589 : kg_env=kg_env, &
590 : matrix_s=matrix_s, &
591 : para_env=para_env, &
592 : particle_set=particle_set, &
593 : qs_kind_set=qs_kind_set, &
594 : rho=rho, &
595 0 : sab_orb=sab_orb)
596 :
597 0 : CPASSERT(ASSOCIATED(kg_env))
598 0 : molecule_set => kg_env%molecule_set
599 0 : atom_to_molecule => kg_env%atom_to_molecule
600 :
601 0 : nspins = dft_control%nspins
602 0 : do_rtp = qs_env%run_rtp
603 :
604 0 : nmom = MIN(max_moment, current_maxl)
605 0 : nm = (6 + 11*nmom + 6*nmom**2 + nmom**3)/6 - 1
606 0 : nmol = SIZE(molecule_set)
607 :
608 0 : NULLIFY (rho_ao, rho_ao_im)
609 0 : CALL qs_rho_get(rho, rho_ao=rho_ao)
610 0 : IF (do_rtp) CALL qs_rho_get(rho, rho_ao_im=rho_ao_im)
611 :
612 0 : ALLOCATE (rmom(nm + 1, 3), rlab(nm + 1))
613 0 : IF (magnetic) ALLOCATE (mmom(3))
614 0 : IF (vel_reprs) THEN
615 0 : IF (com_nl .AND. (nmom >= 2)) THEN
616 : rmom_vel_size = 21
617 : ELSE
618 0 : rmom_vel_size = MAX(nm, 9)
619 : END IF
620 0 : ALLOCATE (rmom_vel(rmom_vel_size))
621 : END IF
622 :
623 0 : DO imol = 1, nmol
624 0 : molecule_kind => molecule_set(imol)%molecule_kind
625 0 : first_atom = molecule_set(imol)%first_atom
626 0 : last_atom = molecule_set(imol)%last_atom
627 0 : CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom, name=mol_name)
628 0 : NULLIFY (ref_point)
629 : CALL get_reference_point(rcc, qs_env=qs_env, reference=use_mom_ref_com, &
630 0 : ref_point=ref_point, ifirst=first_atom, ilast=last_atom)
631 :
632 0 : rmom = 0.0_dp
633 0 : rlab = ""
634 :
635 0 : NULLIFY (moments)
636 0 : CALL dbcsr_allocate_matrix_set(moments, nm)
637 0 : DO i = 1, nm
638 0 : ALLOCATE (moments(i)%matrix)
639 0 : IF (vel_reprs .AND. (nmom >= 2)) THEN
640 : CALL dbcsr_create(moments(i)%matrix, template=matrix_s(1)%matrix, &
641 0 : matrix_type=dbcsr_type_symmetric)
642 0 : CALL cp_dbcsr_alloc_block_from_nbl(moments(i)%matrix, sab_orb)
643 : ELSE
644 0 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
645 : END IF
646 0 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
647 : END DO
648 0 : IF (vel_reprs .AND. (nmom >= 2)) THEN
649 0 : NULLIFY (moments_der)
650 0 : CALL dbcsr_allocate_matrix_set(moments_der, 3, 3)
651 0 : DO i = 1, 3
652 0 : DO idir = 1, 3
653 0 : CALL dbcsr_init_p(moments_der(i, idir)%matrix)
654 : CALL dbcsr_create(moments_der(i, idir)%matrix, template=matrix_s(1)%matrix, &
655 0 : matrix_type=dbcsr_type_antisymmetric)
656 0 : CALL cp_dbcsr_alloc_block_from_nbl(moments_der(i, idir)%matrix, sab_orb)
657 0 : CALL dbcsr_set(moments_der(i, idir)%matrix, 0.0_dp)
658 : END DO
659 : END DO
660 : CALL build_local_moments_der_matrix(qs_env, moments_der, 1, 2, ref_point=rcc, &
661 0 : moments=moments)
662 : ELSE
663 0 : CALL build_local_moment_matrix(qs_env, moments, nmom, ref_point=rcc)
664 : END IF
665 :
666 0 : IF (magnetic) THEN
667 0 : NULLIFY (magmom)
668 0 : CALL dbcsr_allocate_matrix_set(magmom, 3)
669 0 : DO i = 1, 3
670 0 : CALL dbcsr_init_p(magmom(i)%matrix)
671 : CALL dbcsr_create(magmom(i)%matrix, template=matrix_s(1)%matrix, &
672 0 : matrix_type=dbcsr_type_antisymmetric)
673 0 : CALL cp_dbcsr_alloc_block_from_nbl(magmom(i)%matrix, sab_orb)
674 0 : CALL dbcsr_set(magmom(i)%matrix, 0.0_dp)
675 : END DO
676 0 : CALL build_local_magmom_matrix(qs_env, magmom, nmom, ref_point=rcc)
677 0 : mmom = 0.0_dp
678 : END IF
679 :
680 0 : IF (vel_reprs) THEN
681 0 : NULLIFY (momentum)
682 0 : CALL dbcsr_allocate_matrix_set(momentum, 3)
683 0 : DO i = 1, 3
684 0 : CALL dbcsr_init_p(momentum(i)%matrix)
685 : CALL dbcsr_create(momentum(i)%matrix, template=matrix_s(1)%matrix, &
686 0 : matrix_type=dbcsr_type_antisymmetric)
687 0 : CALL cp_dbcsr_alloc_block_from_nbl(momentum(i)%matrix, sab_orb)
688 0 : CALL dbcsr_set(momentum(i)%matrix, 0.0_dp)
689 : END DO
690 0 : CALL build_lin_mom_matrix(qs_env, momentum)
691 0 : rmom_vel = 0.0_dp
692 : END IF
693 :
694 0 : IF (vel_reprs .AND. (nmom >= 2)) THEN
695 0 : ALLOCATE (qupole_der(9))
696 0 : qupole_der = 0.0_dp
697 : END IF
698 :
699 0 : DO ispin = 1, nspins
700 0 : CALL dbcsr_iterator_readonly_start(iter, rho_ao(ispin)%matrix)
701 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
702 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
703 0 : IF (SIZE(pblock) == 0) CYCLE
704 0 : IF (atom_to_molecule(iatom) /= imol) CYCLE
705 0 : IF (atom_to_molecule(jatom) /= imol) CYCLE
706 0 : factor = MERGE(2.0_dp, 1.0_dp, iatom /= jatom)
707 0 : CALL dbcsr_get_readonly_block_p(matrix_s(1)%matrix, iatom, jatom, sblock, found)
708 0 : IF (found) rmom(1, 1) = rmom(1, 1) + factor*SUM(pblock*sblock)
709 0 : DO i = 1, nm
710 0 : CALL dbcsr_get_readonly_block_p(moments(i)%matrix, iatom, jatom, oblock, found)
711 0 : IF (found) rmom(i + 1, 1) = rmom(i + 1, 1) + factor*SUM(pblock*oblock)
712 : END DO
713 : END DO
714 0 : CALL dbcsr_iterator_stop(iter)
715 :
716 0 : IF (do_rtp .AND. (magnetic .OR. vel_reprs)) THEN
717 0 : CALL dbcsr_iterator_readonly_start(iter, rho_ao_im(ispin)%matrix)
718 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
719 0 : CALL dbcsr_iterator_next_block(iter, iatom, jatom, pblock)
720 0 : IF (SIZE(pblock) == 0) CYCLE
721 0 : IF (atom_to_molecule(iatom) /= imol) CYCLE
722 0 : IF (atom_to_molecule(jatom) /= imol) CYCLE
723 0 : factor = MERGE(2.0_dp, 1.0_dp, iatom /= jatom)
724 0 : IF (magnetic) THEN
725 0 : DO i = 1, 3
726 0 : CALL dbcsr_get_readonly_block_p(magmom(i)%matrix, iatom, jatom, oblock, found)
727 0 : IF (found) mmom(i) = mmom(i) + factor*SUM(pblock*oblock)
728 : END DO
729 : END IF
730 0 : IF (vel_reprs) THEN
731 0 : DO i = 1, 3
732 0 : CALL dbcsr_get_readonly_block_p(momentum(i)%matrix, iatom, jatom, oblock, found)
733 0 : IF (found) rmom_vel(i) = rmom_vel(i) + factor*SUM(pblock*oblock)
734 : END DO
735 0 : IF (nmom >= 2) THEN
736 0 : DO i = 1, 3
737 0 : DO idir = 1, 3
738 : CALL dbcsr_get_readonly_block_p(moments_der(i, idir)%matrix, &
739 0 : iatom, jatom, oblock, found)
740 0 : IF (found) &
741 : qupole_der((i - 1)*3 + idir) = &
742 0 : qupole_der((i - 1)*3 + idir) - factor*SUM(pblock*oblock)
743 : END DO
744 : END DO
745 : END IF
746 : END IF
747 : END DO
748 0 : CALL dbcsr_iterator_stop(iter)
749 : END IF
750 : END DO
751 :
752 0 : IF (vel_reprs .AND. (nmom >= 2)) THEN
753 0 : rmom_vel(4) = -2.0_dp*qupole_der(1) - rmom(1, 1)
754 0 : rmom_vel(5) = -qupole_der(2) - qupole_der(4)
755 0 : rmom_vel(6) = -qupole_der(3) - qupole_der(7)
756 0 : rmom_vel(7) = -2.0_dp*qupole_der(5) - rmom(1, 1)
757 0 : rmom_vel(8) = -qupole_der(6) - qupole_der(8)
758 0 : rmom_vel(9) = -2.0_dp*qupole_der(9) - rmom(1, 1)
759 0 : DEALLOCATE (qupole_der)
760 : END IF
761 :
762 0 : CALL para_env%sum(rmom(:, 1))
763 0 : IF (magnetic) CALL para_env%sum(mmom)
764 0 : IF (vel_reprs) CALL para_env%sum(rmom_vel)
765 :
766 0 : DO iatom = first_atom, first_atom + natom - 1
767 0 : ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
768 0 : ria = ria - rcc
769 0 : atomic_kind => particle_set(iatom)%atomic_kind
770 0 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
771 : CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating, &
772 0 : core_charge=charge)
773 0 : IF (ghost .OR. floating) CYCLE
774 0 : rmom(1, 2) = rmom(1, 2) - charge
775 0 : DO l = 1, nm
776 0 : ix = indco(1, l + 1)
777 0 : iy = indco(2, l + 1)
778 0 : iz = indco(3, l + 1)
779 0 : dd = 1.0_dp
780 0 : IF (ix > 0) dd = dd*ria(1)**ix
781 0 : IF (iy > 0) dd = dd*ria(2)**iy
782 0 : IF (iz > 0) dd = dd*ria(3)**iz
783 0 : rmom(l + 1, 2) = rmom(l + 1, 2) - charge*dd
784 0 : CALL set_label(rlab(l + 1), ix, iy, iz)
785 : END DO
786 : END DO
787 0 : rmom(:, :) = -rmom(:, :)
788 0 : rmom(:, 3) = rmom(:, 1) + rmom(:, 2)
789 :
790 0 : IF (com_nl) THEN
791 0 : IF ((nmom >= 1) .AND. vel_reprs) THEN
792 0 : ALLOCATE (nlcom_rv(3))
793 0 : nlcom_rv(:) = 0.0_dp
794 : END IF
795 0 : IF ((nmom >= 2) .AND. vel_reprs) THEN
796 0 : ALLOCATE (nlcom_rrv(6), nlcom_rvr(6), nlcom_rrv_vrr(6))
797 0 : nlcom_rrv(:) = 0.0_dp
798 0 : nlcom_rvr(:) = 0.0_dp
799 0 : nlcom_rrv_vrr(:) = 0.0_dp
800 : END IF
801 0 : IF (magnetic) THEN
802 0 : ALLOCATE (nlcom_rxrv(3))
803 0 : nlcom_rxrv(:) = 0.0_dp
804 : END IF
805 : CALL calculate_commutator_nl_terms(qs_env, nlcom_rv, nlcom_rxrv, nlcom_rrv, &
806 0 : nlcom_rvr, nlcom_rrv_vrr, rcc)
807 : END IF
808 :
809 0 : IF (unit_nr > 0) THEN
810 : WRITE (UNIT=unit_nr, FMT='(/,T3,A,I6,2X,A)') &
811 0 : "# Molecule:", imol, TRIM(mol_name)
812 0 : IF (magnetic .AND. vel_reprs) THEN
813 : CALL print_moments(unit_nr, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., &
814 0 : mmom=mmom, rmom_vel=rmom_vel)
815 0 : ELSE IF (magnetic) THEN
816 : CALL print_moments(unit_nr, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., &
817 0 : mmom=mmom)
818 0 : ELSE IF (vel_reprs) THEN
819 : CALL print_moments(unit_nr, nmom, rmom, rlab, rcc, cell, periodic=.FALSE., &
820 0 : rmom_vel=rmom_vel)
821 : ELSE
822 0 : CALL print_moments(unit_nr, nmom, rmom, rlab, rcc, cell, periodic=.FALSE.)
823 : END IF
824 :
825 0 : IF (com_nl) THEN
826 0 : IF (magnetic) mmom(:) = nlcom_rxrv(:)
827 0 : IF (vel_reprs .AND. (nmom >= 1)) rmom_vel(1:3) = nlcom_rv
828 0 : IF (vel_reprs .AND. (nmom >= 2)) THEN
829 0 : rmom_vel(4:9) = nlcom_rrv
830 0 : rmom_vel(10:15) = nlcom_rvr
831 0 : rmom_vel(16:21) = nlcom_rrv_vrr
832 : END IF
833 0 : IF (magnetic .AND. vel_reprs) THEN
834 0 : CALL print_moments_nl(unit_nr, nmom, rlab, mmom=mmom, rmom_vel=rmom_vel)
835 0 : ELSE IF (magnetic) THEN
836 0 : CALL print_moments_nl(unit_nr, nmom, rlab, mmom=mmom)
837 0 : ELSE IF (vel_reprs) THEN
838 0 : CALL print_moments_nl(unit_nr, nmom, rlab, rmom_vel=rmom_vel)
839 : END IF
840 : END IF
841 : END IF
842 :
843 0 : IF (com_nl) THEN
844 0 : IF ((nmom >= 1) .AND. vel_reprs) DEALLOCATE (nlcom_rv)
845 0 : IF ((nmom >= 2) .AND. vel_reprs) DEALLOCATE (nlcom_rrv, nlcom_rvr, nlcom_rrv_vrr)
846 0 : IF (magnetic) DEALLOCATE (nlcom_rxrv)
847 : END IF
848 0 : CALL dbcsr_deallocate_matrix_set(moments)
849 0 : IF (vel_reprs .AND. (nmom >= 2)) CALL dbcsr_deallocate_matrix_set(moments_der)
850 0 : IF (magnetic) CALL dbcsr_deallocate_matrix_set(magmom)
851 0 : IF (vel_reprs) CALL dbcsr_deallocate_matrix_set(momentum)
852 : END DO
853 :
854 0 : DEALLOCATE (rmom, rlab)
855 0 : IF (magnetic) DEALLOCATE (mmom)
856 0 : IF (vel_reprs) DEALLOCATE (rmom_vel)
857 :
858 0 : CALL timestop(handle)
859 :
860 0 : END SUBROUTINE calculate_kg_moments
861 :
862 : END MODULE localized_moments
|
