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 given the response wavefunctions obtained by the application
10 : !> of the (rxp), p, and ((dk-dl)xp) operators,
11 : !> here the current density vector (jx, jy, jz)
12 : !> is computed for the 3 directions of the magnetic field (Bx, By, Bz)
13 : !> \par History
14 : !> created 02-2006 [MI]
15 : !> \author MI
16 : ! **************************************************************************************************
17 : MODULE qs_linres_current
18 : USE ao_util, ONLY: exp_radius_very_extended
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_p_type,&
21 : gto_basis_set_type
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_array_utils, ONLY: cp_2d_i_p_type,&
25 : cp_2d_r_p_type
26 : USE cp_control_types, ONLY: dft_control_type
27 : USE cp_dbcsr_api, ONLY: &
28 : dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
29 : dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, dbcsr_p_type, dbcsr_put_block, &
30 : dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
31 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
32 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
33 : cp_dbcsr_sm_fm_multiply,&
34 : dbcsr_allocate_matrix_set,&
35 : dbcsr_deallocate_matrix_set
36 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
37 : cp_fm_trace
38 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
39 : cp_fm_struct_release,&
40 : cp_fm_struct_type
41 : USE cp_fm_types, ONLY: cp_fm_create,&
42 : cp_fm_release,&
43 : cp_fm_set_all,&
44 : cp_fm_to_fm,&
45 : cp_fm_type
46 : USE cp_log_handling, ONLY: cp_get_default_logger,&
47 : cp_logger_get_default_io_unit,&
48 : cp_logger_type,&
49 : cp_to_string
50 : USE cp_output_handling, ONLY: cp_p_file,&
51 : cp_print_key_finished_output,&
52 : cp_print_key_should_output,&
53 : cp_print_key_unit_nr
54 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
55 : USE cube_utils, ONLY: compute_cube_center,&
56 : cube_info_type,&
57 : return_cube
58 : USE gaussian_gridlevels, ONLY: gridlevel_info_type
59 : USE grid_api, ONLY: &
60 : GRID_FUNC_AB, GRID_FUNC_ADBmDAB_X, GRID_FUNC_ADBmDAB_Y, GRID_FUNC_ADBmDAB_Z, &
61 : GRID_FUNC_ARDBmDARB_XX, GRID_FUNC_ARDBmDARB_XY, GRID_FUNC_ARDBmDARB_XZ, &
62 : GRID_FUNC_ARDBmDARB_YX, GRID_FUNC_ARDBmDARB_YY, GRID_FUNC_ARDBmDARB_YZ, &
63 : GRID_FUNC_ARDBmDARB_ZX, GRID_FUNC_ARDBmDARB_ZY, GRID_FUNC_ARDBmDARB_ZZ, &
64 : collocate_pgf_product
65 : USE input_constants, ONLY: current_gauge_atom
66 : USE input_section_types, ONLY: section_get_ivals,&
67 : section_get_lval,&
68 : section_vals_get_subs_vals,&
69 : section_vals_type
70 : USE kinds, ONLY: default_path_length,&
71 : default_string_length,&
72 : dp,&
73 : int_8
74 : USE mathconstants, ONLY: twopi
75 : USE memory_utilities, ONLY: reallocate
76 : USE message_passing, ONLY: mp_para_env_type
77 : USE orbital_pointers, ONLY: ncoset
78 : USE particle_list_types, ONLY: particle_list_type
79 : USE particle_methods, ONLY: get_particle_set
80 : USE particle_types, ONLY: particle_type
81 : USE pw_env_types, ONLY: pw_env_get,&
82 : pw_env_type
83 : USE pw_methods, ONLY: pw_axpy,&
84 : pw_integrate_function,&
85 : pw_scale,&
86 : pw_zero
87 : USE pw_pool_types, ONLY: pw_pool_type
88 : USE pw_types, ONLY: pw_c1d_gs_type,&
89 : pw_r3d_rs_type
90 : USE qs_environment_types, ONLY: get_qs_env,&
91 : qs_environment_type
92 : USE qs_kind_types, ONLY: get_qs_kind,&
93 : get_qs_kind_set,&
94 : qs_kind_type
95 : USE qs_linres_atom_current, ONLY: calculate_jrho_atom,&
96 : calculate_jrho_atom_coeff,&
97 : calculate_jrho_atom_rad
98 : USE qs_linres_op, ONLY: fac_vecp,&
99 : fm_scale_by_pbc_AC,&
100 : ind_m2,&
101 : set_vecp,&
102 : set_vecp_rev
103 : USE qs_linres_types, ONLY: current_env_type,&
104 : get_current_env
105 : USE qs_matrix_pools, ONLY: qs_matrix_pools_type
106 : USE qs_mo_types, ONLY: get_mo_set,&
107 : mo_set_type
108 : USE qs_moments, ONLY: build_local_moments_der_matrix
109 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
110 : neighbor_list_iterate,&
111 : neighbor_list_iterator_create,&
112 : neighbor_list_iterator_p_type,&
113 : neighbor_list_iterator_release,&
114 : neighbor_list_set_p_type
115 : USE qs_operators_ao, ONLY: build_lin_mom_matrix
116 : USE qs_rho_types, ONLY: qs_rho_get
117 : USE qs_subsys_types, ONLY: qs_subsys_get,&
118 : qs_subsys_type
119 : USE realspace_grid_types, ONLY: realspace_grid_desc_p_type,&
120 : realspace_grid_desc_type,&
121 : realspace_grid_type,&
122 : rs_grid_create,&
123 : rs_grid_mult_and_add,&
124 : rs_grid_release,&
125 : rs_grid_zero
126 : USE rs_pw_interface, ONLY: density_rs2pw
127 : USE task_list_methods, ONLY: distribute_tasks,&
128 : rs_distribute_matrix,&
129 : task_list_inner_loop
130 : USE task_list_types, ONLY: atom_pair_type,&
131 : reallocate_tasks,&
132 : task_type
133 : #include "./base/base_uses.f90"
134 :
135 : IMPLICIT NONE
136 :
137 : PRIVATE
138 :
139 : ! *** Public subroutines ***
140 : PUBLIC :: current_build_current, current_build_chi, calculate_jrho_resp
141 :
142 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_current'
143 :
144 : TYPE box_type
145 : INTEGER :: n = -1
146 : REAL(dp), POINTER, DIMENSION(:, :) :: r => NULL()
147 : END TYPE box_type
148 : REAL(dp), DIMENSION(3, 3, 3), PARAMETER :: Levi_Civita = RESHAPE([ &
149 : 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
150 : 0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
151 : 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp], [3, 3, 3])
152 :
153 : CONTAINS
154 :
155 : ! **************************************************************************************************
156 : !> \brief First calculate the density matrixes, for each component of the current
157 : !> they are 3 because of the r dependent terms
158 : !> Next it collocates on the grid to have J(r)
159 : !> In the GAPW case one need to collocate on the PW grid only the soft part
160 : !> while the rest goes on Lebedev grids
161 : !> The contributions to the shift and to the susceptibility will be
162 : !> calculated separately and added only at the end
163 : !> The calculation of the shift tensor is performed on the position of the atoms
164 : !> and on other selected points in real space summing up the contributions
165 : !> from the PW grid current density and the local densities
166 : !> Spline interpolation is used
167 : !> \param current_env ...
168 : !> \param qs_env ...
169 : !> \param iB ...
170 : !> \author MI
171 : !> \note
172 : !> The susceptibility is needed to compute the G=0 term of the shift
173 : !> in reciprocal space. \chi_{ij} = \int (r x Jj)_i
174 : !> (where Jj id the current density generated by the field in direction j)
175 : !> To calculate the susceptibility on the PW grids it is necessary to apply
176 : !> the position operator yet another time.
177 : !> This cannot be done on directly on the full J(r) because it is not localized
178 : !> Therefore it is done state by state (see linres_nmr_shift)
179 : ! **************************************************************************************************
180 492 : SUBROUTINE current_build_current(current_env, qs_env, iB)
181 : !
182 : TYPE(current_env_type) :: current_env
183 : TYPE(qs_environment_type), POINTER :: qs_env
184 : INTEGER, INTENT(IN) :: iB
185 :
186 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_current'
187 :
188 : CHARACTER(LEN=default_path_length) :: ext, filename, my_pos
189 : INTEGER :: handle, idir, iiB, iiiB, ispin, istate, &
190 : j, jstate, nao, natom, nmo, nspins, &
191 : nstates(2), output_unit, unit_nr
192 492 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
193 492 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
194 : LOGICAL :: append_cube, gapw, mpi_io
195 : REAL(dp) :: dk(3), jrho_tot_G(3, 3), &
196 : jrho_tot_R(3, 3), maxocc, scale_fac
197 492 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ddk
198 : REAL(dp), EXTERNAL :: DDOT
199 : TYPE(cell_type), POINTER :: cell
200 492 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
201 492 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: p_psi1, psi1
202 492 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
203 492 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_D, psi1_p, psi1_rxp
204 : TYPE(cp_fm_type), POINTER :: mo_coeff
205 : TYPE(cp_logger_type), POINTER :: logger
206 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
207 492 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: density_matrix0, density_matrix_a, &
208 492 : density_matrix_ii, density_matrix_iii
209 : TYPE(dft_control_type), POINTER :: dft_control
210 492 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
211 : TYPE(mp_para_env_type), POINTER :: para_env
212 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
213 492 : POINTER :: sab_all
214 : TYPE(particle_list_type), POINTER :: particles
215 492 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
216 492 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: jrho1_g
217 : TYPE(pw_env_type), POINTER :: pw_env
218 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
219 : TYPE(pw_r3d_rs_type) :: wf_r
220 492 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: jrho1_r
221 492 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
222 : TYPE(qs_matrix_pools_type), POINTER :: mpools
223 : TYPE(qs_subsys_type), POINTER :: subsys
224 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
225 : TYPE(section_vals_type), POINTER :: current_section
226 :
227 492 : CALL timeset(routineN, handle)
228 : !
229 492 : NULLIFY (logger, current_section, density_matrix0, density_matrix_a, &
230 492 : density_matrix_ii, density_matrix_iii, cell, dft_control, mos, &
231 492 : particle_set, pw_env, auxbas_rs_desc, auxbas_pw_pool, &
232 492 : para_env, center_list, mo_coeff, jrho1_r, jrho1_g, &
233 492 : psi1_p, psi1_D, psi1_rxp, sab_all, qs_kind_set)
234 :
235 492 : logger => cp_get_default_logger()
236 492 : output_unit = cp_logger_get_default_io_unit(logger)
237 : !
238 : !
239 : CALL get_current_env(current_env=current_env, &
240 : center_list=center_list, &
241 : psi1_rxp=psi1_rxp, &
242 : psi1_D=psi1_D, &
243 : psi1_p=psi1_p, &
244 : psi0_order=psi0_order, &
245 : nstates=nstates, &
246 492 : nao=nao)
247 : !
248 : !
249 : CALL get_qs_env(qs_env=qs_env, &
250 : cell=cell, &
251 : dft_control=dft_control, &
252 : mos=mos, &
253 : mpools=mpools, &
254 : pw_env=pw_env, &
255 : para_env=para_env, &
256 : subsys=subsys, &
257 : sab_all=sab_all, &
258 : particle_set=particle_set, &
259 : qs_kind_set=qs_kind_set, &
260 492 : dbcsr_dist=dbcsr_dist)
261 :
262 492 : CALL qs_subsys_get(subsys, particles=particles)
263 :
264 492 : gapw = dft_control%qs_control%gapw
265 492 : nspins = dft_control%nspins
266 492 : natom = SIZE(particle_set, 1)
267 : !
268 : ! allocate temporary arrays
269 3396 : ALLOCATE (psi1(nspins), p_psi1(nspins))
270 1206 : DO ispin = 1, nspins
271 714 : CALL cp_fm_create(psi1(ispin), psi0_order(ispin)%matrix_struct)
272 714 : CALL cp_fm_create(p_psi1(ispin), psi0_order(ispin)%matrix_struct)
273 714 : CALL cp_fm_set_all(psi1(ispin), 0.0_dp)
274 1206 : CALL cp_fm_set_all(p_psi1(ispin), 0.0_dp)
275 : END DO
276 : !
277 : !
278 492 : CALL dbcsr_allocate_matrix_set(density_matrix0, nspins)
279 492 : CALL dbcsr_allocate_matrix_set(density_matrix_a, nspins)
280 492 : CALL dbcsr_allocate_matrix_set(density_matrix_ii, nspins)
281 492 : CALL dbcsr_allocate_matrix_set(density_matrix_iii, nspins)
282 : !
283 : ! prepare for allocation
284 1476 : ALLOCATE (first_sgf(natom))
285 984 : ALLOCATE (last_sgf(natom))
286 : CALL get_particle_set(particle_set, qs_kind_set, &
287 : first_sgf=first_sgf, &
288 492 : last_sgf=last_sgf)
289 984 : ALLOCATE (row_blk_sizes(natom))
290 492 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
291 492 : DEALLOCATE (first_sgf)
292 492 : DEALLOCATE (last_sgf)
293 : !
294 : !
295 1206 : DO ispin = 1, nspins
296 : !
297 : !density_matrix0
298 714 : ALLOCATE (density_matrix0(ispin)%matrix)
299 : CALL dbcsr_create(matrix=density_matrix0(ispin)%matrix, &
300 : name="density_matrix0", &
301 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
302 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
303 714 : mutable_work=.TRUE.)
304 714 : CALL cp_dbcsr_alloc_block_from_nbl(density_matrix0(ispin)%matrix, sab_all)
305 : !
306 : !density_matrix_a
307 714 : ALLOCATE (density_matrix_a(ispin)%matrix)
308 : CALL dbcsr_copy(density_matrix_a(ispin)%matrix, density_matrix0(ispin)%matrix, &
309 714 : name="density_matrix_a")
310 : !
311 : !density_matrix_ii
312 714 : ALLOCATE (density_matrix_ii(ispin)%matrix)
313 : CALL dbcsr_copy(density_matrix_ii(ispin)%matrix, density_matrix0(ispin)%matrix, &
314 714 : name="density_matrix_ii")
315 : !
316 : !density_matrix_iii
317 714 : ALLOCATE (density_matrix_iii(ispin)%matrix)
318 : CALL dbcsr_copy(density_matrix_iii(ispin)%matrix, density_matrix0(ispin)%matrix, &
319 1206 : name="density_matrix_iii")
320 : END DO
321 : !
322 492 : DEALLOCATE (row_blk_sizes)
323 : !
324 : !
325 492 : current_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%CURRENT")
326 : !
327 : !
328 492 : jrho_tot_G = 0.0_dp
329 492 : jrho_tot_R = 0.0_dp
330 : !
331 : ! Lets go!
332 492 : CALL set_vecp(iB, iiB, iiiB)
333 1206 : DO ispin = 1, nspins
334 714 : nmo = nstates(ispin)
335 714 : mo_coeff => psi0_order(ispin)
336 : !maxocc = max_occ(ispin)
337 : !
338 714 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
339 : !
340 : !
341 : ! Build the first density matrix
342 714 : CALL dbcsr_set(density_matrix0(ispin)%matrix, 0.0_dp)
343 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix0(ispin)%matrix, &
344 : matrix_v=mo_coeff, matrix_g=mo_coeff, &
345 714 : ncol=nmo, alpha=maxocc)
346 : !
347 : ! Allocate buffer vectors
348 2142 : ALLOCATE (ddk(3, nmo))
349 : !
350 : ! Construct the 3 density matrices for the field in direction iB
351 : !
352 : ! First the full matrix psi_a_iB
353 : ASSOCIATE (psi_a_iB => psi1(ispin), psi_buf => p_psi1(ispin))
354 714 : CALL cp_fm_set_all(psi_a_iB, 0.0_dp)
355 714 : CALL cp_fm_set_all(psi_buf, 0.0_dp)
356 : ! psi_a_iB = - (R_\nu-dk)_ii psi1_piiiB + (R_\nu-dk)_iii psi1_piiB
357 : !
358 : ! contributions from the response psi1_p_ii and psi1_p_iii
359 4296 : DO istate = 1, current_env%nbr_center(ispin)
360 14328 : dk(1:3) = current_env%centers_set(ispin)%array(1:3, istate)
361 : !
362 : ! Copy the vector in the full matrix psi1
363 : !nstate_loc = center_list(ispin)%array(1,icenter+1)-center_list(ispin)%array(1,icenter)
364 8922 : DO j = center_list(ispin)%array(1, istate), center_list(ispin)%array(1, istate + 1) - 1
365 4626 : jstate = center_list(ispin)%array(2, j)
366 4626 : CALL cp_fm_to_fm(psi1_p(ispin, iiB), psi_a_iB, 1, jstate, jstate)
367 4626 : CALL cp_fm_to_fm(psi1_p(ispin, iiiB), psi_buf, 1, jstate, jstate)
368 22086 : ddk(:, jstate) = dk(1:3)
369 : END DO
370 : END DO ! istate
371 714 : CALL fm_scale_by_pbc_AC(psi_a_iB, current_env%basisfun_center, ddk, cell, iiiB)
372 714 : CALL fm_scale_by_pbc_AC(psi_buf, current_env%basisfun_center, ddk, cell, iiB)
373 714 : CALL cp_fm_scale_and_add(-1.0_dp, psi_a_iB, 1.0_dp, psi_buf)
374 : !
375 : !psi_a_iB = psi_a_iB + psi1_rxp
376 : !
377 : ! contribution from the response psi1_rxp
378 714 : CALL cp_fm_scale_and_add(-1.0_dp, psi_a_iB, 1.0_dp, psi1_rxp(ispin, iB))
379 : !
380 : !psi_a_iB = psi_a_iB - psi1_D
381 714 : IF (current_env%full) THEN
382 : !
383 : ! contribution from the response psi1_D
384 582 : CALL cp_fm_scale_and_add(1.0_dp, psi_a_iB, -1.0_dp, psi1_D(ispin, iB))
385 : END IF
386 : !
387 : ! Multiply by the occupation number for the density matrix
388 : !
389 : ! Build the first density matrix
390 714 : CALL dbcsr_set(density_matrix_a(ispin)%matrix, 0.0_dp)
391 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_a(ispin)%matrix, &
392 : matrix_v=mo_coeff, matrix_g=psi_a_iB, &
393 1428 : ncol=nmo, alpha=maxocc)
394 : END ASSOCIATE
395 : !
396 : ! Build the second density matrix
397 714 : CALL dbcsr_set(density_matrix_iii(ispin)%matrix, 0.0_dp)
398 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_iii(ispin)%matrix, &
399 : matrix_v=mo_coeff, matrix_g=psi1_p(ispin, iiiB), &
400 714 : ncol=nmo, alpha=maxocc)
401 : !
402 : ! Build the third density matrix
403 714 : CALL dbcsr_set(density_matrix_ii(ispin)%matrix, 0.0_dp)
404 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_ii(ispin)%matrix, &
405 : matrix_v=mo_coeff, matrix_g=psi1_p(ispin, iiB), &
406 714 : ncol=nmo, alpha=maxocc)
407 2856 : DO idir = 1, 3
408 : !
409 : ! Calculate the current density on the pw grid (only soft if GAPW)
410 : ! idir is the cartesian component of the response current density
411 : ! generated by the magnetic field pointing in cartesian direction iB
412 : ! Use the qs_rho_type already used for rho during the scf
413 2142 : CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_r=jrho1_r)
414 2142 : CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_g=jrho1_g)
415 : ASSOCIATE (jrho_rspace => jrho1_r(ispin), jrho_gspace => jrho1_g(ispin))
416 2142 : CALL pw_zero(jrho_rspace)
417 2142 : CALL pw_zero(jrho_gspace)
418 : CALL calculate_jrho_resp(density_matrix0(ispin)%matrix, &
419 : density_matrix_a(ispin)%matrix, &
420 : density_matrix_ii(ispin)%matrix, &
421 : density_matrix_iii(ispin)%matrix, &
422 : iB, idir, jrho_rspace, jrho_gspace, qs_env, &
423 2142 : current_env, gapw)
424 :
425 2142 : scale_fac = cell%deth/twopi
426 2142 : CALL pw_scale(jrho_rspace, scale_fac)
427 2142 : CALL pw_scale(jrho_gspace, scale_fac)
428 :
429 2142 : jrho_tot_G(idir, iB) = pw_integrate_function(jrho_gspace, isign=-1)
430 4284 : jrho_tot_R(idir, iB) = pw_integrate_function(jrho_rspace, isign=-1)
431 : END ASSOCIATE
432 :
433 2856 : IF (output_unit > 0) THEN
434 : WRITE (output_unit, '(T2,2(A,E24.16))') 'Integrated j_'&
435 1071 : &//ACHAR(idir + 119)//ACHAR(iB + 119)//'(r): G-space=', &
436 2142 : jrho_tot_G(idir, iB), ' R-space=', jrho_tot_R(idir, iB)
437 : END IF
438 : !
439 : END DO ! idir
440 : !
441 : ! Deallocate buffer vectors
442 1920 : DEALLOCATE (ddk)
443 : !
444 : END DO ! ispin
445 :
446 492 : IF (gapw) THEN
447 1152 : DO idir = 1, 3
448 : !
449 : ! compute the atomic response current densities on the spherical grids
450 : ! First the sparse matrices are multiplied by the expansion coefficients
451 : ! this is the equivalent of the CPC for the charge density
452 : CALL calculate_jrho_atom_coeff(qs_env, current_env, &
453 : density_matrix0, &
454 : density_matrix_a, &
455 : density_matrix_ii, &
456 : density_matrix_iii, &
457 864 : iB, idir)
458 : !
459 : ! Then the radial parts are computed on the local radial grid, atom by atom
460 : ! 8 functions are computed for each atom, per grid point
461 : ! and per LM angular momentum. The multiplication by the Clebsh-Gordon
462 : ! coefficients or they correspondent for the derivatives, is also done here
463 864 : CALL calculate_jrho_atom_rad(qs_env, current_env, idir)
464 : !
465 : ! The current on the atomic grids
466 1152 : CALL calculate_jrho_atom(current_env, qs_env, iB, idir)
467 : END DO ! idir
468 : END IF
469 : !
470 : ! Cube files
471 492 : IF (BTEST(cp_print_key_should_output(logger%iter_info, current_section,&
472 : & "PRINT%CURRENT_CUBES"), cp_p_file)) THEN
473 0 : append_cube = section_get_lval(current_section, "PRINT%CURRENT_CUBES%APPEND")
474 0 : my_pos = "REWIND"
475 0 : IF (append_cube) THEN
476 0 : my_pos = "APPEND"
477 : END IF
478 : !
479 : CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, &
480 0 : auxbas_pw_pool=auxbas_pw_pool)
481 : !
482 0 : CALL auxbas_pw_pool%create_pw(wf_r)
483 : !
484 0 : DO idir = 1, 3
485 0 : CALL pw_zero(wf_r)
486 0 : CALL qs_rho_get(current_env%jrho1_set(idir)%rho, rho_r=jrho1_r)
487 0 : DO ispin = 1, nspins
488 0 : CALL pw_axpy(jrho1_r(ispin), wf_r, 1.0_dp)
489 : END DO
490 : !
491 0 : IF (gapw) THEN
492 : ! Add the local hard and soft contributions
493 : ! This can be done atom by atom by a spline extrapolation of the values
494 : ! on the spherical grid to the grid points.
495 0 : CPABORT("GAPW needs to be finalized")
496 : END IF
497 0 : filename = "jresp"
498 0 : mpi_io = .TRUE.
499 0 : WRITE (ext, '(a2,I1,a2,I1,a5)') "iB", iB, "_d", idir, ".cube"
500 0 : WRITE (ext, '(a2,a1,a2,a1,a5)') "iB", ACHAR(iB + 119), "_d", ACHAR(idir + 119), ".cube"
501 : unit_nr = cp_print_key_unit_nr(logger, current_section, "PRINT%CURRENT_CUBES", &
502 : extension=TRIM(ext), middle_name=TRIM(filename), &
503 : log_filename=.FALSE., file_position=my_pos, &
504 0 : mpi_io=mpi_io)
505 :
506 : CALL cp_pw_to_cube(wf_r, unit_nr, "RESPONSE CURRENT DENSITY ", &
507 : particles=particles, &
508 : stride=section_get_ivals(current_section, "PRINT%CURRENT_CUBES%STRIDE"), &
509 0 : mpi_io=mpi_io)
510 : CALL cp_print_key_finished_output(unit_nr, logger, current_section,&
511 0 : & "PRINT%CURRENT_CUBES", mpi_io=mpi_io)
512 : END DO
513 : !
514 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
515 : END IF ! current cube
516 : !
517 : ! Integrated current response checksum
518 492 : IF (output_unit > 0) THEN
519 246 : WRITE (output_unit, '(T2,A,E24.16)') 'CheckSum R-integrated j=', &
520 492 : SQRT(DDOT(9, jrho_tot_R(1, 1), 1, jrho_tot_R(1, 1), 1))
521 : END IF
522 : !
523 : !
524 : ! Dellocate grids for the calculation of jrho and the shift
525 492 : CALL cp_fm_release(psi1)
526 492 : CALL cp_fm_release(p_psi1)
527 :
528 492 : CALL dbcsr_deallocate_matrix_set(density_matrix0)
529 492 : CALL dbcsr_deallocate_matrix_set(density_matrix_a)
530 492 : CALL dbcsr_deallocate_matrix_set(density_matrix_ii)
531 492 : CALL dbcsr_deallocate_matrix_set(density_matrix_iii)
532 : !
533 : ! Finalize
534 492 : CALL timestop(handle)
535 : !
536 1476 : END SUBROUTINE current_build_current
537 :
538 : ! **************************************************************************************************
539 : !> \brief Calculation of the idir component of the response current density
540 : !> in the presence of a constant magnetic field in direction iB
541 : !> the current density is collocated on the pw grid in real space
542 : !> \param mat_d0 ...
543 : !> \param mat_jp ...
544 : !> \param mat_jp_rii ...
545 : !> \param mat_jp_riii ...
546 : !> \param iB ...
547 : !> \param idir ...
548 : !> \param current_rs ...
549 : !> \param current_gs ...
550 : !> \param qs_env ...
551 : !> \param current_env ...
552 : !> \param soft_valid ...
553 : !> \param retain_rsgrid ...
554 : !> \note
555 : !> The collocate is done in three parts, one for each density matrix
556 : !> In all cases the density matrices and therefore the collocation
557 : !> are not symmetric, that means that all the pairs (ab and ba) have
558 : !> to be considered separately
559 : !>
560 : !> mat_jp_{\mu\nu} is multiplied by
561 : !> f_{\mu\nu} = \phi_{\mu} (d\phi_{\nu}/dr)_{idir} -
562 : !> (d\phi_{\mu}/dr)_{idir} \phi_{\nu}
563 : !>
564 : !> mat_jp_rii_{\mu\nu} is multiplied by
565 : !> f_{\mu\nu} = \phi_{\mu} (r - R_{\nu})_{iiiB} (d\phi_{\nu}/dr)_{idir} -
566 : !> (d\phi_{\mu}/dr)_{idir} (r - R_{\nu})_{iiiB} \phi_{\nu} +
567 : !> \phi_{\mu} \phi_{\nu} (last term only if iiiB=idir)
568 : !>
569 : !> mat_jp_riii_{\mu\nu} is multiplied by
570 : !> (be careful: change in sign with respect to previous)
571 : !> f_{\mu\nu} = -\phi_{\mu} (r - R_{\nu})_{iiB} (d\phi_{\nu}/dr)_{idir} +
572 : !> (d\phi_{\mu}/dr)_{idir} (r - R_{\nu})_{iiB} \phi_{\nu} -
573 : !> \phi_{\mu} \phi_{\nu} (last term only if iiB=idir)
574 : !>
575 : !> All the terms sum up to the same grid
576 : ! **************************************************************************************************
577 2286 : SUBROUTINE calculate_jrho_resp(mat_d0, mat_jp, mat_jp_rii, mat_jp_riii, iB, idir, &
578 : current_rs, current_gs, qs_env, current_env, soft_valid, retain_rsgrid)
579 :
580 : TYPE(dbcsr_type), POINTER :: mat_d0, mat_jp, mat_jp_rii, mat_jp_riii
581 : INTEGER, INTENT(IN) :: iB, idir
582 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: current_rs
583 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: current_gs
584 : TYPE(qs_environment_type), POINTER :: qs_env
585 : TYPE(current_env_type) :: current_env
586 : LOGICAL, INTENT(IN), OPTIONAL :: soft_valid, retain_rsgrid
587 :
588 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_jrho_resp'
589 : INTEGER, PARAMETER :: max_tasks = 2000
590 :
591 : CHARACTER(LEN=default_string_length) :: basis_type
592 : INTEGER :: adbmdab_func, bcol, brow, cindex, curr_tasks, handle, i, iatom, iatom_old, idir2, &
593 : igrid_level, iiB, iiiB, ikind, ikind_old, ipgf, iset, iset_old, itask, ithread, jatom, &
594 : jatom_old, jkind, jkind_old, jpgf, jset, jset_old, maxco, maxpgf, maxset, maxsgf, &
595 : maxsgf_set, na1, na2, natom, nb1, nb2, ncoa, ncob, nimages, nkind, nseta, nsetb, ntasks, &
596 : nthread, sgfa, sgfb
597 2286 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
598 2286 : npgfb, nsgfa, nsgfb
599 2286 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
600 : LOGICAL :: atom_pair_changed, den_found, &
601 : den_found_a, distributed_rs_grids, &
602 : do_igaim, my_retain_rsgrid, my_soft
603 2286 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_current, my_gauge, my_rho
604 : REAL(KIND=dp) :: eps_rho_rspace, f, kind_radius_a, &
605 : kind_radius_b, Lxo2, Lyo2, Lzo2, &
606 : prefactor, radius, scale, scale2, zetp
607 : REAL(KIND=dp), DIMENSION(3) :: ra, rab, rb, rp
608 2286 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
609 2286 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: jp_block_a, jp_block_b, jp_block_c, jp_block_d, &
610 2286 : jpab_a, jpab_b, jpab_c, jpab_d, rpgfa, rpgfb, sphi_a, sphi_b, work, zeta, zetb
611 2286 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: jpabt_a, jpabt_b, jpabt_c, jpabt_d, workt
612 2286 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
613 : TYPE(cell_type), POINTER :: cell
614 2286 : TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
615 2286 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltajp_a, deltajp_b, deltajp_c, &
616 2286 : deltajp_d
617 : TYPE(dbcsr_type), POINTER :: mat_a, mat_b, mat_c, mat_d
618 : TYPE(dft_control_type), POINTER :: dft_control
619 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
620 2286 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
621 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, orb_basis_set
622 : TYPE(mp_para_env_type), POINTER :: para_env
623 : TYPE(neighbor_list_iterator_p_type), &
624 2286 : DIMENSION(:), POINTER :: nl_iterator
625 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
626 2286 : POINTER :: sab_orb
627 2286 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
628 : TYPE(pw_env_type), POINTER :: pw_env
629 2286 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
630 : TYPE(qs_kind_type), POINTER :: qs_kind
631 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
632 2286 : POINTER :: rs_descs
633 2286 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_current, rs_rho
634 : TYPE(realspace_grid_type), DIMENSION(:, :), &
635 2286 : POINTER :: rs_gauge
636 2286 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
637 :
638 2286 : NULLIFY (qs_kind, cell, dft_control, orb_basis_set, rs_rho, &
639 2286 : qs_kind_set, sab_orb, particle_set, rs_current, pw_env, &
640 2286 : rs_descs, para_env, set_radius_a, set_radius_b, la_max, la_min, &
641 2286 : lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, rpgfa, rpgfb, &
642 2286 : sphi_a, sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, &
643 2286 : workt, mat_a, mat_b, mat_c, mat_d, rs_gauge)
644 2286 : NULLIFY (deltajp_a, deltajp_b, deltajp_c, deltajp_d)
645 2286 : NULLIFY (jp_block_a, jp_block_b, jp_block_c, jp_block_d)
646 2286 : NULLIFY (jpabt_a, jpabt_b, jpabt_c, jpabt_d)
647 2286 : NULLIFY (atom_pair_send, atom_pair_recv)
648 :
649 2286 : CALL timeset(routineN, handle)
650 :
651 : !
652 : ! Set pointers for the different gauge
653 : ! If do_igaim is False the current_env is never needed
654 2286 : do_igaim = current_env%gauge == current_gauge_atom
655 :
656 2286 : mat_a => mat_jp
657 2286 : mat_b => mat_jp_rii
658 2286 : mat_c => mat_jp_riii
659 2286 : IF (do_igaim) mat_d => mat_d0
660 :
661 2286 : my_retain_rsgrid = .FALSE.
662 2286 : IF (PRESENT(retain_rsgrid)) my_retain_rsgrid = retain_rsgrid
663 :
664 : CALL get_qs_env(qs_env=qs_env, &
665 : qs_kind_set=qs_kind_set, &
666 : cell=cell, &
667 : dft_control=dft_control, &
668 : particle_set=particle_set, &
669 : sab_all=sab_orb, &
670 : para_env=para_env, &
671 2286 : pw_env=pw_env)
672 :
673 2286 : IF (do_igaim) CALL get_current_env(current_env=current_env, rs_gauge=rs_gauge)
674 :
675 : ! Component of appearing in the vector product rxp, iiB and iiiB
676 2286 : CALL set_vecp(iB, iiB, iiiB)
677 : !
678 : !
679 2286 : scale2 = 0.0_dp
680 2286 : idir2 = 1
681 2286 : IF (idir /= iB) THEN
682 1428 : CALL set_vecp_rev(idir, iB, idir2)
683 1428 : scale2 = fac_vecp(idir, iB, idir2)
684 : END IF
685 : !
686 : ! *** assign from pw_env
687 2286 : gridlevel_info => pw_env%gridlevel_info
688 2286 : cube_info => pw_env%cube_info
689 :
690 : ! Check that the neighbor list with all the pairs is associated
691 2286 : CPASSERT(ASSOCIATED(sab_orb))
692 : ! *** set up the pw multi-grids
693 2286 : CPASSERT(ASSOCIATED(pw_env))
694 2286 : CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
695 :
696 2286 : distributed_rs_grids = .FALSE.
697 11394 : DO igrid_level = 1, gridlevel_info%ngrid_levels
698 38718 : IF (.NOT. ALL(rs_descs(igrid_level)%rs_desc%perd == 1)) THEN
699 0 : distributed_rs_grids = .TRUE.
700 : END IF
701 : END DO
702 2286 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
703 2286 : nthread = 1
704 :
705 : ! *** Allocate work storage ***
706 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
707 : maxco=maxco, &
708 : maxsgf=maxsgf, &
709 2286 : maxsgf_set=maxsgf_set)
710 :
711 9144 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
712 9144 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
713 9144 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
714 :
715 2286 : my_soft = .FALSE.
716 2286 : IF (PRESENT(soft_valid)) my_soft = soft_valid
717 2142 : IF (my_soft) THEN
718 1260 : basis_type = "ORB_SOFT"
719 : ELSE
720 1026 : basis_type = "ORB"
721 : END IF
722 :
723 2286 : nkind = SIZE(qs_kind_set)
724 :
725 2286 : CALL reallocate(jpabt_a, 1, maxco, 1, maxco, 0, nthread - 1)
726 2286 : CALL reallocate(jpabt_b, 1, maxco, 1, maxco, 0, nthread - 1)
727 2286 : CALL reallocate(jpabt_c, 1, maxco, 1, maxco, 0, nthread - 1)
728 2286 : CALL reallocate(jpabt_d, 1, maxco, 1, maxco, 0, nthread - 1)
729 2286 : CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
730 2286 : CALL reallocate_tasks(tasks, max_tasks)
731 :
732 2286 : ntasks = 0
733 2286 : curr_tasks = SIZE(tasks)
734 :
735 : ! get maximum numbers
736 2286 : natom = SIZE(particle_set)
737 2286 : maxset = 0
738 2286 : maxpgf = 0
739 :
740 : ! hard code matrix index (no kpoints)
741 2286 : nimages = dft_control%nimages
742 2286 : CPASSERT(nimages == 1)
743 2286 : cindex = 1
744 :
745 6126 : DO ikind = 1, nkind
746 3840 : qs_kind => qs_kind_set(ikind)
747 :
748 3840 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
749 :
750 3840 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
751 :
752 3840 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, npgf=npgfa, nset=nseta)
753 3840 : maxset = MAX(nseta, maxset)
754 15696 : maxpgf = MAX(MAXVAL(npgfa), maxpgf)
755 : END DO
756 :
757 : ! *** Initialize working density matrix ***
758 :
759 : ! distributed rs grids require a matrix that will be changed (distribute_tasks)
760 : ! whereas this is not the case for replicated grids
761 11430 : ALLOCATE (deltajp_a(1), deltajp_b(1), deltajp_c(1), deltajp_d(1))
762 2286 : IF (distributed_rs_grids) THEN
763 0 : ALLOCATE (deltajp_a(1)%matrix, deltajp_b(1)%matrix, deltajp_c(1)%matrix)
764 0 : IF (do_igaim) THEN
765 0 : ALLOCATE (deltajp_d(1)%matrix)
766 : END IF
767 :
768 0 : CALL dbcsr_create(deltajp_a(1)%matrix, template=mat_a, name='deltajp_a')
769 0 : CALL dbcsr_create(deltajp_b(1)%matrix, template=mat_a, name='deltajp_b')
770 0 : CALL dbcsr_create(deltajp_c(1)%matrix, template=mat_a, name='deltajp_c')
771 0 : IF (do_igaim) CALL dbcsr_create(deltajp_d(1)%matrix, template=mat_a, name='deltajp_d')
772 : ELSE
773 2286 : deltajp_a(1)%matrix => mat_a !mat_jp
774 2286 : deltajp_b(1)%matrix => mat_b !mat_jp_rii
775 2286 : deltajp_c(1)%matrix => mat_c !mat_jp_riii
776 2286 : IF (do_igaim) deltajp_d(1)%matrix => mat_d !mat_d0
777 : END IF
778 :
779 10698 : ALLOCATE (basis_set_list(nkind))
780 6126 : DO ikind = 1, nkind
781 3840 : qs_kind => qs_kind_set(ikind)
782 3840 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
783 6126 : IF (ASSOCIATED(basis_set_a)) THEN
784 3840 : basis_set_list(ikind)%gto_basis_set => basis_set_a
785 : ELSE
786 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
787 : END IF
788 : END DO
789 2286 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
790 182301 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
791 180015 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
792 180015 : basis_set_a => basis_set_list(ikind)%gto_basis_set
793 180015 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
794 180015 : basis_set_b => basis_set_list(jkind)%gto_basis_set
795 180015 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
796 180015 : ra(:) = pbc(particle_set(iatom)%r, cell)
797 : ! basis ikind
798 180015 : first_sgfa => basis_set_a%first_sgf
799 180015 : la_max => basis_set_a%lmax
800 180015 : la_min => basis_set_a%lmin
801 180015 : npgfa => basis_set_a%npgf
802 180015 : nseta = basis_set_a%nset
803 180015 : nsgfa => basis_set_a%nsgf_set
804 180015 : rpgfa => basis_set_a%pgf_radius
805 180015 : set_radius_a => basis_set_a%set_radius
806 180015 : kind_radius_a = basis_set_a%kind_radius
807 180015 : sphi_a => basis_set_a%sphi
808 180015 : zeta => basis_set_a%zet
809 : ! basis jkind
810 180015 : first_sgfb => basis_set_b%first_sgf
811 180015 : lb_max => basis_set_b%lmax
812 180015 : lb_min => basis_set_b%lmin
813 180015 : npgfb => basis_set_b%npgf
814 180015 : nsetb = basis_set_b%nset
815 180015 : nsgfb => basis_set_b%nsgf_set
816 180015 : rpgfb => basis_set_b%pgf_radius
817 180015 : set_radius_b => basis_set_b%set_radius
818 180015 : kind_radius_b = basis_set_b%kind_radius
819 180015 : sphi_b => basis_set_b%sphi
820 180015 : zetb => basis_set_b%zet
821 :
822 180015 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) THEN
823 : CYCLE
824 : END IF
825 :
826 11178 : brow = iatom
827 11178 : bcol = jatom
828 :
829 : CALL dbcsr_get_block_p(matrix=mat_a, row=brow, col=bcol, &
830 11178 : block=jp_block_a, found=den_found_a)
831 : CALL dbcsr_get_block_p(matrix=mat_b, row=brow, col=bcol, &
832 11178 : block=jp_block_b, found=den_found)
833 : CALL dbcsr_get_block_p(matrix=mat_c, row=brow, col=bcol, &
834 11178 : block=jp_block_c, found=den_found)
835 11178 : IF (do_igaim) CALL dbcsr_get_block_p(matrix=mat_d, row=brow, col=bcol, &
836 2961 : block=jp_block_d, found=den_found)
837 :
838 11178 : IF (.NOT. ASSOCIATED(jp_block_a)) CYCLE
839 :
840 11061 : IF (distributed_rs_grids) THEN
841 0 : CALL dbcsr_put_block(deltajp_a(1)%matrix, brow, bcol, jp_block_a)
842 0 : CALL dbcsr_put_block(deltajp_b(1)%matrix, brow, bcol, jp_block_b)
843 0 : CALL dbcsr_put_block(deltajp_c(1)%matrix, brow, bcol, jp_block_c)
844 0 : IF (do_igaim) THEN
845 0 : CALL dbcsr_put_block(deltajp_d(1)%matrix, brow, bcol, jp_block_d)
846 : END IF
847 : END IF
848 :
849 : CALL task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, &
850 : dft_control, cube_info, gridlevel_info, cindex, &
851 : iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, &
852 : set_radius_a, set_radius_b, ra, rab, &
853 180015 : la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
854 :
855 : END DO
856 2286 : CALL neighbor_list_iterator_release(nl_iterator)
857 :
858 2286 : DEALLOCATE (basis_set_list)
859 :
860 2286 : IF (distributed_rs_grids) THEN
861 0 : CALL dbcsr_finalize(deltajp_a(1)%matrix)
862 0 : CALL dbcsr_finalize(deltajp_b(1)%matrix)
863 0 : CALL dbcsr_finalize(deltajp_c(1)%matrix)
864 0 : IF (do_igaim) CALL dbcsr_finalize(deltajp_d(1)%matrix)
865 : END IF
866 :
867 : ! sorts / redistributes the task list
868 : CALL distribute_tasks(rs_descs=rs_descs, ntasks=ntasks, natoms=natom, tasks=tasks, &
869 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
870 : symmetric=.FALSE., reorder_rs_grid_ranks=.TRUE., &
871 2286 : skip_load_balance_distributed=.FALSE.)
872 :
873 50256 : ALLOCATE (rs_current(gridlevel_info%ngrid_levels))
874 :
875 11394 : DO igrid_level = 1, gridlevel_info%ngrid_levels
876 : ! Here we need to reallocate the distributed rs_grids, which may have been reordered
877 : ! by distribute_tasks
878 9108 : IF (rs_descs(igrid_level)%rs_desc%distributed .AND. .NOT. my_retain_rsgrid) THEN
879 0 : CALL rs_grid_release(rs_rho(igrid_level))
880 0 : CALL rs_grid_create(rs_rho(igrid_level), rs_descs(igrid_level)%rs_desc)
881 : END IF
882 9108 : CALL rs_grid_zero(rs_rho(igrid_level))
883 9108 : CALL rs_grid_create(rs_current(igrid_level), rs_descs(igrid_level)%rs_desc)
884 11394 : CALL rs_grid_zero(rs_current(igrid_level))
885 : END DO
886 :
887 : !
888 : ! we need to build the gauge here
889 2286 : IF (.NOT. current_env%gauge_init .AND. do_igaim) THEN
890 28 : CALL current_set_gauge(current_env, qs_env)
891 28 : current_env%gauge_init = .TRUE.
892 : END IF
893 : !
894 : ! for any case double check the bounds !
895 360 : IF (do_igaim) THEN
896 1764 : DO igrid_level = 1, gridlevel_info%ngrid_levels
897 1404 : my_rho => rs_rho(igrid_level)%r
898 1404 : my_current => rs_current(igrid_level)%r
899 : IF (LBOUND(my_rho, 3) /= LBOUND(my_current, 3) .OR. &
900 : LBOUND(my_rho, 2) /= LBOUND(my_current, 2) .OR. &
901 : LBOUND(my_rho, 1) /= LBOUND(my_current, 1) .OR. &
902 : UBOUND(my_rho, 3) /= UBOUND(my_current, 3) .OR. &
903 18252 : UBOUND(my_rho, 2) /= UBOUND(my_current, 2) .OR. &
904 : UBOUND(my_rho, 1) /= UBOUND(my_current, 1)) THEN
905 0 : WRITE (*, *) 'LBOUND(my_rho,3),LBOUND(my_current,3)', LBOUND(my_rho, 3), LBOUND(my_current, 3)
906 0 : WRITE (*, *) 'LBOUND(my_rho,2),LBOUND(my_current,2)', LBOUND(my_rho, 2), LBOUND(my_current, 2)
907 0 : WRITE (*, *) 'LBOUND(my_rho,1),LBOUND(my_current,1)', LBOUND(my_rho, 1), LBOUND(my_current, 1)
908 0 : WRITE (*, *) 'UBOUND(my_rho,3),UBOUND(my_current,3)', UBOUND(my_rho, 3), UBOUND(my_current, 3)
909 0 : WRITE (*, *) 'UBOUND(my_rho,2),UBOUND(my_current,2)', UBOUND(my_rho, 2), UBOUND(my_current, 2)
910 0 : WRITE (*, *) 'UBOUND(my_rho,1),UBOUND(my_current,1)', UBOUND(my_rho, 1), UBOUND(my_current, 1)
911 0 : CPABORT("Bug")
912 : END IF
913 :
914 1404 : my_gauge => rs_gauge(1, igrid_level)%r
915 : IF (LBOUND(my_rho, 3) /= LBOUND(my_gauge, 3) .OR. &
916 : LBOUND(my_rho, 2) /= LBOUND(my_gauge, 2) .OR. &
917 : LBOUND(my_rho, 1) /= LBOUND(my_gauge, 1) .OR. &
918 : UBOUND(my_rho, 3) /= UBOUND(my_gauge, 3) .OR. &
919 16848 : UBOUND(my_rho, 2) /= UBOUND(my_gauge, 2) .OR. &
920 360 : UBOUND(my_rho, 1) /= UBOUND(my_gauge, 1)) THEN
921 0 : WRITE (*, *) 'LBOUND(my_rho,3),LBOUND(my_gauge,3)', LBOUND(my_rho, 3), LBOUND(my_gauge, 3)
922 0 : WRITE (*, *) 'LBOUND(my_rho,2),LBOUND(my_gauge,2)', LBOUND(my_rho, 2), LBOUND(my_gauge, 2)
923 0 : WRITE (*, *) 'LBOUND(my_rho,1),LBOUND(my_gauge,1)', LBOUND(my_rho, 1), LBOUND(my_gauge, 1)
924 0 : WRITE (*, *) 'UBOUND(my_rho,3),UbOUND(my_gauge,3)', UBOUND(my_rho, 3), UBOUND(my_gauge, 3)
925 0 : WRITE (*, *) 'UBOUND(my_rho,2),UBOUND(my_gauge,2)', UBOUND(my_rho, 2), UBOUND(my_gauge, 2)
926 0 : WRITE (*, *) 'UBOUND(my_rho,1),UBOUND(my_gauge,1)', UBOUND(my_rho, 1), UBOUND(my_gauge, 1)
927 0 : CPABORT("Bug")
928 : END IF
929 : END DO
930 : END IF
931 : !
932 : !-------------------------------------------------------------
933 :
934 2286 : IF (distributed_rs_grids) THEN
935 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_a, &
936 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
937 0 : nimages=nimages, scatter=.TRUE.)
938 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_b, &
939 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
940 0 : nimages=nimages, scatter=.TRUE.)
941 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_c, &
942 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
943 0 : nimages=nimages, scatter=.TRUE.)
944 0 : IF (do_igaim) CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_d, &
945 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
946 0 : nimages=nimages, scatter=.TRUE.)
947 : END IF
948 :
949 2286 : ithread = 0
950 2286 : jpab_a => jpabt_a(:, :, ithread)
951 2286 : jpab_b => jpabt_b(:, :, ithread)
952 2286 : jpab_c => jpabt_c(:, :, ithread)
953 2286 : IF (do_igaim) jpab_d => jpabt_d(:, :, ithread)
954 2286 : work => workt(:, :, ithread)
955 :
956 2286 : iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
957 2286 : ikind_old = -1; jkind_old = -1
958 :
959 156348 : loop_tasks: DO itask = 1, ntasks
960 154062 : igrid_level = tasks(itask)%grid_level
961 154062 : cindex = tasks(itask)%image
962 154062 : iatom = tasks(itask)%iatom
963 154062 : jatom = tasks(itask)%jatom
964 154062 : iset = tasks(itask)%iset
965 154062 : jset = tasks(itask)%jset
966 154062 : ipgf = tasks(itask)%ipgf
967 154062 : jpgf = tasks(itask)%jpgf
968 :
969 : ! apparently generalised collocation not implemented correctly yet
970 154062 : CPASSERT(tasks(itask)%dist_type /= 2)
971 :
972 154062 : IF (iatom /= iatom_old .OR. jatom /= jatom_old) THEN
973 :
974 23058 : ikind = particle_set(iatom)%atomic_kind%kind_number
975 23058 : jkind = particle_set(jatom)%atomic_kind%kind_number
976 :
977 23058 : IF (iatom /= iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
978 :
979 23058 : brow = iatom
980 23058 : bcol = jatom
981 :
982 23058 : IF (ikind /= ikind_old) THEN
983 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
984 3342 : basis_type=basis_type)
985 :
986 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
987 : first_sgf=first_sgfa, &
988 : lmax=la_max, &
989 : lmin=la_min, &
990 : npgf=npgfa, &
991 : nset=nseta, &
992 : nsgf_set=nsgfa, &
993 : pgf_radius=rpgfa, &
994 : set_radius=set_radius_a, &
995 : sphi=sphi_a, &
996 3342 : zet=zeta)
997 : END IF
998 :
999 23058 : IF (jkind /= jkind_old) THEN
1000 :
1001 : CALL get_qs_kind(qs_kind_set(jkind), &
1002 12252 : basis_set=orb_basis_set, basis_type=basis_type)
1003 :
1004 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1005 : first_sgf=first_sgfb, &
1006 : kind_radius=kind_radius_b, &
1007 : lmax=lb_max, &
1008 : lmin=lb_min, &
1009 : npgf=npgfb, &
1010 : nset=nsetb, &
1011 : nsgf_set=nsgfb, &
1012 : pgf_radius=rpgfb, &
1013 : set_radius=set_radius_b, &
1014 : sphi=sphi_b, &
1015 12252 : zet=zetb)
1016 :
1017 : END IF
1018 :
1019 : CALL dbcsr_get_block_p(matrix=deltajp_a(1)%matrix, row=brow, col=bcol, &
1020 23058 : block=jp_block_a, found=den_found)
1021 : CALL dbcsr_get_block_p(matrix=deltajp_b(1)%matrix, row=brow, col=bcol, &
1022 23058 : block=jp_block_b, found=den_found)
1023 : CALL dbcsr_get_block_p(matrix=deltajp_c(1)%matrix, row=brow, col=bcol, &
1024 23058 : block=jp_block_c, found=den_found)
1025 23058 : IF (do_igaim) CALL dbcsr_get_block_p(matrix=deltajp_d(1)%matrix, row=brow, col=bcol, &
1026 5229 : block=jp_block_d, found=den_found)
1027 :
1028 23058 : IF (.NOT. ASSOCIATED(jp_block_a)) THEN
1029 0 : CPABORT("p_block not associated in deltap")
1030 : END IF
1031 :
1032 : iatom_old = iatom
1033 : jatom_old = jatom
1034 : ikind_old = ikind
1035 : jkind_old = jkind
1036 :
1037 : atom_pair_changed = .TRUE.
1038 :
1039 : ELSE
1040 :
1041 : atom_pair_changed = .FALSE.
1042 :
1043 : END IF
1044 :
1045 154062 : IF (atom_pair_changed .OR. iset_old /= iset .OR. jset_old /= jset) THEN
1046 :
1047 55695 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1048 55695 : sgfa = first_sgfa(1, iset)
1049 55695 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1050 55695 : sgfb = first_sgfb(1, jset)
1051 : ! Decontraction step for the selected blocks of the 3 density matrices
1052 :
1053 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1054 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1055 : jp_block_a(sgfa, sgfb), SIZE(jp_block_a, 1), &
1056 55695 : 0.0_dp, work(1, 1), maxco)
1057 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1058 : 1.0_dp, work(1, 1), maxco, &
1059 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1060 55695 : 0.0_dp, jpab_a(1, 1), maxco)
1061 :
1062 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1063 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1064 : jp_block_b(sgfa, sgfb), SIZE(jp_block_b, 1), &
1065 55695 : 0.0_dp, work(1, 1), maxco)
1066 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1067 : 1.0_dp, work(1, 1), maxco, &
1068 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1069 55695 : 0.0_dp, jpab_b(1, 1), maxco)
1070 :
1071 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1072 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1073 : jp_block_c(sgfa, sgfb), SIZE(jp_block_c, 1), &
1074 55695 : 0.0_dp, work(1, 1), maxco)
1075 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1076 : 1.0_dp, work(1, 1), maxco, &
1077 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1078 55695 : 0.0_dp, jpab_c(1, 1), maxco)
1079 :
1080 55695 : IF (do_igaim) THEN
1081 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1082 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1083 : jp_block_d(sgfa, sgfb), SIZE(jp_block_d, 1), &
1084 12105 : 0.0_dp, work(1, 1), maxco)
1085 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1086 : 1.0_dp, work(1, 1), maxco, &
1087 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1088 12105 : 0.0_dp, jpab_d(1, 1), maxco)
1089 : END IF
1090 :
1091 : iset_old = iset
1092 : jset_old = jset
1093 :
1094 : END IF
1095 :
1096 51354 : SELECT CASE (idir)
1097 : CASE (1)
1098 51354 : adbmdab_func = GRID_FUNC_ADBmDAB_X
1099 : CASE (2)
1100 51354 : adbmdab_func = GRID_FUNC_ADBmDAB_Y
1101 : CASE (3)
1102 51354 : adbmdab_func = GRID_FUNC_ADBmDAB_Z
1103 : CASE DEFAULT
1104 154062 : CPABORT("invalid idir")
1105 : END SELECT
1106 :
1107 616248 : rab(:) = tasks(itask)%rab
1108 616248 : rb(:) = ra(:) + rab(:)
1109 154062 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
1110 154062 : f = zetb(jpgf, jset)/zetp
1111 616248 : prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
1112 616248 : rp(:) = ra(:) + f*rab(:)
1113 :
1114 154062 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
1115 154062 : na2 = ipgf*ncoset(la_max(iset))
1116 154062 : nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
1117 154062 : nb2 = jpgf*ncoset(lb_max(jset))
1118 :
1119 : ! Four calls to the general collocate density, to multply the correct function
1120 : ! to each density matrix
1121 :
1122 : !
1123 : ! here the decontracted mat_jp_{ab} is multiplied by
1124 : ! f_{ab} = g_{a} (dg_{b}/dr)_{idir} - (dg_{a}/dr)_{idir} g_{b}
1125 154062 : scale = 1.0_dp
1126 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1127 : lb_min=lb_min(jset), lb_max=lb_max(jset), &
1128 : ra=ra, rb=rb, rp=rp, zetp=zetp, eps=eps_rho_rspace, &
1129 154062 : prefactor=prefactor, cutoff=1.0_dp)
1130 :
1131 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1132 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1133 : ra, rab, scale, jpab_a, na1 - 1, nb1 - 1, &
1134 : rs_current(igrid_level), &
1135 154062 : radius=radius, ga_gb_function=adbmdab_func)
1136 156348 : IF (do_igaim) THEN
1137 : ! here the decontracted mat_jb_{ab} is multiplied by
1138 : ! f_{ab} = g_{a} * g_{b} ! THIS GOES OUTSIDE THE LOOP !
1139 20574 : IF (scale2 /= 0.0_dp) THEN
1140 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1141 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1142 : ra, rab, scale2, jpab_d, na1 - 1, nb1 - 1, &
1143 : rs_rho(igrid_level), &
1144 13716 : radius=radius, ga_gb_function=GRID_FUNC_AB)
1145 : END IF !rm
1146 : ! here the decontracted mat_jp_rii{ab} is multiplied by
1147 : ! f_{ab} = g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} -
1148 : ! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b}
1149 : scale = 1.0_dp
1150 512067753 : current_env%rs_buf(igrid_level)%r(:, :, :) = 0.0_dp
1151 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1152 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1153 : ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
1154 : radius=radius, &
1155 : ga_gb_function=adbmdab_func, &
1156 20574 : rsgrid=current_env%rs_buf(igrid_level))
1157 : CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level), &
1158 : rsbuf=rs_current(igrid_level), &
1159 : rsgauge=rs_gauge(iiiB, igrid_level), &
1160 : cube_info=cube_info(igrid_level), radius=radius, &
1161 : zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
1162 20574 : ra=ra, rab=rab, ir=iiiB)
1163 :
1164 : ! here the decontracted mat_jp_riii{ab} is multiplied by
1165 : ! f_{ab} = -g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} +
1166 : ! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b}
1167 20574 : scale = -1.0_dp
1168 512067753 : current_env%rs_buf(igrid_level)%r(:, :, :) = 0.0_dp
1169 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1170 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1171 : ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
1172 : radius=radius, &
1173 : ga_gb_function=adbmdab_func, &
1174 20574 : rsgrid=current_env%rs_buf(igrid_level))
1175 : CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level), &
1176 : rsbuf=rs_current(igrid_level), &
1177 : rsgauge=rs_gauge(iiB, igrid_level), &
1178 : cube_info=cube_info(igrid_level), radius=radius, &
1179 : zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
1180 20574 : ra=ra, rab=rab, ir=iiB)
1181 : ELSE
1182 : ! here the decontracted mat_jp_rii{ab} is multiplied by
1183 : ! f_{ab} = g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} -
1184 : ! (dg_{a}/dr)_{idir} (r - R_{b})_{iiB} g_{b}
1185 : scale = 1.0_dp
1186 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1187 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1188 : ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
1189 : rs_current(igrid_level), &
1190 : radius=radius, &
1191 133488 : ga_gb_function=encode_ardbmdarb_func(idir=idir, ir=iiiB))
1192 : ! here the decontracted mat_jp_riii{ab} is multiplied by
1193 : ! f_{ab} = -g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} +
1194 : ! (dg_{a}/dr)_{idir} (r - R_{b})_{iiB} g_{b}
1195 133488 : scale = -1.0_dp
1196 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1197 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1198 : ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
1199 : rs_current(igrid_level), &
1200 : radius=radius, &
1201 133488 : ga_gb_function=encode_ardbmdarb_func(idir=idir, ir=iiB))
1202 : END IF
1203 :
1204 : END DO loop_tasks
1205 : !
1206 : ! Scale the density with the gauge rho * ( r - d(r) ) if needed
1207 2286 : IF (do_igaim) THEN
1208 1764 : DO igrid_level = 1, gridlevel_info%ngrid_levels
1209 : CALL rs_grid_mult_and_add(rs_current(igrid_level), rs_rho(igrid_level), &
1210 1764 : rs_gauge(idir2, igrid_level), 1.0_dp)
1211 : END DO
1212 : END IF
1213 : ! *** Release work storage ***
1214 :
1215 2286 : IF (distributed_rs_grids) THEN
1216 0 : CALL dbcsr_deallocate_matrix(deltajp_a(1)%matrix)
1217 0 : CALL dbcsr_deallocate_matrix(deltajp_b(1)%matrix)
1218 0 : CALL dbcsr_deallocate_matrix(deltajp_c(1)%matrix)
1219 0 : IF (do_igaim) CALL dbcsr_deallocate_matrix(deltajp_d(1)%matrix)
1220 : END IF
1221 2286 : DEALLOCATE (deltajp_a, deltajp_b, deltajp_c, deltajp_d)
1222 :
1223 2286 : DEALLOCATE (jpabt_a, jpabt_b, jpabt_c, jpabt_d, workt, tasks)
1224 :
1225 2286 : IF (ASSOCIATED(atom_pair_send)) DEALLOCATE (atom_pair_send)
1226 2286 : IF (ASSOCIATED(atom_pair_recv)) DEALLOCATE (atom_pair_recv)
1227 :
1228 2286 : CALL density_rs2pw(pw_env, rs_current, current_rs, current_gs)
1229 11394 : DO i = 1, SIZE(rs_current)
1230 11394 : CALL rs_grid_release(rs_current(i))
1231 : END DO
1232 :
1233 11394 : DO i = 1, SIZE(rs_rho)
1234 11394 : IF (rs_descs(i)%rs_desc%distributed .AND. .NOT. my_retain_rsgrid) THEN
1235 0 : CALL rs_grid_release(rs_rho(i))
1236 : END IF
1237 : END DO
1238 :
1239 : ! Free the array of grids (grids themselves are released in density_rs2pw)
1240 2286 : DEALLOCATE (rs_current)
1241 :
1242 2286 : CALL timestop(handle)
1243 :
1244 6858 : END SUBROUTINE calculate_jrho_resp
1245 :
1246 : ! **************************************************************************************************
1247 : !> \brief ...
1248 : !> \param idir ...
1249 : !> \param ir ...
1250 : !> \return ...
1251 : ! **************************************************************************************************
1252 266976 : FUNCTION encode_ardbmdarb_func(idir, ir) RESULT(func)
1253 : INTEGER, INTENT(IN) :: idir, ir
1254 : INTEGER :: func
1255 :
1256 266976 : CPASSERT(1 <= idir .AND. idir <= 3 .AND. 1 <= ir .AND. ir <= 3)
1257 266976 : SELECT CASE (10*idir + ir)
1258 : CASE (11)
1259 30426 : func = GRID_FUNC_ARDBmDARB_XX
1260 : CASE (12)
1261 30426 : func = GRID_FUNC_ARDBmDARB_XY
1262 : CASE (13)
1263 30426 : func = GRID_FUNC_ARDBmDARB_XZ
1264 : CASE (21)
1265 30426 : func = GRID_FUNC_ARDBmDARB_YX
1266 : CASE (22)
1267 28140 : func = GRID_FUNC_ARDBmDARB_YY
1268 : CASE (23)
1269 30426 : func = GRID_FUNC_ARDBmDARB_YZ
1270 : CASE (31)
1271 30426 : func = GRID_FUNC_ARDBmDARB_ZX
1272 : CASE (32)
1273 30426 : func = GRID_FUNC_ARDBmDARB_ZY
1274 : CASE (33)
1275 28140 : func = GRID_FUNC_ARDBmDARB_ZZ
1276 : CASE DEFAULT
1277 266976 : CPABORT("invalid idir or iiiB")
1278 : END SELECT
1279 266976 : END FUNCTION encode_ardbmdarb_func
1280 :
1281 : ! **************************************************************************************************
1282 : !> \brief ...
1283 : !> \param rsgrid ...
1284 : !> \param rsbuf ...
1285 : !> \param rsgauge ...
1286 : !> \param cube_info ...
1287 : !> \param radius ...
1288 : !> \param ra ...
1289 : !> \param rab ...
1290 : !> \param zeta ...
1291 : !> \param zetb ...
1292 : !> \param ir ...
1293 : ! **************************************************************************************************
1294 41148 : SUBROUTINE collocate_gauge_ortho(rsgrid, rsbuf, rsgauge, cube_info, radius, ra, rab, zeta, zetb, ir)
1295 : TYPE(realspace_grid_type) :: rsgrid, rsbuf, rsgauge
1296 : TYPE(cube_info_type), INTENT(IN) :: cube_info
1297 : REAL(KIND=dp), INTENT(IN) :: radius
1298 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab
1299 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb
1300 : INTEGER, INTENT(IN) :: ir
1301 :
1302 : INTEGER :: cmax, i, ig, igmax, igmin, j, j2, jg, &
1303 : jg2, jgmin, k, k2, kg, kg2, kgmin, &
1304 : length, offset, sci, start
1305 41148 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map
1306 : INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, ub_cube
1307 41148 : INTEGER, DIMENSION(:), POINTER :: sphere_bounds
1308 : REAL(KIND=dp) :: f, point(3, 4), res(4), x, y, y2, z, z2, &
1309 : zetp
1310 : REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp
1311 41148 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gauge, grid, grid_buf
1312 :
1313 0 : CPASSERT(rsgrid%desc%orthorhombic)
1314 41148 : NULLIFY (sphere_bounds)
1315 :
1316 41148 : grid => rsgrid%r(:, :, :)
1317 41148 : grid_buf => rsbuf%r(:, :, :)
1318 41148 : gauge => rsgauge%r(:, :, :)
1319 :
1320 : ! *** center of gaussians and their product
1321 41148 : zetp = zeta + zetb
1322 41148 : f = zetb/zetp
1323 164592 : rap(:) = f*rab(:)
1324 164592 : rbp(:) = rap(:) - rab(:)
1325 164592 : rp(:) = ra(:) + rap(:)
1326 164592 : rb(:) = ra(:) + rab(:)
1327 :
1328 : ! *** properties of the grid ***
1329 164592 : ng(:) = rsgrid%desc%npts(:)
1330 41148 : dr(1) = rsgrid%desc%dh(1, 1)
1331 41148 : dr(2) = rsgrid%desc%dh(2, 2)
1332 41148 : dr(3) = rsgrid%desc%dh(3, 3)
1333 :
1334 : ! *** get the sub grid properties for the given radius ***
1335 41148 : CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
1336 164592 : cmax = MAXVAL(ub_cube)
1337 :
1338 : ! *** position of the gaussian product
1339 : !
1340 : ! this is the actual definition of the position on the grid
1341 : ! i.e. a point rp(:) gets here grid coordinates
1342 : ! MODULO(rp(:)/dr(:),ng(:))+1
1343 : ! hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid ((1,1,1) on grid)
1344 : !
1345 164592 : ALLOCATE (map(-cmax:cmax, 3))
1346 3221964 : map(:, :) = -1
1347 41148 : CALL compute_cube_center(cubecenter, rsgrid%desc, zeta, zetb, ra, rab)
1348 164592 : roffset(:) = rp(:) - REAL(cubecenter(:), dp)*dr(:)
1349 :
1350 : ! *** a mapping so that the ig corresponds to the right grid point
1351 164592 : DO i = 1, 3
1352 164592 : IF (rsgrid%desc%perd(i) == 1) THEN
1353 123444 : start = lb_cube(i)
1354 123354 : DO
1355 246798 : offset = MODULO(cubecenter(i) + start, ng(i)) + 1 - start
1356 246798 : length = MIN(ub_cube(i), ng(i) - offset) - start
1357 3180726 : DO ig = start, start + length
1358 3180726 : map(ig, i) = ig + offset
1359 : END DO
1360 246798 : IF (start + length >= ub_cube(i)) EXIT
1361 123354 : start = start + length + 1
1362 : END DO
1363 : ELSE
1364 : ! this takes partial grid + border regions into account
1365 0 : offset = MODULO(cubecenter(i) + lb_cube(i) + rsgrid%desc%lb(i) - rsgrid%lb_local(i), ng(i)) + 1 - lb_cube(i)
1366 : ! check for out of bounds
1367 0 : IF (ub_cube(i) + offset > UBOUND(grid, i) .OR. lb_cube(i) + offset < LBOUND(grid, i)) THEN
1368 0 : CPABORT("Cube bound out of grid bounds")
1369 : END IF
1370 0 : DO ig = lb_cube(i), ub_cube(i)
1371 0 : map(ig, i) = ig + offset
1372 : END DO
1373 : END IF
1374 : END DO
1375 :
1376 : ! *** actually loop over the grid
1377 41148 : sci = 1
1378 41148 : kgmin = sphere_bounds(sci)
1379 41148 : sci = sci + 1
1380 530136 : DO kg = kgmin, 0
1381 488988 : kg2 = 1 - kg
1382 488988 : k = map(kg, 3)
1383 488988 : k2 = map(kg2, 3)
1384 488988 : jgmin = sphere_bounds(sci)
1385 488988 : sci = sci + 1
1386 488988 : z = (REAL(kg, dp) + REAL(cubecenter(3), dp))*dr(3)
1387 488988 : z2 = (REAL(kg2, dp) + REAL(cubecenter(3), dp))*dr(3)
1388 5604174 : DO jg = jgmin, 0
1389 5074038 : jg2 = 1 - jg
1390 5074038 : j = map(jg, 2)
1391 5074038 : j2 = map(jg2, 2)
1392 5074038 : igmin = sphere_bounds(sci)
1393 5074038 : sci = sci + 1
1394 5074038 : igmax = 1 - igmin
1395 5074038 : y = (REAL(jg, dp) + REAL(cubecenter(2), dp))*dr(2)
1396 5074038 : y2 = (REAL(jg2, dp) + REAL(cubecenter(2), dp))*dr(2)
1397 118256958 : DO ig = igmin, igmax
1398 112693932 : i = map(ig, 1)
1399 112693932 : x = (REAL(ig, dp) + REAL(cubecenter(1), dp))*dr(1)
1400 112693932 : point(1, 1) = x; point(2, 1) = y; point(3, 1) = z
1401 112693932 : point(1, 2) = x; point(2, 2) = y2; point(3, 2) = z
1402 112693932 : point(1, 3) = x; point(2, 3) = y; point(3, 3) = z2
1403 112693932 : point(1, 4) = x; point(2, 4) = y2; point(3, 4) = z2
1404 : !
1405 112693932 : res(1) = (point(ir, 1) - rb(ir)) - gauge(i, j, k)
1406 112693932 : res(2) = (point(ir, 2) - rb(ir)) - gauge(i, j2, k)
1407 112693932 : res(3) = (point(ir, 3) - rb(ir)) - gauge(i, j, k2)
1408 112693932 : res(4) = (point(ir, 4) - rb(ir)) - gauge(i, j2, k2)
1409 : !
1410 112693932 : grid_buf(i, j, k) = grid_buf(i, j, k) + grid(i, j, k)*res(1)
1411 112693932 : grid_buf(i, j2, k) = grid_buf(i, j2, k) + grid(i, j2, k)*res(2)
1412 112693932 : grid_buf(i, j, k2) = grid_buf(i, j, k2) + grid(i, j, k2)*res(3)
1413 117767970 : grid_buf(i, j2, k2) = grid_buf(i, j2, k2) + grid(i, j2, k2)*res(4)
1414 : END DO
1415 : END DO
1416 : END DO
1417 82296 : END SUBROUTINE collocate_gauge_ortho
1418 :
1419 : ! **************************************************************************************************
1420 : !> \brief ...
1421 : !> \param current_env ...
1422 : !> \param qs_env ...
1423 : ! **************************************************************************************************
1424 28 : SUBROUTINE current_set_gauge(current_env, qs_env)
1425 : !
1426 : TYPE(current_env_type) :: current_env
1427 : TYPE(qs_environment_type), POINTER :: qs_env
1428 :
1429 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_set_gauge'
1430 :
1431 : INTEGER :: idir
1432 : REAL(dp) :: dbox(3)
1433 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: box_data
1434 : INTEGER :: handle, igrid_level, nbox(3), gauge
1435 28 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: box_ptr
1436 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1437 28 : POINTER :: rs_descs
1438 : TYPE(pw_env_type), POINTER :: pw_env
1439 28 : TYPE(realspace_grid_type), DIMENSION(:, :), POINTER :: rs_gauge
1440 :
1441 28 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
1442 : LOGICAL :: use_old_gauge_atom
1443 :
1444 28 : NULLIFY (rs_gauge, box)
1445 :
1446 28 : CALL timeset(routineN, handle)
1447 :
1448 : CALL get_current_env(current_env=current_env, &
1449 : use_old_gauge_atom=use_old_gauge_atom, &
1450 : rs_gauge=rs_gauge, &
1451 28 : gauge=gauge)
1452 :
1453 28 : IF (gauge == current_gauge_atom) THEN
1454 : CALL get_qs_env(qs_env=qs_env, &
1455 28 : pw_env=pw_env)
1456 : CALL pw_env_get(pw_env=pw_env, &
1457 28 : rs_descs=rs_descs)
1458 : !
1459 : ! box the atoms
1460 28 : IF (use_old_gauge_atom) THEN
1461 22 : CALL box_atoms(qs_env)
1462 : ELSE
1463 6 : CALL box_atoms_new(current_env, qs_env, box)
1464 : END IF
1465 : !
1466 : ! allocate and build the gauge
1467 136 : DO igrid_level = pw_env%gridlevel_info%ngrid_levels, 1, -1
1468 :
1469 432 : DO idir = 1, 3
1470 432 : CALL rs_grid_create(rs_gauge(idir, igrid_level), rs_descs(igrid_level)%rs_desc)
1471 : END DO
1472 :
1473 136 : IF (use_old_gauge_atom) THEN
1474 : CALL collocate_gauge(current_env, qs_env, &
1475 : rs_gauge(1, igrid_level), &
1476 : rs_gauge(2, igrid_level), &
1477 88 : rs_gauge(3, igrid_level))
1478 : ELSE
1479 : CALL collocate_gauge_new(current_env, qs_env, &
1480 : rs_gauge(1, igrid_level), &
1481 : rs_gauge(2, igrid_level), &
1482 : rs_gauge(3, igrid_level), &
1483 20 : box)
1484 : END IF
1485 : END DO
1486 : !
1487 : ! allocate the buf
1488 724 : ALLOCATE (current_env%rs_buf(pw_env%gridlevel_info%ngrid_levels))
1489 136 : DO igrid_level = 1, pw_env%gridlevel_info%ngrid_levels
1490 136 : CALL rs_grid_create(current_env%rs_buf(igrid_level), rs_descs(igrid_level)%rs_desc)
1491 : END DO
1492 : !
1493 28 : DEALLOCATE (box_ptr, box_data)
1494 28 : CALL deallocate_box(box)
1495 : END IF
1496 :
1497 56 : CALL timestop(handle)
1498 :
1499 : CONTAINS
1500 :
1501 : ! **************************************************************************************************
1502 : !> \brief ...
1503 : !> \param qs_env ...
1504 : ! **************************************************************************************************
1505 22 : SUBROUTINE box_atoms(qs_env)
1506 : TYPE(qs_environment_type), POINTER :: qs_env
1507 :
1508 : REAL(kind=dp), PARAMETER :: box_size_guess = 5.0_dp
1509 :
1510 : INTEGER :: i, iatom, ibox, ii, jbox, kbox, natms
1511 : REAL(dp) :: offset(3)
1512 22 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom
1513 : TYPE(cell_type), POINTER :: cell
1514 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1515 22 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1516 :
1517 : CALL get_qs_env(qs_env=qs_env, &
1518 : qs_kind_set=qs_kind_set, &
1519 : cell=cell, &
1520 22 : particle_set=particle_set)
1521 :
1522 22 : natms = SIZE(particle_set, 1)
1523 66 : ALLOCATE (ratom(3, natms))
1524 : !
1525 : ! box the atoms
1526 22 : nbox(1) = CEILING(cell%hmat(1, 1)/box_size_guess)
1527 22 : nbox(2) = CEILING(cell%hmat(2, 2)/box_size_guess)
1528 22 : nbox(3) = CEILING(cell%hmat(3, 3)/box_size_guess)
1529 : !write(*,*) 'nbox',nbox
1530 22 : dbox(1) = cell%hmat(1, 1)/REAL(nbox(1), dp)
1531 22 : dbox(2) = cell%hmat(2, 2)/REAL(nbox(2), dp)
1532 22 : dbox(3) = cell%hmat(3, 3)/REAL(nbox(3), dp)
1533 : !write(*,*) 'dbox',dbox
1534 132 : ALLOCATE (box_ptr(0:nbox(1), 0:nbox(2) - 1, 0:nbox(3) - 1), box_data(3, natms))
1535 390 : box_data(:, :) = HUGE(0.0_dp)
1536 418 : box_ptr(:, :, :) = HUGE(0)
1537 : !
1538 22 : offset(1) = cell%hmat(1, 1)*0.5_dp
1539 22 : offset(2) = cell%hmat(2, 2)*0.5_dp
1540 22 : offset(3) = cell%hmat(3, 3)*0.5_dp
1541 114 : DO iatom = 1, natms
1542 390 : ratom(:, iatom) = pbc(particle_set(iatom)%r(:), cell) + offset(:)
1543 : END DO
1544 : !
1545 22 : i = 1
1546 66 : DO kbox = 0, nbox(3) - 1
1547 154 : DO jbox = 0, nbox(2) - 1
1548 88 : box_ptr(0, jbox, kbox) = i
1549 308 : DO ibox = 0, nbox(1) - 1
1550 : ii = 0
1551 912 : DO iatom = 1, natms
1552 : IF (INT(ratom(1, iatom)/dbox(1)) == ibox .AND. &
1553 736 : INT(ratom(2, iatom)/dbox(2)) == jbox .AND. &
1554 176 : INT(ratom(3, iatom)/dbox(3)) == kbox) THEN
1555 368 : box_data(:, i) = ratom(:, iatom) - offset(:)
1556 92 : i = i + 1
1557 92 : ii = ii + 1
1558 : END IF
1559 : END DO
1560 264 : box_ptr(ibox + 1, jbox, kbox) = box_ptr(ibox, jbox, kbox) + ii
1561 : END DO
1562 : END DO
1563 : END DO
1564 : !
1565 : IF (.FALSE.) THEN
1566 : DO kbox = 0, nbox(3) - 1
1567 : DO jbox = 0, nbox(2) - 1
1568 : DO ibox = 0, nbox(1) - 1
1569 : WRITE (*, *) 'box=', ibox, jbox, kbox
1570 : WRITE (*, *) 'nbr atom=', box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox)
1571 : DO iatom = box_ptr(ibox, jbox, kbox), box_ptr(ibox + 1, jbox, kbox) - 1
1572 : WRITE (*, *) 'iatom=', iatom
1573 : WRITE (*, '(A,3E14.6)') 'coor=', box_data(:, iatom)
1574 : END DO
1575 : END DO
1576 : END DO
1577 : END DO
1578 : END IF
1579 22 : DEALLOCATE (ratom)
1580 22 : END SUBROUTINE box_atoms
1581 :
1582 : ! **************************************************************************************************
1583 : !> \brief ...
1584 : !> \param current_env ...
1585 : !> \param qs_env ...
1586 : !> \param rs_grid_x ...
1587 : !> \param rs_grid_y ...
1588 : !> \param rs_grid_z ...
1589 : ! **************************************************************************************************
1590 88 : SUBROUTINE collocate_gauge(current_env, qs_env, rs_grid_x, rs_grid_y, rs_grid_z)
1591 : !
1592 : TYPE(current_env_type) :: current_env
1593 : TYPE(qs_environment_type), POINTER :: qs_env
1594 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid_x, rs_grid_y, rs_grid_z
1595 :
1596 : INTEGER :: i, iatom, ibeg, ibox, iend, imax, imin, &
1597 : j, jatom, jbox, jmax, jmin, k, kbox, &
1598 : kmax, kmin, lb(3), lb_local(3), natms, &
1599 : natms_local, ng(3)
1600 : REAL(KIND=dp) :: ab, buf_tmp, dist, dr(3), &
1601 : gauge_atom_radius, offset(3), pa, pb, &
1602 : point(3), pra(3), r(3), res(3), summe, &
1603 : tmp, x, y, z
1604 88 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, nrm_atms_pnt
1605 88 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atms_pnt, ratom
1606 88 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid_x, grid_y, grid_z
1607 : TYPE(cell_type), POINTER :: cell
1608 88 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1609 88 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1610 :
1611 : !
1612 :
1613 : CALL get_current_env(current_env=current_env, &
1614 88 : gauge_atom_radius=gauge_atom_radius)
1615 : !
1616 : CALL get_qs_env(qs_env=qs_env, &
1617 : qs_kind_set=qs_kind_set, &
1618 : cell=cell, &
1619 88 : particle_set=particle_set)
1620 : !
1621 88 : natms = SIZE(particle_set, 1)
1622 88 : dr(1) = rs_grid_x%desc%dh(1, 1)
1623 88 : dr(2) = rs_grid_x%desc%dh(2, 2)
1624 88 : dr(3) = rs_grid_x%desc%dh(3, 3)
1625 352 : lb(:) = rs_grid_x%desc%lb(:)
1626 352 : lb_local(:) = rs_grid_x%lb_local(:)
1627 88 : grid_x => rs_grid_x%r(:, :, :)
1628 88 : grid_y => rs_grid_y%r(:, :, :)
1629 88 : grid_z => rs_grid_z%r(:, :, :)
1630 352 : ng(:) = UBOUND(grid_x)
1631 88 : offset(1) = cell%hmat(1, 1)*0.5_dp
1632 88 : offset(2) = cell%hmat(2, 2)*0.5_dp
1633 88 : offset(3) = cell%hmat(3, 3)*0.5_dp
1634 616 : ALLOCATE (buf(natms), ratom(3, natms), atms_pnt(3, natms), nrm_atms_pnt(natms))
1635 : !
1636 : ! go over the grid
1637 1350 : DO k = 1, ng(3)
1638 26188 : DO j = 1, ng(2)
1639 618230 : DO i = 1, ng(1)
1640 : !
1641 592130 : point(3) = REAL(k - 1 + lb_local(3) - lb(3), dp)*dr(3)
1642 592130 : point(2) = REAL(j - 1 + lb_local(2) - lb(2), dp)*dr(2)
1643 592130 : point(1) = REAL(i - 1 + lb_local(1) - lb(1), dp)*dr(1)
1644 2368520 : point = pbc(point, cell)
1645 : !
1646 : ! run over the overlaping boxes
1647 592130 : natms_local = 0
1648 592130 : kmin = INT((point(3) + offset(3) - gauge_atom_radius)/dbox(3))
1649 592130 : kmax = INT((point(3) + offset(3) + gauge_atom_radius)/dbox(3))
1650 592130 : IF (kmax - kmin + 1 > nbox(3)) THEN
1651 527618 : kmin = 0
1652 527618 : kmax = nbox(3) - 1
1653 : END IF
1654 1735750 : DO kbox = kmin, kmax
1655 1143620 : jmin = INT((point(2) + offset(2) - gauge_atom_radius)/dbox(2))
1656 1143620 : jmax = INT((point(2) + offset(2) + gauge_atom_radius)/dbox(2))
1657 1143620 : IF (jmax - jmin + 1 > nbox(2)) THEN
1658 1055236 : jmin = 0
1659 1055236 : jmax = nbox(2) - 1
1660 : END IF
1661 3967326 : DO jbox = jmin, jmax
1662 2231576 : imin = INT((point(1) + offset(1) - gauge_atom_radius)/dbox(1))
1663 2231576 : imax = INT((point(1) + offset(1) + gauge_atom_radius)/dbox(1))
1664 2231576 : IF (imax - imin + 1 > nbox(1)) THEN
1665 2110472 : imin = 0
1666 2110472 : imax = nbox(1) - 1
1667 : END IF
1668 7762096 : DO ibox = imin, imax
1669 4386900 : ibeg = box_ptr(MODULO(ibox, nbox(1)), MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3)))
1670 4386900 : iend = box_ptr(MODULO(ibox, nbox(1)) + 1, MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3))) - 1
1671 9170920 : DO iatom = ibeg, iend
1672 20419552 : r(:) = pbc(box_data(:, iatom) - point(:), cell) + point(:)
1673 2552444 : dist = (r(1) - point(1))**2 + (r(2) - point(2))**2 + (r(3) - point(3))**2
1674 6939344 : IF (dist < gauge_atom_radius**2) THEN
1675 2505134 : natms_local = natms_local + 1
1676 10020536 : ratom(:, natms_local) = r(:)
1677 : !
1678 : ! compute the distance atoms-point
1679 2505134 : x = point(1) - r(1)
1680 2505134 : y = point(2) - r(2)
1681 2505134 : z = point(3) - r(3)
1682 2505134 : atms_pnt(1, natms_local) = x
1683 2505134 : atms_pnt(2, natms_local) = y
1684 2505134 : atms_pnt(3, natms_local) = z
1685 2505134 : nrm_atms_pnt(natms_local) = SQRT(x*x + y*y + z*z)
1686 : END IF
1687 : END DO
1688 : END DO
1689 : END DO
1690 : END DO
1691 : !
1692 616968 : IF (natms_local > 0) THEN
1693 : !
1694 : !
1695 3034004 : DO iatom = 1, natms_local
1696 2505134 : buf_tmp = 1.0_dp
1697 2505134 : pra(1) = atms_pnt(1, iatom)
1698 2505134 : pra(2) = atms_pnt(2, iatom)
1699 2505134 : pra(3) = atms_pnt(3, iatom)
1700 2505134 : pa = nrm_atms_pnt(iatom)
1701 15477460 : DO jatom = 1, natms_local
1702 12972326 : IF (iatom == jatom) CYCLE
1703 10467192 : pb = nrm_atms_pnt(jatom)
1704 10467192 : x = pra(1) - atms_pnt(1, jatom)
1705 10467192 : y = pra(2) - atms_pnt(2, jatom)
1706 10467192 : z = pra(3) - atms_pnt(3, jatom)
1707 10467192 : ab = SQRT(x*x + y*y + z*z)
1708 : !
1709 10467192 : tmp = (pa - pb)/ab
1710 10467192 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1711 10467192 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1712 10467192 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1713 15477460 : buf_tmp = buf_tmp*0.5_dp*(1.0_dp - tmp)
1714 : END DO
1715 3034004 : buf(iatom) = buf_tmp
1716 : END DO
1717 : res(1) = 0.0_dp
1718 : res(2) = 0.0_dp
1719 : res(3) = 0.0_dp
1720 : summe = 0.0_dp
1721 3034004 : DO iatom = 1, natms_local
1722 2505134 : res(1) = res(1) + ratom(1, iatom)*buf(iatom)
1723 2505134 : res(2) = res(2) + ratom(2, iatom)*buf(iatom)
1724 2505134 : res(3) = res(3) + ratom(3, iatom)*buf(iatom)
1725 3034004 : summe = summe + buf(iatom)
1726 : END DO
1727 528870 : res(1) = res(1)/summe
1728 528870 : res(2) = res(2)/summe
1729 528870 : res(3) = res(3)/summe
1730 528870 : grid_x(i, j, k) = point(1) - res(1)
1731 528870 : grid_y(i, j, k) = point(2) - res(2)
1732 528870 : grid_z(i, j, k) = point(3) - res(3)
1733 : ELSE
1734 63260 : grid_x(i, j, k) = 0.0_dp
1735 63260 : grid_y(i, j, k) = 0.0_dp
1736 63260 : grid_z(i, j, k) = 0.0_dp
1737 : END IF
1738 : END DO
1739 : END DO
1740 : END DO
1741 :
1742 88 : DEALLOCATE (buf, ratom, atms_pnt, nrm_atms_pnt)
1743 :
1744 88 : END SUBROUTINE collocate_gauge
1745 :
1746 : ! **************************************************************************************************
1747 : !> \brief ...
1748 : !> \param current_env ...
1749 : !> \param qs_env ...
1750 : !> \param box ...
1751 : ! **************************************************************************************************
1752 6 : SUBROUTINE box_atoms_new(current_env, qs_env, box)
1753 : TYPE(current_env_type) :: current_env
1754 : TYPE(qs_environment_type), POINTER :: qs_env
1755 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
1756 :
1757 : CHARACTER(LEN=*), PARAMETER :: routineN = 'box_atoms_new'
1758 :
1759 : INTEGER :: handle, i, iatom, ibeg, ibox, iend, &
1760 : ifind, ii, imax, imin, j, jatom, jbox, &
1761 : jmax, jmin, k, kbox, kmax, kmin, &
1762 : natms, natms_local
1763 : REAL(dp) :: gauge_atom_radius, offset(3), scale
1764 6 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom
1765 6 : REAL(dp), DIMENSION(:, :), POINTER :: r_ptr
1766 : REAL(kind=dp) :: box_center(3), box_center_wrap(3), &
1767 : box_size_guess, r(3)
1768 : TYPE(cell_type), POINTER :: cell
1769 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1770 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1771 :
1772 6 : CALL timeset(routineN, handle)
1773 :
1774 : CALL get_qs_env(qs_env=qs_env, &
1775 : qs_kind_set=qs_kind_set, &
1776 : cell=cell, &
1777 6 : particle_set=particle_set)
1778 :
1779 : CALL get_current_env(current_env=current_env, &
1780 6 : gauge_atom_radius=gauge_atom_radius)
1781 :
1782 6 : scale = 2.0_dp
1783 :
1784 6 : box_size_guess = gauge_atom_radius/scale
1785 :
1786 6 : natms = SIZE(particle_set, 1)
1787 18 : ALLOCATE (ratom(3, natms))
1788 :
1789 : !
1790 : ! box the atoms
1791 6 : nbox(1) = CEILING(cell%hmat(1, 1)/box_size_guess)
1792 6 : nbox(2) = CEILING(cell%hmat(2, 2)/box_size_guess)
1793 6 : nbox(3) = CEILING(cell%hmat(3, 3)/box_size_guess)
1794 6 : dbox(1) = cell%hmat(1, 1)/REAL(nbox(1), dp)
1795 6 : dbox(2) = cell%hmat(2, 2)/REAL(nbox(2), dp)
1796 6 : dbox(3) = cell%hmat(3, 3)/REAL(nbox(3), dp)
1797 36 : ALLOCATE (box_ptr(0:nbox(1), 0:nbox(2) - 1, 0:nbox(3) - 1), box_data(3, natms))
1798 62 : box_data(:, :) = HUGE(0.0_dp)
1799 27054 : box_ptr(:, :, :) = HUGE(0)
1800 : !
1801 6 : offset(1) = cell%hmat(1, 1)*0.5_dp
1802 6 : offset(2) = cell%hmat(2, 2)*0.5_dp
1803 6 : offset(3) = cell%hmat(3, 3)*0.5_dp
1804 20 : DO iatom = 1, natms
1805 20 : ratom(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
1806 : END DO
1807 : !
1808 6 : i = 1
1809 88 : DO kbox = 0, nbox(3) - 1
1810 1430 : DO jbox = 0, nbox(2) - 1
1811 1342 : box_ptr(0, jbox, kbox) = i
1812 25706 : DO ibox = 0, nbox(1) - 1
1813 24282 : ii = 0
1814 88846 : DO iatom = 1, natms
1815 : IF (MODULO(FLOOR(ratom(1, iatom)/dbox(1)), nbox(1)) == ibox .AND. &
1816 64564 : MODULO(FLOOR(ratom(2, iatom)/dbox(2)), nbox(2)) == jbox .AND. &
1817 24282 : MODULO(FLOOR(ratom(3, iatom)/dbox(3)), nbox(3)) == kbox) THEN
1818 56 : box_data(:, i) = ratom(:, iatom)
1819 14 : i = i + 1
1820 14 : ii = ii + 1
1821 : END IF
1822 : END DO
1823 25624 : box_ptr(ibox + 1, jbox, kbox) = box_ptr(ibox, jbox, kbox) + ii
1824 : END DO
1825 : END DO
1826 : END DO
1827 : !
1828 : IF (.FALSE.) THEN
1829 : DO kbox = 0, nbox(3) - 1
1830 : DO jbox = 0, nbox(2) - 1
1831 : DO ibox = 0, nbox(1) - 1
1832 : IF (box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox) > 0) THEN
1833 : WRITE (*, *) 'box=', ibox, jbox, kbox
1834 : WRITE (*, *) 'nbr atom=', box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox)
1835 : DO iatom = box_ptr(ibox, jbox, kbox), box_ptr(ibox + 1, jbox, kbox) - 1
1836 : WRITE (*, '(A,I3,3E14.6)') 'coor=', iatom, box_data(:, iatom)
1837 : END DO
1838 : END IF
1839 : END DO
1840 : END DO
1841 : END DO
1842 : END IF
1843 : !
1844 6 : NULLIFY (box)
1845 25736 : ALLOCATE (box(0:nbox(1) - 1, 0:nbox(2) - 1, 0:nbox(3) - 1))
1846 : !
1847 : ! build the list
1848 88 : DO k = 0, nbox(3) - 1
1849 1430 : DO j = 0, nbox(2) - 1
1850 25706 : DO i = 0, nbox(1) - 1
1851 : !
1852 24282 : box_center(1) = (REAL(i, dp) + 0.5_dp)*dbox(1)
1853 24282 : box_center(2) = (REAL(j, dp) + 0.5_dp)*dbox(2)
1854 24282 : box_center(3) = (REAL(k, dp) + 0.5_dp)*dbox(3)
1855 24282 : box_center_wrap = pbc(box_center, cell)
1856 : !
1857 : ! find the atoms that are in the overlaping boxes
1858 24282 : natms_local = 0
1859 24282 : kmin = FLOOR((box_center(3) - gauge_atom_radius)/dbox(3))
1860 24282 : kmax = FLOOR((box_center(3) + gauge_atom_radius)/dbox(3))
1861 24282 : IF (kmax - kmin + 1 > nbox(3)) THEN
1862 0 : kmin = 0
1863 0 : kmax = nbox(3) - 1
1864 : END IF
1865 145692 : DO kbox = kmin, kmax
1866 121410 : jmin = FLOOR((box_center(2) - gauge_atom_radius)/dbox(2))
1867 121410 : jmax = FLOOR((box_center(2) + gauge_atom_radius)/dbox(2))
1868 121410 : IF (jmax - jmin + 1 > nbox(2)) THEN
1869 450 : jmin = 0
1870 450 : jmax = nbox(2) - 1
1871 : END IF
1872 751842 : DO jbox = jmin, jmax
1873 606150 : imin = FLOOR((box_center(1) - gauge_atom_radius)/dbox(1))
1874 606150 : imax = FLOOR((box_center(1) + gauge_atom_radius)/dbox(1))
1875 606150 : IF (imax - imin + 1 > nbox(1)) THEN
1876 1350 : imin = 0
1877 1350 : imax = nbox(1) - 1
1878 : END IF
1879 3755610 : DO ibox = imin, imax
1880 3028050 : ibeg = box_ptr(MODULO(ibox, nbox(1)), MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3)))
1881 3028050 : iend = box_ptr(MODULO(ibox, nbox(1)) + 1, MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3))) - 1
1882 3635630 : DO iatom = ibeg, iend
1883 5720 : r = pbc(box_center_wrap(:) - box_data(:, iatom), cell)
1884 : IF (ABS(r(1)) <= (scale + 0.5_dp)*dbox(1) .AND. &
1885 1430 : ABS(r(2)) <= (scale + 0.5_dp)*dbox(2) .AND. &
1886 3028050 : ABS(r(3)) <= (scale + 0.5_dp)*dbox(3)) THEN
1887 1430 : natms_local = natms_local + 1
1888 5720 : ratom(:, natms_local) = box_data(:, iatom)
1889 : END IF
1890 : END DO
1891 : END DO ! box
1892 : END DO
1893 : END DO
1894 : !
1895 : ! set the list
1896 24282 : box(i, j, k)%n = natms_local
1897 24282 : NULLIFY (box(i, j, k)%r)
1898 25624 : IF (natms_local > 0) THEN
1899 3840 : ALLOCATE (box(i, j, k)%r(3, natms_local))
1900 1280 : r_ptr => box(i, j, k)%r
1901 1280 : CALL dcopy(3*natms_local, ratom(1, 1), 1, r_ptr(1, 1), 1)
1902 : END IF
1903 : END DO ! list
1904 : END DO
1905 : END DO
1906 :
1907 : IF (.FALSE.) THEN
1908 : DO k = 0, nbox(3) - 1
1909 : DO j = 0, nbox(2) - 1
1910 : DO i = 0, nbox(1) - 1
1911 : IF (box(i, j, k)%n > 0) THEN
1912 : WRITE (*, *)
1913 : WRITE (*, *) 'box=', i, j, k
1914 : box_center(1) = (REAL(i, dp) + 0.5_dp)*dbox(1)
1915 : box_center(2) = (REAL(j, dp) + 0.5_dp)*dbox(2)
1916 : box_center(3) = (REAL(k, dp) + 0.5_dp)*dbox(3)
1917 : box_center = pbc(box_center, cell)
1918 : WRITE (*, '(A,3E14.6)') 'box_center=', box_center
1919 : WRITE (*, *) 'nbr atom=', box(i, j, k)%n
1920 : r_ptr => box(i, j, k)%r
1921 : DO iatom = 1, box(i, j, k)%n
1922 : WRITE (*, '(A,I3,3E14.6)') 'coor=', iatom, r_ptr(:, iatom)
1923 : r(:) = pbc(box_center(:) - r_ptr(:, iatom), cell)
1924 : IF (ABS(r(1)) > (scale + 0.5_dp)*dbox(1) .OR. &
1925 : ABS(r(2)) > (scale + 0.5_dp)*dbox(2) .OR. &
1926 : ABS(r(3)) > (scale + 0.5_dp)*dbox(3)) THEN
1927 : WRITE (*, *) 'error too many atoms'
1928 : WRITE (*, *) 'dist=', ABS(r(:))
1929 : WRITE (*, *) 'large_dist=', (scale + 0.5_dp)*dbox
1930 : CPABORT("Some atoms are too far from box center")
1931 : END IF
1932 : END DO
1933 : END IF
1934 : END DO ! list
1935 : END DO
1936 : END DO
1937 : END IF
1938 :
1939 : IF (.TRUE.) THEN
1940 88 : DO k = 0, nbox(3) - 1
1941 1430 : DO j = 0, nbox(2) - 1
1942 25706 : DO i = 0, nbox(1) - 1
1943 24282 : box_center(1) = (REAL(i, dp) + 0.5_dp)*dbox(1)
1944 24282 : box_center(2) = (REAL(j, dp) + 0.5_dp)*dbox(2)
1945 24282 : box_center(3) = (REAL(k, dp) + 0.5_dp)*dbox(3)
1946 97128 : box_center = pbc(box_center, cell)
1947 24282 : r_ptr => box(i, j, k)%r
1948 90188 : DO iatom = 1, natms
1949 258256 : r(:) = pbc(box_center(:) - ratom(:, iatom), cell)
1950 64564 : ifind = 0
1951 68174 : DO jatom = 1, box(i, j, k)%n
1952 79004 : IF (SUM(ABS(ratom(:, iatom) - r_ptr(:, jatom))) < 1E-10_dp) ifind = 1
1953 : END DO
1954 :
1955 88846 : IF (ifind == 0) THEN
1956 : ! SQRT(DOT_PRODUCT(r, r)) < gauge_atom_radius
1957 252536 : IF (DOT_PRODUCT(r, r) < (gauge_atom_radius*gauge_atom_radius)) THEN
1958 0 : WRITE (*, *) 'error atom too close'
1959 0 : WRITE (*, *) 'iatom', iatom
1960 0 : WRITE (*, *) 'box_center', box_center
1961 0 : WRITE (*, *) 'ratom', ratom(:, iatom)
1962 0 : WRITE (*, *) 'gauge_atom_radius', gauge_atom_radius
1963 0 : CPABORT("Some atoms are too close within gauge_atom_radius")
1964 : END IF
1965 : END IF
1966 : END DO
1967 : END DO ! list
1968 : END DO
1969 : END DO
1970 : END IF
1971 :
1972 6 : DEALLOCATE (ratom)
1973 :
1974 6 : CALL timestop(handle)
1975 :
1976 12 : END SUBROUTINE box_atoms_new
1977 :
1978 : ! **************************************************************************************************
1979 : !> \brief ...
1980 : !> \param current_env ...
1981 : !> \param qs_env ...
1982 : !> \param rs_grid_x ...
1983 : !> \param rs_grid_y ...
1984 : !> \param rs_grid_z ...
1985 : !> \param box ...
1986 : ! **************************************************************************************************
1987 20 : SUBROUTINE collocate_gauge_new(current_env, qs_env, rs_grid_x, rs_grid_y, rs_grid_z, box)
1988 : !
1989 : TYPE(current_env_type) :: current_env
1990 : TYPE(qs_environment_type), POINTER :: qs_env
1991 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid_x, rs_grid_y, rs_grid_z
1992 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
1993 :
1994 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_gauge_new'
1995 :
1996 : INTEGER :: box_end(3), box_start(3), delta_lb(3), grid_end(3), grid_start(3), handle, i, &
1997 : iatom, ib, ibox, j, jatom, jb, jbox, k, kb, kbox, lb(3), lb_local(3), natms, &
1998 : natms_local0, natms_local1, ng(3)
1999 : INTEGER(KIND=int_8) :: box_lower(3), box_upper(3), &
2000 : grid_lower(3), grid_upper(3), npts(3)
2001 20 : REAL(dp), DIMENSION(:, :), POINTER :: r_ptr
2002 : REAL(KIND=dp) :: ab, box_center(3), buf_tmp, dist, dr(3), &
2003 : gauge_atom_radius, offset(3), pa, pb, &
2004 : point(3), pra(3), r(3), res(3), summe, &
2005 : tmp, x, y, z
2006 20 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, nrm_atms_pnt
2007 20 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atms_pnt, ratom
2008 20 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid_x, grid_y, grid_z
2009 : TYPE(cell_type), POINTER :: cell
2010 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2011 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2012 :
2013 20 : CALL timeset(routineN, handle)
2014 :
2015 : !
2016 : CALL get_current_env(current_env=current_env, &
2017 20 : gauge_atom_radius=gauge_atom_radius)
2018 : !
2019 : CALL get_qs_env(qs_env=qs_env, &
2020 : qs_kind_set=qs_kind_set, &
2021 : cell=cell, &
2022 20 : particle_set=particle_set)
2023 : !
2024 20 : natms = SIZE(particle_set, 1)
2025 20 : dr(1) = rs_grid_x%desc%dh(1, 1)
2026 20 : dr(2) = rs_grid_x%desc%dh(2, 2)
2027 20 : dr(3) = rs_grid_x%desc%dh(3, 3)
2028 80 : lb(:) = rs_grid_x%desc%lb(:)
2029 80 : lb_local(:) = rs_grid_x%lb_local(:)
2030 20 : grid_x => rs_grid_x%r(:, :, :)
2031 20 : grid_y => rs_grid_y%r(:, :, :)
2032 20 : grid_z => rs_grid_z%r(:, :, :)
2033 80 : ng(:) = UBOUND(grid_x)
2034 80 : npts(:) = INT(rs_grid_x%desc%npts(:), KIND=int_8)
2035 80 : delta_lb(:) = lb_local(:) - lb(:)
2036 20 : offset(1) = cell%hmat(1, 1)*0.5_dp
2037 20 : offset(2) = cell%hmat(2, 2)*0.5_dp
2038 20 : offset(3) = cell%hmat(3, 3)*0.5_dp
2039 140 : ALLOCATE (buf(natms), ratom(3, natms), atms_pnt(3, natms), nrm_atms_pnt(natms))
2040 : !
2041 : ! Both the grid and the boxes partition the same cell. Cancel the cell length:
2042 : ! g*dr >= b*dbox is equivalent to g*nbox >= b*npts.
2043 : ! For a positive divisor d, floor(a/d) = (a - MODULO(a,d))/d,
2044 : ! including negative indices in the grid's halo. Use 64-bit products.
2045 80 : grid_lower = INT(delta_lb, KIND=int_8)*nbox
2046 80 : grid_upper = (INT(ng, KIND=int_8) - 1_int_8 + delta_lb)*nbox
2047 80 : box_start = INT((grid_lower - MODULO(grid_lower, npts))/npts)
2048 80 : box_end = INT((grid_upper - MODULO(grid_upper, npts))/npts)
2049 : !
2050 : ! go over the box-list
2051 314 : DO kb = box_start(3), box_end(3)
2052 5148 : DO jb = box_start(2), box_end(2)
2053 89870 : DO ib = box_start(1), box_end(1)
2054 84742 : ibox = MODULO(ib, nbox(1))
2055 84742 : jbox = MODULO(jb, nbox(2))
2056 84742 : kbox = MODULO(kb, nbox(3))
2057 : !
2058 : ! The existing box intervals include both endpoints. Compute ceil(b*npts/nbox)
2059 : ! and floor((b+1)*npts/nbox) exactly, then clip to this rank's local grid.
2060 : ! ceil(a/d) = (a + MODULO(-a,d))/d also handles negative box indices.
2061 338968 : box_lower = INT([ib, jb, kb], KIND=int_8)*npts
2062 338968 : box_upper = (INT([ib, jb, kb], KIND=int_8) + 1_int_8)*npts
2063 338968 : grid_start = INT((box_lower + MODULO(-box_lower, INT(nbox, KIND=int_8)))/nbox)
2064 338968 : grid_end = INT((box_upper - MODULO(box_upper, INT(nbox, KIND=int_8)))/nbox)
2065 338968 : grid_start = MAX(grid_start, delta_lb) - delta_lb + 1
2066 338968 : grid_end = MIN(grid_end, ng - 1 + delta_lb) - delta_lb + 1
2067 : !
2068 : ! Check the box boundaries in the same exact integer coordinates.
2069 338968 : grid_lower = (INT(grid_start, KIND=int_8) - 1_int_8 + delta_lb)*nbox
2070 338968 : grid_upper = (INT(grid_end, KIND=int_8) - 1_int_8 + delta_lb)*nbox
2071 338968 : CPASSERT(ALL(grid_lower >= box_lower))
2072 338968 : CPASSERT(ALL(grid_upper <= box_upper))
2073 : !
2074 : ! the center of the box
2075 84742 : box_center(1) = (REAL(ibox, dp) + 0.5_dp)*dbox(1)
2076 84742 : box_center(2) = (REAL(jbox, dp) + 0.5_dp)*dbox(2)
2077 84742 : box_center(3) = (REAL(kbox, dp) + 0.5_dp)*dbox(3)
2078 : !
2079 : ! find the atoms that are in the overlaping boxes
2080 84742 : natms_local0 = box(ibox, jbox, kbox)%n
2081 84742 : r_ptr => box(ibox, jbox, kbox)%r
2082 : !
2083 : ! go over the grid inside the box
2084 89576 : IF (natms_local0 > 0) THEN
2085 : !
2086 : ! here there are some atoms...
2087 10132 : DO k = grid_start(3), grid_end(3)
2088 32420 : DO j = grid_start(2), grid_end(2)
2089 159542 : DO i = grid_start(1), grid_end(1)
2090 130750 : point(1) = REAL(i - 1 + delta_lb(1), dp)*dr(1)
2091 130750 : point(2) = REAL(j - 1 + delta_lb(2), dp)*dr(2)
2092 130750 : point(3) = REAL(k - 1 + delta_lb(3), dp)*dr(3)
2093 523000 : point = pbc(point, cell)
2094 : !
2095 : ! compute atom-point distances
2096 130750 : natms_local1 = 0
2097 372536 : DO iatom = 1, natms_local0
2098 1934288 : r(:) = pbc(r_ptr(:, iatom) - point(:), cell) + point(:) !needed?
2099 241786 : dist = (r(1) - point(1))**2 + (r(2) - point(2))**2 + (r(3) - point(3))**2
2100 372536 : IF (dist < gauge_atom_radius**2) THEN
2101 154830 : natms_local1 = natms_local1 + 1
2102 619320 : ratom(:, natms_local1) = r(:)
2103 : !
2104 : ! compute the distance atoms-point
2105 154830 : x = point(1) - r(1)
2106 154830 : y = point(2) - r(2)
2107 154830 : z = point(3) - r(3)
2108 154830 : atms_pnt(1, natms_local1) = x
2109 154830 : atms_pnt(2, natms_local1) = y
2110 154830 : atms_pnt(3, natms_local1) = z
2111 154830 : nrm_atms_pnt(natms_local1) = SQRT(x*x + y*y + z*z)
2112 : END IF
2113 : END DO
2114 : !
2115 : !
2116 153038 : IF (natms_local1 > 0) THEN
2117 : !
2118 : ! build the step
2119 240028 : DO iatom = 1, natms_local1
2120 154830 : buf_tmp = 1.0_dp
2121 154830 : pra(1) = atms_pnt(1, iatom)
2122 154830 : pra(2) = atms_pnt(2, iatom)
2123 154830 : pra(3) = atms_pnt(3, iatom)
2124 154830 : pa = nrm_atms_pnt(iatom)
2125 448924 : DO jatom = 1, natms_local1
2126 294094 : IF (iatom == jatom) CYCLE
2127 139264 : pb = nrm_atms_pnt(jatom)
2128 139264 : x = pra(1) - atms_pnt(1, jatom)
2129 139264 : y = pra(2) - atms_pnt(2, jatom)
2130 139264 : z = pra(3) - atms_pnt(3, jatom)
2131 139264 : ab = SQRT(x*x + y*y + z*z)
2132 : !
2133 139264 : tmp = (pa - pb)/ab
2134 139264 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2135 139264 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2136 139264 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2137 448924 : buf_tmp = buf_tmp*0.5_dp*(1.0_dp - tmp)
2138 : END DO
2139 240028 : buf(iatom) = buf_tmp
2140 : END DO
2141 : res(1) = 0.0_dp
2142 : res(2) = 0.0_dp
2143 : res(3) = 0.0_dp
2144 : summe = 0.0_dp
2145 240028 : DO iatom = 1, natms_local1
2146 154830 : res(1) = res(1) + ratom(1, iatom)*buf(iatom)
2147 154830 : res(2) = res(2) + ratom(2, iatom)*buf(iatom)
2148 154830 : res(3) = res(3) + ratom(3, iatom)*buf(iatom)
2149 240028 : summe = summe + buf(iatom)
2150 : END DO
2151 85198 : res(1) = res(1)/summe
2152 85198 : res(2) = res(2)/summe
2153 85198 : res(3) = res(3)/summe
2154 85198 : grid_x(i, j, k) = point(1) - res(1)
2155 85198 : grid_y(i, j, k) = point(2) - res(2)
2156 85198 : grid_z(i, j, k) = point(3) - res(3)
2157 : ELSE
2158 45552 : grid_x(i, j, k) = 0.0_dp
2159 45552 : grid_y(i, j, k) = 0.0_dp
2160 45552 : grid_z(i, j, k) = 0.0_dp
2161 : END IF
2162 : END DO ! grid
2163 : END DO
2164 : END DO
2165 : !
2166 : ELSE
2167 : !
2168 : ! here there is no atom
2169 192148 : DO k = grid_start(3), grid_end(3)
2170 399546 : DO j = grid_start(2), grid_end(2)
2171 803654 : DO i = grid_start(1), grid_end(1)
2172 485222 : grid_x(i, j, k) = 0.0_dp
2173 485222 : grid_y(i, j, k) = 0.0_dp
2174 692620 : grid_z(i, j, k) = 0.0_dp
2175 : END DO ! grid
2176 : END DO
2177 : END DO
2178 : !
2179 : END IF
2180 : !
2181 : END DO ! list
2182 : END DO
2183 : END DO
2184 :
2185 20 : DEALLOCATE (buf, ratom, atms_pnt, nrm_atms_pnt)
2186 :
2187 20 : CALL timestop(handle)
2188 :
2189 20 : END SUBROUTINE collocate_gauge_new
2190 :
2191 : ! **************************************************************************************************
2192 : !> \brief ...
2193 : !> \param box ...
2194 : ! **************************************************************************************************
2195 28 : SUBROUTINE deallocate_box(box)
2196 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
2197 :
2198 : INTEGER :: i, j, k
2199 :
2200 28 : IF (ASSOCIATED(box)) THEN
2201 100 : DO k = LBOUND(box, 3), UBOUND(box, 3)
2202 1594 : DO j = LBOUND(box, 2), UBOUND(box, 2)
2203 28390 : DO i = LBOUND(box, 1), UBOUND(box, 1)
2204 25624 : IF (ASSOCIATED(box(i, j, k)%r)) THEN
2205 1280 : DEALLOCATE (box(i, j, k)%r)
2206 : END IF
2207 : END DO
2208 : END DO
2209 : END DO
2210 6 : DEALLOCATE (box)
2211 : END IF
2212 28 : END SUBROUTINE deallocate_box
2213 : END SUBROUTINE current_set_gauge
2214 :
2215 : ! **************************************************************************************************
2216 : !> \brief ...
2217 : !> \param current_env ...
2218 : !> \param qs_env ...
2219 : !> \param iB ...
2220 : ! **************************************************************************************************
2221 492 : SUBROUTINE current_build_chi(current_env, qs_env, iB)
2222 : !
2223 : TYPE(current_env_type) :: current_env
2224 : TYPE(qs_environment_type), POINTER :: qs_env
2225 : INTEGER, INTENT(IN) :: iB
2226 :
2227 492 : IF (current_env%full) THEN
2228 396 : CALL current_build_chi_many_centers(current_env, qs_env, iB)
2229 96 : ELSE IF (current_env%nbr_center(1) > 1) THEN
2230 0 : CALL current_build_chi_many_centers(current_env, qs_env, iB)
2231 : ELSE
2232 96 : CALL current_build_chi_one_center(current_env, qs_env, iB)
2233 : END IF
2234 :
2235 492 : END SUBROUTINE current_build_chi
2236 :
2237 : ! **************************************************************************************************
2238 : !> \brief ...
2239 : !> \param current_env ...
2240 : !> \param qs_env ...
2241 : !> \param iB ...
2242 : ! **************************************************************************************************
2243 396 : SUBROUTINE current_build_chi_many_centers(current_env, qs_env, iB)
2244 : !
2245 : TYPE(current_env_type) :: current_env
2246 : TYPE(qs_environment_type), POINTER :: qs_env
2247 : INTEGER, INTENT(IN) :: iB
2248 :
2249 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_chi_many_centers'
2250 :
2251 : CHARACTER(LEN=default_string_length) :: basis_type
2252 : INTEGER :: handle, icenter, idir, idir2, ii, iiB, iii, iiiB, ispin, istate, j, jstate, &
2253 : max_states, nao, natom, nbr_center(2), nmo, nspins, nstate_loc, nstates(2), output_unit
2254 396 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
2255 396 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
2256 : LOGICAL :: chi_pbc, gapw
2257 : REAL(dp) :: chi(3), chi_tmp, contrib, contrib2, &
2258 : dk(3), int_current(3), &
2259 : int_current_tmp, maxocc
2260 : TYPE(cell_type), POINTER :: cell
2261 396 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
2262 396 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
2263 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
2264 : TYPE(cp_fm_type) :: psi0, psi_D, psi_p1, psi_p2, psi_rxp
2265 4356 : TYPE(cp_fm_type), DIMENSION(3) :: p_rxp, r_p1, r_p2
2266 36036 : TYPE(cp_fm_type), DIMENSION(9, 3) :: rr_p1, rr_p2, rr_rxp
2267 396 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
2268 396 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_D, psi1_p, psi1_rxp
2269 : TYPE(cp_fm_type), POINTER :: mo_coeff
2270 : TYPE(cp_logger_type), POINTER :: logger
2271 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
2272 396 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_mom_ao, op_p_ao
2273 396 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_mom_der_ao
2274 : TYPE(dft_control_type), POINTER :: dft_control
2275 396 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2276 : TYPE(mp_para_env_type), POINTER :: para_env
2277 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2278 396 : POINTER :: sab_all, sab_orb
2279 396 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2280 396 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2281 :
2282 : !
2283 :
2284 396 : CALL timeset(routineN, handle)
2285 : !
2286 396 : NULLIFY (dft_control, mos, para_env, mo_coeff, op_mom_ao, &
2287 396 : op_mom_der_ao, center_list, centers_set, &
2288 396 : op_p_ao, psi1_p, psi1_rxp, psi1_D, &
2289 396 : cell, particle_set, qs_kind_set)
2290 :
2291 396 : logger => cp_get_default_logger()
2292 396 : output_unit = cp_logger_get_default_io_unit(logger)
2293 :
2294 : CALL get_qs_env(qs_env=qs_env, &
2295 : dft_control=dft_control, &
2296 : mos=mos, &
2297 : para_env=para_env, &
2298 : cell=cell, &
2299 : dbcsr_dist=dbcsr_dist, &
2300 : particle_set=particle_set, &
2301 : qs_kind_set=qs_kind_set, &
2302 : sab_all=sab_all, &
2303 396 : sab_orb=sab_orb)
2304 :
2305 396 : nspins = dft_control%nspins
2306 396 : gapw = dft_control%qs_control%gapw
2307 396 : basis_type = "ORB"
2308 396 : IF (gapw) basis_type = "ORB_SOFT"
2309 :
2310 : CALL get_current_env(current_env=current_env, &
2311 : chi_pbc=chi_pbc, &
2312 : nao=nao, &
2313 : nbr_center=nbr_center, &
2314 : center_list=center_list, &
2315 : centers_set=centers_set, &
2316 : psi1_p=psi1_p, &
2317 : psi1_rxp=psi1_rxp, &
2318 : psi1_D=psi1_D, &
2319 : nstates=nstates, &
2320 396 : psi0_order=psi0_order)
2321 : !
2322 : ! get max nbr of states per center
2323 396 : max_states = 0
2324 978 : DO ispin = 1, nspins
2325 4428 : DO icenter = 1, nbr_center(ispin)
2326 : max_states = MAX(max_states, center_list(ispin)%array(1, icenter + 1)&
2327 4032 : & - center_list(ispin)%array(1, icenter))
2328 : END DO
2329 : END DO
2330 : !
2331 : ! Allocate sparse matrices for dipole, quadrupole and their derivatives => 9x3
2332 : ! Remember the derivatives are antisymmetric
2333 396 : CALL dbcsr_allocate_matrix_set(op_mom_ao, 9)
2334 396 : CALL dbcsr_allocate_matrix_set(op_mom_der_ao, 9, 3)
2335 : !
2336 : ! prepare for allocation
2337 396 : natom = SIZE(particle_set, 1)
2338 1188 : ALLOCATE (first_sgf(natom))
2339 792 : ALLOCATE (last_sgf(natom))
2340 : CALL get_particle_set(particle_set, qs_kind_set, &
2341 : first_sgf=first_sgf, &
2342 396 : last_sgf=last_sgf)
2343 792 : ALLOCATE (row_blk_sizes(natom))
2344 396 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
2345 396 : DEALLOCATE (first_sgf)
2346 396 : DEALLOCATE (last_sgf)
2347 : !
2348 : !
2349 396 : ALLOCATE (op_mom_ao(1)%matrix)
2350 : CALL dbcsr_create(matrix=op_mom_ao(1)%matrix, &
2351 : name="op_mom", &
2352 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
2353 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2354 396 : mutable_work=.TRUE.)
2355 396 : CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix, sab_all)
2356 :
2357 1584 : DO idir2 = 1, 3
2358 1188 : ALLOCATE (op_mom_der_ao(1, idir2)%matrix)
2359 : CALL dbcsr_copy(op_mom_der_ao(1, idir2)%matrix, op_mom_ao(1)%matrix, &
2360 1584 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir2))))
2361 : END DO
2362 :
2363 3564 : DO idir = 2, SIZE(op_mom_ao, 1)
2364 3168 : ALLOCATE (op_mom_ao(idir)%matrix)
2365 : CALL dbcsr_copy(op_mom_ao(idir)%matrix, op_mom_ao(1)%matrix, &
2366 3168 : "op_mom_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2367 13068 : DO idir2 = 1, 3
2368 9504 : ALLOCATE (op_mom_der_ao(idir, idir2)%matrix)
2369 : CALL dbcsr_copy(op_mom_der_ao(idir, idir2)%matrix, op_mom_ao(1)%matrix, &
2370 12672 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir*idir2))))
2371 : END DO
2372 : END DO
2373 : !
2374 396 : CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
2375 396 : ALLOCATE (op_p_ao(1)%matrix)
2376 : CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
2377 : name="op_p_ao", &
2378 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
2379 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2380 396 : mutable_work=.TRUE.)
2381 396 : CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
2382 :
2383 1188 : DO idir = 2, 3
2384 792 : ALLOCATE (op_p_ao(idir)%matrix)
2385 : CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
2386 1188 : "op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2387 : END DO
2388 : !
2389 : !
2390 396 : DEALLOCATE (row_blk_sizes)
2391 : !
2392 : !
2393 : ! Allocate full matrices for only one vector
2394 396 : mo_coeff => psi0_order(1)
2395 396 : NULLIFY (tmp_fm_struct)
2396 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
2397 : ncol_global=max_states, para_env=para_env, &
2398 396 : context=mo_coeff%matrix_struct%context)
2399 396 : CALL cp_fm_create(psi0, tmp_fm_struct)
2400 396 : CALL cp_fm_create(psi_D, tmp_fm_struct)
2401 396 : CALL cp_fm_create(psi_rxp, tmp_fm_struct)
2402 396 : CALL cp_fm_create(psi_p1, tmp_fm_struct)
2403 396 : CALL cp_fm_create(psi_p2, tmp_fm_struct)
2404 396 : CALL cp_fm_struct_release(tmp_fm_struct)
2405 : !
2406 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
2407 : ncol_global=max_states, para_env=para_env, &
2408 396 : context=mo_coeff%matrix_struct%context)
2409 1584 : DO idir = 1, 3
2410 1188 : CALL cp_fm_create(p_rxp(idir), tmp_fm_struct, set_zero=.TRUE.)
2411 1188 : CALL cp_fm_create(r_p1(idir), tmp_fm_struct, set_zero=.TRUE.)
2412 1188 : CALL cp_fm_create(r_p2(idir), tmp_fm_struct, set_zero=.TRUE.)
2413 12276 : DO idir2 = 1, 9
2414 10692 : CALL cp_fm_create(rr_rxp(idir2, idir), tmp_fm_struct, set_zero=.TRUE.)
2415 10692 : CALL cp_fm_create(rr_p1(idir2, idir), tmp_fm_struct, set_zero=.TRUE.)
2416 11880 : CALL cp_fm_create(rr_p2(idir2, idir), tmp_fm_struct, set_zero=.TRUE.)
2417 : END DO
2418 : END DO
2419 396 : CALL cp_fm_struct_release(tmp_fm_struct)
2420 : !
2421 : !
2422 : !
2423 : ! recompute the linear momentum matrices
2424 396 : CALL build_lin_mom_matrix(qs_env, op_p_ao)
2425 : !
2426 : !
2427 : ! get iiB and iiiB
2428 396 : CALL set_vecp(iB, iiB, iiiB)
2429 978 : DO ispin = 1, nspins
2430 : !
2431 : ! get ground state MOS
2432 582 : nmo = nstates(ispin)
2433 582 : mo_coeff => psi0_order(ispin)
2434 582 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
2435 : !
2436 : ! Initialize the temporary vector chi
2437 582 : chi = 0.0_dp
2438 : int_current = 0.0_dp
2439 : !
2440 : ! Start loop over the occupied states
2441 4032 : DO icenter = 1, nbr_center(ispin)
2442 : !
2443 : ! Get the Wannier center of the istate-th ground state orbital
2444 13800 : dk(1:3) = centers_set(ispin)%array(1:3, icenter)
2445 : !
2446 : ! Compute the multipole integrals for the state istate,
2447 : ! using as reference center the corresponding Wannier center
2448 34500 : DO idir = 1, 9
2449 31050 : CALL dbcsr_set(op_mom_ao(idir)%matrix, 0.0_dp)
2450 127650 : DO idir2 = 1, 3
2451 124200 : CALL dbcsr_set(op_mom_der_ao(idir, idir2)%matrix, 0.0_dp)
2452 : END DO
2453 : END DO
2454 : CALL build_local_moments_der_matrix(qs_env, op_mom_der_ao, 2, 2, ref_point=dk, &
2455 : moments=op_mom_ao, basis_type=basis_type, &
2456 3450 : minimum_image=.FALSE., ordered=.TRUE.)
2457 : !
2458 : ! collecte the states that belong to a given center
2459 3450 : CALL cp_fm_set_all(psi0, 0.0_dp)
2460 3450 : CALL cp_fm_set_all(psi_rxp, 0.0_dp)
2461 3450 : CALL cp_fm_set_all(psi_D, 0.0_dp)
2462 3450 : CALL cp_fm_set_all(psi_p1, 0.0_dp)
2463 3450 : CALL cp_fm_set_all(psi_p2, 0.0_dp)
2464 3450 : nstate_loc = center_list(ispin)%array(1, icenter + 1) - center_list(ispin)%array(1, icenter)
2465 3450 : jstate = 1
2466 7278 : DO j = center_list(ispin)%array(1, icenter), center_list(ispin)%array(1, icenter + 1) - 1
2467 3828 : istate = center_list(ispin)%array(2, j)
2468 : !
2469 : ! block the states that belong to this center
2470 3828 : CALL cp_fm_to_fm(mo_coeff, psi0, 1, istate, jstate)
2471 : !
2472 3828 : CALL cp_fm_to_fm(psi1_rxp(ispin, iB), psi_rxp, 1, istate, jstate)
2473 3828 : IF (current_env%full) CALL cp_fm_to_fm(psi1_D(ispin, iB), psi_D, 1, istate, jstate)
2474 : !
2475 : ! psi1_p_iiB_istate and psi1_p_iiiB_istate
2476 3828 : CALL cp_fm_to_fm(psi1_p(ispin, iiB), psi_p1, 1, istate, jstate)
2477 3828 : CALL cp_fm_to_fm(psi1_p(ispin, iiiB), psi_p2, 1, istate, jstate)
2478 : !
2479 7278 : jstate = jstate + 1
2480 : END DO ! istate
2481 : !
2482 : ! scale the ordered mos
2483 3450 : IF (current_env%full) CALL cp_fm_scale_and_add(1.0_dp, psi_rxp, -1.0_dp, psi_D)
2484 : !
2485 13800 : DO idir = 1, 3
2486 10350 : CALL set_vecp(idir, ii, iii)
2487 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix, psi_rxp, &
2488 10350 : p_rxp(idir), ncol=nstate_loc, alpha=1.e0_dp)
2489 10350 : IF (iiiB == iii .OR. iiiB == ii) THEN
2490 : CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, psi_p1, &
2491 6900 : r_p1(idir), ncol=nstate_loc, alpha=1.e0_dp)
2492 : END IF
2493 10350 : IF (iiB == iii .OR. iiB == ii) THEN
2494 : CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, psi_p2, &
2495 6900 : r_p2(idir), ncol=nstate_loc, alpha=1.e0_dp)
2496 : END IF
2497 117300 : DO idir2 = 1, 9
2498 93150 : IF (idir2 == ii .OR. idir2 == iii) THEN
2499 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_rxp, &
2500 20700 : rr_rxp(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2501 : END IF
2502 : !
2503 93150 : IF (idir2 == ind_m2(ii, iiiB) .OR. idir2 == ind_m2(iii, iiiB) .OR. idir2 == iiiB) THEN
2504 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_p1, &
2505 31050 : rr_p1(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2506 : END IF
2507 : !
2508 103500 : IF (idir2 == ind_m2(ii, iiB) .OR. idir2 == ind_m2(iii, iiB) .OR. idir2 == iiB) THEN
2509 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_p2, &
2510 31050 : rr_p2(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2511 : END IF
2512 : END DO
2513 : END DO
2514 : !
2515 : ! Multuply left and right by the appropriate coefficients and sum into the
2516 : ! correct component of the chi tensor using the appropriate multiplicative factor
2517 : ! (don't forget the occupation number)
2518 : ! Loop over the cartesian components of the tensor
2519 : ! The loop over the components of the external field is external, thereby
2520 : ! only one column of the chi tensor is computed here
2521 14382 : DO idir = 1, 3
2522 10350 : chi_tmp = 0.0_dp
2523 10350 : int_current_tmp = 0.0_dp
2524 : !
2525 : ! get ii and iii
2526 10350 : CALL set_vecp(idir, ii, iii)
2527 : !
2528 : ! term: 2[C0| (r-dk)_ii |d_iii(C1(rxp-D))]-2[C0| (r-dk)_iii |d_ii(C1(rxp-D))]
2529 : ! the factor 2 should be already included in the matrix elements
2530 : contrib = 0.0_dp
2531 10350 : CALL cp_fm_trace(psi0, rr_rxp(ii, iii), contrib)
2532 10350 : chi_tmp = chi_tmp + 2.0_dp*contrib
2533 : !
2534 : contrib = 0.0_dp
2535 10350 : CALL cp_fm_trace(psi0, rr_rxp(iii, ii), contrib)
2536 10350 : chi_tmp = chi_tmp - 2.0_dp*contrib
2537 : !
2538 : ! correction: dk_ii*2[C0| d_iii(C1(rxp-D))] - dk_iii*2[C0| d_ii(C1(rxp-D))]
2539 : ! factor 2 not included in the matrix elements
2540 : contrib = 0.0_dp
2541 10350 : CALL cp_fm_trace(psi0, p_rxp(iii), contrib)
2542 10350 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(ii)*contrib
2543 10350 : int_current_tmp = int_current_tmp + 2.0_dp*contrib
2544 : !
2545 : contrib2 = 0.0_dp
2546 10350 : CALL cp_fm_trace(psi0, p_rxp(ii), contrib2)
2547 10350 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(iii)*contrib2
2548 : !
2549 : ! term: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] \
2550 : ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))]
2551 : ! the factor 2 should be already included in the matrix elements
2552 : contrib = 0.0_dp
2553 10350 : idir2 = ind_m2(ii, iiB)
2554 10350 : CALL cp_fm_trace(psi0, rr_p2(idir2, iii), contrib)
2555 10350 : chi_tmp = chi_tmp - 2.0_dp*contrib
2556 : contrib2 = 0.0_dp
2557 10350 : IF (iiB == iii) THEN
2558 3450 : CALL cp_fm_trace(psi0, r_p2(ii), contrib2)
2559 3450 : chi_tmp = chi_tmp - contrib2
2560 : END IF
2561 : !
2562 : contrib = 0.0_dp
2563 10350 : idir2 = ind_m2(iii, iiB)
2564 10350 : CALL cp_fm_trace(psi0, rr_p2(idir2, ii), contrib)
2565 10350 : chi_tmp = chi_tmp + 2.0_dp*contrib
2566 : contrib2 = 0.0_dp
2567 10350 : IF (iiB == ii) THEN
2568 3450 : CALL cp_fm_trace(psi0, r_p2(iii), contrib2)
2569 3450 : chi_tmp = chi_tmp + contrib2
2570 : END IF
2571 : !
2572 : ! correction: -dk_ii * 2[C0|(r-dk)_iiB | d_iii(C1(piiiB))] \
2573 : ! +dk_iii * 2[C0|(r-dk)_iiB | d_ii(C1(piiiB))]
2574 : ! the factor 2 should be already included in the matrix elements
2575 : ! no additional correction terms because of the orthogonality between C0 and C1
2576 : contrib = 0.0_dp
2577 10350 : CALL cp_fm_trace(psi0, rr_p2(iiB, iii), contrib)
2578 10350 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(ii)*contrib
2579 10350 : int_current_tmp = int_current_tmp - 2.0_dp*contrib
2580 : !
2581 : contrib2 = 0.0_dp
2582 10350 : CALL cp_fm_trace(psi0, rr_p2(iiB, ii), contrib2)
2583 10350 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(iii)*contrib2
2584 : !
2585 : ! term: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] \
2586 : ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))]
2587 : ! the factor 2 should be already included in the matrix elements
2588 : contrib = 0.0_dp
2589 10350 : idir2 = ind_m2(ii, iiiB)
2590 10350 : CALL cp_fm_trace(psi0, rr_p1(idir2, iii), contrib)
2591 10350 : chi_tmp = chi_tmp + 2.0_dp*contrib
2592 : contrib2 = 0.0_dp
2593 10350 : IF (iiiB == iii) THEN
2594 3450 : CALL cp_fm_trace(psi0, r_p1(ii), contrib2)
2595 3450 : chi_tmp = chi_tmp + contrib2
2596 : END IF
2597 : !
2598 : contrib = 0.0_dp
2599 10350 : idir2 = ind_m2(iii, iiiB)
2600 10350 : CALL cp_fm_trace(psi0, rr_p1(idir2, ii), contrib)
2601 10350 : chi_tmp = chi_tmp - 2.0_dp*contrib
2602 : contrib2 = 0.0_dp
2603 10350 : IF (iiiB == ii) THEN
2604 3450 : CALL cp_fm_trace(psi0, r_p1(iii), contrib2)
2605 3450 : chi_tmp = chi_tmp - contrib2
2606 : END IF
2607 : !
2608 : ! correction: +dk_ii * 2[C0|(r-dk)_iiiB | d_iii(C1(piiB))] +\
2609 : ! -dk_iii * 2[C0|(r-dk)_iiiB | d_ii(C1(piiB))]
2610 : ! the factor 2 should be already included in the matrix elements
2611 : contrib = 0.0_dp
2612 10350 : CALL cp_fm_trace(psi0, rr_p1(iiiB, iii), contrib)
2613 10350 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(ii)*contrib
2614 10350 : int_current_tmp = int_current_tmp + 2.0_dp*contrib
2615 : !
2616 : contrib2 = 0.0_dp
2617 10350 : CALL cp_fm_trace(psi0, rr_p1(iiiB, ii), contrib2)
2618 10350 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(iii)*contrib2
2619 : !
2620 : ! accumulate
2621 10350 : chi(idir) = chi(idir) + maxocc*chi_tmp
2622 106950 : int_current(iii) = int_current(iii) + int_current_tmp
2623 : END DO ! idir
2624 :
2625 : END DO ! icenter
2626 : !
2627 3306 : DO idir = 1, 3
2628 : current_env%chi_tensor(idir, iB, ispin) = current_env%chi_tensor(idir, iB, ispin) + &
2629 1746 : chi(idir)
2630 582 : IF (output_unit > 0) THEN
2631 : !WRITE(output_unit,'(A,E12.6)') ' chi_'//ACHAR(119+idir)//ACHAR(119+iB)//&
2632 : ! & ' = ',chi(idir)
2633 : !WRITE(output_unit,'(A,E12.6)') ' analytic \int j_'//ACHAR(119+idir)//ACHAR(119+iB)//&
2634 : ! & '(r) d^3r = ',int_current(idir)
2635 : END IF
2636 : END DO
2637 : !
2638 : END DO ! ispin
2639 : !
2640 : ! deallocate the sparse matrices
2641 396 : CALL dbcsr_deallocate_matrix_set(op_mom_ao)
2642 396 : CALL dbcsr_deallocate_matrix_set(op_mom_der_ao)
2643 396 : CALL dbcsr_deallocate_matrix_set(op_p_ao)
2644 :
2645 396 : CALL cp_fm_release(psi0)
2646 396 : CALL cp_fm_release(psi_rxp)
2647 396 : CALL cp_fm_release(psi_D)
2648 396 : CALL cp_fm_release(psi_p1)
2649 396 : CALL cp_fm_release(psi_p2)
2650 1584 : DO idir = 1, 3
2651 1188 : CALL cp_fm_release(p_rxp(idir))
2652 1188 : CALL cp_fm_release(r_p1(idir))
2653 1188 : CALL cp_fm_release(r_p2(idir))
2654 12276 : DO idir2 = 1, 9
2655 10692 : CALL cp_fm_release(rr_rxp(idir2, idir))
2656 10692 : CALL cp_fm_release(rr_p1(idir2, idir))
2657 11880 : CALL cp_fm_release(rr_p2(idir2, idir))
2658 : END DO
2659 : END DO
2660 :
2661 396 : CALL timestop(handle)
2662 :
2663 2376 : END SUBROUTINE current_build_chi_many_centers
2664 :
2665 : ! **************************************************************************************************
2666 : !> \brief ...
2667 : !> \param current_env ...
2668 : !> \param qs_env ...
2669 : !> \param iB ...
2670 : ! **************************************************************************************************
2671 96 : SUBROUTINE current_build_chi_one_center(current_env, qs_env, iB)
2672 : !
2673 : TYPE(current_env_type) :: current_env
2674 : TYPE(qs_environment_type), POINTER :: qs_env
2675 : INTEGER, INTENT(IN) :: iB
2676 :
2677 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_chi_one_center'
2678 :
2679 : CHARACTER(LEN=default_string_length) :: basis_type
2680 : INTEGER :: handle, idir, idir2, iiB, iiiB, ispin, jdir, jjdir, kdir, max_states, nao, natom, &
2681 : nbr_center(2), nmo, nspins, nstates(2), output_unit
2682 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
2683 96 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
2684 : LOGICAL :: chi_pbc, gapw
2685 : REAL(dp) :: chi(3), contrib, dk(3), int_current(3), &
2686 : maxocc
2687 : TYPE(cell_type), POINTER :: cell
2688 96 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
2689 96 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
2690 : TYPE(cp_fm_type) :: buf
2691 96 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
2692 96 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_p, psi1_rxp
2693 : TYPE(cp_fm_type), POINTER :: mo_coeff
2694 : TYPE(cp_logger_type), POINTER :: logger
2695 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
2696 96 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_mom_ao, op_p_ao
2697 96 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_mom_der_ao
2698 : TYPE(dft_control_type), POINTER :: dft_control
2699 96 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2700 : TYPE(mp_para_env_type), POINTER :: para_env
2701 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2702 96 : POINTER :: sab_all, sab_orb
2703 96 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2704 96 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2705 :
2706 : !
2707 :
2708 96 : CALL timeset(routineN, handle)
2709 : !
2710 96 : NULLIFY (dft_control, mos, para_env, mo_coeff, op_mom_ao, &
2711 96 : op_mom_der_ao, center_list, centers_set, &
2712 96 : op_p_ao, psi1_p, psi1_rxp, cell, psi0_order)
2713 :
2714 96 : logger => cp_get_default_logger()
2715 96 : output_unit = cp_logger_get_default_io_unit(logger)
2716 :
2717 : CALL get_qs_env(qs_env=qs_env, &
2718 : dft_control=dft_control, &
2719 : mos=mos, &
2720 : para_env=para_env, &
2721 : cell=cell, &
2722 : particle_set=particle_set, &
2723 : qs_kind_set=qs_kind_set, &
2724 : sab_all=sab_all, &
2725 : sab_orb=sab_orb, &
2726 96 : dbcsr_dist=dbcsr_dist)
2727 :
2728 96 : nspins = dft_control%nspins
2729 96 : gapw = dft_control%qs_control%gapw
2730 96 : basis_type = "ORB"
2731 96 : IF (gapw) basis_type = "ORB_SOFT"
2732 :
2733 : CALL get_current_env(current_env=current_env, &
2734 : chi_pbc=chi_pbc, &
2735 : nao=nao, &
2736 : nbr_center=nbr_center, &
2737 : center_list=center_list, &
2738 : centers_set=centers_set, &
2739 : psi1_p=psi1_p, &
2740 : psi1_rxp=psi1_rxp, &
2741 : nstates=nstates, &
2742 96 : psi0_order=psi0_order)
2743 : !
2744 228 : max_states = MAXVAL(nstates(1:nspins))
2745 : !
2746 : ! Allocate sparse matrices for dipole, quadrupole and their derivatives => 9x3
2747 : ! Remember the derivatives are antisymmetric
2748 96 : CALL dbcsr_allocate_matrix_set(op_mom_ao, 9)
2749 96 : CALL dbcsr_allocate_matrix_set(op_mom_der_ao, 9, 3)
2750 : !
2751 : ! prepare for allocation
2752 96 : natom = SIZE(particle_set, 1)
2753 288 : ALLOCATE (first_sgf(natom))
2754 192 : ALLOCATE (last_sgf(natom))
2755 : CALL get_particle_set(particle_set, qs_kind_set, &
2756 : first_sgf=first_sgf, &
2757 96 : last_sgf=last_sgf)
2758 192 : ALLOCATE (row_blk_sizes(natom))
2759 96 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
2760 96 : DEALLOCATE (first_sgf)
2761 96 : DEALLOCATE (last_sgf)
2762 : !
2763 : !
2764 96 : ALLOCATE (op_mom_ao(1)%matrix)
2765 : CALL dbcsr_create(matrix=op_mom_ao(1)%matrix, &
2766 : name="op_mom", &
2767 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
2768 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2769 96 : mutable_work=.TRUE.)
2770 96 : CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix, sab_all)
2771 :
2772 384 : DO idir2 = 1, 3
2773 288 : ALLOCATE (op_mom_der_ao(1, idir2)%matrix)
2774 : CALL dbcsr_copy(op_mom_der_ao(1, idir2)%matrix, op_mom_ao(1)%matrix, &
2775 384 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir2))))
2776 : END DO
2777 :
2778 864 : DO idir = 2, SIZE(op_mom_ao, 1)
2779 768 : ALLOCATE (op_mom_ao(idir)%matrix)
2780 : CALL dbcsr_copy(op_mom_ao(idir)%matrix, op_mom_ao(1)%matrix, &
2781 768 : "op_mom_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2782 3168 : DO idir2 = 1, 3
2783 2304 : ALLOCATE (op_mom_der_ao(idir, idir2)%matrix)
2784 : CALL dbcsr_copy(op_mom_der_ao(idir, idir2)%matrix, op_mom_ao(1)%matrix, &
2785 3072 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir*idir2))))
2786 : END DO
2787 : END DO
2788 : !
2789 96 : CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
2790 96 : ALLOCATE (op_p_ao(1)%matrix)
2791 : CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
2792 : name="op_p_ao", &
2793 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
2794 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2795 96 : mutable_work=.TRUE.)
2796 96 : CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
2797 :
2798 288 : DO idir = 2, 3
2799 192 : ALLOCATE (op_p_ao(idir)%matrix)
2800 : CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
2801 288 : "op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2802 : END DO
2803 : !
2804 : !
2805 96 : DEALLOCATE (row_blk_sizes)
2806 : !
2807 : ! recompute the linear momentum matrices
2808 96 : CALL build_lin_mom_matrix(qs_env, op_p_ao)
2809 : !
2810 : !
2811 : ! get iiB and iiiB
2812 96 : CALL set_vecp(iB, iiB, iiiB)
2813 228 : DO ispin = 1, nspins
2814 : !
2815 132 : CPASSERT(nbr_center(ispin) == 1)
2816 : !
2817 : ! get ground state MOS
2818 132 : nmo = nstates(ispin)
2819 132 : mo_coeff => psi0_order(ispin)
2820 132 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
2821 : !
2822 : ! Create buffer matrix
2823 132 : CALL cp_fm_create(buf, mo_coeff%matrix_struct)
2824 : !
2825 : ! Initialize the temporary vector chi
2826 132 : chi = 0.0_dp
2827 : int_current = 0.0_dp
2828 : !
2829 : !
2830 : ! Get the Wannier center of the istate-th ground state orbital
2831 528 : dk(1:3) = centers_set(ispin)%array(1:3, 1)
2832 : !
2833 : ! Compute the multipole integrals for the state istate,
2834 : ! using as reference center the corresponding Wannier center
2835 1320 : DO idir = 1, 9
2836 1188 : CALL dbcsr_set(op_mom_ao(idir)%matrix, 0.0_dp)
2837 4884 : DO idir2 = 1, 3
2838 4752 : CALL dbcsr_set(op_mom_der_ao(idir, idir2)%matrix, 0.0_dp)
2839 : END DO
2840 : END DO
2841 : CALL build_local_moments_der_matrix(qs_env, op_mom_der_ao, 2, 2, ref_point=dk, &
2842 : moments=op_mom_ao, basis_type=basis_type, &
2843 132 : minimum_image=.FALSE., ordered=.TRUE.)
2844 : !
2845 : !
2846 : ! Multuply left and right by the appropriate coefficients and sum into the
2847 : ! correct component of the chi tensor using the appropriate multiplicative factor
2848 : ! (don't forget the occupation number)
2849 : ! Loop over the cartesian components of the tensor
2850 : ! The loop over the components of the external field is external, thereby
2851 : ! only one column of the chi tensor is computed here
2852 528 : DO idir = 1, 3
2853 : !
2854 : !
2855 : !
2856 : ! term: dk_ii*2[C0| d_iii(C1(rxp-D))] - dk_iii*2[C0| d_ii(C1(rxp-D))]
2857 396 : IF (.NOT. chi_pbc) THEN
2858 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix, mo_coeff, &
2859 396 : buf, ncol=nmo, alpha=1.e0_dp)
2860 1584 : DO jdir = 1, 3
2861 5148 : DO kdir = 1, 3
2862 3564 : IF (Levi_Civita(kdir, jdir, idir) == 0.0_dp) CYCLE
2863 792 : CALL cp_fm_trace(buf, psi1_rxp(ispin, iB), contrib)
2864 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jdir, idir)*2.0_dp*dk(jdir)*contrib
2865 : END DO
2866 : END DO
2867 : END IF
2868 : !
2869 : !
2870 : !
2871 : ! term: 2[C0| (r-dk)_ii |d_iii(C1(rxp-D))]-2[C0| (r-dk)_iii |d_ii(C1(rxp-D))]
2872 : ! and
2873 : ! term: -dk_ii * 2[C0|(r-dk)_iiB | d_iii(C1(piiiB))] +
2874 : ! +dk_iii * 2[C0|(r-dk)_iiB | d_ii(C1(piiiB))]
2875 : ! and
2876 : ! term: +dk_ii * 2[C0|(r-dk)_iiiB | d_iii(C1(piiB))] +
2877 : ! -dk_iii * 2[C0|(r-dk)_iiiB | d_ii(C1(piiB))]
2878 1584 : DO jdir = 1, 3
2879 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(jdir, idir)%matrix, mo_coeff, &
2880 1188 : buf, ncol=nmo, alpha=1.e0_dp)
2881 4752 : DO kdir = 1, 3
2882 3564 : IF (Levi_Civita(kdir, jdir, idir) == 0.0_dp) CYCLE
2883 792 : CALL cp_fm_trace(buf, psi1_rxp(ispin, iB), contrib)
2884 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jdir, idir)*2.0_dp*contrib
2885 : END DO
2886 : !
2887 1584 : IF (.NOT. chi_pbc) THEN
2888 1188 : IF (jdir == iiB) THEN
2889 1584 : DO jjdir = 1, 3
2890 5148 : DO kdir = 1, 3
2891 3564 : IF (Levi_Civita(kdir, jjdir, idir) == 0.0_dp) CYCLE
2892 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiiB), contrib)
2893 4752 : chi(kdir) = chi(kdir) + Levi_Civita(kdir, jjdir, idir)*2.0_dp*dk(jjdir)*contrib
2894 : END DO
2895 : END DO
2896 : END IF
2897 : !
2898 1188 : IF (jdir == iiiB) THEN
2899 1584 : DO jjdir = 1, 3
2900 5148 : DO kdir = 1, 3
2901 3564 : IF (Levi_Civita(kdir, jjdir, idir) == 0.0_dp) CYCLE
2902 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiB), contrib)
2903 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jjdir, idir)*2.0_dp*dk(jjdir)*contrib
2904 : END DO
2905 : END DO
2906 : END IF
2907 : END IF
2908 : END DO
2909 : !
2910 : !
2911 : !
2912 : ! term1: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] +
2913 : ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))]
2914 : ! and
2915 : ! term1: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] +
2916 : ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))]
2917 : ! HERE THERE IS ONE EXTRA MULTIPLY
2918 1584 : DO jdir = 1, 3
2919 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir, iiB), idir)%matrix, mo_coeff, &
2920 1188 : buf, ncol=nmo, alpha=1.e0_dp)
2921 4752 : DO kdir = 1, 3
2922 3564 : IF (Levi_Civita(kdir, jdir, idir) == 0.0_dp) CYCLE
2923 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiiB), contrib)
2924 4752 : chi(kdir) = chi(kdir) + Levi_Civita(kdir, jdir, idir)*2.0_dp*contrib
2925 : END DO
2926 : !
2927 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir, iiiB), idir)%matrix, mo_coeff, &
2928 1188 : buf, ncol=nmo, alpha=1.e0_dp)
2929 5148 : DO kdir = 1, 3
2930 3564 : IF (Levi_Civita(kdir, jdir, idir) == 0.0_dp) CYCLE
2931 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiB), contrib)
2932 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jdir, idir)*2.0_dp*contrib
2933 : END DO
2934 : END DO
2935 : !
2936 : !
2937 : !
2938 : ! term2: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] +
2939 : ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))]
2940 : ! and
2941 : ! term2: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] +
2942 : ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))]
2943 : CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, mo_coeff, &
2944 396 : buf, ncol=nmo, alpha=1.e0_dp)
2945 1584 : DO jdir = 1, 3
2946 5148 : DO kdir = 1, 3
2947 3564 : IF (Levi_Civita(kdir, idir, jdir) == 0.0_dp) CYCLE
2948 1980 : IF (iiB == jdir) THEN
2949 264 : CALL cp_fm_trace(buf, psi1_p(ispin, iiiB), contrib)
2950 264 : chi(kdir) = chi(kdir) + Levi_Civita(kdir, idir, jdir)*contrib
2951 : END IF
2952 : END DO
2953 : END DO
2954 : !
2955 1716 : DO jdir = 1, 3
2956 5148 : DO kdir = 1, 3
2957 3564 : IF (Levi_Civita(kdir, idir, jdir) == 0.0_dp) CYCLE
2958 1980 : IF (iiiB == jdir) THEN
2959 264 : CALL cp_fm_trace(buf, psi1_p(ispin, iiB), contrib)
2960 264 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, idir, jdir)*contrib
2961 : END IF
2962 : !
2963 : END DO
2964 : END DO
2965 : !
2966 : !
2967 : !
2968 : !
2969 : END DO ! idir
2970 : !
2971 528 : DO idir = 1, 3
2972 : current_env%chi_tensor(idir, iB, ispin) = current_env%chi_tensor(idir, iB, ispin) + &
2973 396 : maxocc*chi(idir)
2974 132 : IF (output_unit > 0) THEN
2975 : !WRITE(output_unit,'(A,E12.6)') ' chi_'//ACHAR(119+idir)//ACHAR(119+iB)//&
2976 : ! & ' = ',maxocc * chi(idir)
2977 : END IF
2978 : END DO
2979 : !
2980 360 : CALL cp_fm_release(buf)
2981 : END DO ! ispin
2982 : !
2983 : ! deallocate the sparse matrices
2984 96 : CALL dbcsr_deallocate_matrix_set(op_mom_ao)
2985 96 : CALL dbcsr_deallocate_matrix_set(op_mom_der_ao)
2986 96 : CALL dbcsr_deallocate_matrix_set(op_p_ao)
2987 :
2988 96 : CALL timestop(handle)
2989 :
2990 192 : END SUBROUTINE current_build_chi_one_center
2991 :
2992 0 : END MODULE qs_linres_current
|
