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