Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 : ! **************************************************************************************************
18 : MODULE 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_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
31 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply,&
32 : dbcsr_allocate_matrix_set,&
33 : dbcsr_deallocate_matrix_set
34 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add
35 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36 : cp_fm_struct_release,&
37 : cp_fm_struct_type
38 : USE cp_fm_types, ONLY: &
39 : cp_fm_create, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_indxl2g, cp_fm_release, &
40 : cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type
41 : USE cp_log_handling, ONLY: cp_get_default_logger,&
42 : cp_logger_type,&
43 : cp_to_string
44 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
45 : cp_print_key_unit_nr
46 : USE dbcsr_api, ONLY: &
47 : dbcsr_checksum, dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, &
48 : dbcsr_deallocate_matrix, dbcsr_distribution_type, dbcsr_get_block_p, &
49 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
50 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_set, dbcsr_type, &
51 : dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
52 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
53 : section_vals_type
54 : USE kinds, ONLY: dp
55 : USE mathconstants, ONLY: twopi
56 : USE message_passing, ONLY: mp_para_env_type
57 : USE orbital_pointers, ONLY: coset
58 : USE parallel_gemm_api, ONLY: parallel_gemm
59 : USE particle_methods, ONLY: get_particle_set
60 : USE particle_types, ONLY: particle_type
61 : USE qs_elec_field, ONLY: build_efg_matrix
62 : USE qs_environment_types, ONLY: get_qs_env,&
63 : qs_environment_type
64 : USE qs_fermi_contact, ONLY: build_fermi_contact_matrix
65 : USE qs_kind_types, ONLY: get_qs_kind_set,&
66 : qs_kind_type
67 : USE qs_linres_types, ONLY: current_env_type,&
68 : get_current_env,&
69 : get_issc_env,&
70 : get_polar_env,&
71 : issc_env_type,&
72 : linres_control_type,&
73 : polar_env_type
74 : USE qs_mo_types, ONLY: get_mo_set,&
75 : mo_set_type
76 : USE qs_moments, ONLY: build_berry_moment_matrix,&
77 : build_local_moment_matrix
78 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
79 : USE qs_operators_ao, ONLY: build_ang_mom_matrix,&
80 : build_lin_mom_matrix,&
81 : rRc_xyz_ao
82 : USE qs_spin_orbit, ONLY: build_pso_matrix
83 : #include "./base/base_uses.f90"
84 :
85 : IMPLICIT NONE
86 :
87 : PRIVATE
88 : PUBLIC :: current_operators, issc_operators, fac_vecp, ind_m2, set_vecp, set_vecp_rev, &
89 : fm_scale_by_pbc_AC, polar_operators
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_op'
92 :
93 : ! **************************************************************************************************
94 :
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief Calculate the first order hamiltonian applied to the ao
99 : !> and then apply them to the ground state orbitals,
100 : !> the h1_psi1 full matrices are then ready to solve the
101 : !> non-homogeneous linear equations that give the psi1
102 : !> linear response orbitals.
103 : !> \param current_env ...
104 : !> \param qs_env ...
105 : !> \par History
106 : !> 07.2005 created [MI]
107 : !> \author MI
108 : !> \note
109 : !> For the operators rxp and D the h1 depends on the psi0 to which
110 : !> is applied, or better the center of charge of the psi0 is
111 : !> used to define the position operator
112 : !> The centers of the orbitals result form the orbital localization procedure
113 : !> that typically uses the berry phase operator to define the Wannier centers.
114 : ! **************************************************************************************************
115 174 : SUBROUTINE current_operators(current_env, qs_env)
116 :
117 : TYPE(current_env_type) :: current_env
118 : TYPE(qs_environment_type), POINTER :: qs_env
119 :
120 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_operators'
121 :
122 : INTEGER :: handle, iao, icenter, idir, ii, iii, &
123 : ispin, istate, j, nao, natom, &
124 : nbr_center(2), nmo, nsgf, nspins, &
125 : nstates(2), output_unit
126 348 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
127 174 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
128 : REAL(dp) :: chk(3), ck(3), ckdk(3), dk(3)
129 348 : REAL(dp), DIMENSION(:, :), POINTER :: basisfun_center, vecbuf_c0
130 : TYPE(cell_type), POINTER :: cell
131 174 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
132 696 : TYPE(cp_2d_r_p_type), DIMENSION(3) :: vecbuf_RmdC0
133 174 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
134 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
135 : TYPE(cp_fm_type) :: fm_work1
136 696 : TYPE(cp_fm_type), DIMENSION(3) :: fm_Rmd_mos
137 174 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
138 174 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: p_psi0, rxp_psi0
139 : TYPE(cp_fm_type), POINTER :: mo_coeff
140 : TYPE(cp_logger_type), POINTER :: logger
141 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
142 174 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_ao
143 : TYPE(dft_control_type), POINTER :: dft_control
144 : TYPE(linres_control_type), POINTER :: linres_control
145 : TYPE(mp_para_env_type), POINTER :: para_env
146 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
147 174 : POINTER :: sab_all, sab_orb
148 174 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
149 174 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
150 : TYPE(section_vals_type), POINTER :: lr_section
151 :
152 174 : CALL timeset(routineN, handle)
153 :
154 174 : NULLIFY (qs_kind_set, cell, dft_control, linres_control, &
155 174 : logger, particle_set, lr_section, &
156 174 : basisfun_center, centers_set, center_list, p_psi0, &
157 174 : rxp_psi0, vecbuf_c0, psi0_order, &
158 174 : mo_coeff, op_ao, sab_all)
159 :
160 174 : logger => cp_get_default_logger()
161 : lr_section => section_vals_get_subs_vals(qs_env%input, &
162 174 : "PROPERTIES%LINRES")
163 :
164 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
165 174 : extension=".linresLog")
166 174 : IF (output_unit > 0) THEN
167 : WRITE (output_unit, FMT="(T2,A,/)") &
168 87 : "CURRENT| Calculation of the p and (r-d)xp operators applied to psi0"
169 : END IF
170 :
171 : CALL get_qs_env(qs_env=qs_env, &
172 : qs_kind_set=qs_kind_set, &
173 : cell=cell, &
174 : dft_control=dft_control, &
175 : linres_control=linres_control, &
176 : para_env=para_env, &
177 : particle_set=particle_set, &
178 : sab_all=sab_all, &
179 : sab_orb=sab_orb, &
180 174 : dbcsr_dist=dbcsr_dist)
181 :
182 174 : nspins = dft_control%nspins
183 :
184 : CALL get_current_env(current_env=current_env, nao=nao, centers_set=centers_set, &
185 : center_list=center_list, basisfun_center=basisfun_center, &
186 : nbr_center=nbr_center, p_psi0=p_psi0, rxp_psi0=rxp_psi0, &
187 : psi0_order=psi0_order, &
188 174 : nstates=nstates)
189 :
190 522 : ALLOCATE (vecbuf_c0(1, nao))
191 696 : DO idir = 1, 3
192 522 : NULLIFY (vecbuf_Rmdc0(idir)%array)
193 1740 : ALLOCATE (vecbuf_Rmdc0(idir)%array(1, nao))
194 : END DO
195 :
196 174 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, nsgf=nsgf)
197 :
198 174 : natom = SIZE(particle_set, 1)
199 522 : ALLOCATE (first_sgf(natom))
200 522 : ALLOCATE (last_sgf(natom))
201 :
202 : CALL get_particle_set(particle_set, qs_kind_set, &
203 : first_sgf=first_sgf, &
204 174 : last_sgf=last_sgf)
205 :
206 : ! Calculate the (r - dk)xp operator applied to psi0k
207 : ! One possible way to go is to use the distributive property of the vector product and calculatr
208 : ! (r-c)xp + (c-d)xp
209 : ! where c depends on the contracted functions and not on the states
210 : ! d is the center of a specific state and a loop over states is needed
211 : ! the second term can be added in a second moment as a correction
212 : ! notice: (r-c) and p are operators, whereas (c-d) is a multiplicative factor
213 :
214 : ! !First term: operator matrix elements
215 : ! CALL rmc_x_p_xyz_ao(op_rmd_ao,qs_env,minimum_image=.FALSE.)
216 : !************************************************************
217 : !
218 : ! Since many psi0 vector can have the same center, depending on how the center is selected,
219 : ! the (r - dk)xp operator matrix is computed Ncenter times,
220 : ! where Ncenter is the total number of different centers
221 : ! and each time it is multiplied by all the psi0 with center dk to get the rxp_psi0 matrix
222 :
223 : !
224 : ! prepare for allocation
225 522 : ALLOCATE (row_blk_sizes(natom))
226 174 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
227 : !
228 : !
229 174 : CALL dbcsr_allocate_matrix_set(op_ao, 3)
230 174 : ALLOCATE (op_ao(1)%matrix, op_ao(2)%matrix, op_ao(3)%matrix)
231 :
232 : CALL dbcsr_create(matrix=op_ao(1)%matrix, &
233 : name="op_ao", &
234 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
235 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
236 174 : nze=0, mutable_work=.TRUE.)
237 174 : CALL cp_dbcsr_alloc_block_from_nbl(op_ao(1)%matrix, sab_all)
238 174 : CALL dbcsr_set(op_ao(1)%matrix, 0.0_dp)
239 :
240 522 : DO idir = 2, 3
241 : CALL dbcsr_copy(op_ao(idir)%matrix, op_ao(1)%matrix, &
242 348 : "op_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
243 522 : CALL dbcsr_set(op_ao(idir)%matrix, 0.0_dp)
244 : END DO
245 :
246 174 : chk(:) = 0.0_dp
247 424 : DO ispin = 1, nspins
248 250 : mo_coeff => psi0_order(ispin)
249 250 : nmo = nstates(ispin)
250 250 : CALL cp_fm_set_all(p_psi0(ispin, 1), 0.0_dp)
251 250 : CALL cp_fm_set_all(p_psi0(ispin, 2), 0.0_dp)
252 250 : CALL cp_fm_set_all(p_psi0(ispin, 3), 0.0_dp)
253 1500 : DO icenter = 1, nbr_center(ispin)
254 1250 : CALL dbcsr_set(op_ao(1)%matrix, 0.0_dp)
255 1250 : CALL dbcsr_set(op_ao(2)%matrix, 0.0_dp)
256 1250 : CALL dbcsr_set(op_ao(3)%matrix, 0.0_dp)
257 : !CALL rmc_x_p_xyz_ao(op_ao,qs_env,minimum_image=.FALSE.,&
258 : ! & wancen=centers_set(ispin)%array(1:3,icenter))
259 : ! &
260 1250 : CALL build_ang_mom_matrix(qs_env, op_ao, centers_set(ispin)%array(1:3, icenter))
261 : !
262 : ! accumulate checksums
263 1250 : chk(1) = chk(1) + dbcsr_checksum(op_ao(1)%matrix)
264 1250 : chk(2) = chk(2) + dbcsr_checksum(op_ao(2)%matrix)
265 1250 : chk(3) = chk(3) + dbcsr_checksum(op_ao(3)%matrix)
266 5250 : DO idir = 1, 3
267 3750 : CALL cp_fm_set_all(rxp_psi0(ispin, idir), 0.0_dp)
268 : CALL cp_dbcsr_sm_fm_multiply(op_ao(idir)%matrix, mo_coeff, &
269 : rxp_psi0(ispin, idir), ncol=nmo, &
270 3750 : alpha=-1.0_dp)
271 9794 : DO j = center_list(ispin)%array(1, icenter), center_list(ispin)%array(1, icenter + 1) - 1
272 4794 : istate = center_list(ispin)%array(2, j)
273 : ! the p_psi0 fm is used as temporary matrix to store the results for the psi0 centered in dk
274 : CALL cp_fm_to_fm(rxp_psi0(ispin, idir), &
275 8544 : p_psi0(ispin, idir), 1, istate, istate)
276 : END DO
277 : END DO
278 : END DO
279 250 : CALL cp_fm_to_fm(p_psi0(ispin, 1), rxp_psi0(ispin, 1))
280 250 : CALL cp_fm_to_fm(p_psi0(ispin, 2), rxp_psi0(ispin, 2))
281 424 : CALL cp_fm_to_fm(p_psi0(ispin, 3), rxp_psi0(ispin, 3))
282 : END DO
283 : !
284 174 : CALL dbcsr_deallocate_matrix_set(op_ao)
285 : !
286 : ! print checksums
287 174 : IF (output_unit > 0) THEN
288 87 : WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum L_x =', chk(1)
289 87 : WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum L_y =', chk(2)
290 87 : WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum L_z =', chk(3)
291 : END IF
292 : !
293 : ! Calculate the px py pz operators
294 174 : CALL dbcsr_allocate_matrix_set(op_ao, 3)
295 174 : ALLOCATE (op_ao(1)%matrix, op_ao(2)%matrix, op_ao(3)%matrix)
296 :
297 : CALL dbcsr_create(matrix=op_ao(1)%matrix, &
298 : name="op_ao", &
299 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
300 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
301 174 : nze=0, mutable_work=.TRUE.)
302 174 : CALL cp_dbcsr_alloc_block_from_nbl(op_ao(1)%matrix, sab_orb)
303 174 : CALL dbcsr_set(op_ao(1)%matrix, 0.0_dp)
304 :
305 522 : DO idir = 2, 3
306 : CALL dbcsr_copy(op_ao(idir)%matrix, op_ao(1)%matrix, &
307 348 : "op_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
308 522 : CALL dbcsr_set(op_ao(idir)%matrix, 0.0_dp)
309 : END DO
310 : !
311 174 : CALL build_lin_mom_matrix(qs_env, op_ao)
312 : !CALL p_xyz_ao(op_ao,qs_env,minimum_image=.FALSE.)
313 : !
314 : ! print checksums
315 174 : chk(1) = dbcsr_checksum(op_ao(1)%matrix)
316 174 : chk(2) = dbcsr_checksum(op_ao(2)%matrix)
317 174 : chk(3) = dbcsr_checksum(op_ao(3)%matrix)
318 174 : IF (output_unit > 0) THEN
319 87 : WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum P_x =', chk(1)
320 87 : WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum P_y =', chk(2)
321 87 : WRITE (output_unit, '(T2,A,E23.16)') 'CURRENT| current_operators: CheckSum P_z =', chk(3)
322 : END IF
323 : ! Apply the p operator to the psi0
324 696 : DO idir = 1, 3
325 1446 : DO ispin = 1, nspins
326 750 : mo_coeff => psi0_order(ispin)
327 750 : nmo = nstates(ispin)
328 750 : CALL cp_fm_set_all(p_psi0(ispin, idir), 0.0_dp)
329 : CALL cp_dbcsr_sm_fm_multiply(op_ao(idir)%matrix, mo_coeff, &
330 : p_psi0(ispin, idir), ncol=nmo, &
331 1272 : alpha=-1.0_dp)
332 : END DO
333 : END DO
334 : !
335 174 : CALL dbcsr_deallocate_matrix_set(op_ao)
336 : !
337 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
338 174 : "PRINT%PROGRAM_RUN_INFO")
339 :
340 : ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
341 : ! This part is not necessary with the present implementation
342 : ! the angular momentum operator is computed directly for each dk independently
343 : ! and multiplied by the proper psi0 (i.e. those centered in dk)
344 : ! If Wannier centers are used, and no grouping of states with close centers is applied
345 : ! the (r-dk)xp operator is computed Nstate times and each time applied to only one vector psi0
346 : !
347 : ! Apply the (r-c)xp operator to the psi0
348 : !DO ispin = 1,nspins
349 : ! CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, nmo=nmo, homo=homo)
350 : ! DO idir = 1,3
351 : ! CALL cp_fm_set_all(rxp_psi0(ispin,idir),0.0_dp)
352 : ! CALL cp_sm_fm_multiply(op_rmd_ao(idir)%matrix,mo_coeff,&
353 : ! rxp_psi0(ispin,idir),ncol=nmo,alpha=-1.0_dp)
354 : ! END DO
355 : !END DO
356 :
357 : !Calculate the second term of the operator state by state
358 : !!!! what follows is a way to avoid calculating the L matrix for each centers.
359 : !!!! not tested
360 : IF (.FALSE.) THEN
361 : DO ispin = 1, nspins
362 : ! Allocate full matrices as working storage in the calculation
363 : ! of the rxp operator matrix. 3 matrices for the 3 Cartesian direction
364 : ! plus one to apply the momentum oprator to the modified mos fm
365 : mo_coeff => psi0_order(ispin)
366 : nmo = nstates(ispin)
367 : NULLIFY (tmp_fm_struct)
368 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
369 : ncol_global=nmo, para_env=para_env, &
370 : context=mo_coeff%matrix_struct%context)
371 : DO idir = 1, 3
372 : CALL cp_fm_create(fm_Rmd_mos(idir), tmp_fm_struct)
373 : END DO
374 : CALL cp_fm_create(fm_work1, tmp_fm_struct)
375 : CALL cp_fm_struct_release(tmp_fm_struct)
376 :
377 : ! This part should be done better, using the full matrix distribution
378 : DO istate = 1, nmo
379 : CALL cp_fm_get_submatrix(psi0_order(ispin), vecbuf_c0, 1, istate, nao, 1, &
380 : transpose=.TRUE.)
381 : !center of the localized psi0 state istate
382 : dk(1:3) = centers_set(ispin)%array(1:3, istate)
383 : DO idir = 1, 3
384 : ! This loop should be distributed over the processors
385 : DO iao = 1, nao
386 : ck(1:3) = basisfun_center(1:3, iao)
387 : ckdk = pbc(dk, ck, cell)
388 : vecbuf_Rmdc0(idir)%array(1, iao) = vecbuf_c0(1, iao)*ckdk(idir)
389 : END DO ! iao
390 : CALL cp_fm_set_submatrix(fm_Rmd_mos(idir), vecbuf_Rmdc0(idir)%array, &
391 : 1, istate, nao, 1, transpose=.TRUE.)
392 : END DO ! idir
393 : END DO ! istate
394 :
395 : DO idir = 1, 3
396 : CALL set_vecp(idir, ii, iii)
397 :
398 : !Add the second term to the idir component
399 : CALL cp_fm_set_all(fm_work1, 0.0_dp)
400 : CALL cp_dbcsr_sm_fm_multiply(op_ao(iii)%matrix, fm_Rmd_mos(ii), &
401 : fm_work1, ncol=nmo, alpha=-1.0_dp)
402 : CALL cp_fm_scale_and_add(1.0_dp, rxp_psi0(ispin, idir), &
403 : 1.0_dp, fm_work1)
404 :
405 : CALL cp_fm_set_all(fm_work1, 0.0_dp)
406 : CALL cp_dbcsr_sm_fm_multiply(op_ao(ii)%matrix, fm_Rmd_mos(iii), &
407 : fm_work1, ncol=nmo, alpha=-1.0_dp)
408 : CALL cp_fm_scale_and_add(1.0_dp, rxp_psi0(ispin, idir), &
409 : -1.0_dp, fm_work1)
410 :
411 : END DO ! idir
412 :
413 : DO idir = 1, 3
414 : CALL cp_fm_release(fm_Rmd_mos(idir))
415 : END DO
416 : CALL cp_fm_release(fm_work1)
417 :
418 : END DO ! ispin
419 : END IF
420 :
421 174 : DEALLOCATE (row_blk_sizes)
422 :
423 174 : DEALLOCATE (first_sgf, last_sgf)
424 :
425 174 : DEALLOCATE (vecbuf_c0)
426 696 : DO idir = 1, 3
427 696 : DEALLOCATE (vecbuf_Rmdc0(idir)%array)
428 : END DO
429 :
430 174 : CALL timestop(handle)
431 :
432 696 : END SUBROUTINE current_operators
433 :
434 : ! **************************************************************************************************
435 : !> \brief ...
436 : !> \param issc_env ...
437 : !> \param qs_env ...
438 : !> \param iatom ...
439 : ! **************************************************************************************************
440 44 : SUBROUTINE issc_operators(issc_env, qs_env, iatom)
441 :
442 : TYPE(issc_env_type) :: issc_env
443 : TYPE(qs_environment_type), POINTER :: qs_env
444 : INTEGER, INTENT(IN) :: iatom
445 :
446 : CHARACTER(LEN=*), PARAMETER :: routineN = 'issc_operators'
447 :
448 : INTEGER :: handle, idir, ispin, nmo, nspins, &
449 : output_unit
450 : LOGICAL :: do_dso, do_fc, do_pso, do_sd
451 : REAL(dp) :: chk(20), r_i(3)
452 : TYPE(cell_type), POINTER :: cell
453 44 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fc_psi0
454 44 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dso_psi0, efg_psi0, pso_psi0
455 : TYPE(cp_fm_type), POINTER :: mo_coeff
456 : TYPE(cp_logger_type), POINTER :: logger
457 44 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_dso, matrix_efg, matrix_fc, &
458 44 : matrix_pso
459 : TYPE(dft_control_type), POINTER :: dft_control
460 : TYPE(linres_control_type), POINTER :: linres_control
461 44 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
462 : TYPE(mp_para_env_type), POINTER :: para_env
463 44 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
464 44 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
465 : TYPE(section_vals_type), POINTER :: lr_section
466 :
467 44 : CALL timeset(routineN, handle)
468 :
469 44 : NULLIFY (matrix_fc, matrix_pso, matrix_efg)
470 44 : NULLIFY (efg_psi0, pso_psi0, fc_psi0)
471 :
472 44 : logger => cp_get_default_logger()
473 : lr_section => section_vals_get_subs_vals(qs_env%input, &
474 44 : "PROPERTIES%LINRES")
475 :
476 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
477 44 : extension=".linresLog")
478 :
479 : CALL get_qs_env(qs_env=qs_env, &
480 : qs_kind_set=qs_kind_set, &
481 : cell=cell, &
482 : dft_control=dft_control, &
483 : linres_control=linres_control, &
484 : para_env=para_env, &
485 : mos=mos, &
486 44 : particle_set=particle_set)
487 :
488 44 : nspins = dft_control%nspins
489 :
490 : CALL get_issc_env(issc_env=issc_env, &
491 : matrix_efg=matrix_efg, & !this is used only here alloc/dealloc here???
492 : matrix_pso=matrix_pso, & !this is used only here alloc/dealloc here???
493 : matrix_fc=matrix_fc, & !this is used only here alloc/dealloc here???
494 : matrix_dso=matrix_dso, & !this is used only here alloc/dealloc here???
495 : efg_psi0=efg_psi0, &
496 : pso_psi0=pso_psi0, &
497 : dso_psi0=dso_psi0, &
498 : fc_psi0=fc_psi0, &
499 : do_fc=do_fc, &
500 : do_sd=do_sd, &
501 : do_pso=do_pso, &
502 44 : do_dso=do_dso)
503 : !
504 : !
505 176 : r_i = particle_set(iatom)%r !pbc(particle_set(iatom)%r,cell)
506 : !write(*,*) 'issc_operators iatom=',iatom,' r_i=',r_i
507 44 : chk = 0.0_dp
508 : !
509 : !
510 : !
511 : ! Fermi contact integral
512 : !IF(do_fc) THEN
513 : IF (.TRUE.) THEN ! for the moment we build it (regs)
514 44 : CALL dbcsr_set(matrix_fc(1)%matrix, 0.0_dp)
515 44 : CALL build_fermi_contact_matrix(qs_env, matrix_fc, r_i)
516 :
517 44 : chk(1) = dbcsr_checksum(matrix_fc(1)%matrix)
518 :
519 44 : IF (output_unit > 0) THEN
520 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| fermi_contact: CheckSum =', chk(1)
521 : END IF
522 : END IF
523 : !
524 : ! spin-orbit integral
525 : !IF(do_pso) THEN
526 : IF (.TRUE.) THEN ! for the moment we build it (regs)
527 44 : CALL dbcsr_set(matrix_pso(1)%matrix, 0.0_dp)
528 44 : CALL dbcsr_set(matrix_pso(2)%matrix, 0.0_dp)
529 44 : CALL dbcsr_set(matrix_pso(3)%matrix, 0.0_dp)
530 44 : CALL build_pso_matrix(qs_env, matrix_pso, r_i)
531 :
532 44 : chk(2) = dbcsr_checksum(matrix_pso(1)%matrix)
533 44 : chk(3) = dbcsr_checksum(matrix_pso(2)%matrix)
534 44 : chk(4) = dbcsr_checksum(matrix_pso(3)%matrix)
535 :
536 44 : IF (output_unit > 0) THEN
537 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| pso_x: CheckSum =', chk(2)
538 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| pso_y: CheckSum =', chk(3)
539 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| pso_z: CheckSum =', chk(4)
540 : END IF
541 : END IF
542 : !
543 : ! electric field integral
544 : !IF(do_sd) THEN
545 : IF (.TRUE.) THEN ! for the moment we build it (regs)
546 44 : CALL dbcsr_set(matrix_efg(1)%matrix, 0.0_dp)
547 44 : CALL dbcsr_set(matrix_efg(2)%matrix, 0.0_dp)
548 44 : CALL dbcsr_set(matrix_efg(3)%matrix, 0.0_dp)
549 44 : CALL dbcsr_set(matrix_efg(4)%matrix, 0.0_dp)
550 44 : CALL dbcsr_set(matrix_efg(5)%matrix, 0.0_dp)
551 44 : CALL dbcsr_set(matrix_efg(6)%matrix, 0.0_dp)
552 44 : CALL build_efg_matrix(qs_env, matrix_efg, r_i)
553 :
554 44 : chk(5) = dbcsr_checksum(matrix_efg(1)%matrix)
555 44 : chk(6) = dbcsr_checksum(matrix_efg(2)%matrix)
556 44 : chk(7) = dbcsr_checksum(matrix_efg(3)%matrix)
557 44 : chk(8) = dbcsr_checksum(matrix_efg(4)%matrix)
558 44 : chk(9) = dbcsr_checksum(matrix_efg(5)%matrix)
559 44 : chk(10) = dbcsr_checksum(matrix_efg(6)%matrix)
560 :
561 44 : IF (output_unit > 0) THEN
562 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg (3xx-rr)/3: CheckSum =', chk(5)
563 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg (3yy-rr)/3: CheckSum =', chk(6)
564 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg (3zz-rr)/3: CheckSum =', chk(7)
565 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg xy: CheckSum =', chk(8)
566 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg xz: CheckSum =', chk(9)
567 22 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| efg yz: CheckSum =', chk(10)
568 : END IF
569 : END IF
570 : !
571 : !
572 44 : IF (output_unit > 0) THEN
573 242 : WRITE (output_unit, '(T2,A,E23.16)') 'ISSC| all operator: CheckSum =', SUM(chk(1:10))
574 : END IF
575 : !
576 : !>>> debugging only here we build the dipole matrix... debugging the kernel...
577 44 : IF (do_dso) THEN
578 2 : CALL dbcsr_set(matrix_dso(1)%matrix, 0.0_dp)
579 2 : CALL dbcsr_set(matrix_dso(2)%matrix, 0.0_dp)
580 2 : CALL dbcsr_set(matrix_dso(3)%matrix, 0.0_dp)
581 2 : CALL rRc_xyz_ao(matrix_dso, qs_env, (/0.0_dp, 0.0_dp, 0.0_dp/), 1)
582 : END IF
583 : !
584 : ! multiply by the mos
585 92 : DO ispin = 1, nspins
586 : !
587 48 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
588 48 : CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
589 : !
590 : ! EFG
591 48 : IF (do_sd) THEN
592 0 : DO idir = 1, 6
593 : CALL cp_dbcsr_sm_fm_multiply(matrix_efg(idir)%matrix, mo_coeff, &
594 : efg_psi0(ispin, idir), ncol=nmo, &
595 0 : alpha=1.0_dp)
596 : END DO
597 : END IF
598 : !
599 : ! PSO
600 48 : IF (do_pso) THEN
601 152 : DO idir = 1, 3
602 : CALL cp_dbcsr_sm_fm_multiply(matrix_pso(idir)%matrix, mo_coeff, &
603 : pso_psi0(ispin, idir), ncol=nmo, &
604 152 : alpha=-1.0_dp)
605 : END DO
606 : END IF
607 : !
608 : ! FC
609 48 : IF (do_fc) THEN
610 : CALL cp_dbcsr_sm_fm_multiply(matrix_fc(1)%matrix, mo_coeff, &
611 : fc_psi0(ispin), ncol=nmo, &
612 0 : alpha=1.0_dp)
613 : END IF
614 : !
615 : !>>> for debugging only
616 140 : IF (do_dso) THEN
617 8 : DO idir = 1, 3
618 : CALL cp_dbcsr_sm_fm_multiply(matrix_dso(idir)%matrix, mo_coeff, &
619 : dso_psi0(ispin, idir), ncol=nmo, &
620 8 : alpha=-1.0_dp)
621 : END DO
622 : END IF
623 : !<<< for debugging only
624 : END DO
625 :
626 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
627 44 : "PRINT%PROGRAM_RUN_INFO")
628 :
629 44 : CALL timestop(handle)
630 :
631 44 : END SUBROUTINE issc_operators
632 :
633 : ! **************************************************************************************************
634 : !> \brief Calculate the dipole operator in the AO basis and its derivative wrt to MOs
635 : !>
636 : !> \param qs_env ...
637 : ! **************************************************************************************************
638 108 : SUBROUTINE polar_operators(qs_env)
639 :
640 : TYPE(qs_environment_type), POINTER :: qs_env
641 :
642 : LOGICAL :: do_periodic
643 : TYPE(dft_control_type), POINTER :: dft_control
644 : TYPE(polar_env_type), POINTER :: polar_env
645 :
646 108 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, polar_env=polar_env)
647 108 : CALL get_polar_env(polar_env=polar_env, do_periodic=do_periodic)
648 108 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
649 8 : IF (do_periodic) THEN
650 4 : CALL polar_tb_operators_berry(qs_env)
651 : ELSE
652 4 : CALL polar_tb_operators_local(qs_env)
653 : END IF
654 : ELSE
655 100 : IF (do_periodic) THEN
656 14 : CALL polar_operators_berry(qs_env)
657 : ELSE
658 86 : CALL polar_operators_local(qs_env)
659 : END IF
660 : END IF
661 :
662 108 : END SUBROUTINE polar_operators
663 :
664 : ! **************************************************************************************************
665 : !> \brief Calculate the Berry phase operator in the AO basis and
666 : !> then the derivative of the Berry phase operator with respect to
667 : !> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
668 : !> afterwards multiply with the ground state MO coefficients
669 : !>
670 : !> \param qs_env ...
671 : !> \par History
672 : !> 01.2013 created [SL]
673 : !> 06.2018 polar_env integrated into qs_env (MK)
674 : !> \author SL
675 : ! **************************************************************************************************
676 :
677 14 : SUBROUTINE polar_operators_berry(qs_env)
678 :
679 : TYPE(qs_environment_type), POINTER :: qs_env
680 :
681 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_operators_berry'
682 : COMPLEX(KIND=dp), PARAMETER :: one = (1.0_dp, 0.0_dp), &
683 : zero = (0.0_dp, 0.0_dp)
684 :
685 : COMPLEX(DP) :: zdet, zdeta
686 : INTEGER :: handle, i, idim, ispin, nao, nmo, &
687 : nspins, tmp_dim, z
688 : LOGICAL :: do_raman
689 : REAL(dp) :: kvec(3), maxocc
690 : TYPE(cell_type), POINTER :: cell
691 14 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: eigrmat
692 14 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:, :) :: inv_mat
693 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
694 14 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: op_fm_set, opvec
695 14 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :, :) :: inv_work
696 14 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0
697 : TYPE(cp_fm_type), POINTER :: mo_coeff
698 14 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
699 : TYPE(dbcsr_type), POINTER :: cosmat, sinmat
700 : TYPE(dft_control_type), POINTER :: dft_control
701 14 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
702 : TYPE(mp_para_env_type), POINTER :: para_env
703 : TYPE(polar_env_type), POINTER :: polar_env
704 :
705 14 : CALL timeset(routineN, handle)
706 :
707 14 : NULLIFY (dBerry_psi0, sinmat, cosmat)
708 14 : NULLIFY (polar_env)
709 :
710 14 : NULLIFY (cell, dft_control, mos, matrix_s)
711 : CALL get_qs_env(qs_env=qs_env, &
712 : cell=cell, &
713 : dft_control=dft_control, &
714 : para_env=para_env, &
715 : polar_env=polar_env, &
716 : mos=mos, &
717 14 : matrix_s=matrix_s)
718 :
719 14 : nspins = dft_control%nspins
720 :
721 : CALL get_polar_env(polar_env=polar_env, &
722 : do_raman=do_raman, &
723 14 : dBerry_psi0=dBerry_psi0)
724 : !SL calculate dipole berry phase
725 14 : IF (do_raman) THEN
726 :
727 56 : DO i = 1, 3
728 98 : DO ispin = 1, nspins
729 84 : CALL cp_fm_set_all(dBerry_psi0(i, ispin), 0.0_dp)
730 : END DO
731 : END DO
732 :
733 : ! initialize all work matrices needed
734 84 : ALLOCATE (opvec(2, dft_control%nspins))
735 84 : ALLOCATE (op_fm_set(2, dft_control%nspins))
736 56 : ALLOCATE (eigrmat(dft_control%nspins))
737 98 : ALLOCATE (inv_mat(3, dft_control%nspins))
738 182 : ALLOCATE (inv_work(2, 3, dft_control%nspins))
739 :
740 : ! A bit to allocate for the wavefunction
741 28 : DO ispin = 1, dft_control%nspins
742 14 : NULLIFY (tmp_fm_struct, mo_coeff)
743 14 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
744 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmo, &
745 14 : ncol_global=nmo, para_env=para_env, context=mo_coeff%matrix_struct%context)
746 42 : DO i = 1, SIZE(op_fm_set, 1)
747 28 : CALL cp_fm_create(opvec(i, ispin), mo_coeff%matrix_struct)
748 42 : CALL cp_fm_create(op_fm_set(i, ispin), tmp_fm_struct)
749 : END DO
750 14 : CALL cp_cfm_create(eigrmat(ispin), op_fm_set(1, ispin)%matrix_struct)
751 14 : CALL cp_fm_struct_release(tmp_fm_struct)
752 84 : DO i = 1, 3
753 42 : CALL cp_cfm_create(inv_mat(i, ispin), op_fm_set(1, ispin)%matrix_struct)
754 42 : CALL cp_fm_create(inv_work(2, i, ispin), op_fm_set(2, ispin)%matrix_struct)
755 56 : CALL cp_fm_create(inv_work(1, i, ispin), op_fm_set(1, ispin)%matrix_struct)
756 : END DO
757 : END DO
758 :
759 14 : NULLIFY (cosmat, sinmat)
760 14 : ALLOCATE (cosmat, sinmat)
761 14 : CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
762 14 : CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
763 :
764 56 : DO i = 1, 3
765 168 : kvec(:) = twopi*cell%h_inv(i, :)
766 42 : CALL dbcsr_set(cosmat, 0.0_dp)
767 42 : CALL dbcsr_set(sinmat, 0.0_dp)
768 42 : CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
769 :
770 84 : DO ispin = 1, dft_control%nspins ! spin
771 42 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, mo_coeff=mo_coeff, nmo=nmo)
772 :
773 42 : CALL cp_dbcsr_sm_fm_multiply(cosmat, mo_coeff, opvec(1, ispin), ncol=nmo)
774 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(1, ispin), 0.0_dp, &
775 42 : op_fm_set(1, ispin))
776 42 : CALL cp_dbcsr_sm_fm_multiply(sinmat, mo_coeff, opvec(2, ispin), ncol=nmo)
777 : CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, mo_coeff, opvec(2, ispin), 0.0_dp, &
778 126 : op_fm_set(2, ispin))
779 :
780 : END DO
781 :
782 : ! Second step invert C^T S_berry C
783 42 : zdet = one
784 84 : DO ispin = 1, dft_control%nspins
785 42 : CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
786 210 : DO idim = 1, tmp_dim
787 : eigrmat(ispin)%local_data(:, idim) = &
788 : CMPLX(op_fm_set(1, ispin)%local_data(:, idim), &
789 546 : -op_fm_set(2, ispin)%local_data(:, idim), dp)
790 : END DO
791 42 : CALL cp_cfm_set_all(inv_mat(i, ispin), zero, one)
792 126 : CALL cp_cfm_solve(eigrmat(ispin), inv_mat(i, ispin), zdeta)
793 : END DO
794 :
795 : ! Compute the derivative and add the result to mo_derivatives
796 98 : DO ispin = 1, dft_control%nspins
797 42 : CALL cp_cfm_get_info(eigrmat(ispin), ncol_local=tmp_dim)
798 42 : CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo, maxocc=maxocc)
799 210 : DO z = 1, tmp_dim
800 504 : inv_work(1, i, ispin)%local_data(:, z) = REAL(inv_mat(i, ispin)%local_data(:, z), dp)
801 546 : inv_work(2, i, ispin)%local_data(:, z) = AIMAG(inv_mat(i, ispin)%local_data(:, z))
802 : END DO
803 : CALL parallel_gemm("N", "N", nao, nmo, nmo, -1.0_dp, opvec(1, ispin), inv_work(2, i, ispin), &
804 42 : 0.0_dp, dBerry_psi0(i, ispin))
805 : CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, opvec(2, ispin), inv_work(1, i, ispin), &
806 126 : 1.0_dp, dBerry_psi0(i, ispin))
807 : END DO
808 : END DO !x/y/z-direction
809 : !SL we omit here the multiplication with hmat (this scaling back done at the end of the response calc)
810 :
811 28 : DO ispin = 1, dft_control%nspins
812 14 : CALL cp_cfm_release(eigrmat(ispin))
813 70 : DO i = 1, 3
814 56 : CALL cp_cfm_release(inv_mat(i, ispin))
815 : END DO
816 : END DO
817 14 : DEALLOCATE (inv_mat)
818 14 : DEALLOCATE (eigrmat)
819 :
820 14 : CALL cp_fm_release(inv_work)
821 14 : CALL cp_fm_release(opvec)
822 14 : CALL cp_fm_release(op_fm_set)
823 :
824 14 : CALL dbcsr_deallocate_matrix(cosmat)
825 14 : CALL dbcsr_deallocate_matrix(sinmat)
826 :
827 : END IF ! do_raman
828 :
829 14 : CALL timestop(handle)
830 :
831 28 : END SUBROUTINE polar_operators_berry
832 :
833 : ! **************************************************************************************************
834 : !> \brief Calculate the Berry phase operator in the AO basis and
835 : !> then the derivative of the Berry phase operator with respect to
836 : !> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
837 : !> afterwards multiply with the ground state MO coefficients
838 : !>
839 : !> \param qs_env ...
840 : !> \par History
841 : !> 01.2013 created [SL]
842 : !> 06.2018 polar_env integrated into qs_env (MK)
843 : !> 08.2020 adapt for xTB/DFTB (JHU)
844 : !> \author SL
845 : ! **************************************************************************************************
846 :
847 4 : SUBROUTINE polar_tb_operators_berry(qs_env)
848 :
849 : TYPE(qs_environment_type), POINTER :: qs_env
850 :
851 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_tb_operators_berry'
852 :
853 : COMPLEX(dp) :: zdeta
854 : INTEGER :: blk, handle, i, icol, idir, irow, ispin, &
855 : nmo, nspins
856 : LOGICAL :: do_raman, found
857 : REAL(dp) :: dd, fdir
858 : REAL(dp), DIMENSION(3) :: kvec, ria, rib
859 : REAL(dp), DIMENSION(3, 3) :: hmat
860 4 : REAL(dp), DIMENSION(:, :), POINTER :: d_block, s_block
861 : TYPE(cell_type), POINTER :: cell
862 4 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0
863 : TYPE(cp_fm_type), POINTER :: mo_coeff
864 : TYPE(dbcsr_iterator_type) :: iter
865 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
866 : TYPE(dft_control_type), POINTER :: dft_control
867 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
868 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
869 : TYPE(polar_env_type), POINTER :: polar_env
870 :
871 4 : CALL timeset(routineN, handle)
872 :
873 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
874 : cell=cell, particle_set=particle_set, &
875 4 : polar_env=polar_env, mos=mos, matrix_s=matrix_s)
876 :
877 4 : nspins = dft_control%nspins
878 :
879 : CALL get_polar_env(polar_env=polar_env, &
880 : do_raman=do_raman, &
881 4 : dBerry_psi0=dBerry_psi0)
882 :
883 4 : IF (do_raman) THEN
884 :
885 16 : ALLOCATE (dipmat(3))
886 16 : DO i = 1, 3
887 12 : ALLOCATE (dipmat(i)%matrix)
888 12 : CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
889 16 : CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
890 : END DO
891 :
892 52 : hmat = cell%hmat(:, :)/twopi
893 :
894 4 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
895 20 : DO WHILE (dbcsr_iterator_blocks_left(iter))
896 16 : NULLIFY (s_block, d_block)
897 16 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
898 64 : ria = particle_set(irow)%r
899 64 : rib = particle_set(icol)%r
900 68 : DO idir = 1, 3
901 192 : kvec(:) = twopi*cell%h_inv(idir, :)
902 192 : dd = SUM(kvec(:)*ria(:))
903 48 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
904 48 : fdir = AIMAG(LOG(zdeta))
905 192 : dd = SUM(kvec(:)*rib(:))
906 48 : zdeta = CMPLX(COS(dd), SIN(dd), KIND=dp)
907 48 : fdir = fdir + AIMAG(LOG(zdeta))
908 : CALL dbcsr_get_block_p(matrix=dipmat(idir)%matrix, &
909 48 : row=irow, col=icol, BLOCK=d_block, found=found)
910 48 : CPASSERT(found)
911 1168 : d_block = d_block + 0.5_dp*fdir*s_block
912 : END DO
913 : END DO
914 4 : CALL dbcsr_iterator_stop(iter)
915 :
916 : ! Compute the derivative and add the result to mo_derivatives
917 8 : DO ispin = 1, dft_control%nspins ! spin
918 4 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
919 20 : DO i = 1, 3
920 : CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
921 16 : dBerry_psi0(i, ispin), ncol=nmo)
922 : END DO !x/y/z-direction
923 : END DO
924 :
925 16 : DO i = 1, 3
926 16 : CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
927 : END DO
928 4 : DEALLOCATE (dipmat)
929 :
930 : END IF ! do_raman
931 :
932 4 : CALL timestop(handle)
933 4 : END SUBROUTINE polar_tb_operators_berry
934 :
935 : ! **************************************************************************************************
936 : !> \brief Calculate the Berry phase operator in the AO basis and
937 : !> then the derivative of the Berry phase operator with respect to
938 : !> the ground state wave function (see paper Putrino et al., JCP, 13, 7102) for the AOs;
939 : !> afterwards multiply with the ground state MO coefficients
940 : !>
941 : !> \param qs_env ...
942 : !> \par History
943 : !> 01.2013 created [SL]
944 : !> 06.2018 polar_env integrated into qs_env (MK)
945 : !> \author SL
946 : ! **************************************************************************************************
947 86 : SUBROUTINE polar_operators_local(qs_env)
948 :
949 : TYPE(qs_environment_type), POINTER :: qs_env
950 :
951 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_operators_local'
952 :
953 : INTEGER :: handle, i, ispin, nmo, nspins
954 : LOGICAL :: do_raman
955 86 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0
956 : TYPE(cp_fm_type), POINTER :: mo_coeff
957 86 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
958 : TYPE(dft_control_type), POINTER :: dft_control
959 86 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
960 : TYPE(polar_env_type), POINTER :: polar_env
961 :
962 86 : CALL timeset(routineN, handle)
963 :
964 : CALL get_qs_env(qs_env=qs_env, &
965 : dft_control=dft_control, &
966 : polar_env=polar_env, &
967 : mos=mos, &
968 86 : matrix_s=matrix_s)
969 :
970 86 : nspins = dft_control%nspins
971 :
972 : CALL get_polar_env(polar_env=polar_env, &
973 : do_raman=do_raman, &
974 86 : dBerry_psi0=dBerry_psi0)
975 :
976 : !SL calculate dipole berry phase
977 86 : IF (do_raman) THEN
978 :
979 344 : ALLOCATE (dipmat(3))
980 344 : DO i = 1, 3
981 258 : ALLOCATE (dipmat(i)%matrix)
982 258 : CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
983 344 : CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
984 : END DO
985 86 : CALL build_local_moment_matrix(qs_env, dipmat, 1)
986 :
987 : ! Compute the derivative and add the result to mo_derivatives
988 178 : DO ispin = 1, dft_control%nspins ! spin
989 92 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
990 454 : DO i = 1, 3
991 : CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
992 368 : dBerry_psi0(i, ispin), ncol=nmo)
993 : END DO !x/y/z-direction
994 : END DO
995 :
996 344 : DO i = 1, 3
997 344 : CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
998 : END DO
999 86 : DEALLOCATE (dipmat)
1000 :
1001 : END IF ! do_raman
1002 :
1003 86 : CALL timestop(handle)
1004 :
1005 86 : END SUBROUTINE polar_operators_local
1006 :
1007 : ! **************************************************************************************************
1008 : !> \brief Calculate the local dipole operator in the AO basis
1009 : !> afterwards multiply with the ground state MO coefficients
1010 : !>
1011 : !> \param qs_env ...
1012 : !> \par History
1013 : !> 01.2013 created [SL]
1014 : !> 06.2018 polar_env integrated into qs_env (MK)
1015 : !> 08.2020 TB version (JHU)
1016 : !> \author SL
1017 : ! **************************************************************************************************
1018 4 : SUBROUTINE polar_tb_operators_local(qs_env)
1019 :
1020 : TYPE(qs_environment_type), POINTER :: qs_env
1021 :
1022 : CHARACTER(LEN=*), PARAMETER :: routineN = 'polar_tb_operators_local'
1023 :
1024 : INTEGER :: blk, handle, i, icol, irow, ispin, nmo, &
1025 : nspins
1026 : LOGICAL :: do_raman, found
1027 : REAL(dp) :: fdir
1028 : REAL(dp), DIMENSION(3) :: ria, rib
1029 4 : REAL(dp), DIMENSION(:, :), POINTER :: d_block, s_block
1030 : TYPE(cell_type), POINTER :: cell
1031 4 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: dBerry_psi0
1032 : TYPE(cp_fm_type), POINTER :: mo_coeff
1033 : TYPE(dbcsr_iterator_type) :: iter
1034 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat, matrix_s
1035 : TYPE(dft_control_type), POINTER :: dft_control
1036 4 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1037 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1038 : TYPE(polar_env_type), POINTER :: polar_env
1039 :
1040 4 : CALL timeset(routineN, handle)
1041 :
1042 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
1043 : cell=cell, particle_set=particle_set, &
1044 4 : polar_env=polar_env, mos=mos, matrix_s=matrix_s)
1045 :
1046 4 : nspins = dft_control%nspins
1047 :
1048 : CALL get_polar_env(polar_env=polar_env, &
1049 : do_raman=do_raman, &
1050 4 : dBerry_psi0=dBerry_psi0)
1051 :
1052 4 : IF (do_raman) THEN
1053 :
1054 16 : ALLOCATE (dipmat(3))
1055 16 : DO i = 1, 3
1056 12 : ALLOCATE (dipmat(i)%matrix)
1057 16 : CALL dbcsr_copy(dipmat(i)%matrix, matrix_s(1)%matrix, 'dipole')
1058 : END DO
1059 :
1060 4 : CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
1061 20 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1062 16 : NULLIFY (s_block, d_block)
1063 16 : CALL dbcsr_iterator_next_block(iter, irow, icol, s_block, blk)
1064 64 : ria = particle_set(irow)%r
1065 64 : ria = pbc(ria, cell)
1066 64 : rib = particle_set(icol)%r
1067 64 : rib = pbc(rib, cell)
1068 68 : DO i = 1, 3
1069 : CALL dbcsr_get_block_p(matrix=dipmat(i)%matrix, &
1070 48 : row=irow, col=icol, BLOCK=d_block, found=found)
1071 48 : CPASSERT(found)
1072 48 : fdir = 0.5_dp*(ria(i) + rib(i))
1073 1168 : d_block = s_block*fdir
1074 : END DO
1075 : END DO
1076 4 : CALL dbcsr_iterator_stop(iter)
1077 :
1078 : ! Compute the derivative and add the result to mo_derivatives
1079 8 : DO ispin = 1, dft_control%nspins ! spin
1080 4 : CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
1081 20 : DO i = 1, 3
1082 : CALL cp_dbcsr_sm_fm_multiply(dipmat(i)%matrix, mo_coeff, &
1083 16 : dBerry_psi0(i, ispin), ncol=nmo)
1084 : END DO !x/y/z-direction
1085 : END DO
1086 :
1087 16 : DO i = 1, 3
1088 16 : CALL dbcsr_deallocate_matrix(dipmat(i)%matrix)
1089 : END DO
1090 4 : DEALLOCATE (dipmat)
1091 :
1092 : END IF ! do_raman
1093 :
1094 4 : CALL timestop(handle)
1095 :
1096 4 : END SUBROUTINE polar_tb_operators_local
1097 :
1098 : ! **************************************************************************************************
1099 : !> \brief ...
1100 : !> \param a ...
1101 : !> \param b ...
1102 : !> \param c ...
1103 : !> \return ...
1104 : ! **************************************************************************************************
1105 7656 : FUNCTION fac_vecp(a, b, c) RESULT(factor)
1106 :
1107 : INTEGER :: a, b, c
1108 : REAL(dp) :: factor
1109 :
1110 7656 : factor = 0.0_dp
1111 :
1112 7656 : IF ((b .EQ. a + 1 .OR. b .EQ. a - 2) .AND. (c .EQ. b + 1 .OR. c .EQ. b - 2)) THEN
1113 : factor = 1.0_dp
1114 4215 : ELSEIF ((b .EQ. a - 1 .OR. b .EQ. a + 2) .AND. (c .EQ. b - 1 .OR. c .EQ. b + 2)) THEN
1115 4215 : factor = -1.0_dp
1116 : END IF
1117 :
1118 7656 : END FUNCTION fac_vecp
1119 :
1120 : ! **************************************************************************************************
1121 : !> \brief ...
1122 : !> \param ii ...
1123 : !> \param iii ...
1124 : !> \return ...
1125 : ! **************************************************************************************************
1126 414828 : FUNCTION ind_m2(ii, iii) RESULT(i)
1127 :
1128 : INTEGER :: ii, iii, i
1129 :
1130 : INTEGER :: l(3)
1131 :
1132 414828 : i = 0
1133 414828 : l(1:3) = 0
1134 414828 : IF (ii == 0) THEN
1135 0 : l(iii) = 1
1136 414828 : ELSEIF (iii == 0) THEN
1137 0 : l(ii) = 1
1138 414828 : ELSEIF (ii == iii) THEN
1139 138276 : l(ii) = 2
1140 138276 : i = coset(l(1), l(2), l(3)) - 1
1141 : ELSE
1142 276552 : l(ii) = 1
1143 276552 : l(iii) = 1
1144 : END IF
1145 414828 : i = coset(l(1), l(2), l(3)) - 1
1146 414828 : END FUNCTION ind_m2
1147 :
1148 : ! **************************************************************************************************
1149 : !> \brief ...
1150 : !> \param i1 ...
1151 : !> \param i2 ...
1152 : !> \param i3 ...
1153 : ! **************************************************************************************************
1154 44593 : SUBROUTINE set_vecp(i1, i2, i3)
1155 :
1156 : INTEGER, INTENT(IN) :: i1
1157 : INTEGER, INTENT(OUT) :: i2, i3
1158 :
1159 44593 : IF (i1 == 1) THEN
1160 14031 : i2 = 2
1161 14031 : i3 = 3
1162 30562 : ELSEIF (i1 == 2) THEN
1163 15281 : i2 = 3
1164 15281 : i3 = 1
1165 15281 : ELSEIF (i1 == 3) THEN
1166 15281 : i2 = 1
1167 15281 : i3 = 2
1168 : ELSE
1169 : END IF
1170 :
1171 44593 : END SUBROUTINE set_vecp
1172 : ! **************************************************************************************************
1173 : !> \brief ...
1174 : !> \param i1 ...
1175 : !> \param i2 ...
1176 : !> \param i3 ...
1177 : ! **************************************************************************************************
1178 7458 : SUBROUTINE set_vecp_rev(i1, i2, i3)
1179 :
1180 : INTEGER, INTENT(IN) :: i1, i2
1181 : INTEGER, INTENT(OUT) :: i3
1182 :
1183 7458 : IF ((i1 + i2) == 3) THEN
1184 2486 : i3 = 3
1185 4972 : ELSEIF ((i1 + i2) == 4) THEN
1186 2486 : i3 = 2
1187 2486 : ELSEIF ((i1 + i2) == 5) THEN
1188 2486 : i3 = 1
1189 : ELSE
1190 : END IF
1191 :
1192 7458 : END SUBROUTINE set_vecp_rev
1193 :
1194 : ! **************************************************************************************************
1195 : !> \brief scale a matrix as a_ij = a_ij * pbc(rc(:,j),ra(:,i))(ixyz)
1196 : !> \param matrix ...
1197 : !> \param ra ...
1198 : !> \param rc ...
1199 : !> \param cell ...
1200 : !> \param ixyz ...
1201 : !> \author vw
1202 : ! **************************************************************************************************
1203 1500 : SUBROUTINE fm_scale_by_pbc_AC(matrix, ra, rc, cell, ixyz)
1204 : TYPE(cp_fm_type), INTENT(IN) :: matrix
1205 : REAL(KIND=dp), DIMENSION(:, :), INTENT(in) :: ra, rc
1206 : TYPE(cell_type), POINTER :: cell
1207 : INTEGER, INTENT(IN) :: ixyz
1208 :
1209 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fm_scale_by_pbc_AC'
1210 :
1211 : INTEGER :: handle, icol_global, icol_local, irow_global, irow_local, m, mypcol, myprow, n, &
1212 : ncol_block, ncol_global, ncol_local, npcol, nprow, nrow_block, nrow_global, nrow_local
1213 : REAL(KIND=dp) :: dist(3), rra(3), rrc(3)
1214 1500 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: a
1215 :
1216 1500 : CALL timeset(routineN, handle)
1217 :
1218 1500 : myprow = matrix%matrix_struct%context%mepos(1)
1219 1500 : mypcol = matrix%matrix_struct%context%mepos(2)
1220 1500 : nprow = matrix%matrix_struct%context%num_pe(1)
1221 1500 : npcol = matrix%matrix_struct%context%num_pe(2)
1222 :
1223 1500 : nrow_block = matrix%matrix_struct%nrow_block
1224 1500 : ncol_block = matrix%matrix_struct%ncol_block
1225 1500 : nrow_global = matrix%matrix_struct%nrow_global
1226 1500 : ncol_global = matrix%matrix_struct%ncol_global
1227 1500 : nrow_local = matrix%matrix_struct%nrow_locals(myprow)
1228 1500 : ncol_local = matrix%matrix_struct%ncol_locals(mypcol)
1229 :
1230 1500 : n = SIZE(rc, 2)
1231 1500 : m = SIZE(ra, 2)
1232 :
1233 1500 : a => matrix%local_data
1234 11088 : DO icol_local = 1, ncol_local
1235 : icol_global = cp_fm_indxl2g(icol_local, ncol_block, mypcol, &
1236 9588 : matrix%matrix_struct%first_p_pos(2), npcol)
1237 9588 : IF (icol_global .GT. n) CYCLE
1238 38352 : rrc = rc(:, icol_global)
1239 131400 : DO irow_local = 1, nrow_local
1240 : irow_global = cp_fm_indxl2g(irow_local, nrow_block, myprow, &
1241 120312 : matrix%matrix_struct%first_p_pos(1), nprow)
1242 120312 : IF (irow_global .GT. m) CYCLE
1243 481248 : rra = ra(:, irow_global)
1244 120312 : dist = pbc(rrc, rra, cell)
1245 129900 : a(irow_local, icol_local) = a(irow_local, icol_local)*dist(ixyz)
1246 : END DO
1247 : END DO
1248 :
1249 1500 : CALL timestop(handle)
1250 :
1251 1500 : END SUBROUTINE fm_scale_by_pbc_AC
1252 :
1253 : END MODULE qs_linres_op
1254 :
|
