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