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