(git:6d276e9)
Loading...
Searching...
No Matches
qs_linres_op.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the operators p rxp and D needed in the optimization
10!> of the different contribution of the firs order response orbitals
11!> in a epr calculation
12!> \note
13!> The interactions are considered only within the minimum image convention
14!> \par History
15!> created 07-2005 [MI]
16!> \author MI
17! **************************************************************************************************
18MODULE qs_linres_op
19 USE cell_types, ONLY: cell_type,&
20 pbc
21 USE cp_array_utils, ONLY: cp_2d_i_p_type,&
22 cp_2d_r_p_type
23 USE cp_cfm_basic_linalg, ONLY: cp_cfm_solve
24 USE cp_cfm_types, ONLY: cp_cfm_create,&
25 cp_cfm_get_info,&
26 cp_cfm_release,&
27 cp_cfm_set_all,&
28 cp_cfm_type
29 USE cp_control_types, ONLY: dft_control_type
30 USE cp_dbcsr_api, ONLY: &
31 dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
32 dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
33 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
34 dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
35 USE cp_dbcsr_contrib, ONLY: dbcsr_checksum
36 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
37 USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
38 dbcsr_allocate_matrix_set,&
39 dbcsr_deallocate_matrix_set
40 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
41 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
42 cp_fm_struct_release,&
43 cp_fm_struct_type
44 USE cp_fm_types, ONLY: cp_fm_create,&
45 cp_fm_get_info,&
46 cp_fm_get_submatrix,&
47 cp_fm_release,&
48 cp_fm_set_all,&
49 cp_fm_set_submatrix,&
50 cp_fm_to_fm,&
51 cp_fm_type
52 USE cp_log_handling, ONLY: cp_get_default_logger,&
53 cp_logger_type,&
54 cp_to_string
55 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
56 cp_print_key_unit_nr
57 USE input_section_types, ONLY: section_vals_get_subs_vals,&
58 section_vals_type
59 USE kinds, ONLY: dp
60 USE mathconstants, ONLY: twopi
61 USE message_passing, ONLY: mp_para_env_type
62 USE molecule_types, ONLY: molecule_of_atom,&
63 molecule_type
64 USE orbital_pointers, ONLY: coset
65 USE parallel_gemm_api, ONLY: parallel_gemm
66 USE particle_methods, ONLY: get_particle_set
67 USE particle_types, ONLY: particle_type
68 USE qs_dcdr_utils, ONLY: multiply_localization,&
69 shift_wannier_into_cell
70 USE qs_elec_field, ONLY: build_efg_matrix
71 USE qs_environment_types, ONLY: get_qs_env,&
72 qs_environment_type
73 USE qs_fermi_contact, ONLY: build_fermi_contact_matrix
74 USE qs_kind_types, ONLY: get_qs_kind_set,&
75 qs_kind_type
76 USE qs_linres_types, ONLY: current_env_type,&
77 dcdr_env_type,&
78 get_current_env,&
79 get_issc_env,&
80 get_polar_env,&
81 issc_env_type,&
82 linres_control_type,&
83 polar_env_type
84 USE qs_mo_types, ONLY: get_mo_set,&
85 mo_set_type
86 USE qs_moments, ONLY: build_berry_moment_matrix,&
87 build_local_moment_matrix
88 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
89 USE qs_operators_ao, ONLY: build_ang_mom_matrix,&
90 build_lin_mom_matrix
91 USE qs_spin_orbit, ONLY: build_pso_matrix
92#include "./base/base_uses.f90"
93
94 IMPLICIT NONE
95
96 PRIVATE
97 PUBLIC :: current_operators, issc_operators, fac_vecp, ind_m2, set_vecp, set_vecp_rev, &
98 fm_scale_by_pbc_ac, polar_operators, polar_operators_local, &
99 polar_operators_local_wannier, polar_operators_berry
100
101 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_op'
102
103! **************************************************************************************************
104
105CONTAINS
106
107! **************************************************************************************************
108!> \brief Calculate the first order hamiltonian applied to the ao
109!> and then apply them to the ground state orbitals,
110!> the h1_psi1 full matrices are then ready to solve the
111!> non-homogeneous linear equations that give the psi1
112!> linear response orbitals.
113!> \param current_env ...
114!> \param qs_env ...
115!> \par History
116!> 07.2005 created [MI]
117!> \author MI
118!> \note
119!> For the operators rxp and D the h1 depends on the psi0 to which
120!> is applied, or better the center of charge of the psi0 is
121!> used to define the position operator
122!> The centers of the orbitals result form the orbital localization procedure
123!> that typically uses the berry phase operator to define the Wannier centers.
124! **************************************************************************************************
125 SUBROUTINE current_operators(current_env, qs_env)
126
127 TYPE(current_env_type) :: current_env
128 TYPE(qs_environment_type), POINTER :: qs_env
129
130 CHARACTER(LEN=*), PARAMETER :: routinen = 'current_operators'
131
132 INTEGER :: handle, iao, icenter, idir, ii, iii, &
133 ispin, istate, j, nao, natom, &
134 nbr_center(2), nmo, nsgf, nspins, &
135 nstates(2), output_unit
136 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
137 INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
138 REAL(dp) :: chk(3), ck(3), ckdk(3), dk(3)
139 REAL(dp), DIMENSION(:, :), POINTER :: basisfun_center, vecbuf_c0
140 TYPE(cell_type), POINTER :: cell
141 TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
142 TYPE(cp_2d_r_p_type), DIMENSION(3) :: vecbuf_rmdc0
143 TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
144 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
145 TYPE(cp_fm_type) :: fm_work1
146 TYPE(cp_fm_type), DIMENSION(3) :: fm_rmd_mos
147 TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
148 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: p_psi0, rxp_psi0
149 TYPE(cp_fm_type), POINTER :: mo_coeff
150 TYPE(cp_logger_type), POINTER :: logger
151 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
152 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_ao
153 TYPE(dft_control_type), POINTER :: dft_control
154 TYPE(linres_control_type), POINTER :: linres_control
155 TYPE(mp_para_env_type), POINTER :: para_env
156 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
157 POINTER :: sab_all, sab_orb
158 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
159 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
160 TYPE(section_vals_type), POINTER :: lr_section
161
162 CALL timeset(routinen, handle)
163
164 NULLIFY (qs_kind_set, cell, dft_control, linres_control, &
165 logger, particle_set, lr_section, &
166 basisfun_center, centers_set, center_list, p_psi0, &
167 rxp_psi0, vecbuf_c0, psi0_order, &
168 mo_coeff, op_ao, sab_all)
169
170 logger => cp_get_default_logger()
171 lr_section => section_vals_get_subs_vals(qs_env%input, &
172 "PROPERTIES%LINRES")
173
174 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
175 extension=".linresLog")
176 IF (output_unit > 0) THEN
177 WRITE (output_unit, fmt="(T2,A,/)") &
178 "CURRENT| Calculation of the p and (r-d)xp operators applied to psi0"
179 END IF
180
181 CALL get_qs_env(qs_env=qs_env, &
182 qs_kind_set=qs_kind_set, &
183 cell=cell, &
184 dft_control=dft_control, &
185 linres_control=linres_control, &
186 para_env=para_env, &
187 particle_set=particle_set, &
188 sab_all=sab_all, &
189 sab_orb=sab_orb, &
190 dbcsr_dist=dbcsr_dist)
191
192 nspins = dft_control%nspins
193
194 CALL get_current_env(current_env=current_env, nao=nao, centers_set=centers_set, &
195 center_list=center_list, basisfun_center=basisfun_center, &
196 nbr_center=nbr_center, p_psi0=p_psi0, rxp_psi0=rxp_psi0, &
197 psi0_order=psi0_order, &
198 nstates=nstates)
199
200 ALLOCATE (vecbuf_c0(1, nao))
201 DO idir = 1, 3
202 NULLIFY (vecbuf_rmdc0(idir)%array)
203 ALLOCATE (vecbuf_rmdc0(idir)%array(1, nao))
204 END DO
205
206 CALL get_qs_kind_set(qs_kind_set=qs_kind_set, nsgf=nsgf)
207
208 natom = SIZE(particle_set, 1)
209 ALLOCATE (first_sgf(natom))
210 ALLOCATE (last_sgf(natom))
211
212 CALL get_particle_set(particle_set, qs_kind_set, &
213 first_sgf=first_sgf, &
214 last_sgf=last_sgf)
215
216 ! Calculate the (r - dk)xp operator applied to psi0k
217 ! One possible way to go is to use the distributive property of the vector product and calculatr
218 ! (r-c)xp + (c-d)xp
219 ! where c depends on the contracted functions and not on the states
220 ! d is the center of a specific state and a loop over states is needed
221 ! the second term can be added in a second moment as a correction
222 ! notice: (r-c) and p are operators, whereas (c-d) is a multiplicative factor
223
224 ! !First term: operator matrix elements
225 ! CALL rmc_x_p_xyz_ao(op_rmd_ao,qs_env,minimum_image=.FALSE.)
226 !************************************************************
227 !
228 ! Since many psi0 vector can have the same center, depending on how the center is selected,
229 ! the (r - dk)xp operator matrix is computed Ncenter times,
230 ! where Ncenter is the total number of different centers
231 ! and each time it is multiplied by all the psi0 with center dk to get the rxp_psi0 matrix
232
233 !
234 ! prepare for allocation
235 ALLOCATE (row_blk_sizes(natom))
236 CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
237 !
238 !
239 CALL dbcsr_allocate_matrix_set(op_ao, 3)
240 ALLOCATE (op_ao(1)%matrix, op_ao(2)%matrix, op_ao(3)%matrix)
241
242 CALL dbcsr_create(matrix=op_ao(1)%matrix, &
243 name="op_ao", &
244 dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
245 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
246 mutable_work=.true.)
247 CALL cp_dbcsr_alloc_block_from_nbl(op_ao(1)%matrix, sab_all)
248 CALL dbcsr_set(op_ao(1)%matrix, 0.0_dp)
249
250 DO idir = 2, 3
251 CALL dbcsr_copy(op_ao(idir)%matrix, op_ao(1)%matrix, &
252 "op_ao"//"-"//trim(adjustl(cp_to_string(idir))))
253 CALL dbcsr_set(op_ao(idir)%matrix, 0.0_dp)
254 END DO
255
256 chk(:) = 0.0_dp
257 DO ispin = 1, nspins
258 mo_coeff => psi0_order(ispin)
259 nmo = nstates(ispin)
260 CALL cp_fm_set_all(p_psi0(ispin, 1), 0.0_dp)
261 CALL cp_fm_set_all(p_psi0(ispin, 2), 0.0_dp)
262 CALL cp_fm_set_all(p_psi0(ispin, 3), 0.0_dp)
263 DO icenter = 1, nbr_center(ispin)
264 CALL dbcsr_set(op_ao(1)%matrix, 0.0_dp)
265 CALL dbcsr_set(op_ao(2)%matrix, 0.0_dp)
266 CALL dbcsr_set(op_ao(3)%matrix, 0.0_dp)
267 !CALL rmc_x_p_xyz_ao(op_ao,qs_env,minimum_image=.FALSE.,&
268 ! & wancen=centers_set(ispin)%array(1:3,icenter))
269 ! &
270 CALL build_ang_mom_matrix(qs_env, op_ao, centers_set(ispin)%array(1:3, icenter))
271 !
272 ! accumulate checksums
273 chk(1) = chk(1) + dbcsr_checksum(op_ao(1)%matrix)
274 chk(2) = chk(2) + dbcsr_checksum(op_ao(2)%matrix)
275 chk(3) = chk(3) + dbcsr_checksum(op_ao(3)%matrix)
276 DO idir = 1, 3
277 CALL cp_fm_set_all(rxp_psi0(ispin, idir), 0.0_dp)
278 CALL cp_dbcsr_sm_fm_multiply(op_ao(idir)%matrix, mo_coeff, &
279 rxp_psi0(ispin, idir), ncol=nmo, &
280 alpha=-1.0_dp)
281 DO j = center_list(ispin)%array(1, icenter), center_list(ispin)%array(1, icenter + 1) - 1
282 istate = center_list(ispin)%array(2, j)
283 ! the p_psi0 fm is used as temporary matrix to store the results for the psi0 centered in dk
284 CALL cp_fm_to_fm(rxp_psi0(ispin, idir), &
285 p_psi0(ispin, idir), 1, istate, istate)
286 END DO
287 END DO
288 END DO
289 CALL cp_fm_to_fm(p_psi0(ispin, 1), rxp_psi0(ispin, 1))
290 CALL cp_fm_to_fm(p_psi0(ispin, 2), rxp_psi0(ispin, 2))
291 CALL cp_fm_to_fm(p_psi0(ispin, 3), rxp_psi0(ispin, 3))
292 END DO
293 !
294 CALL dbcsr_deallocate_matrix_set(op_ao)
295 !
296 ! print checksums
297 IF (output_unit > 0) THEN
298 WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum L_x =', chk(1)
299 WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum L_y =', chk(2)
300 WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum L_z =', chk(3)
301 END IF
302 !
303 ! Calculate the px py pz operators
304 CALL dbcsr_allocate_matrix_set(op_ao, 3)
305 ALLOCATE (op_ao(1)%matrix, op_ao(2)%matrix, op_ao(3)%matrix)
306
307 CALL dbcsr_create(matrix=op_ao(1)%matrix, &
308 name="op_ao", &
309 dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
310 row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
311 mutable_work=.true.)
312 CALL cp_dbcsr_alloc_block_from_nbl(op_ao(1)%matrix, sab_orb)
313 CALL dbcsr_set(op_ao(1)%matrix, 0.0_dp)
314
315 DO idir = 2, 3
316 CALL dbcsr_copy(op_ao(idir)%matrix, op_ao(1)%matrix, &
317 "op_ao"//"-"//trim(adjustl(cp_to_string(idir))))
318 CALL dbcsr_set(op_ao(idir)%matrix, 0.0_dp)
319 END DO
320 !
321 CALL build_lin_mom_matrix(qs_env, op_ao)
322 !
323 ! print checksums
324 chk(1) = dbcsr_checksum(op_ao(1)%matrix)
325 chk(2) = dbcsr_checksum(op_ao(2)%matrix)
326 chk(3) = dbcsr_checksum(op_ao(3)%matrix)
327 IF (output_unit > 0) THEN
328 WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum P_x =', chk(1)
329 WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum P_y =', chk(2)
330 WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum P_z =', chk(3)
331 END IF
332 ! Apply the p operator to the psi0
333 DO idir = 1, 3
334 DO ispin = 1, nspins
335 mo_coeff => psi0_order(ispin)
336 nmo = nstates(ispin)
337 CALL cp_fm_set_all(p_psi0(ispin, idir), 0.0_dp)
338 CALL cp_dbcsr_sm_fm_multiply(op_ao(idir)%matrix, mo_coeff, &
339 p_psi0(ispin, idir), ncol=nmo, &
340 alpha=-1.0_dp)
341 END DO
342 END DO
343 !
344 CALL dbcsr_deallocate_matrix_set(op_ao)
345 !
346 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
347 "PRINT%PROGRAM_RUN_INFO")
348
349! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
350 ! This part is not necessary with the present implementation
351 ! the angular momentum operator is computed directly for each dk independently
352 ! and multiplied by the proper psi0 (i.e. those centered in dk)
353 ! If Wannier centers are used, and no grouping of states with close centers is applied
354 ! the (r-dk)xp operator is computed Nstate times and each time applied to only one vector psi0
355 !
356 ! Apply the (r-c)xp operator to the psi0
357 !DO ispin = 1,nspins
358 ! CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo, homo=homo)
359 ! DO idir = 1,3
360 ! CALL cp_fm_set_all(rxp_psi0(ispin,idir),0.0_dp)
361 ! CALL cp_sm_fm_multiply(op_rmd_ao(idir)%matrix,mo_coeff,&
362 ! rxp_psi0(ispin,idir),ncol=nmo,alpha=-1.0_dp)
363 ! END DO
364 !END DO
365
366 !Calculate the second term of the operator state by state
367 !!!! what follows is a way to avoid calculating the L matrix for each centers.
368 !!!! not tested
369 IF (.false.) THEN
370 DO ispin = 1, nspins
371 ! Allocate full matrices as working storage in the calculation
372 ! of the rxp operator matrix. 3 matrices for the 3 Cartesian direction
373 ! plus one to apply the momentum oprator to the modified mos fm
374 mo_coeff => psi0_order(ispin)
375 nmo = nstates(ispin)
376 NULLIFY (tmp_fm_struct)
377 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
378 ncol_global=nmo, para_env=para_env, &
379 context=mo_coeff%matrix_struct%context)
380 DO idir = 1, 3
381 CALL cp_fm_create(fm_rmd_mos(idir), tmp_fm_struct)
382 END DO
383 CALL cp_fm_create(fm_work1, tmp_fm_struct)
384 CALL cp_fm_struct_release(tmp_fm_struct)
385
386 ! This part should be done better, using the full matrix distribution
387 DO istate = 1, nmo
388 CALL cp_fm_get_submatrix(psi0_order(ispin), vecbuf_c0, 1, istate, nao, 1, &
389 transpose=.true.)
390 !center of the localized psi0 state istate
391 dk(1:3) = centers_set(ispin)%array(1:3, istate)
392 DO idir = 1, 3
393 ! This loop should be distributed over the processors
394 DO iao = 1, nao
395 ck(1:3) = basisfun_center(1:3, iao)
396 ckdk = pbc(dk, ck, cell)
397 vecbuf_rmdc0(idir)%array(1, iao) = vecbuf_c0(1, iao)*ckdk(idir)
398 END DO ! iao
399 CALL cp_fm_set_submatrix(fm_rmd_mos(idir), vecbuf_rmdc0(idir)%array, &
400 1, istate, nao, 1, transpose=.true.)
401 END DO ! idir
402 END DO ! istate
403
404 DO idir = 1, 3
405 CALL set_vecp(idir, ii, iii)
406
407 !Add the second term to the idir component
408 CALL cp_fm_set_all(fm_work1, 0.0_dp)
409 CALL cp_dbcsr_sm_fm_multiply(op_ao(iii)%matrix, fm_rmd_mos(ii), &
410 fm_work1, ncol=nmo, alpha=-1.0_dp)
411 CALL cp_fm_scale_and_add(1.0_dp, rxp_psi0(ispin, idir), &
412 1.0_dp, fm_work1)
413
414 CALL cp_fm_set_all(fm_work1, 0.0_dp)
415 CALL cp_dbcsr_sm_fm_multiply(op_ao(ii)%matrix, fm_rmd_mos(iii), &
416 fm_work1, ncol=nmo, alpha=-1.0_dp)
417 CALL cp_fm_scale_and_add(1.0_dp, rxp_psi0(ispin, idir), &
418 -1.0_dp, fm_work1)
419
420 END DO ! idir
421
422 DO idir = 1, 3
423 CALL cp_fm_release(fm_rmd_mos(idir))
424 END DO
425 CALL cp_fm_release(fm_work1)
426
427 END DO ! ispin
428 END IF
429
430 DEALLOCATE (row_blk_sizes)
431
432 DEALLOCATE (first_sgf, last_sgf)
433
434 DEALLOCATE (vecbuf_c0)
435 DO idir = 1, 3
436 DEALLOCATE (vecbuf_rmdc0(idir)%array)
437 END DO
438
439 CALL timestop(handle)
440
441 END SUBROUTINE current_operators
442
443! **************************************************************************************************
444!> \brief ...
445!> \param issc_env ...
446!> \param qs_env ...
447!> \param iatom ...
448! **************************************************************************************************
449 SUBROUTINE issc_operators(issc_env, qs_env, iatom)
450
451 TYPE(issc_env_type) :: issc_env
452 TYPE(qs_environment_type), POINTER :: qs_env
453 INTEGER, INTENT(IN) :: iatom
454
455 CHARACTER(LEN=*), PARAMETER :: routinen = 'issc_operators'
456
457 INTEGER :: handle, idir, ispin, nmo, nspins, &
458 output_unit
459 LOGICAL :: do_dso, do_fc, do_pso, do_sd
460 REAL(dp) :: chk(20), r_i(3)
461 TYPE(cell_type), POINTER :: cell
462 TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0
463 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dso_psi0, efg_psi0, pso_psi0
464 TYPE(cp_fm_type), POINTER :: mo_coeff
465 TYPE(cp_logger_type), POINTER :: logger
466 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dso, matrix_efg, matrix_fc, &
467 matrix_pso
468 TYPE(dft_control_type), POINTER :: dft_control
469 TYPE(linres_control_type), POINTER :: linres_control
470 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
471 TYPE(mp_para_env_type), POINTER :: para_env
472 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
473 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
474 TYPE(section_vals_type), POINTER :: lr_section
475
476 CALL timeset(routinen, handle)
477
478 NULLIFY (matrix_fc, matrix_pso, matrix_efg)
479 NULLIFY (efg_psi0, pso_psi0, fc_psi0)
480
481 logger => cp_get_default_logger()
482 lr_section => section_vals_get_subs_vals(qs_env%input, &
483 "PROPERTIES%LINRES")
484
485 output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
486 extension=".linresLog")
487
488 CALL get_qs_env(qs_env=qs_env, &
489 qs_kind_set=qs_kind_set, &
490 cell=cell, &
491 dft_control=dft_control, &
492 linres_control=linres_control, &
493 para_env=para_env, &
494 mos=mos, &
495 particle_set=particle_set)
496
497 nspins = dft_control%nspins
498
499 CALL get_issc_env(issc_env=issc_env, &
500 matrix_efg=matrix_efg, & !this is used only here alloc/dealloc here???
501 matrix_pso=matrix_pso, & !this is used only here alloc/dealloc here???
502 matrix_fc=matrix_fc, & !this is used only here alloc/dealloc here???
503 matrix_dso=matrix_dso, & !this is used only here alloc/dealloc here???
504 efg_psi0=efg_psi0, &
505 pso_psi0=pso_psi0, &
506 dso_psi0=dso_psi0, &
507 fc_psi0=fc_psi0, &
508 do_fc=do_fc, &
509 do_sd=do_sd, &
510 do_pso=do_pso, &
511 do_dso=do_dso)
512 !
513 !
514 r_i = particle_set(iatom)%r !pbc(particle_set(iatom)%r,cell)
515 !write(*,*) 'issc_operators iatom=',iatom,' r_i=',r_i
516 chk = 0.0_dp
517 !
518 !
519 !
520 ! Fermi contact integral
521 !IF(do_fc) THEN
522 IF (.true.) THEN ! for the moment we build it (regs)
523 CALL dbcsr_set(matrix_fc(1)%matrix, 0.0_dp)
524 CALL build_fermi_contact_matrix(qs_env, matrix_fc, r_i)
525
526 chk(1) = dbcsr_checksum(matrix_fc(1)%matrix)
527
528 IF (output_unit > 0) THEN
529 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| fermi_contact: CheckSum =', chk(1)
530 END IF
531 END IF
532 !
533 ! spin-orbit integral
534 !IF(do_pso) THEN
535 IF (.true.) THEN ! for the moment we build it (regs)
536 CALL dbcsr_set(matrix_pso(1)%matrix, 0.0_dp)
537 CALL dbcsr_set(matrix_pso(2)%matrix, 0.0_dp)
538 CALL dbcsr_set(matrix_pso(3)%matrix, 0.0_dp)
539 CALL build_pso_matrix(qs_env, matrix_pso, r_i)
540
541 chk(2) = dbcsr_checksum(matrix_pso(1)%matrix)
542 chk(3) = dbcsr_checksum(matrix_pso(2)%matrix)
543 chk(4) = dbcsr_checksum(matrix_pso(3)%matrix)
544
545 IF (output_unit > 0) THEN
546 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| pso_x: CheckSum =', chk(2)
547 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| pso_y: CheckSum =', chk(3)
548 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| pso_z: CheckSum =', chk(4)
549 END IF
550 END IF
551 !
552 ! electric field integral
553 !IF(do_sd) THEN
554 IF (.true.) THEN ! for the moment we build it (regs)
555 CALL dbcsr_set(matrix_efg(1)%matrix, 0.0_dp)
556 CALL dbcsr_set(matrix_efg(2)%matrix, 0.0_dp)
557 CALL dbcsr_set(matrix_efg(3)%matrix, 0.0_dp)
558 CALL dbcsr_set(matrix_efg(4)%matrix, 0.0_dp)
559 CALL dbcsr_set(matrix_efg(5)%matrix, 0.0_dp)
560 CALL dbcsr_set(matrix_efg(6)%matrix, 0.0_dp)
561 CALL build_efg_matrix(qs_env, matrix_efg, r_i)
562
563 chk(5) = dbcsr_checksum(matrix_efg(1)%matrix)
564 chk(6) = dbcsr_checksum(matrix_efg(2)%matrix)
565 chk(7) = dbcsr_checksum(matrix_efg(3)%matrix)
566 chk(8) = dbcsr_checksum(matrix_efg(4)%matrix)
567 chk(9) = dbcsr_checksum(matrix_efg(5)%matrix)
568 chk(10) = dbcsr_checksum(matrix_efg(6)%matrix)
569
570 IF (output_unit > 0) THEN
571 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg (3xx-rr)/3: CheckSum =', chk(5)
572 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg (3yy-rr)/3: CheckSum =', chk(6)
573 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg (3zz-rr)/3: CheckSum =', chk(7)
574 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg xy: CheckSum =', chk(8)
575 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg xz: CheckSum =', chk(9)
576 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg yz: CheckSum =', chk(10)
577 END IF
578 END IF
579 !
580 !
581 IF (output_unit > 0) THEN
582 WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| all operator: CheckSum =', sum(chk(1:10))
583 END IF
584 !
585 !>>> debugging only here we build the dipole matrix... debugging the kernel...
586 IF (do_dso) THEN
587 CALL dbcsr_set(matrix_dso(1)%matrix, 0.0_dp)
588 CALL dbcsr_set(matrix_dso(2)%matrix, 0.0_dp)
589 CALL dbcsr_set(matrix_dso(3)%matrix, 0.0_dp)
590 CALL build_local_moment_matrix(qs_env, matrix_dso, 1, &
591 ref_point=[0.0_dp, 0.0_dp, 0.0_dp], all_images=.true.)
592 END IF
593 !
594 ! multiply by the mos
595 DO ispin = 1, nspins
596 !
597 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
598 CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
599 !
600 ! EFG
601 IF (do_sd) THEN
602 DO idir = 1, 6
603 CALL cp_dbcsr_sm_fm_multiply(matrix_efg(idir)%matrix, mo_coeff, &
604 efg_psi0(ispin, idir), ncol=nmo, &
605 alpha=1.0_dp)
606 END DO
607 END IF
608 !
609 ! PSO
610 IF (do_pso) THEN
611 DO idir = 1, 3
612 CALL cp_dbcsr_sm_fm_multiply(matrix_pso(idir)%matrix, mo_coeff, &
613 pso_psi0(ispin, idir), ncol=nmo, &
614 alpha=-1.0_dp)
615 END DO
616 END IF
617 !
618 ! FC
619 IF (do_fc) THEN
620 CALL cp_dbcsr_sm_fm_multiply(matrix_fc(1)%matrix, mo_coeff, &
621 fc_psi0(ispin), ncol=nmo, &
622 alpha=1.0_dp)
623 END IF
624 !
625 !>>> for debugging only
626 IF (do_dso) THEN
627 DO idir = 1, 3
628 CALL cp_dbcsr_sm_fm_multiply(matrix_dso(idir)%matrix, mo_coeff, &
629 dso_psi0(ispin, idir), ncol=nmo, &
630 alpha=-1.0_dp)
631 END DO
632 END IF
633 !<<< for debugging only
634 END DO
635
636 CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
637 "PRINT%PROGRAM_RUN_INFO")
638
639 CALL timestop(handle)
640
641 END SUBROUTINE issc_operators
642
643! **************************************************************************************************
644!> \brief Calculate the dipole operator in the AO basis and its derivative wrt to MOs
645!>
646!> \param qs_env ...
647! **************************************************************************************************
648 SUBROUTINE polar_operators(qs_env)
649
650 TYPE(qs_environment_type), POINTER :: qs_env
651
652 LOGICAL :: do_periodic
653 TYPE(dft_control_type), POINTER :: dft_control
654 TYPE(polar_env_type), POINTER :: polar_env
655
656 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, polar_env=polar_env)
657 CALL get_polar_env(polar_env=polar_env, do_periodic=do_periodic)
658 IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
659 IF (do_periodic) THEN
660 CALL polar_tb_operators_berry(qs_env)
661 ELSE
662 CALL polar_tb_operators_local(qs_env)
663 END IF
664 ELSE
665 IF (do_periodic) THEN
666 CALL polar_operators_berry(qs_env)
667 ELSE
668 CALL polar_operators_local(qs_env)
669 END IF
670 END IF
671
672 END SUBROUTINE polar_operators
673
674! **************************************************************************************************
675!> \brief Calculate the Berry phase operator in the AO basis and
676!> then the derivative of the Berry phase operator with respect to
677!> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
678!> afterwards multiply with the ground state MO coefficients
679!>
680!> \param qs_env ...
681!> \par History
682!> 01.2013 created [SL]
683!> 06.2018 polar_env integrated into qs_env (MK)
684!> \author SL
685! **************************************************************************************************
686
687 SUBROUTINE polar_operators_berry(qs_env)
688
689 TYPE(qs_environment_type), POINTER :: qs_env
690
691 CHARACTER(LEN=*), PARAMETER :: routinen = 'polar_operators_berry'
692 COMPLEX(KIND=dp), PARAMETER :: one = (1.0_dp, 0.0_dp), &
693 zero = (0.0_dp, 0.0_dp)
694
695 COMPLEX(DP) :: zdet, zdeta
696 INTEGER :: handle, i, idim, ispin, nao, nmo, &
697 nspins, tmp_dim, z
698 LOGICAL :: do_raman
699 REAL(dp) :: kvec(3), maxocc
700 TYPE(cell_type), POINTER :: cell
701 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
702 TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_mat
703 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
704 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set, opvec
705 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: inv_work
706 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0
707 TYPE(cp_fm_type), POINTER :: mo_coeff
708 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
709 TYPE(dbcsr_type), POINTER :: cosmat, sinmat
710 TYPE(dft_control_type), POINTER :: dft_control
711 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
712 TYPE(mp_para_env_type), POINTER :: para_env
713 TYPE(polar_env_type), POINTER :: polar_env
714
715 CALL timeset(routinen, handle)
716
717 NULLIFY (dberry_psi0, sinmat, cosmat)
718 NULLIFY (polar_env)
719
720 NULLIFY (cell, dft_control, mos, matrix_s)
721 CALL get_qs_env(qs_env=qs_env, &
722 cell=cell, &
723 dft_control=dft_control, &
724 para_env=para_env, &
725 polar_env=polar_env, &
726 mos=mos, &
727 matrix_s=matrix_s)
728
729 nspins = dft_control%nspins
730
731 CALL get_polar_env(polar_env=polar_env, &
732 do_raman=do_raman, &
733 dberry_psi0=dberry_psi0)
734 !SL calculate dipole berry phase
735 IF (do_raman) THEN
736
737 DO i = 1, 3
738 DO ispin = 1, nspins
739 CALL cp_fm_set_all(dberry_psi0(i, ispin), 0.0_dp)
740 END DO
741 END DO
742
743 ! initialize all work matrices needed
744 ALLOCATE (opvec(2, dft_control%nspins))
745 ALLOCATE (op_fm_set(2, dft_control%nspins))
746 ALLOCATE (eigrmat(dft_control%nspins))
747 ALLOCATE (inv_mat(3, dft_control%nspins))
748 ALLOCATE (inv_work(2, 3, dft_control%nspins))
749
750 ! A bit to allocate for the wavefunction
751 DO ispin = 1, dft_control%nspins
752 NULLIFY (tmp_fm_struct, mo_coeff)
753 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
754 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
755 ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
756 DO i = 1, SIZE(op_fm_set, 1)
757 CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
758 CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
759 END DO
760 CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
761 CALL cp_fm_struct_release(tmp_fm_struct)
762 DO i = 1, 3
763 CALL cp_cfm_create(inv_mat(i, ispin), op_fm_set(1, ispin)%matrix_struct)
764 CALL cp_fm_create(inv_work(2, i, ispin), op_fm_set(2, ispin)%matrix_struct)
765 CALL cp_fm_create(inv_work(1, i, ispin), op_fm_set(1, ispin)%matrix_struct)
766 END DO
767 END DO
768
769 NULLIFY (cosmat, sinmat)
770 ALLOCATE (cosmat, sinmat)
771 CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
772 CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
773
774 DO i = 1, 3
775 kvec(:) = twopi*cell%h_inv(i, :)
776 CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
777
778 DO ispin = 1, dft_control%nspins ! spin
779 CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
780
781 CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(1, ispin), ncol=nmo)
782 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(1, ispin), 0.0_dp, &
783 op_fm_set(1, ispin))
784 CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(2, ispin), ncol=nmo)
785 CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(2, ispin), 0.0_dp, &
786 op_fm_set(2, ispin))
787
788 END DO
789
790 ! Second step invert C^T S_berry C
791 zdet = one
792 DO ispin = 1, dft_control%nspins
793 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
794 DO idim = 1, tmp_dim
795 eigrmat(ispin)%local_data(:, idim) = &
796 cmplx(op_fm_set(1, ispin)%local_data(:, idim), &
797 -op_fm_set(2, ispin)%local_data(:, idim), dp)
798 END DO
799 CALL cp_cfm_set_all(inv_mat(i, ispin), zero, one)
800 CALL cp_cfm_solve(eigrmat(ispin), inv_mat(i, ispin), zdeta)
801 END DO
802
803 ! Compute the derivative and add the result to mo_derivatives
804 DO ispin = 1, dft_control%nspins
805 CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
806 CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo, maxocc=maxocc)
807 DO z = 1, tmp_dim
808 inv_work(1, i, ispin)%local_data(:, z) = real(inv_mat(i, ispin)%local_data(:, z), dp)
809 inv_work(2, i, ispin)%local_data(:, z) = aimag(inv_mat(i, ispin)%local_data(:, z))
810 END DO
811 CALL parallel_gemm("N", "N", nao, nmo, nmo, -1.0_dp, opvec(1, ispin), inv_work(2, i, ispin), &
812 0.0_dp, dberry_psi0(i, ispin))
813 CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, opvec(2, ispin), inv_work(1, i, ispin), &
814 1.0_dp, dberry_psi0(i, ispin))
815 END DO
816 END DO !x/y/z-direction
817 !SL we omit here the multiplication with hmat (this scaling back done at the end of the response calc)
818
819 DO ispin = 1, dft_control%nspins
820 CALL cp_cfm_release(eigrmat(ispin))
821 DO i = 1, 3
822 CALL cp_cfm_release(inv_mat(i, ispin))
823 END DO
824 END DO
825 DEALLOCATE (inv_mat)
826 DEALLOCATE (eigrmat)
827
828 CALL cp_fm_release(inv_work)
829 CALL cp_fm_release(opvec)
830 CALL cp_fm_release(op_fm_set)
831
832 CALL dbcsr_deallocate_matrix(cosmat)
833 CALL dbcsr_deallocate_matrix(sinmat)
834
835 END IF ! do_raman
836
837 CALL timestop(handle)
838
839 END SUBROUTINE polar_operators_berry
840
841! **************************************************************************************************
842!> \brief Calculate the Berry phase operator in the AO basis and
843!> then the derivative of the Berry phase operator with respect to
844!> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
845!> afterwards multiply with the ground state MO coefficients
846!>
847!> \param qs_env ...
848!> \par History
849!> 01.2013 created [SL]
850!> 06.2018 polar_env integrated into qs_env (MK)
851!> 08.2020 adapt for xTB/DFTB (JHU)
852!> \author SL
853! **************************************************************************************************
854
855 SUBROUTINE polar_tb_operators_berry(qs_env)
856
857 TYPE(qs_environment_type), POINTER :: qs_env
858
859 CHARACTER(LEN=*), PARAMETER :: routinen = 'polar_tb_operators_berry'
860
861 COMPLEX(dp) :: zdeta
862 INTEGER :: handle, i, icol, idir, irow, ispin, nmo, &
863 nspins
864 LOGICAL :: do_raman, found
865 REAL(dp) :: dd, fdir
866 REAL(dp), DIMENSION(3) :: kvec, ria, rib
867 REAL(dp), DIMENSION(3, 3) :: hmat
868 REAL(dp), DIMENSION(:, :), POINTER :: d_block, s_block
869 TYPE(cell_type), POINTER :: cell
870 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0
871 TYPE(cp_fm_type), POINTER :: mo_coeff
872 TYPE(dbcsr_iterator_type) :: iter
873 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
874 TYPE(dft_control_type), POINTER :: dft_control
875 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
876 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
877 TYPE(polar_env_type), POINTER :: polar_env
878
879 CALL timeset(routinen, handle)
880
881 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
882 cell=cell, particle_set=particle_set, &
883 polar_env=polar_env, mos=mos, matrix_s=matrix_s)
884
885 nspins = dft_control%nspins
886
887 CALL get_polar_env(polar_env=polar_env, &
888 do_raman=do_raman, &
889 dberry_psi0=dberry_psi0)
890
891 IF (do_raman) THEN
892
893 ALLOCATE (dipmat(3))
894 DO i = 1, 3
895 ALLOCATE (dipmat(i)%matrix)
896 CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
897 CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
898 END DO
899
900 hmat = cell%hmat(:, :)/twopi
901
902 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
903 DO WHILE (dbcsr_iterator_blocks_left(iter))
904 NULLIFY (s_block, d_block)
905 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
906 ria = particle_set(irow)%r
907 rib = particle_set(icol)%r
908 DO idir = 1, 3
909 kvec(:) = twopi*cell%h_inv(idir, :)
910 dd = sum(kvec(:)*ria(:))
911 zdeta = cmplx(cos(dd), sin(dd), kind=dp)
912 fdir = aimag(log(zdeta))
913 dd = sum(kvec(:)*rib(:))
914 zdeta = cmplx(cos(dd), sin(dd), kind=dp)
915 fdir = fdir + aimag(log(zdeta))
916 CALL dbcsr_get_block_p(matrix=dipmat(idir)%matrix, &
917 row=irow, col=icol, block=d_block, found=found)
918 cpassert(found)
919 d_block = d_block + 0.5_dp*fdir*s_block
920 END DO
921 END DO
922 CALL dbcsr_iterator_stop(iter)
923
924 ! Compute the derivative and add the result to mo_derivatives
925 DO ispin = 1, dft_control%nspins ! spin
926 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
927 DO i = 1, 3
928 CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
929 dberry_psi0(i, ispin), ncol=nmo)
930 END DO !x/y/z-direction
931 END DO
932
933 DO i = 1, 3
934 CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
935 END DO
936 DEALLOCATE (dipmat)
937
938 END IF ! do_raman
939
940 CALL timestop(handle)
941 END SUBROUTINE polar_tb_operators_berry
942
943! **************************************************************************************************
944!> \brief Calculate the Berry phase operator in the AO basis and
945!> then the derivative of the Berry phase operator with respect to
946!> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
947!> afterwards multiply with the ground state MO coefficients
948!>
949!> \param qs_env ...
950!> \par History
951!> 01.2013 created [SL]
952!> 06.2018 polar_env integrated into qs_env (MK)
953!> \author SL
954! **************************************************************************************************
955 SUBROUTINE polar_operators_local(qs_env)
956
957 TYPE(qs_environment_type), POINTER :: qs_env
958
959 CHARACTER(LEN=*), PARAMETER :: routinen = 'polar_operators_local'
960
961 INTEGER :: handle, i, ispin, nmo, nspins
962 LOGICAL :: do_raman
963 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0
964 TYPE(cp_fm_type), POINTER :: mo_coeff
965 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
966 TYPE(dft_control_type), POINTER :: dft_control
967 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
968 TYPE(polar_env_type), POINTER :: polar_env
969
970 CALL timeset(routinen, handle)
971
972 CALL get_qs_env(qs_env=qs_env, &
973 dft_control=dft_control, &
974 polar_env=polar_env, &
975 mos=mos, &
976 matrix_s=matrix_s)
977
978 nspins = dft_control%nspins
979
980 CALL get_polar_env(polar_env=polar_env, &
981 do_raman=do_raman, &
982 dberry_psi0=dberry_psi0)
983
984 !SL calculate dipole berry phase
985 IF (do_raman) THEN
986
987 ALLOCATE (dipmat(3))
988 DO i = 1, 3
989 ALLOCATE (dipmat(i)%matrix)
990 CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
991 CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
992 END DO
993 CALL build_local_moment_matrix(qs_env, dipmat, 1)
994
995 ! Compute the derivative and add the result to mo_derivatives
996 DO ispin = 1, dft_control%nspins ! spin
997 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
998 DO i = 1, 3
999 CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
1000 dberry_psi0(i, ispin), ncol=nmo)
1001 END DO !x/y/z-direction
1002 END DO
1003
1004 DO i = 1, 3
1005 CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
1006 END DO
1007 DEALLOCATE (dipmat)
1008
1009 END IF ! do_raman
1010
1011 CALL timestop(handle)
1012
1013 END SUBROUTINE polar_operators_local
1014
1015 ! **************************************************************************************************
1016!> \brief Calculate the dipole operator referenced at the Wannier centers in the MO basis
1017!> \param qs_env ...
1018!> \param dcdr_env ...
1019!> \par History
1020!> 01.2013 created [SL]
1021!> 06.2018 polar_env integrated into qs_env (MK)
1022!> \authors Ravi Kumar
1023!> Rangsiman Ketkaew
1024! **************************************************************************************************
1025 SUBROUTINE polar_operators_local_wannier(qs_env, dcdr_env)
1026 TYPE(qs_environment_type), POINTER :: qs_env
1027 TYPE(dcdr_env_type) :: dcdr_env
1028
1029 CHARACTER(LEN=*), PARAMETER :: routinen = 'polar_operators_local_wannier'
1030
1031 INTEGER :: alpha, handle, i, icenter, ispin, &
1032 map_atom, map_molecule, &
1033 max_nbr_center, nao, natom, nmo, &
1034 nsubset
1035 INTEGER, ALLOCATABLE, DIMENSION(:) :: mapping_atom_molecule
1036 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: mapping_wannier_atom
1037 REAL(dp) :: f_spin, smallest_r, tmp_r
1038 REAL(dp), DIMENSION(3) :: distance, r_shifted
1039 REAL(dp), DIMENSION(:, :, :), POINTER :: apt_el, apt_nuc
1040 REAL(dp), DIMENSION(:, :, :, :), POINTER :: apt_center, apt_subset
1041 TYPE(cell_type), POINTER :: cell
1042 TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
1043 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0
1044 TYPE(cp_fm_type), POINTER :: mo_coeff, overlap1_mo, tmp_fm, &
1045 tmp_fm_like_mos, tmp_fm_momo
1046 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
1047 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1048 TYPE(polar_env_type), POINTER :: polar_env
1049 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1050
1051 CALL timeset(routinen, handle)
1052
1053 NULLIFY (qs_kind_set, particle_set, molecule_set, cell)
1054
1055 CALL get_qs_env(qs_env=qs_env, &
1056 qs_kind_set=qs_kind_set, &
1057 particle_set=particle_set, &
1058 molecule_set=molecule_set, &
1059 polar_env=polar_env, &
1060 cell=cell)
1061
1062 CALL get_polar_env(polar_env=polar_env, dberry_psi0=dberry_psi0)
1063
1064 nsubset = SIZE(molecule_set)
1065 natom = SIZE(particle_set)
1066 apt_el => dcdr_env%apt_el_dcdr
1067 apt_nuc => dcdr_env%apt_nuc_dcdr
1068 apt_subset => dcdr_env%apt_el_dcdr_per_subset
1069 apt_center => dcdr_env%apt_el_dcdr_per_center
1070
1071 ! Map wannier functions to atoms
1072 IF (dcdr_env%nspins == 1) THEN
1073 max_nbr_center = dcdr_env%nbr_center(1)
1074 ELSE
1075 max_nbr_center = max(dcdr_env%nbr_center(1), dcdr_env%nbr_center(2))
1076 END IF
1077 ALLOCATE (mapping_wannier_atom(max_nbr_center, dcdr_env%nspins))
1078 ALLOCATE (mapping_atom_molecule(natom))
1079 centers_set => dcdr_env%centers_set
1080 DO ispin = 1, dcdr_env%nspins
1081 DO icenter = 1, dcdr_env%nbr_center(ispin)
1082 ! For every center we check which atom is closest
1083 CALL shift_wannier_into_cell(r=centers_set(ispin)%array(1:3, icenter), &
1084 cell=cell, &
1085 r_shifted=r_shifted)
1086
1087 smallest_r = huge(0._dp)
1088 DO i = 1, natom
1089 distance = pbc(r_shifted, particle_set(i)%r(1:3), cell)
1090 tmp_r = sum(distance**2)
1091 IF (tmp_r < smallest_r) THEN
1092 mapping_wannier_atom(icenter, ispin) = i
1093 smallest_r = tmp_r
1094 END IF
1095 END DO
1096 END DO
1097
1098 ! Map atoms to molecules
1099 CALL molecule_of_atom(molecule_set, atom_to_mol=mapping_atom_molecule)
1100 IF (dcdr_env%lambda == 1 .AND. dcdr_env%beta == 1) THEN
1101 DO icenter = 1, dcdr_env%nbr_center(ispin)
1102 map_atom = mapping_wannier_atom(icenter, ispin)
1103 map_molecule = mapping_atom_molecule(map_atom)
1104 END DO
1105 END IF
1106 END DO !ispin
1107
1108 nao = dcdr_env%nao
1109 f_spin = 2._dp/dcdr_env%nspins
1110
1111 DO ispin = 1, dcdr_env%nspins
1112 ! Compute S^(1,R)_(ij)
1113
1114 ALLOCATE (tmp_fm_like_mos)
1115 ALLOCATE (overlap1_mo)
1116 CALL cp_fm_create(tmp_fm_like_mos, dcdr_env%likemos_fm_struct(ispin)%struct)
1117 CALL cp_fm_create(overlap1_mo, dcdr_env%momo_fm_struct(ispin)%struct)
1118 nmo = dcdr_env%nmo(ispin)
1119 mo_coeff => dcdr_env%mo_coeff(ispin)
1120 CALL cp_fm_set_all(tmp_fm_like_mos, 0.0_dp)
1121 CALL cp_fm_scale_and_add(0._dp, dcdr_env%dCR_prime(ispin), 1._dp, dcdr_env%dCR(ispin))
1122 ! CALL cp_dbcsr_sm_fm_multiply(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, mo_coeff, &
1123 ! tmp_fm_like_mos, ncol=nmo)
1124 CALL parallel_gemm("T", "N", nmo, nmo, nao, &
1125 1.0_dp, mo_coeff, tmp_fm_like_mos, &
1126 0.0_dp, overlap1_mo)
1127
1128 ! C^1 <- -dCR - 0.5 * mo_coeff @ S1_ij
1129 ! We get the negative of the coefficients out of the linres solver
1130 ! And apply the constant correction due to the overlap derivative.
1131 CALL parallel_gemm("N", "N", nao, nmo, nmo, &
1132 -0.5_dp, mo_coeff, overlap1_mo, &
1133 -1.0_dp, dcdr_env%dCR_prime(ispin))
1134 CALL cp_fm_release(overlap1_mo)
1135
1136 ! Allocate temporary matrices
1137 ALLOCATE (tmp_fm)
1138 ALLOCATE (tmp_fm_momo)
1139 CALL cp_fm_create(tmp_fm, dcdr_env%likemos_fm_struct(ispin)%struct)
1140 CALL cp_fm_create(tmp_fm_momo, dcdr_env%momo_fm_struct(ispin)%struct)
1141
1142 ! this_factor = -2._dp*f_spin
1143 DO alpha = 1, 3
1144 DO icenter = 1, dcdr_env%nbr_center(ispin)
1145 CALL dbcsr_set(dcdr_env%moments(alpha)%matrix, 0.0_dp)
1146 CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, &
1147 ref_point=centers_set(ispin)%array(1:3, icenter))
1148 CALL multiply_localization(ao_matrix=dcdr_env%moments(alpha)%matrix, &
1149 mo_coeff=mo_coeff, work=tmp_fm, nmo=nmo, &
1150 icenter=icenter, &
1151 res=dberry_psi0(alpha, ispin))
1152 END DO
1153
1154 END DO
1155
1156 CALL cp_fm_release(tmp_fm)
1157 CALL cp_fm_release(tmp_fm_like_mos)
1158 CALL cp_fm_release(tmp_fm_momo)
1159 DEALLOCATE (overlap1_mo)
1160 DEALLOCATE (tmp_fm)
1161 DEALLOCATE (tmp_fm_like_mos)
1162 DEALLOCATE (tmp_fm_momo)
1163 END DO !ispin
1164
1165 ! And deallocate all the things!
1166
1167 CALL timestop(handle)
1168 END SUBROUTINE polar_operators_local_wannier
1169
1170! **************************************************************************************************
1171!> \brief Calculate the local dipole operator in the AO basis
1172!> afterwards multiply with the ground state MO coefficients
1173!>
1174!> \param qs_env ...
1175!> \par History
1176!> 01.2013 created [SL]
1177!> 06.2018 polar_env integrated into qs_env (MK)
1178!> 08.2020 TB version (JHU)
1179!> \author SL
1180! **************************************************************************************************
1181 SUBROUTINE polar_tb_operators_local(qs_env)
1182
1183 TYPE(qs_environment_type), POINTER :: qs_env
1184
1185 CHARACTER(LEN=*), PARAMETER :: routinen = 'polar_tb_operators_local'
1186
1187 INTEGER :: handle, i, icol, irow, ispin, nmo, nspins
1188 LOGICAL :: do_raman, found
1189 REAL(dp) :: fdir
1190 REAL(dp), DIMENSION(3) :: ria, rib
1191 REAL(dp), DIMENSION(:, :), POINTER :: d_block, s_block
1192 TYPE(cell_type), POINTER :: cell
1193 TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dberry_psi0
1194 TYPE(cp_fm_type), POINTER :: mo_coeff
1195 TYPE(dbcsr_iterator_type) :: iter
1196 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
1197 TYPE(dft_control_type), POINTER :: dft_control
1198 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1199 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1200 TYPE(polar_env_type), POINTER :: polar_env
1201
1202 CALL timeset(routinen, handle)
1203
1204 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
1205 cell=cell, particle_set=particle_set, &
1206 polar_env=polar_env, mos=mos, matrix_s=matrix_s)
1207
1208 nspins = dft_control%nspins
1209
1210 CALL get_polar_env(polar_env=polar_env, &
1211 do_raman=do_raman, &
1212 dberry_psi0=dberry_psi0)
1213
1214 IF (do_raman) THEN
1215
1216 ALLOCATE (dipmat(3))
1217 DO i = 1, 3
1218 ALLOCATE (dipmat(i)%matrix)
1219 CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
1220 END DO
1221
1222 CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1223 DO WHILE (dbcsr_iterator_blocks_left(iter))
1224 NULLIFY (s_block, d_block)
1225 CALL dbcsr_iterator_next_block(iter, irow, icol, s_block)
1226 ria = particle_set(irow)%r
1227 ria = pbc(ria, cell)
1228 rib = particle_set(icol)%r
1229 rib = pbc(rib, cell)
1230 DO i = 1, 3
1231 CALL dbcsr_get_block_p(matrix=dipmat(i)%matrix, &
1232 row=irow, col=icol, block=d_block, found=found)
1233 cpassert(found)
1234 fdir = 0.5_dp*(ria(i) + rib(i))
1235 d_block = s_block*fdir
1236 END DO
1237 END DO
1238 CALL dbcsr_iterator_stop(iter)
1239
1240 ! Compute the derivative and add the result to mo_derivatives
1241 DO ispin = 1, dft_control%nspins ! spin
1242 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1243 DO i = 1, 3
1244 CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
1245 dberry_psi0(i, ispin), ncol=nmo)
1246 END DO !x/y/z-direction
1247 END DO
1248
1249 DO i = 1, 3
1250 CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
1251 END DO
1252 DEALLOCATE (dipmat)
1253
1254 END IF ! do_raman
1255
1256 CALL timestop(handle)
1257
1258 END SUBROUTINE polar_tb_operators_local
1259
1260! **************************************************************************************************
1261!> \brief ...
1262!> \param a ...
1263!> \param b ...
1264!> \param c ...
1265!> \return ...
1266! **************************************************************************************************
1267 FUNCTION fac_vecp(a, b, c) RESULT(factor)
1268
1269 INTEGER :: a, b, c
1270 REAL(dp) :: factor
1271
1272 factor = 0.0_dp
1273
1274 IF ((b == a + 1 .OR. b == a - 2) .AND. (c == b + 1 .OR. c == b - 2)) THEN
1275 factor = 1.0_dp
1276 ELSE IF ((b == a - 1 .OR. b == a + 2) .AND. (c == b - 1 .OR. c == b + 2)) THEN
1277 factor = -1.0_dp
1278 END IF
1279
1280 END FUNCTION fac_vecp
1281
1282! **************************************************************************************************
1283!> \brief ...
1284!> \param ii ...
1285!> \param iii ...
1286!> \return ...
1287! **************************************************************************************************
1288 FUNCTION ind_m2(ii, iii) RESULT(i)
1289
1290 INTEGER :: ii, iii, i
1291
1292 INTEGER :: l(3)
1293
1294 i = 0
1295 l(1:3) = 0
1296 IF (ii == 0) THEN
1297 l(iii) = 1
1298 ELSE IF (iii == 0) THEN
1299 l(ii) = 1
1300 ELSE IF (ii == iii) THEN
1301 l(ii) = 2
1302 i = coset(l(1), l(2), l(3)) - 1
1303 ELSE
1304 l(ii) = 1
1305 l(iii) = 1
1306 END IF
1307 i = coset(l(1), l(2), l(3)) - 1
1308 END FUNCTION ind_m2
1309
1310! **************************************************************************************************
1311!> \brief ...
1312!> \param i1 ...
1313!> \param i2 ...
1314!> \param i3 ...
1315! **************************************************************************************************
1316 SUBROUTINE set_vecp(i1, i2, i3)
1317
1318 INTEGER, INTENT(IN) :: i1
1319 INTEGER, INTENT(OUT) :: i2, i3
1320
1321 IF (i1 == 1) THEN
1322 i2 = 2
1323 i3 = 3
1324 ELSE IF (i1 == 2) THEN
1325 i2 = 3
1326 i3 = 1
1327 ELSE IF (i1 == 3) THEN
1328 i2 = 1
1329 i3 = 2
1330 ELSE
1331 END IF
1332
1333 END SUBROUTINE set_vecp
1334! **************************************************************************************************
1335!> \brief ...
1336!> \param i1 ...
1337!> \param i2 ...
1338!> \param i3 ...
1339! **************************************************************************************************
1340 SUBROUTINE set_vecp_rev(i1, i2, i3)
1341
1342 INTEGER, INTENT(IN) :: i1, i2
1343 INTEGER, INTENT(OUT) :: i3
1344
1345 IF ((i1 + i2) == 3) THEN
1346 i3 = 3
1347 ELSE IF ((i1 + i2) == 4) THEN
1348 i3 = 2
1349 ELSE IF ((i1 + i2) == 5) THEN
1350 i3 = 1
1351 ELSE
1352 END IF
1353
1354 END SUBROUTINE set_vecp_rev
1355
1356! **************************************************************************************************
1357!> \brief scale a matrix as a_ij = a_ij * pbc(rc(:,j),ra(:,i))(ixyz)
1358!> \param matrix ...
1359!> \param ra ...
1360!> \param rc ...
1361!> \param cell ...
1362!> \param ixyz ...
1363!> \author vw
1364! **************************************************************************************************
1365 SUBROUTINE fm_scale_by_pbc_ac(matrix, ra, rc, cell, ixyz)
1366 TYPE(cp_fm_type), INTENT(IN) :: matrix
1367 REAL(kind=dp), DIMENSION(:, :), INTENT(in) :: ra, rc
1368 TYPE(cell_type), POINTER :: cell
1369 INTEGER, INTENT(IN) :: ixyz
1370
1371 CHARACTER(LEN=*), PARAMETER :: routinen = 'fm_scale_by_pbc_AC'
1372
1373 INTEGER :: handle, icol_global, icol_local, &
1374 irow_global, irow_local, m, mypcol, &
1375 myprow, n, ncol_local, nrow_local
1376 REAL(kind=dp) :: dist(3), rra(3), rrc(3)
1377 REAL(kind=dp), DIMENSION(:, :), POINTER :: a
1378
1379 CALL timeset(routinen, handle)
1380
1381 myprow = matrix%matrix_struct%context%mepos(1)
1382 mypcol = matrix%matrix_struct%context%mepos(2)
1383
1384 nrow_local = matrix%matrix_struct%nrow_locals(myprow)
1385 ncol_local = matrix%matrix_struct%ncol_locals(mypcol)
1386
1387 n = SIZE(rc, 2)
1388 m = SIZE(ra, 2)
1389
1390 a => matrix%local_data
1391 DO icol_local = 1, ncol_local
1392 icol_global = matrix%matrix_struct%col_indices(icol_local)
1393 IF (icol_global > n) cycle
1394 rrc = rc(:, icol_global)
1395 DO irow_local = 1, nrow_local
1396 irow_global = matrix%matrix_struct%row_indices(irow_local)
1397 IF (irow_global > m) cycle
1398 rra = ra(:, irow_global)
1399 dist = pbc(rrc, rra, cell)
1400 a(irow_local, icol_local) = a(irow_local, icol_local)*dist(ixyz)
1401 END DO
1402 END DO
1403
1404 CALL timestop(handle)
1405
1406 END SUBROUTINE fm_scale_by_pbc_ac
1407
1408END MODULE qs_linres_op
cell_types::pbc
Definition cell_types.F:104
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:39
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
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_solve
subroutine, public cp_cfm_solve(matrix_a, general_a, determinant)
Solve the system of linear equations A*b=A_general using LU decomposition. Pay attention that both ma...
Definition cp_cfm_basic_linalg.F:731
cp_cfm_types
Represents a complex full matrix distributed on many processors.
Definition cp_cfm_types.F:12
cp_cfm_types::cp_cfm_release
subroutine, public cp_cfm_release(matrix)
Releases a full matrix.
Definition cp_cfm_types.F:182
cp_cfm_types::cp_cfm_create
subroutine, public cp_cfm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
Creates a new full matrix with the given structure.
Definition cp_cfm_types.F:125
cp_cfm_types::cp_cfm_get_info
subroutine, public cp_cfm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, matrix_struct, para_env)
Returns information about a full matrix.
Definition cp_cfm_types.F:631
cp_cfm_types::cp_cfm_set_all
subroutine, public cp_cfm_set_all(matrix, alpha, beta)
Set all elements of the full matrix to alpha. Besides, set all diagonal matrix elements to beta (if g...
Definition cp_cfm_types.F:202
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:944
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1050
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_get_block_p
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
Definition cp_dbcsr_api.F:693
cp_dbcsr_api::dbcsr_iterator_next_block
subroutine, public dbcsr_iterator_next_block(iterator, row, column, block, block_number_argument_has_been_removed, row_size, col_size, row_offset, col_offset, transposed)
...
Definition cp_dbcsr_api.F:971
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1012
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_checksum
real(kind=dp) function, public dbcsr_checksum(matrix, pos)
Calculates the checksum of a DBCSR matrix.
Definition cp_dbcsr_contrib.F:628
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:548
cp_fm_basic_linalg
Basic linear algebra operations for full matrices.
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_scale_and_add
subroutine, public cp_fm_scale_and_add(alpha, matrix_a, beta, matrix_b)
calc A <- alpha*A + beta*B optimized for alpha == 1.0 (just add beta*B) and beta == 0....
Definition cp_fm_basic_linalg.F:170
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::one
real(kind=dp), parameter, public one
Definition mathconstants.F:123
mathconstants::twopi
real(kind=dp), parameter, public twopi
Definition mathconstants.F:112
mathconstants::zero
real(kind=dp), parameter, public zero
Definition mathconstants.F:123
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::molecule_of_atom
subroutine, public molecule_of_atom(molecule_set, atom_to_mol)
finds for each atom the molecule it belongs to
Definition molecule_types.F:478
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::coset
integer, dimension(:, :, :), allocatable, public coset
Definition orbital_pointers.F:45
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_dcdr_utils
Calculate the derivatives of the MO coefficients wrt nuclear coordinates.
Definition qs_dcdr_utils.F:13
qs_dcdr_utils::multiply_localization
subroutine, public multiply_localization(ao_matrix, mo_coeff, work, nmo, icenter, res)
Multiply (ao_matrix @ mo_coeff) and store the column icenter in res.
Definition qs_dcdr_utils.F:113
qs_dcdr_utils::shift_wannier_into_cell
subroutine, public shift_wannier_into_cell(r, cell, r_shifted)
...
Definition qs_dcdr_utils.F:528
qs_elec_field
Distribution of the electric field gradient integral matrix.
Definition qs_elec_field.F:13
qs_elec_field::build_efg_matrix
subroutine, public build_efg_matrix(qs_env, matrix_efg, rc)
Calculation of the electric field gradient matrix over Cartesian Gaussian functions.
Definition qs_elec_field.F:75
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_fermi_contact
Distribution of the Fermi contact integral matrix.
Definition qs_fermi_contact.F:13
qs_fermi_contact::build_fermi_contact_matrix
subroutine, public build_fermi_contact_matrix(qs_env, matrix_fc, rc)
Calculation of the Fermi contact matrix over Cartesian Gaussian functions.
Definition qs_fermi_contact.F:75
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
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, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_linres_op
Calculate the operators p rxp and D needed in the optimization of the different contribution of the f...
Definition qs_linres_op.F:18
qs_linres_op::polar_operators_berry
subroutine, public polar_operators_berry(qs_env)
Calculate the Berry phase operator in the AO basis and then the derivative of the Berry phase operato...
Definition qs_linres_op.F:688
qs_linres_op::ind_m2
integer function, public ind_m2(ii, iii)
...
Definition qs_linres_op.F:1289
qs_linres_op::polar_operators_local_wannier
subroutine, public polar_operators_local_wannier(qs_env, dcdr_env)
Calculate the dipole operator referenced at the Wannier centers in the MO basis.
Definition qs_linres_op.F:1026
qs_linres_op::current_operators
subroutine, public current_operators(current_env, qs_env)
Calculate the first order hamiltonian applied to the ao and then apply them to the ground state orbit...
Definition qs_linres_op.F:126
qs_linres_op::set_vecp
subroutine, public set_vecp(i1, i2, i3)
...
Definition qs_linres_op.F:1317
qs_linres_op::fac_vecp
real(dp) function, public fac_vecp(a, b, c)
...
Definition qs_linres_op.F:1268
qs_linres_op::polar_operators
subroutine, public polar_operators(qs_env)
Calculate the dipole operator in the AO basis and its derivative wrt to MOs.
Definition qs_linres_op.F:649
qs_linres_op::polar_operators_local
subroutine, public polar_operators_local(qs_env)
Calculate the Berry phase operator in the AO basis and then the derivative of the Berry phase operato...
Definition qs_linres_op.F:956
qs_linres_op::fm_scale_by_pbc_ac
subroutine, public fm_scale_by_pbc_ac(matrix, ra, rc, cell, ixyz)
scale a matrix as a_ij = a_ij * pbc(rc(:,j),ra(:,i))(ixyz)
Definition qs_linres_op.F:1366
qs_linres_op::set_vecp_rev
subroutine, public set_vecp_rev(i1, i2, i3)
...
Definition qs_linres_op.F:1341
qs_linres_op::issc_operators
subroutine, public issc_operators(issc_env, qs_env, iatom)
...
Definition qs_linres_op.F:450
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_linres_types::get_current_env
subroutine, public get_current_env(current_env, simple_done, simple_converged, full_done, nao, nstates, gauge, list_cubes, statetrueindex, gauge_name, basisfun_center, nbr_center, center_list, centers_set, psi1_p, psi1_rxp, psi1_d, p_psi0, rxp_psi0, jrho1_atom_set, jrho1_set, chi_tensor, chi_tensor_loc, gauge_atom_radius, rs_gauge, use_old_gauge_atom, chi_pbc, psi0_order)
...
Definition qs_linres_types.F:429
qs_linres_types::get_issc_env
subroutine, public get_issc_env(issc_env, issc_on_atom_list, issc_gapw_radius, issc_loc, do_fc, do_sd, do_pso, do_dso, issc, interpolate_issc, psi1_efg, psi1_pso, psi1_dso, psi1_fc, efg_psi0, pso_psi0, dso_psi0, fc_psi0, matrix_efg, matrix_pso, matrix_dso, matrix_fc)
...
Definition qs_linres_types.F:561
qs_linres_types::get_polar_env
subroutine, public get_polar_env(polar_env, do_raman, do_periodic, dberry_psi0, polar, psi1_dberry, run_stopped)
...
Definition qs_linres_types.F:1242
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition qs_mo_types.F:401
qs_moments
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition qs_moments.F:14
qs_moments::build_local_moment_matrix
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type, all_images, minimum_image, neighbor_image, first_component)
...
Definition qs_moments.F:160
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:1045
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_operators_ao
Definition qs_operators_ao.F:13
qs_operators_ao::build_ang_mom_matrix
subroutine, public build_ang_mom_matrix(qs_env, matrix, rc)
Calculation of the angular momentum matrix over Cartesian Gaussian functions.
Definition qs_operators_ao.F:513
qs_operators_ao::build_lin_mom_matrix
subroutine, public build_lin_mom_matrix(qs_env, matrix, minimum_image)
Calculation of the linear momentum matrix <mu|∂|nu> over Cartesian Gaussian functions.
Definition qs_operators_ao.F:275
qs_spin_orbit
Distribution of the spin orbit integral matrix.
Definition qs_spin_orbit.F:13
qs_spin_orbit::build_pso_matrix
subroutine, public build_pso_matrix(qs_env, matrix_so, rc)
Calculation of the paramagnetic spin orbit matrix over Cartesian Gaussian functions.
Definition qs_spin_orbit.F:74
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_array_utils::cp_2d_i_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:231
cp_array_utils::cp_2d_r_p_type
represent a pointer to a 2d array
Definition cp_array_utils.F:109
cp_cfm_types::cp_cfm_type
Represent a complex full matrix.
Definition cp_cfm_types.F:69
cp_control_types::dft_control_type
Definition cp_control_types.F:758
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_iterator_type
Definition cp_dbcsr_api.F:188
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
molecule_types::molecule_type
Definition molecule_types.F:159
particle_types::particle_type
Definition particle_types.F:35
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_linres_types::current_env_type
Definition qs_linres_types.F:94
qs_linres_types::dcdr_env_type
Definition qs_linres_types.F:247
qs_linres_types::issc_env_type
Definition qs_linres_types.F:146
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_linres_types::polar_env_type
Definition qs_linres_types.F:136
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67