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