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