Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief given the response wavefunctions obtained by the application
10 : !> of the (rxp), p, and ((dk-dl)xp) operators,
11 : !> here the current density vector (jx, jy, jz)
12 : !> is computed for the 3 directions of the magnetic field (Bx, By, Bz)
13 : !> \par History
14 : !> created 02-2006 [MI]
15 : !> \author MI
16 : ! **************************************************************************************************
17 : MODULE qs_linres_current
18 : USE ao_util, ONLY: exp_radius_very_extended
19 : USE basis_set_types, ONLY: get_gto_basis_set,&
20 : gto_basis_set_p_type,&
21 : gto_basis_set_type
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_array_utils, ONLY: cp_2d_i_p_type,&
25 : cp_2d_r_p_type
26 : USE cp_control_types, ONLY: dft_control_type
27 : USE cp_dbcsr_api, ONLY: &
28 : dbcsr_convert_offsets_to_sizes, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, &
29 : dbcsr_distribution_type, dbcsr_finalize, dbcsr_get_block_p, dbcsr_p_type, dbcsr_put_block, &
30 : dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, dbcsr_type_no_symmetry
31 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
32 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
33 : cp_dbcsr_sm_fm_multiply,&
34 : dbcsr_allocate_matrix_set,&
35 : dbcsr_deallocate_matrix_set
36 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
37 : cp_fm_trace
38 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
39 : cp_fm_struct_release,&
40 : cp_fm_struct_type
41 : USE cp_fm_types, ONLY: cp_fm_create,&
42 : cp_fm_release,&
43 : cp_fm_set_all,&
44 : cp_fm_to_fm,&
45 : cp_fm_type
46 : USE cp_log_handling, ONLY: cp_get_default_logger,&
47 : cp_logger_get_default_io_unit,&
48 : cp_logger_type,&
49 : cp_to_string
50 : USE cp_output_handling, ONLY: cp_p_file,&
51 : cp_print_key_finished_output,&
52 : cp_print_key_should_output,&
53 : cp_print_key_unit_nr
54 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
55 : USE cube_utils, ONLY: compute_cube_center,&
56 : cube_info_type,&
57 : return_cube
58 : USE gaussian_gridlevels, ONLY: gridlevel_info_type
59 : USE grid_api, ONLY: &
60 : GRID_FUNC_AB, GRID_FUNC_ADBmDAB_X, GRID_FUNC_ADBmDAB_Y, GRID_FUNC_ADBmDAB_Z, &
61 : GRID_FUNC_ARDBmDARB_XX, GRID_FUNC_ARDBmDARB_XY, GRID_FUNC_ARDBmDARB_XZ, &
62 : GRID_FUNC_ARDBmDARB_YX, GRID_FUNC_ARDBmDARB_YY, GRID_FUNC_ARDBmDARB_YZ, &
63 : GRID_FUNC_ARDBmDARB_ZX, GRID_FUNC_ARDBmDARB_ZY, GRID_FUNC_ARDBmDARB_ZZ, &
64 : collocate_pgf_product
65 : USE input_constants, ONLY: current_gauge_atom
66 : USE input_section_types, ONLY: section_get_ivals,&
67 : section_get_lval,&
68 : section_vals_get_subs_vals,&
69 : section_vals_type
70 : USE kinds, ONLY: default_path_length,&
71 : default_string_length,&
72 : dp
73 : 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 : 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, rpgfa, rpgfb, sphi_a, sphi_b, work, zeta, zetb
610 2394 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: jpabt_a, jpabt_b, jpabt_c, jpabt_d, workt
611 2394 : TYPE(atom_pair_type), DIMENSION(:), POINTER :: atom_pair_recv, atom_pair_send
612 : TYPE(cell_type), POINTER :: cell
613 2394 : TYPE(cube_info_type), DIMENSION(:), POINTER :: cube_info
614 2394 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: deltajp_a, deltajp_b, deltajp_c, &
615 2394 : deltajp_d
616 : TYPE(dbcsr_type), POINTER :: mat_a, mat_b, mat_c, mat_d
617 : TYPE(dft_control_type), POINTER :: dft_control
618 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
619 2394 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
620 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b, orb_basis_set
621 : TYPE(mp_para_env_type), POINTER :: para_env
622 : TYPE(neighbor_list_iterator_p_type), &
623 2394 : DIMENSION(:), POINTER :: nl_iterator
624 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
625 2394 : POINTER :: sab_orb
626 2394 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
627 : TYPE(pw_env_type), POINTER :: pw_env
628 2394 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
629 : TYPE(qs_kind_type), POINTER :: qs_kind
630 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
631 2394 : POINTER :: rs_descs
632 2394 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_current, rs_rho
633 : TYPE(realspace_grid_type), DIMENSION(:, :), &
634 2394 : POINTER :: rs_gauge
635 2394 : TYPE(task_type), DIMENSION(:), POINTER :: tasks
636 :
637 2394 : NULLIFY (qs_kind, cell, dft_control, orb_basis_set, rs_rho, &
638 2394 : qs_kind_set, sab_orb, particle_set, rs_current, pw_env, &
639 2394 : rs_descs, para_env, set_radius_a, set_radius_b, la_max, la_min, &
640 2394 : lb_max, lb_min, npgfa, npgfb, nsgfa, nsgfb, rpgfa, rpgfb, &
641 2394 : sphi_a, sphi_b, zeta, zetb, first_sgfa, first_sgfb, tasks, &
642 2394 : workt, mat_a, mat_b, mat_c, mat_d, rs_gauge)
643 2394 : NULLIFY (deltajp_a, deltajp_b, deltajp_c, deltajp_d)
644 2394 : NULLIFY (jp_block_a, jp_block_b, jp_block_c, jp_block_d)
645 2394 : NULLIFY (jpabt_a, jpabt_b, jpabt_c, jpabt_d)
646 2394 : NULLIFY (atom_pair_send, atom_pair_recv)
647 :
648 2394 : CALL timeset(routineN, handle)
649 :
650 : !
651 : ! Set pointers for the different gauge
652 : ! If do_igaim is False the current_env is never needed
653 2394 : do_igaim = current_env%gauge .EQ. current_gauge_atom
654 :
655 2394 : mat_a => mat_jp
656 2394 : mat_b => mat_jp_rii
657 2394 : mat_c => mat_jp_riii
658 2394 : IF (do_igaim) mat_d => mat_d0
659 :
660 2394 : my_retain_rsgrid = .FALSE.
661 2394 : IF (PRESENT(retain_rsgrid)) my_retain_rsgrid = retain_rsgrid
662 :
663 : CALL get_qs_env(qs_env=qs_env, &
664 : qs_kind_set=qs_kind_set, &
665 : cell=cell, &
666 : dft_control=dft_control, &
667 : particle_set=particle_set, &
668 : sab_all=sab_orb, &
669 : para_env=para_env, &
670 2394 : pw_env=pw_env)
671 :
672 2394 : IF (do_igaim) CALL get_current_env(current_env=current_env, rs_gauge=rs_gauge)
673 :
674 : ! Component of appearing in the vector product rxp, iiB and iiiB
675 2394 : CALL set_vecp(iB, iiB, iiiB)
676 : !
677 : !
678 2394 : scale2 = 0.0_dp
679 2394 : idir2 = 1
680 2394 : IF (idir .NE. iB) THEN
681 1500 : CALL set_vecp_rev(idir, iB, idir2)
682 1500 : scale2 = fac_vecp(idir, iB, idir2)
683 : END IF
684 : !
685 : ! *** assign from pw_env
686 2394 : gridlevel_info => pw_env%gridlevel_info
687 2394 : cube_info => pw_env%cube_info
688 :
689 : ! Check that the neighbor list with all the pairs is associated
690 2394 : CPASSERT(ASSOCIATED(sab_orb))
691 : ! *** set up the pw multi-grids
692 2394 : CPASSERT(ASSOCIATED(pw_env))
693 2394 : CALL pw_env_get(pw_env, rs_descs=rs_descs, rs_grids=rs_rho)
694 :
695 2394 : distributed_rs_grids = .FALSE.
696 11934 : DO igrid_level = 1, gridlevel_info%ngrid_levels
697 40554 : IF (.NOT. ALL(rs_descs(igrid_level)%rs_desc%perd == 1)) THEN
698 0 : distributed_rs_grids = .TRUE.
699 : END IF
700 : END DO
701 2394 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
702 2394 : nthread = 1
703 :
704 : ! *** Allocate work storage ***
705 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
706 : maxco=maxco, &
707 : maxsgf=maxsgf, &
708 2394 : maxsgf_set=maxsgf_set)
709 :
710 9576 : Lxo2 = SQRT(SUM(cell%hmat(:, 1)**2))/2.0_dp
711 9576 : Lyo2 = SQRT(SUM(cell%hmat(:, 2)**2))/2.0_dp
712 9576 : Lzo2 = SQRT(SUM(cell%hmat(:, 3)**2))/2.0_dp
713 :
714 2394 : my_soft = .FALSE.
715 2394 : IF (PRESENT(soft_valid)) my_soft = soft_valid
716 2250 : IF (my_soft) THEN
717 1260 : basis_type = "ORB_SOFT"
718 : ELSE
719 1134 : basis_type = "ORB"
720 : END IF
721 :
722 2394 : nkind = SIZE(qs_kind_set)
723 :
724 2394 : CALL reallocate(jpabt_a, 1, maxco, 1, maxco, 0, nthread - 1)
725 2394 : CALL reallocate(jpabt_b, 1, maxco, 1, maxco, 0, nthread - 1)
726 2394 : CALL reallocate(jpabt_c, 1, maxco, 1, maxco, 0, nthread - 1)
727 2394 : CALL reallocate(jpabt_d, 1, maxco, 1, maxco, 0, nthread - 1)
728 2394 : CALL reallocate(workt, 1, maxco, 1, maxsgf_set, 0, nthread - 1)
729 2394 : CALL reallocate_tasks(tasks, max_tasks)
730 :
731 2394 : ntasks = 0
732 2394 : curr_tasks = SIZE(tasks)
733 :
734 : ! get maximum numbers
735 2394 : natom = SIZE(particle_set)
736 2394 : maxset = 0
737 2394 : maxpgf = 0
738 :
739 : ! hard code matrix index (no kpoints)
740 2394 : nimages = dft_control%nimages
741 2394 : CPASSERT(nimages == 1)
742 2394 : cindex = 1
743 :
744 6414 : DO ikind = 1, nkind
745 4020 : qs_kind => qs_kind_set(ikind)
746 :
747 4020 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
748 :
749 4020 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
750 :
751 4020 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, npgf=npgfa, nset=nseta)
752 4020 : maxset = MAX(nseta, maxset)
753 16344 : maxpgf = MAX(MAXVAL(npgfa), maxpgf)
754 : END DO
755 :
756 : ! *** Initialize working density matrix ***
757 :
758 : ! distributed rs grids require a matrix that will be changed (distribute_tasks)
759 : ! whereas this is not the case for replicated grids
760 11970 : ALLOCATE (deltajp_a(1), deltajp_b(1), deltajp_c(1), deltajp_d(1))
761 2394 : IF (distributed_rs_grids) THEN
762 0 : ALLOCATE (deltajp_a(1)%matrix, deltajp_b(1)%matrix, deltajp_c(1)%matrix)
763 0 : IF (do_igaim) THEN
764 0 : ALLOCATE (deltajp_d(1)%matrix)
765 : END IF
766 :
767 0 : CALL dbcsr_create(deltajp_a(1)%matrix, template=mat_a, name='deltajp_a')
768 0 : CALL dbcsr_create(deltajp_b(1)%matrix, template=mat_a, name='deltajp_b')
769 0 : CALL dbcsr_create(deltajp_c(1)%matrix, template=mat_a, name='deltajp_c')
770 0 : IF (do_igaim) CALL dbcsr_create(deltajp_d(1)%matrix, template=mat_a, name='deltajp_d')
771 : ELSE
772 2394 : deltajp_a(1)%matrix => mat_a !mat_jp
773 2394 : deltajp_b(1)%matrix => mat_b !mat_jp_rii
774 2394 : deltajp_c(1)%matrix => mat_c !mat_jp_riii
775 2394 : IF (do_igaim) deltajp_d(1)%matrix => mat_d !mat_d0
776 : END IF
777 :
778 11202 : ALLOCATE (basis_set_list(nkind))
779 6414 : DO ikind = 1, nkind
780 4020 : qs_kind => qs_kind_set(ikind)
781 4020 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set_a, basis_type=basis_type)
782 6414 : IF (ASSOCIATED(basis_set_a)) THEN
783 4020 : basis_set_list(ikind)%gto_basis_set => basis_set_a
784 : ELSE
785 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
786 : END IF
787 : END DO
788 2394 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
789 190833 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
790 188439 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
791 188439 : basis_set_a => basis_set_list(ikind)%gto_basis_set
792 188439 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
793 188439 : basis_set_b => basis_set_list(jkind)%gto_basis_set
794 188439 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
795 188439 : ra(:) = pbc(particle_set(iatom)%r, cell)
796 : ! basis ikind
797 188439 : first_sgfa => basis_set_a%first_sgf
798 188439 : la_max => basis_set_a%lmax
799 188439 : la_min => basis_set_a%lmin
800 188439 : npgfa => basis_set_a%npgf
801 188439 : nseta = basis_set_a%nset
802 188439 : nsgfa => basis_set_a%nsgf_set
803 188439 : rpgfa => basis_set_a%pgf_radius
804 188439 : set_radius_a => basis_set_a%set_radius
805 188439 : kind_radius_a = basis_set_a%kind_radius
806 188439 : sphi_a => basis_set_a%sphi
807 188439 : zeta => basis_set_a%zet
808 : ! basis jkind
809 188439 : first_sgfb => basis_set_b%first_sgf
810 188439 : lb_max => basis_set_b%lmax
811 188439 : lb_min => basis_set_b%lmin
812 188439 : npgfb => basis_set_b%npgf
813 188439 : nsetb = basis_set_b%nset
814 188439 : nsgfb => basis_set_b%nsgf_set
815 188439 : rpgfb => basis_set_b%pgf_radius
816 188439 : set_radius_b => basis_set_b%set_radius
817 188439 : kind_radius_b = basis_set_b%kind_radius
818 188439 : sphi_b => basis_set_b%sphi
819 188439 : zetb => basis_set_b%zet
820 :
821 188439 : IF (ABS(rab(1)) > Lxo2 .OR. ABS(rab(2)) > Lyo2 .OR. ABS(rab(3)) > Lzo2) THEN
822 : CYCLE
823 : END IF
824 :
825 11574 : brow = iatom
826 11574 : bcol = jatom
827 :
828 : CALL dbcsr_get_block_p(matrix=mat_a, row=brow, col=bcol, &
829 11574 : block=jp_block_a, found=den_found_a)
830 : CALL dbcsr_get_block_p(matrix=mat_b, row=brow, col=bcol, &
831 11574 : block=jp_block_b, found=den_found)
832 : CALL dbcsr_get_block_p(matrix=mat_c, row=brow, col=bcol, &
833 11574 : block=jp_block_c, found=den_found)
834 11574 : IF (do_igaim) CALL dbcsr_get_block_p(matrix=mat_d, row=brow, col=bcol, &
835 2961 : block=jp_block_d, found=den_found)
836 :
837 11574 : IF (.NOT. ASSOCIATED(jp_block_a)) CYCLE
838 :
839 11457 : IF (distributed_rs_grids) THEN
840 0 : CALL dbcsr_put_block(deltajp_a(1)%matrix, brow, bcol, jp_block_a)
841 0 : CALL dbcsr_put_block(deltajp_b(1)%matrix, brow, bcol, jp_block_b)
842 0 : CALL dbcsr_put_block(deltajp_c(1)%matrix, brow, bcol, jp_block_c)
843 0 : IF (do_igaim) THEN
844 0 : CALL dbcsr_put_block(deltajp_d(1)%matrix, brow, bcol, jp_block_d)
845 : END IF
846 : END IF
847 :
848 : CALL task_list_inner_loop(tasks, ntasks, curr_tasks, rs_descs, &
849 : dft_control, cube_info, gridlevel_info, cindex, &
850 : iatom, jatom, rpgfa, rpgfb, zeta, zetb, kind_radius_b, &
851 : set_radius_a, set_radius_b, ra, rab, &
852 188439 : la_max, la_min, lb_max, lb_min, npgfa, npgfb, nseta, nsetb)
853 :
854 : END DO
855 2394 : CALL neighbor_list_iterator_release(nl_iterator)
856 :
857 2394 : DEALLOCATE (basis_set_list)
858 :
859 2394 : IF (distributed_rs_grids) THEN
860 0 : CALL dbcsr_finalize(deltajp_a(1)%matrix)
861 0 : CALL dbcsr_finalize(deltajp_b(1)%matrix)
862 0 : CALL dbcsr_finalize(deltajp_c(1)%matrix)
863 0 : IF (do_igaim) CALL dbcsr_finalize(deltajp_d(1)%matrix)
864 : END IF
865 :
866 : ! sorts / redistributes the task list
867 : CALL distribute_tasks(rs_descs=rs_descs, ntasks=ntasks, natoms=natom, tasks=tasks, &
868 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
869 : symmetric=.FALSE., reorder_rs_grid_ranks=.TRUE., &
870 2394 : skip_load_balance_distributed=.FALSE.)
871 :
872 52632 : ALLOCATE (rs_current(gridlevel_info%ngrid_levels))
873 :
874 11934 : DO igrid_level = 1, gridlevel_info%ngrid_levels
875 : ! Here we need to reallocate the distributed rs_grids, which may have been reordered
876 : ! by distribute_tasks
877 9540 : IF (rs_descs(igrid_level)%rs_desc%distributed .AND. .NOT. my_retain_rsgrid) THEN
878 0 : CALL rs_grid_release(rs_rho(igrid_level))
879 0 : CALL rs_grid_create(rs_rho(igrid_level), rs_descs(igrid_level)%rs_desc)
880 : END IF
881 9540 : CALL rs_grid_zero(rs_rho(igrid_level))
882 9540 : CALL rs_grid_create(rs_current(igrid_level), rs_descs(igrid_level)%rs_desc)
883 11934 : CALL rs_grid_zero(rs_current(igrid_level))
884 : END DO
885 :
886 : !
887 : ! we need to build the gauge here
888 2394 : IF (.NOT. current_env%gauge_init .AND. do_igaim) THEN
889 28 : CALL current_set_gauge(current_env, qs_env)
890 28 : current_env%gauge_init = .TRUE.
891 : END IF
892 : !
893 : ! for any case double check the bounds !
894 360 : IF (do_igaim) THEN
895 1764 : DO igrid_level = 1, gridlevel_info%ngrid_levels
896 1404 : my_rho => rs_rho(igrid_level)%r
897 1404 : my_current => rs_current(igrid_level)%r
898 : IF (LBOUND(my_rho, 3) .NE. LBOUND(my_current, 3) .OR. &
899 : LBOUND(my_rho, 2) .NE. LBOUND(my_current, 2) .OR. &
900 : LBOUND(my_rho, 1) .NE. LBOUND(my_current, 1) .OR. &
901 : UBOUND(my_rho, 3) .NE. UBOUND(my_current, 3) .OR. &
902 18252 : UBOUND(my_rho, 2) .NE. UBOUND(my_current, 2) .OR. &
903 : UBOUND(my_rho, 1) .NE. UBOUND(my_current, 1)) THEN
904 0 : WRITE (*, *) 'LBOUND(my_rho,3),LBOUND(my_current,3)', LBOUND(my_rho, 3), LBOUND(my_current, 3)
905 0 : WRITE (*, *) 'LBOUND(my_rho,2),LBOUND(my_current,2)', LBOUND(my_rho, 2), LBOUND(my_current, 2)
906 0 : WRITE (*, *) 'LBOUND(my_rho,1),LBOUND(my_current,1)', LBOUND(my_rho, 1), LBOUND(my_current, 1)
907 0 : WRITE (*, *) 'UBOUND(my_rho,3),UBOUND(my_current,3)', UBOUND(my_rho, 3), UBOUND(my_current, 3)
908 0 : WRITE (*, *) 'UBOUND(my_rho,2),UBOUND(my_current,2)', UBOUND(my_rho, 2), UBOUND(my_current, 2)
909 0 : WRITE (*, *) 'UBOUND(my_rho,1),UBOUND(my_current,1)', UBOUND(my_rho, 1), UBOUND(my_current, 1)
910 0 : CPABORT("Bug")
911 : END IF
912 :
913 1404 : my_gauge => rs_gauge(1, igrid_level)%r
914 : IF (LBOUND(my_rho, 3) .NE. LBOUND(my_gauge, 3) .OR. &
915 : LBOUND(my_rho, 2) .NE. LBOUND(my_gauge, 2) .OR. &
916 : LBOUND(my_rho, 1) .NE. LBOUND(my_gauge, 1) .OR. &
917 : UBOUND(my_rho, 3) .NE. UBOUND(my_gauge, 3) .OR. &
918 16848 : UBOUND(my_rho, 2) .NE. UBOUND(my_gauge, 2) .OR. &
919 360 : UBOUND(my_rho, 1) .NE. UBOUND(my_gauge, 1)) THEN
920 0 : WRITE (*, *) 'LBOUND(my_rho,3),LBOUND(my_gauge,3)', LBOUND(my_rho, 3), LBOUND(my_gauge, 3)
921 0 : WRITE (*, *) 'LBOUND(my_rho,2),LBOUND(my_gauge,2)', LBOUND(my_rho, 2), LBOUND(my_gauge, 2)
922 0 : WRITE (*, *) 'LBOUND(my_rho,1),LBOUND(my_gauge,1)', LBOUND(my_rho, 1), LBOUND(my_gauge, 1)
923 0 : WRITE (*, *) 'UBOUND(my_rho,3),UbOUND(my_gauge,3)', UBOUND(my_rho, 3), UBOUND(my_gauge, 3)
924 0 : WRITE (*, *) 'UBOUND(my_rho,2),UBOUND(my_gauge,2)', UBOUND(my_rho, 2), UBOUND(my_gauge, 2)
925 0 : WRITE (*, *) 'UBOUND(my_rho,1),UBOUND(my_gauge,1)', UBOUND(my_rho, 1), UBOUND(my_gauge, 1)
926 0 : CPABORT("Bug")
927 : END IF
928 : END DO
929 : END IF
930 : !
931 : !-------------------------------------------------------------
932 :
933 2394 : IF (distributed_rs_grids) THEN
934 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_a, &
935 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
936 0 : nimages=nimages, scatter=.TRUE.)
937 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_b, &
938 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
939 0 : nimages=nimages, scatter=.TRUE.)
940 : CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_c, &
941 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
942 0 : nimages=nimages, scatter=.TRUE.)
943 0 : IF (do_igaim) CALL rs_distribute_matrix(rs_descs=rs_descs, pmats=deltajp_d, &
944 : atom_pair_send=atom_pair_send, atom_pair_recv=atom_pair_recv, &
945 0 : nimages=nimages, scatter=.TRUE.)
946 : END IF
947 :
948 2394 : ithread = 0
949 2394 : jpab_a => jpabt_a(:, :, ithread)
950 2394 : jpab_b => jpabt_b(:, :, ithread)
951 2394 : jpab_c => jpabt_c(:, :, ithread)
952 2394 : IF (do_igaim) jpab_d => jpabt_d(:, :, ithread)
953 2394 : work => workt(:, :, ithread)
954 :
955 2394 : iatom_old = -1; jatom_old = -1; iset_old = -1; jset_old = -1
956 2394 : ikind_old = -1; jkind_old = -1
957 :
958 166284 : loop_tasks: DO itask = 1, ntasks
959 163890 : igrid_level = tasks(itask)%grid_level
960 163890 : cindex = tasks(itask)%image
961 163890 : iatom = tasks(itask)%iatom
962 163890 : jatom = tasks(itask)%jatom
963 163890 : iset = tasks(itask)%iset
964 163890 : jset = tasks(itask)%jset
965 163890 : ipgf = tasks(itask)%ipgf
966 163890 : jpgf = tasks(itask)%jpgf
967 :
968 : ! apparently generalised collocation not implemented correctly yet
969 163890 : CPASSERT(tasks(itask)%dist_type .NE. 2)
970 :
971 163890 : IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old) THEN
972 :
973 24138 : ikind = particle_set(iatom)%atomic_kind%kind_number
974 24138 : jkind = particle_set(jatom)%atomic_kind%kind_number
975 :
976 24138 : IF (iatom .NE. iatom_old) ra(:) = pbc(particle_set(iatom)%r, cell)
977 :
978 24138 : brow = iatom
979 24138 : bcol = jatom
980 :
981 24138 : IF (ikind .NE. ikind_old) THEN
982 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
983 3450 : basis_type=basis_type)
984 :
985 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
986 : first_sgf=first_sgfa, &
987 : lmax=la_max, &
988 : lmin=la_min, &
989 : npgf=npgfa, &
990 : nset=nseta, &
991 : nsgf_set=nsgfa, &
992 : pgf_radius=rpgfa, &
993 : set_radius=set_radius_a, &
994 : sphi=sphi_a, &
995 3450 : zet=zeta)
996 : END IF
997 :
998 24138 : IF (jkind .NE. jkind_old) THEN
999 :
1000 : CALL get_qs_kind(qs_kind_set(jkind), &
1001 12900 : basis_set=orb_basis_set, basis_type=basis_type)
1002 :
1003 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1004 : first_sgf=first_sgfb, &
1005 : kind_radius=kind_radius_b, &
1006 : lmax=lb_max, &
1007 : lmin=lb_min, &
1008 : npgf=npgfb, &
1009 : nset=nsetb, &
1010 : nsgf_set=nsgfb, &
1011 : pgf_radius=rpgfb, &
1012 : set_radius=set_radius_b, &
1013 : sphi=sphi_b, &
1014 12900 : zet=zetb)
1015 :
1016 : END IF
1017 :
1018 : CALL dbcsr_get_block_p(matrix=deltajp_a(1)%matrix, row=brow, col=bcol, &
1019 24138 : block=jp_block_a, found=den_found)
1020 : CALL dbcsr_get_block_p(matrix=deltajp_b(1)%matrix, row=brow, col=bcol, &
1021 24138 : block=jp_block_b, found=den_found)
1022 : CALL dbcsr_get_block_p(matrix=deltajp_c(1)%matrix, row=brow, col=bcol, &
1023 24138 : block=jp_block_c, found=den_found)
1024 24138 : IF (do_igaim) CALL dbcsr_get_block_p(matrix=deltajp_d(1)%matrix, row=brow, col=bcol, &
1025 5229 : block=jp_block_d, found=den_found)
1026 :
1027 24138 : IF (.NOT. ASSOCIATED(jp_block_a)) &
1028 0 : CPABORT("p_block not associated in deltap")
1029 :
1030 : iatom_old = iatom
1031 : jatom_old = jatom
1032 : ikind_old = ikind
1033 : jkind_old = jkind
1034 :
1035 : atom_pair_changed = .TRUE.
1036 :
1037 : ELSE
1038 :
1039 : atom_pair_changed = .FALSE.
1040 :
1041 : END IF
1042 :
1043 163890 : IF (atom_pair_changed .OR. iset_old .NE. iset .OR. jset_old .NE. jset) THEN
1044 :
1045 58395 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1046 58395 : sgfa = first_sgfa(1, iset)
1047 58395 : ncob = npgfb(jset)*ncoset(lb_max(jset))
1048 58395 : sgfb = first_sgfb(1, jset)
1049 : ! Decontraction step for the selected blocks of the 3 density matrices
1050 :
1051 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1052 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1053 : jp_block_a(sgfa, sgfb), SIZE(jp_block_a, 1), &
1054 58395 : 0.0_dp, work(1, 1), maxco)
1055 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1056 : 1.0_dp, work(1, 1), maxco, &
1057 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1058 58395 : 0.0_dp, jpab_a(1, 1), maxco)
1059 :
1060 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1061 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1062 : jp_block_b(sgfa, sgfb), SIZE(jp_block_b, 1), &
1063 58395 : 0.0_dp, work(1, 1), maxco)
1064 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1065 : 1.0_dp, work(1, 1), maxco, &
1066 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1067 58395 : 0.0_dp, jpab_b(1, 1), maxco)
1068 :
1069 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1070 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1071 : jp_block_c(sgfa, sgfb), SIZE(jp_block_c, 1), &
1072 58395 : 0.0_dp, work(1, 1), maxco)
1073 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1074 : 1.0_dp, work(1, 1), maxco, &
1075 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1076 58395 : 0.0_dp, jpab_c(1, 1), maxco)
1077 :
1078 58395 : IF (do_igaim) THEN
1079 : CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
1080 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1081 : jp_block_d(sgfa, sgfb), SIZE(jp_block_d, 1), &
1082 12105 : 0.0_dp, work(1, 1), maxco)
1083 : CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
1084 : 1.0_dp, work(1, 1), maxco, &
1085 : sphi_b(1, sgfb), SIZE(sphi_b, 1), &
1086 12105 : 0.0_dp, jpab_d(1, 1), maxco)
1087 : END IF
1088 :
1089 : iset_old = iset
1090 : jset_old = jset
1091 :
1092 : END IF
1093 :
1094 54630 : SELECT CASE (idir)
1095 : CASE (1)
1096 54630 : adbmdab_func = GRID_FUNC_ADBmDAB_X
1097 : CASE (2)
1098 54630 : adbmdab_func = GRID_FUNC_ADBmDAB_Y
1099 : CASE (3)
1100 54630 : adbmdab_func = GRID_FUNC_ADBmDAB_Z
1101 : CASE DEFAULT
1102 163890 : CPABORT("invalid idir")
1103 : END SELECT
1104 :
1105 655560 : rab(:) = tasks(itask)%rab
1106 655560 : rb(:) = ra(:) + rab(:)
1107 163890 : zetp = zeta(ipgf, iset) + zetb(jpgf, jset)
1108 163890 : f = zetb(jpgf, jset)/zetp
1109 655560 : prefactor = EXP(-zeta(ipgf, iset)*f*DOT_PRODUCT(rab, rab))
1110 655560 : rp(:) = ra(:) + f*rab(:)
1111 :
1112 163890 : na1 = (ipgf - 1)*ncoset(la_max(iset)) + 1
1113 163890 : na2 = ipgf*ncoset(la_max(iset))
1114 163890 : nb1 = (jpgf - 1)*ncoset(lb_max(jset)) + 1
1115 163890 : nb2 = jpgf*ncoset(lb_max(jset))
1116 :
1117 : ! Four calls to the general collocate density, to multply the correct function
1118 : ! to each density matrix
1119 :
1120 : !
1121 : ! here the decontracted mat_jp_{ab} is multiplied by
1122 : ! f_{ab} = g_{a} (dg_{b}/dr)_{idir} - (dg_{a}/dr)_{idir} g_{b}
1123 163890 : scale = 1.0_dp
1124 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1125 : lb_min=lb_min(jset), lb_max=lb_max(jset), &
1126 : ra=ra, rb=rb, rp=rp, zetp=zetp, eps=eps_rho_rspace, &
1127 163890 : prefactor=prefactor, cutoff=1.0_dp)
1128 :
1129 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1130 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1131 : ra, rab, scale, jpab_a, na1 - 1, nb1 - 1, &
1132 : rs_current(igrid_level), &
1133 163890 : radius=radius, ga_gb_function=adbmdab_func)
1134 166284 : IF (do_igaim) THEN
1135 : ! here the decontracted mat_jb_{ab} is multiplied by
1136 : ! f_{ab} = g_{a} * g_{b} ! THIS GOES OUTSIDE THE LOOP !
1137 20574 : IF (scale2 .NE. 0.0_dp) THEN
1138 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1139 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1140 : ra, rab, scale2, jpab_d, na1 - 1, nb1 - 1, &
1141 : rs_rho(igrid_level), &
1142 13716 : radius=radius, ga_gb_function=GRID_FUNC_AB)
1143 : END IF !rm
1144 : ! here the decontracted mat_jp_rii{ab} is multiplied by
1145 : ! f_{ab} = g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} -
1146 : ! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b}
1147 : scale = 1.0_dp
1148 512067753 : current_env%rs_buf(igrid_level)%r(:, :, :) = 0.0_dp
1149 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1150 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1151 : ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
1152 : radius=radius, &
1153 : ga_gb_function=adbmdab_func, &
1154 20574 : rsgrid=current_env%rs_buf(igrid_level))
1155 : CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level), &
1156 : rsbuf=rs_current(igrid_level), &
1157 : rsgauge=rs_gauge(iiiB, igrid_level), &
1158 : cube_info=cube_info(igrid_level), radius=radius, &
1159 : zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
1160 20574 : ra=ra, rab=rab, ir=iiiB)
1161 :
1162 : ! here the decontracted mat_jp_riii{ab} is multiplied by
1163 : ! f_{ab} = -g_{a} (d(r) - R_{b})_{iiB} (dg_{b}/dr)_{idir} +
1164 : ! (dg_{a}/dr)_{idir} (d(r) - R_{b})_{iiB} g_{b}
1165 20574 : scale = -1.0_dp
1166 512067753 : current_env%rs_buf(igrid_level)%r(:, :, :) = 0.0_dp
1167 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1168 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1169 : ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
1170 : radius=radius, &
1171 : ga_gb_function=adbmdab_func, &
1172 20574 : rsgrid=current_env%rs_buf(igrid_level))
1173 : CALL collocate_gauge_ortho(rsgrid=current_env%rs_buf(igrid_level), &
1174 : rsbuf=rs_current(igrid_level), &
1175 : rsgauge=rs_gauge(iiB, igrid_level), &
1176 : cube_info=cube_info(igrid_level), radius=radius, &
1177 : zeta=zeta(ipgf, iset), zetb=zetb(jpgf, jset), &
1178 20574 : ra=ra, rab=rab, ir=iiB)
1179 : ELSE
1180 : ! here the decontracted mat_jp_rii{ab} is multiplied by
1181 : ! f_{ab} = g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} -
1182 : ! (dg_{a}/dr)_{idir} (r - R_{b})_{iiB} g_{b}
1183 : scale = 1.0_dp
1184 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1185 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1186 : ra, rab, scale, jpab_b, na1 - 1, nb1 - 1, &
1187 : rs_current(igrid_level), &
1188 : radius=radius, &
1189 143316 : ga_gb_function=encode_ardbmdarb_func(idir=idir, ir=iiiB))
1190 : ! here the decontracted mat_jp_riii{ab} is multiplied by
1191 : ! f_{ab} = -g_{a} (r - R_{b})_{iiB} (dg_{b}/dr)_{idir} +
1192 : ! (dg_{a}/dr)_{idir} (r - R_{b})_{iiB} g_{b}
1193 143316 : scale = -1.0_dp
1194 : CALL collocate_pgf_product(la_max(iset), zeta(ipgf, iset), &
1195 : la_min(iset), lb_max(jset), zetb(jpgf, jset), lb_min(jset), &
1196 : ra, rab, scale, jpab_c, na1 - 1, nb1 - 1, &
1197 : rs_current(igrid_level), &
1198 : radius=radius, &
1199 143316 : ga_gb_function=encode_ardbmdarb_func(idir=idir, ir=iiB))
1200 : END IF
1201 :
1202 : END DO loop_tasks
1203 : !
1204 : ! Scale the density with the gauge rho * ( r - d(r) ) if needed
1205 2394 : IF (do_igaim) THEN
1206 1764 : DO igrid_level = 1, gridlevel_info%ngrid_levels
1207 : CALL rs_grid_mult_and_add(rs_current(igrid_level), rs_rho(igrid_level), &
1208 1764 : rs_gauge(idir2, igrid_level), 1.0_dp)
1209 : END DO
1210 : END IF
1211 : ! *** Release work storage ***
1212 :
1213 2394 : IF (distributed_rs_grids) THEN
1214 0 : CALL dbcsr_deallocate_matrix(deltajp_a(1)%matrix)
1215 0 : CALL dbcsr_deallocate_matrix(deltajp_b(1)%matrix)
1216 0 : CALL dbcsr_deallocate_matrix(deltajp_c(1)%matrix)
1217 0 : IF (do_igaim) CALL dbcsr_deallocate_matrix(deltajp_d(1)%matrix)
1218 : END IF
1219 2394 : DEALLOCATE (deltajp_a, deltajp_b, deltajp_c, deltajp_d)
1220 :
1221 2394 : DEALLOCATE (jpabt_a, jpabt_b, jpabt_c, jpabt_d, workt, tasks)
1222 :
1223 2394 : IF (ASSOCIATED(atom_pair_send)) DEALLOCATE (atom_pair_send)
1224 2394 : IF (ASSOCIATED(atom_pair_recv)) DEALLOCATE (atom_pair_recv)
1225 :
1226 2394 : CALL density_rs2pw(pw_env, rs_current, current_rs, current_gs)
1227 11934 : DO i = 1, SIZE(rs_current)
1228 11934 : CALL rs_grid_release(rs_current(i))
1229 : END DO
1230 :
1231 11934 : DO i = 1, SIZE(rs_rho)
1232 11934 : IF (rs_descs(i)%rs_desc%distributed .AND. .NOT. my_retain_rsgrid) THEN
1233 0 : CALL rs_grid_release(rs_rho(i))
1234 : END IF
1235 : END DO
1236 :
1237 : ! Free the array of grids (grids themselves are released in density_rs2pw)
1238 2394 : DEALLOCATE (rs_current)
1239 :
1240 2394 : CALL timestop(handle)
1241 :
1242 7182 : END SUBROUTINE calculate_jrho_resp
1243 :
1244 : ! **************************************************************************************************
1245 : !> \brief ...
1246 : !> \param idir ...
1247 : !> \param ir ...
1248 : !> \return ...
1249 : ! **************************************************************************************************
1250 286632 : FUNCTION encode_ardbmdarb_func(idir, ir) RESULT(func)
1251 : INTEGER, INTENT(IN) :: idir, ir
1252 : INTEGER :: func
1253 :
1254 286632 : CPASSERT(1 <= idir .AND. idir <= 3 .AND. 1 <= ir .AND. ir <= 3)
1255 286632 : SELECT CASE (10*idir + ir)
1256 : CASE (11)
1257 32610 : func = GRID_FUNC_ARDBmDARB_XX
1258 : CASE (12)
1259 32610 : func = GRID_FUNC_ARDBmDARB_XY
1260 : CASE (13)
1261 32610 : func = GRID_FUNC_ARDBmDARB_XZ
1262 : CASE (21)
1263 32610 : func = GRID_FUNC_ARDBmDARB_YX
1264 : CASE (22)
1265 30324 : func = GRID_FUNC_ARDBmDARB_YY
1266 : CASE (23)
1267 32610 : func = GRID_FUNC_ARDBmDARB_YZ
1268 : CASE (31)
1269 32610 : func = GRID_FUNC_ARDBmDARB_ZX
1270 : CASE (32)
1271 32610 : func = GRID_FUNC_ARDBmDARB_ZY
1272 : CASE (33)
1273 30324 : func = GRID_FUNC_ARDBmDARB_ZZ
1274 : CASE DEFAULT
1275 286632 : CPABORT("invalid idir or iiiB")
1276 : END SELECT
1277 286632 : END FUNCTION encode_ardbmdarb_func
1278 :
1279 : ! **************************************************************************************************
1280 : !> \brief ...
1281 : !> \param rsgrid ...
1282 : !> \param rsbuf ...
1283 : !> \param rsgauge ...
1284 : !> \param cube_info ...
1285 : !> \param radius ...
1286 : !> \param ra ...
1287 : !> \param rab ...
1288 : !> \param zeta ...
1289 : !> \param zetb ...
1290 : !> \param ir ...
1291 : ! **************************************************************************************************
1292 41148 : SUBROUTINE collocate_gauge_ortho(rsgrid, rsbuf, rsgauge, cube_info, radius, ra, rab, zeta, zetb, ir)
1293 : TYPE(realspace_grid_type) :: rsgrid, rsbuf, rsgauge
1294 : TYPE(cube_info_type), INTENT(IN) :: cube_info
1295 : REAL(KIND=dp), INTENT(IN) :: radius
1296 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rab
1297 : REAL(KIND=dp), INTENT(IN) :: zeta, zetb
1298 : INTEGER, INTENT(IN) :: ir
1299 :
1300 : INTEGER :: cmax, i, ig, igmax, igmin, j, j2, jg, &
1301 : jg2, jgmin, k, k2, kg, kg2, kgmin, &
1302 : length, offset, sci, start
1303 41148 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map
1304 : INTEGER, DIMENSION(3) :: cubecenter, lb_cube, ng, ub_cube
1305 41148 : INTEGER, DIMENSION(:), POINTER :: sphere_bounds
1306 : REAL(KIND=dp) :: f, point(3, 4), res(4), x, y, y2, z, z2, &
1307 : zetp
1308 : REAL(KIND=dp), DIMENSION(3) :: dr, rap, rb, rbp, roffset, rp
1309 41148 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gauge, grid, grid_buf
1310 :
1311 0 : CPASSERT(rsgrid%desc%orthorhombic)
1312 41148 : NULLIFY (sphere_bounds)
1313 :
1314 41148 : grid => rsgrid%r(:, :, :)
1315 41148 : grid_buf => rsbuf%r(:, :, :)
1316 41148 : gauge => rsgauge%r(:, :, :)
1317 :
1318 : ! *** center of gaussians and their product
1319 41148 : zetp = zeta + zetb
1320 41148 : f = zetb/zetp
1321 164592 : rap(:) = f*rab(:)
1322 164592 : rbp(:) = rap(:) - rab(:)
1323 164592 : rp(:) = ra(:) + rap(:)
1324 164592 : rb(:) = ra(:) + rab(:)
1325 :
1326 : ! *** properties of the grid ***
1327 164592 : ng(:) = rsgrid%desc%npts(:)
1328 41148 : dr(1) = rsgrid%desc%dh(1, 1)
1329 41148 : dr(2) = rsgrid%desc%dh(2, 2)
1330 41148 : dr(3) = rsgrid%desc%dh(3, 3)
1331 :
1332 : ! *** get the sub grid properties for the given radius ***
1333 41148 : CALL return_cube(cube_info, radius, lb_cube, ub_cube, sphere_bounds)
1334 164592 : cmax = MAXVAL(ub_cube)
1335 :
1336 : ! *** position of the gaussian product
1337 : !
1338 : ! this is the actual definition of the position on the grid
1339 : ! i.e. a point rp(:) gets here grid coordinates
1340 : ! MODULO(rp(:)/dr(:),ng(:))+1
1341 : ! hence (0.0,0.0,0.0) in real space is rsgrid%lb on the rsgrid ((1,1,1) on grid)
1342 : !
1343 164592 : ALLOCATE (map(-cmax:cmax, 3))
1344 3221964 : map(:, :) = -1
1345 41148 : CALL compute_cube_center(cubecenter, rsgrid%desc, zeta, zetb, ra, rab)
1346 164592 : roffset(:) = rp(:) - REAL(cubecenter(:), dp)*dr(:)
1347 :
1348 : ! *** a mapping so that the ig corresponds to the right grid point
1349 164592 : DO i = 1, 3
1350 164592 : IF (rsgrid%desc%perd(i) == 1) THEN
1351 123444 : start = lb_cube(i)
1352 123354 : DO
1353 246798 : offset = MODULO(cubecenter(i) + start, ng(i)) + 1 - start
1354 246798 : length = MIN(ub_cube(i), ng(i) - offset) - start
1355 3180726 : DO ig = start, start + length
1356 3180726 : map(ig, i) = ig + offset
1357 : END DO
1358 246798 : IF (start + length .GE. ub_cube(i)) EXIT
1359 123354 : start = start + length + 1
1360 : END DO
1361 : ELSE
1362 : ! this takes partial grid + border regions into account
1363 0 : offset = MODULO(cubecenter(i) + lb_cube(i) + rsgrid%desc%lb(i) - rsgrid%lb_local(i), ng(i)) + 1 - lb_cube(i)
1364 : ! check for out of bounds
1365 0 : IF (ub_cube(i) + offset > UBOUND(grid, i) .OR. lb_cube(i) + offset < LBOUND(grid, i)) THEN
1366 0 : CPASSERT(.FALSE.)
1367 : END IF
1368 0 : DO ig = lb_cube(i), ub_cube(i)
1369 0 : map(ig, i) = ig + offset
1370 : END DO
1371 : END IF
1372 : END DO
1373 :
1374 : ! *** actually loop over the grid
1375 41148 : sci = 1
1376 41148 : kgmin = sphere_bounds(sci)
1377 41148 : sci = sci + 1
1378 530136 : DO kg = kgmin, 0
1379 488988 : kg2 = 1 - kg
1380 488988 : k = map(kg, 3)
1381 488988 : k2 = map(kg2, 3)
1382 488988 : jgmin = sphere_bounds(sci)
1383 488988 : sci = sci + 1
1384 488988 : z = (REAL(kg, dp) + REAL(cubecenter(3), dp))*dr(3)
1385 488988 : z2 = (REAL(kg2, dp) + REAL(cubecenter(3), dp))*dr(3)
1386 5604174 : DO jg = jgmin, 0
1387 5074038 : jg2 = 1 - jg
1388 5074038 : j = map(jg, 2)
1389 5074038 : j2 = map(jg2, 2)
1390 5074038 : igmin = sphere_bounds(sci)
1391 5074038 : sci = sci + 1
1392 5074038 : igmax = 1 - igmin
1393 5074038 : y = (REAL(jg, dp) + REAL(cubecenter(2), dp))*dr(2)
1394 5074038 : y2 = (REAL(jg2, dp) + REAL(cubecenter(2), dp))*dr(2)
1395 118256958 : DO ig = igmin, igmax
1396 112693932 : i = map(ig, 1)
1397 112693932 : x = (REAL(ig, dp) + REAL(cubecenter(1), dp))*dr(1)
1398 112693932 : point(1, 1) = x; point(2, 1) = y; point(3, 1) = z
1399 112693932 : point(1, 2) = x; point(2, 2) = y2; point(3, 2) = z
1400 112693932 : point(1, 3) = x; point(2, 3) = y; point(3, 3) = z2
1401 112693932 : point(1, 4) = x; point(2, 4) = y2; point(3, 4) = z2
1402 : !
1403 112693932 : res(1) = (point(ir, 1) - rb(ir)) - gauge(i, j, k)
1404 112693932 : res(2) = (point(ir, 2) - rb(ir)) - gauge(i, j2, k)
1405 112693932 : res(3) = (point(ir, 3) - rb(ir)) - gauge(i, j, k2)
1406 112693932 : res(4) = (point(ir, 4) - rb(ir)) - gauge(i, j2, k2)
1407 : !
1408 112693932 : grid_buf(i, j, k) = grid_buf(i, j, k) + grid(i, j, k)*res(1)
1409 112693932 : grid_buf(i, j2, k) = grid_buf(i, j2, k) + grid(i, j2, k)*res(2)
1410 112693932 : grid_buf(i, j, k2) = grid_buf(i, j, k2) + grid(i, j, k2)*res(3)
1411 117767970 : grid_buf(i, j2, k2) = grid_buf(i, j2, k2) + grid(i, j2, k2)*res(4)
1412 : END DO
1413 : END DO
1414 : END DO
1415 82296 : END SUBROUTINE collocate_gauge_ortho
1416 :
1417 : ! **************************************************************************************************
1418 : !> \brief ...
1419 : !> \param current_env ...
1420 : !> \param qs_env ...
1421 : ! **************************************************************************************************
1422 28 : SUBROUTINE current_set_gauge(current_env, qs_env)
1423 : !
1424 : TYPE(current_env_type) :: current_env
1425 : TYPE(qs_environment_type), POINTER :: qs_env
1426 :
1427 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_set_gauge'
1428 :
1429 : INTEGER :: idir
1430 : REAL(dp) :: dbox(3)
1431 28 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: box_data
1432 : INTEGER :: handle, igrid_level, nbox(3), gauge
1433 28 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: box_ptr
1434 : TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
1435 28 : POINTER :: rs_descs
1436 : TYPE(pw_env_type), POINTER :: pw_env
1437 28 : TYPE(realspace_grid_type), DIMENSION(:, :), POINTER :: rs_gauge
1438 :
1439 28 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
1440 : LOGICAL :: use_old_gauge_atom
1441 :
1442 28 : NULLIFY (rs_gauge, box)
1443 :
1444 28 : CALL timeset(routineN, handle)
1445 :
1446 : CALL get_current_env(current_env=current_env, &
1447 : use_old_gauge_atom=use_old_gauge_atom, &
1448 : rs_gauge=rs_gauge, &
1449 28 : gauge=gauge)
1450 :
1451 28 : IF (gauge .EQ. current_gauge_atom) THEN
1452 : CALL get_qs_env(qs_env=qs_env, &
1453 28 : pw_env=pw_env)
1454 : CALL pw_env_get(pw_env=pw_env, &
1455 28 : rs_descs=rs_descs)
1456 : !
1457 : ! box the atoms
1458 28 : IF (use_old_gauge_atom) THEN
1459 22 : CALL box_atoms(qs_env)
1460 : ELSE
1461 6 : CALL box_atoms_new(current_env, qs_env, box)
1462 : END IF
1463 : !
1464 : ! allocate and build the gauge
1465 136 : DO igrid_level = pw_env%gridlevel_info%ngrid_levels, 1, -1
1466 :
1467 432 : DO idir = 1, 3
1468 432 : CALL rs_grid_create(rs_gauge(idir, igrid_level), rs_descs(igrid_level)%rs_desc)
1469 : END DO
1470 :
1471 136 : IF (use_old_gauge_atom) THEN
1472 : CALL collocate_gauge(current_env, qs_env, &
1473 : rs_gauge(1, igrid_level), &
1474 : rs_gauge(2, igrid_level), &
1475 88 : rs_gauge(3, igrid_level))
1476 : ELSE
1477 : CALL collocate_gauge_new(current_env, qs_env, &
1478 : rs_gauge(1, igrid_level), &
1479 : rs_gauge(2, igrid_level), &
1480 : rs_gauge(3, igrid_level), &
1481 20 : box)
1482 : END IF
1483 : END DO
1484 : !
1485 : ! allocate the buf
1486 724 : ALLOCATE (current_env%rs_buf(pw_env%gridlevel_info%ngrid_levels))
1487 136 : DO igrid_level = 1, pw_env%gridlevel_info%ngrid_levels
1488 136 : CALL rs_grid_create(current_env%rs_buf(igrid_level), rs_descs(igrid_level)%rs_desc)
1489 : END DO
1490 : !
1491 28 : DEALLOCATE (box_ptr, box_data)
1492 28 : CALL deallocate_box(box)
1493 : END IF
1494 :
1495 56 : CALL timestop(handle)
1496 :
1497 : CONTAINS
1498 :
1499 : ! **************************************************************************************************
1500 : !> \brief ...
1501 : !> \param qs_env ...
1502 : ! **************************************************************************************************
1503 22 : SUBROUTINE box_atoms(qs_env)
1504 : TYPE(qs_environment_type), POINTER :: qs_env
1505 :
1506 : REAL(kind=dp), PARAMETER :: box_size_guess = 5.0_dp
1507 :
1508 : INTEGER :: i, iatom, ibox, ii, jbox, kbox, natms
1509 : REAL(dp) :: offset(3)
1510 22 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom
1511 : TYPE(cell_type), POINTER :: cell
1512 22 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1513 22 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1514 :
1515 : CALL get_qs_env(qs_env=qs_env, &
1516 : qs_kind_set=qs_kind_set, &
1517 : cell=cell, &
1518 22 : particle_set=particle_set)
1519 :
1520 22 : natms = SIZE(particle_set, 1)
1521 66 : ALLOCATE (ratom(3, natms))
1522 : !
1523 : ! box the atoms
1524 22 : nbox(1) = CEILING(cell%hmat(1, 1)/box_size_guess)
1525 22 : nbox(2) = CEILING(cell%hmat(2, 2)/box_size_guess)
1526 22 : nbox(3) = CEILING(cell%hmat(3, 3)/box_size_guess)
1527 : !write(*,*) 'nbox',nbox
1528 22 : dbox(1) = cell%hmat(1, 1)/REAL(nbox(1), dp)
1529 22 : dbox(2) = cell%hmat(2, 2)/REAL(nbox(2), dp)
1530 22 : dbox(3) = cell%hmat(3, 3)/REAL(nbox(3), dp)
1531 : !write(*,*) 'dbox',dbox
1532 132 : ALLOCATE (box_ptr(0:nbox(1), 0:nbox(2) - 1, 0:nbox(3) - 1), box_data(3, natms))
1533 390 : box_data(:, :) = HUGE(0.0_dp)
1534 418 : box_ptr(:, :, :) = HUGE(0)
1535 : !
1536 22 : offset(1) = cell%hmat(1, 1)*0.5_dp
1537 22 : offset(2) = cell%hmat(2, 2)*0.5_dp
1538 22 : offset(3) = cell%hmat(3, 3)*0.5_dp
1539 114 : DO iatom = 1, natms
1540 390 : ratom(:, iatom) = pbc(particle_set(iatom)%r(:), cell) + offset(:)
1541 : END DO
1542 : !
1543 22 : i = 1
1544 66 : DO kbox = 0, nbox(3) - 1
1545 154 : DO jbox = 0, nbox(2) - 1
1546 88 : box_ptr(0, jbox, kbox) = i
1547 308 : DO ibox = 0, nbox(1) - 1
1548 : ii = 0
1549 912 : DO iatom = 1, natms
1550 : IF (INT(ratom(1, iatom)/dbox(1)) .EQ. ibox .AND. &
1551 736 : INT(ratom(2, iatom)/dbox(2)) .EQ. jbox .AND. &
1552 176 : INT(ratom(3, iatom)/dbox(3)) .EQ. kbox) THEN
1553 368 : box_data(:, i) = ratom(:, iatom) - offset(:)
1554 92 : i = i + 1
1555 92 : ii = ii + 1
1556 : END IF
1557 : END DO
1558 264 : box_ptr(ibox + 1, jbox, kbox) = box_ptr(ibox, jbox, kbox) + ii
1559 : END DO
1560 : END DO
1561 : END DO
1562 : !
1563 : IF (.FALSE.) THEN
1564 : DO kbox = 0, nbox(3) - 1
1565 : DO jbox = 0, nbox(2) - 1
1566 : DO ibox = 0, nbox(1) - 1
1567 : WRITE (*, *) 'box=', ibox, jbox, kbox
1568 : WRITE (*, *) 'nbr atom=', box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox)
1569 : DO iatom = box_ptr(ibox, jbox, kbox), box_ptr(ibox + 1, jbox, kbox) - 1
1570 : WRITE (*, *) 'iatom=', iatom
1571 : WRITE (*, '(A,3E14.6)') 'coor=', box_data(:, iatom)
1572 : END DO
1573 : END DO
1574 : END DO
1575 : END DO
1576 : END IF
1577 22 : DEALLOCATE (ratom)
1578 22 : END SUBROUTINE box_atoms
1579 :
1580 : ! **************************************************************************************************
1581 : !> \brief ...
1582 : !> \param current_env ...
1583 : !> \param qs_env ...
1584 : !> \param rs_grid_x ...
1585 : !> \param rs_grid_y ...
1586 : !> \param rs_grid_z ...
1587 : ! **************************************************************************************************
1588 88 : SUBROUTINE collocate_gauge(current_env, qs_env, rs_grid_x, rs_grid_y, rs_grid_z)
1589 : !
1590 : TYPE(current_env_type) :: current_env
1591 : TYPE(qs_environment_type), POINTER :: qs_env
1592 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid_x, rs_grid_y, rs_grid_z
1593 :
1594 : INTEGER :: i, iatom, ibeg, ibox, iend, imax, imin, &
1595 : j, jatom, jbox, jmax, jmin, k, kbox, &
1596 : kmax, kmin, lb(3), lb_local(3), natms, &
1597 : natms_local, ng(3)
1598 : REAL(KIND=dp) :: ab, buf_tmp, dist, dr(3), &
1599 : gauge_atom_radius, offset(3), pa, pb, &
1600 : point(3), pra(3), r(3), res(3), summe, &
1601 : tmp, x, y, z
1602 88 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, nrm_atms_pnt
1603 88 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atms_pnt, ratom
1604 88 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid_x, grid_y, grid_z
1605 : TYPE(cell_type), POINTER :: cell
1606 88 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1607 88 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1608 :
1609 : !
1610 :
1611 : CALL get_current_env(current_env=current_env, &
1612 88 : gauge_atom_radius=gauge_atom_radius)
1613 : !
1614 : CALL get_qs_env(qs_env=qs_env, &
1615 : qs_kind_set=qs_kind_set, &
1616 : cell=cell, &
1617 88 : particle_set=particle_set)
1618 : !
1619 88 : natms = SIZE(particle_set, 1)
1620 88 : dr(1) = rs_grid_x%desc%dh(1, 1)
1621 88 : dr(2) = rs_grid_x%desc%dh(2, 2)
1622 88 : dr(3) = rs_grid_x%desc%dh(3, 3)
1623 352 : lb(:) = rs_grid_x%desc%lb(:)
1624 352 : lb_local(:) = rs_grid_x%lb_local(:)
1625 88 : grid_x => rs_grid_x%r(:, :, :)
1626 88 : grid_y => rs_grid_y%r(:, :, :)
1627 88 : grid_z => rs_grid_z%r(:, :, :)
1628 352 : ng(:) = UBOUND(grid_x)
1629 88 : offset(1) = cell%hmat(1, 1)*0.5_dp
1630 88 : offset(2) = cell%hmat(2, 2)*0.5_dp
1631 88 : offset(3) = cell%hmat(3, 3)*0.5_dp
1632 616 : ALLOCATE (buf(natms), ratom(3, natms), atms_pnt(3, natms), nrm_atms_pnt(natms))
1633 : !
1634 : ! go over the grid
1635 1350 : DO k = 1, ng(3)
1636 26188 : DO j = 1, ng(2)
1637 618230 : DO i = 1, ng(1)
1638 : !
1639 592130 : point(3) = REAL(k - 1 + lb_local(3) - lb(3), dp)*dr(3)
1640 592130 : point(2) = REAL(j - 1 + lb_local(2) - lb(2), dp)*dr(2)
1641 592130 : point(1) = REAL(i - 1 + lb_local(1) - lb(1), dp)*dr(1)
1642 2368520 : point = pbc(point, cell)
1643 : !
1644 : ! run over the overlaping boxes
1645 592130 : natms_local = 0
1646 592130 : kmin = INT((point(3) + offset(3) - gauge_atom_radius)/dbox(3))
1647 592130 : kmax = INT((point(3) + offset(3) + gauge_atom_radius)/dbox(3))
1648 592130 : IF (kmax - kmin + 1 .GT. nbox(3)) THEN
1649 527618 : kmin = 0
1650 527618 : kmax = nbox(3) - 1
1651 : END IF
1652 1735750 : DO kbox = kmin, kmax
1653 1143620 : jmin = INT((point(2) + offset(2) - gauge_atom_radius)/dbox(2))
1654 1143620 : jmax = INT((point(2) + offset(2) + gauge_atom_radius)/dbox(2))
1655 1143620 : IF (jmax - jmin + 1 .GT. nbox(2)) THEN
1656 1055236 : jmin = 0
1657 1055236 : jmax = nbox(2) - 1
1658 : END IF
1659 3967326 : DO jbox = jmin, jmax
1660 2231576 : imin = INT((point(1) + offset(1) - gauge_atom_radius)/dbox(1))
1661 2231576 : imax = INT((point(1) + offset(1) + gauge_atom_radius)/dbox(1))
1662 2231576 : IF (imax - imin + 1 .GT. nbox(1)) THEN
1663 2110472 : imin = 0
1664 2110472 : imax = nbox(1) - 1
1665 : END IF
1666 7762096 : DO ibox = imin, imax
1667 4386900 : ibeg = box_ptr(MODULO(ibox, nbox(1)), MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3)))
1668 4386900 : iend = box_ptr(MODULO(ibox, nbox(1)) + 1, MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3))) - 1
1669 9170920 : DO iatom = ibeg, iend
1670 20419552 : r(:) = pbc(box_data(:, iatom) - point(:), cell) + point(:)
1671 2552444 : dist = (r(1) - point(1))**2 + (r(2) - point(2))**2 + (r(3) - point(3))**2
1672 6939344 : IF (dist .LT. gauge_atom_radius**2) THEN
1673 2505134 : natms_local = natms_local + 1
1674 10020536 : ratom(:, natms_local) = r(:)
1675 : !
1676 : ! compute the distance atoms-point
1677 2505134 : x = point(1) - r(1)
1678 2505134 : y = point(2) - r(2)
1679 2505134 : z = point(3) - r(3)
1680 2505134 : atms_pnt(1, natms_local) = x
1681 2505134 : atms_pnt(2, natms_local) = y
1682 2505134 : atms_pnt(3, natms_local) = z
1683 2505134 : nrm_atms_pnt(natms_local) = SQRT(x*x + y*y + z*z)
1684 : END IF
1685 : END DO
1686 : END DO
1687 : END DO
1688 : END DO
1689 : !
1690 616968 : IF (natms_local .GT. 0) THEN
1691 : !
1692 : !
1693 3034004 : DO iatom = 1, natms_local
1694 2505134 : buf_tmp = 1.0_dp
1695 2505134 : pra(1) = atms_pnt(1, iatom)
1696 2505134 : pra(2) = atms_pnt(2, iatom)
1697 2505134 : pra(3) = atms_pnt(3, iatom)
1698 2505134 : pa = nrm_atms_pnt(iatom)
1699 15477460 : DO jatom = 1, natms_local
1700 12972326 : IF (iatom .EQ. jatom) CYCLE
1701 10467192 : pb = nrm_atms_pnt(jatom)
1702 10467192 : x = pra(1) - atms_pnt(1, jatom)
1703 10467192 : y = pra(2) - atms_pnt(2, jatom)
1704 10467192 : z = pra(3) - atms_pnt(3, jatom)
1705 10467192 : ab = SQRT(x*x + y*y + z*z)
1706 : !
1707 10467192 : tmp = (pa - pb)/ab
1708 10467192 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1709 10467192 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1710 10467192 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
1711 15477460 : buf_tmp = buf_tmp*0.5_dp*(1.0_dp - tmp)
1712 : END DO
1713 3034004 : buf(iatom) = buf_tmp
1714 : END DO
1715 : res(1) = 0.0_dp
1716 : res(2) = 0.0_dp
1717 : res(3) = 0.0_dp
1718 : summe = 0.0_dp
1719 3034004 : DO iatom = 1, natms_local
1720 2505134 : res(1) = res(1) + ratom(1, iatom)*buf(iatom)
1721 2505134 : res(2) = res(2) + ratom(2, iatom)*buf(iatom)
1722 2505134 : res(3) = res(3) + ratom(3, iatom)*buf(iatom)
1723 3034004 : summe = summe + buf(iatom)
1724 : END DO
1725 528870 : res(1) = res(1)/summe
1726 528870 : res(2) = res(2)/summe
1727 528870 : res(3) = res(3)/summe
1728 528870 : grid_x(i, j, k) = point(1) - res(1)
1729 528870 : grid_y(i, j, k) = point(2) - res(2)
1730 528870 : grid_z(i, j, k) = point(3) - res(3)
1731 : ELSE
1732 63260 : grid_x(i, j, k) = 0.0_dp
1733 63260 : grid_y(i, j, k) = 0.0_dp
1734 63260 : grid_z(i, j, k) = 0.0_dp
1735 : END IF
1736 : END DO
1737 : END DO
1738 : END DO
1739 :
1740 88 : DEALLOCATE (buf, ratom, atms_pnt, nrm_atms_pnt)
1741 :
1742 88 : END SUBROUTINE collocate_gauge
1743 :
1744 : ! **************************************************************************************************
1745 : !> \brief ...
1746 : !> \param current_env ...
1747 : !> \param qs_env ...
1748 : !> \param box ...
1749 : ! **************************************************************************************************
1750 6 : SUBROUTINE box_atoms_new(current_env, qs_env, box)
1751 : TYPE(current_env_type) :: current_env
1752 : TYPE(qs_environment_type), POINTER :: qs_env
1753 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
1754 :
1755 : CHARACTER(LEN=*), PARAMETER :: routineN = 'box_atoms_new'
1756 :
1757 : INTEGER :: handle, i, iatom, ibeg, ibox, iend, &
1758 : ifind, ii, imax, imin, j, jatom, jbox, &
1759 : jmax, jmin, k, kbox, kmax, kmin, &
1760 : natms, natms_local
1761 : REAL(dp) :: gauge_atom_radius, offset(3), scale
1762 6 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: ratom
1763 6 : REAL(dp), DIMENSION(:, :), POINTER :: r_ptr
1764 : REAL(kind=dp) :: box_center(3), box_center_wrap(3), &
1765 : box_size_guess, r(3)
1766 : TYPE(cell_type), POINTER :: cell
1767 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1768 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1769 :
1770 6 : CALL timeset(routineN, handle)
1771 :
1772 : CALL get_qs_env(qs_env=qs_env, &
1773 : qs_kind_set=qs_kind_set, &
1774 : cell=cell, &
1775 6 : particle_set=particle_set)
1776 :
1777 : CALL get_current_env(current_env=current_env, &
1778 6 : gauge_atom_radius=gauge_atom_radius)
1779 :
1780 6 : scale = 2.0_dp
1781 :
1782 6 : box_size_guess = gauge_atom_radius/scale
1783 :
1784 6 : natms = SIZE(particle_set, 1)
1785 18 : ALLOCATE (ratom(3, natms))
1786 :
1787 : !
1788 : ! box the atoms
1789 6 : nbox(1) = CEILING(cell%hmat(1, 1)/box_size_guess)
1790 6 : nbox(2) = CEILING(cell%hmat(2, 2)/box_size_guess)
1791 6 : nbox(3) = CEILING(cell%hmat(3, 3)/box_size_guess)
1792 6 : dbox(1) = cell%hmat(1, 1)/REAL(nbox(1), dp)
1793 6 : dbox(2) = cell%hmat(2, 2)/REAL(nbox(2), dp)
1794 6 : dbox(3) = cell%hmat(3, 3)/REAL(nbox(3), dp)
1795 36 : ALLOCATE (box_ptr(0:nbox(1), 0:nbox(2) - 1, 0:nbox(3) - 1), box_data(3, natms))
1796 62 : box_data(:, :) = HUGE(0.0_dp)
1797 27054 : box_ptr(:, :, :) = HUGE(0)
1798 : !
1799 6 : offset(1) = cell%hmat(1, 1)*0.5_dp
1800 6 : offset(2) = cell%hmat(2, 2)*0.5_dp
1801 6 : offset(3) = cell%hmat(3, 3)*0.5_dp
1802 20 : DO iatom = 1, natms
1803 20 : ratom(:, iatom) = pbc(particle_set(iatom)%r(:), cell)
1804 : END DO
1805 : !
1806 6 : i = 1
1807 88 : DO kbox = 0, nbox(3) - 1
1808 1430 : DO jbox = 0, nbox(2) - 1
1809 1342 : box_ptr(0, jbox, kbox) = i
1810 25706 : DO ibox = 0, nbox(1) - 1
1811 24282 : ii = 0
1812 88846 : DO iatom = 1, natms
1813 : IF (MODULO(FLOOR(ratom(1, iatom)/dbox(1)), nbox(1)) .EQ. ibox .AND. &
1814 64564 : MODULO(FLOOR(ratom(2, iatom)/dbox(2)), nbox(2)) .EQ. jbox .AND. &
1815 24282 : MODULO(FLOOR(ratom(3, iatom)/dbox(3)), nbox(3)) .EQ. kbox) THEN
1816 56 : box_data(:, i) = ratom(:, iatom)
1817 14 : i = i + 1
1818 14 : ii = ii + 1
1819 : END IF
1820 : END DO
1821 25624 : box_ptr(ibox + 1, jbox, kbox) = box_ptr(ibox, jbox, kbox) + ii
1822 : END DO
1823 : END DO
1824 : END DO
1825 : !
1826 : IF (.FALSE.) THEN
1827 : DO kbox = 0, nbox(3) - 1
1828 : DO jbox = 0, nbox(2) - 1
1829 : DO ibox = 0, nbox(1) - 1
1830 : IF (box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox) .GT. 0) THEN
1831 : WRITE (*, *) 'box=', ibox, jbox, kbox
1832 : WRITE (*, *) 'nbr atom=', box_ptr(ibox + 1, jbox, kbox) - box_ptr(ibox, jbox, kbox)
1833 : DO iatom = box_ptr(ibox, jbox, kbox), box_ptr(ibox + 1, jbox, kbox) - 1
1834 : WRITE (*, '(A,I3,3E14.6)') 'coor=', iatom, box_data(:, iatom)
1835 : END DO
1836 : END IF
1837 : END DO
1838 : END DO
1839 : END DO
1840 : END IF
1841 : !
1842 6 : NULLIFY (box)
1843 25736 : ALLOCATE (box(0:nbox(1) - 1, 0:nbox(2) - 1, 0:nbox(3) - 1))
1844 : !
1845 : ! build the list
1846 88 : DO k = 0, nbox(3) - 1
1847 1430 : DO j = 0, nbox(2) - 1
1848 25706 : DO i = 0, nbox(1) - 1
1849 : !
1850 24282 : box_center(1) = (REAL(i, dp) + 0.5_dp)*dbox(1)
1851 24282 : box_center(2) = (REAL(j, dp) + 0.5_dp)*dbox(2)
1852 24282 : box_center(3) = (REAL(k, dp) + 0.5_dp)*dbox(3)
1853 24282 : box_center_wrap = pbc(box_center, cell)
1854 : !
1855 : ! find the atoms that are in the overlaping boxes
1856 24282 : natms_local = 0
1857 24282 : kmin = FLOOR((box_center(3) - gauge_atom_radius)/dbox(3))
1858 24282 : kmax = FLOOR((box_center(3) + gauge_atom_radius)/dbox(3))
1859 24282 : IF (kmax - kmin + 1 .GT. nbox(3)) THEN
1860 0 : kmin = 0
1861 0 : kmax = nbox(3) - 1
1862 : END IF
1863 145692 : DO kbox = kmin, kmax
1864 121410 : jmin = FLOOR((box_center(2) - gauge_atom_radius)/dbox(2))
1865 121410 : jmax = FLOOR((box_center(2) + gauge_atom_radius)/dbox(2))
1866 121410 : IF (jmax - jmin + 1 .GT. nbox(2)) THEN
1867 450 : jmin = 0
1868 450 : jmax = nbox(2) - 1
1869 : END IF
1870 751842 : DO jbox = jmin, jmax
1871 606150 : imin = FLOOR((box_center(1) - gauge_atom_radius)/dbox(1))
1872 606150 : imax = FLOOR((box_center(1) + gauge_atom_radius)/dbox(1))
1873 606150 : IF (imax - imin + 1 .GT. nbox(1)) THEN
1874 1350 : imin = 0
1875 1350 : imax = nbox(1) - 1
1876 : END IF
1877 3755610 : DO ibox = imin, imax
1878 3028050 : ibeg = box_ptr(MODULO(ibox, nbox(1)), MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3)))
1879 3028050 : iend = box_ptr(MODULO(ibox, nbox(1)) + 1, MODULO(jbox, nbox(2)), MODULO(kbox, nbox(3))) - 1
1880 3635630 : DO iatom = ibeg, iend
1881 5720 : r = pbc(box_center_wrap(:) - box_data(:, iatom), cell)
1882 : IF (ABS(r(1)) .LE. (scale + 0.5_dp)*dbox(1) .AND. &
1883 1430 : ABS(r(2)) .LE. (scale + 0.5_dp)*dbox(2) .AND. &
1884 3028050 : ABS(r(3)) .LE. (scale + 0.5_dp)*dbox(3)) THEN
1885 1430 : natms_local = natms_local + 1
1886 5720 : ratom(:, natms_local) = box_data(:, iatom)
1887 : END IF
1888 : END DO
1889 : END DO ! box
1890 : END DO
1891 : END DO
1892 : !
1893 : ! set the list
1894 24282 : box(i, j, k)%n = natms_local
1895 24282 : NULLIFY (box(i, j, k)%r)
1896 25624 : IF (natms_local .GT. 0) THEN
1897 3840 : ALLOCATE (box(i, j, k)%r(3, natms_local))
1898 1280 : r_ptr => box(i, j, k)%r
1899 1280 : CALL dcopy(3*natms_local, ratom(1, 1), 1, r_ptr(1, 1), 1)
1900 : END IF
1901 : END DO ! list
1902 : END DO
1903 : END DO
1904 :
1905 : IF (.FALSE.) THEN
1906 : DO k = 0, nbox(3) - 1
1907 : DO j = 0, nbox(2) - 1
1908 : DO i = 0, nbox(1) - 1
1909 : IF (box(i, j, k)%n .GT. 0) THEN
1910 : WRITE (*, *)
1911 : WRITE (*, *) 'box=', i, j, k
1912 : box_center(1) = (REAL(i, dp) + 0.5_dp)*dbox(1)
1913 : box_center(2) = (REAL(j, dp) + 0.5_dp)*dbox(2)
1914 : box_center(3) = (REAL(k, dp) + 0.5_dp)*dbox(3)
1915 : box_center = pbc(box_center, cell)
1916 : WRITE (*, '(A,3E14.6)') 'box_center=', box_center
1917 : WRITE (*, *) 'nbr atom=', box(i, j, k)%n
1918 : r_ptr => box(i, j, k)%r
1919 : DO iatom = 1, box(i, j, k)%n
1920 : WRITE (*, '(A,I3,3E14.6)') 'coor=', iatom, r_ptr(:, iatom)
1921 : r(:) = pbc(box_center(:) - r_ptr(:, iatom), cell)
1922 : IF (ABS(r(1)) .GT. (scale + 0.5_dp)*dbox(1) .OR. &
1923 : ABS(r(2)) .GT. (scale + 0.5_dp)*dbox(2) .OR. &
1924 : ABS(r(3)) .GT. (scale + 0.5_dp)*dbox(3)) THEN
1925 : WRITE (*, *) 'error too many atoms'
1926 : WRITE (*, *) 'dist=', ABS(r(:))
1927 : WRITE (*, *) 'large_dist=', (scale + 0.5_dp)*dbox
1928 : CPABORT("")
1929 : END IF
1930 : END DO
1931 : END IF
1932 : END DO ! list
1933 : END DO
1934 : END DO
1935 : END IF
1936 :
1937 : IF (.TRUE.) THEN
1938 88 : DO k = 0, nbox(3) - 1
1939 1430 : DO j = 0, nbox(2) - 1
1940 25706 : DO i = 0, nbox(1) - 1
1941 24282 : box_center(1) = (REAL(i, dp) + 0.5_dp)*dbox(1)
1942 24282 : box_center(2) = (REAL(j, dp) + 0.5_dp)*dbox(2)
1943 24282 : box_center(3) = (REAL(k, dp) + 0.5_dp)*dbox(3)
1944 97128 : box_center = pbc(box_center, cell)
1945 24282 : r_ptr => box(i, j, k)%r
1946 90188 : DO iatom = 1, natms
1947 258256 : r(:) = pbc(box_center(:) - ratom(:, iatom), cell)
1948 64564 : ifind = 0
1949 68174 : DO jatom = 1, box(i, j, k)%n
1950 79004 : IF (SUM(ABS(ratom(:, iatom) - r_ptr(:, jatom))) .LT. 1E-10_dp) ifind = 1
1951 : END DO
1952 :
1953 88846 : IF (ifind .EQ. 0) THEN
1954 : ! SQRT(DOT_PRODUCT(r, r)) .LT. gauge_atom_radius
1955 252536 : IF (DOT_PRODUCT(r, r) .LT. (gauge_atom_radius*gauge_atom_radius)) THEN
1956 0 : WRITE (*, *) 'error atom too close'
1957 0 : WRITE (*, *) 'iatom', iatom
1958 0 : WRITE (*, *) 'box_center', box_center
1959 0 : WRITE (*, *) 'ratom', ratom(:, iatom)
1960 0 : WRITE (*, *) 'gauge_atom_radius', gauge_atom_radius
1961 0 : CPABORT("")
1962 : END IF
1963 : END IF
1964 : END DO
1965 : END DO ! list
1966 : END DO
1967 : END DO
1968 : END IF
1969 :
1970 6 : DEALLOCATE (ratom)
1971 :
1972 6 : CALL timestop(handle)
1973 :
1974 12 : END SUBROUTINE box_atoms_new
1975 :
1976 : ! **************************************************************************************************
1977 : !> \brief ...
1978 : !> \param current_env ...
1979 : !> \param qs_env ...
1980 : !> \param rs_grid_x ...
1981 : !> \param rs_grid_y ...
1982 : !> \param rs_grid_z ...
1983 : !> \param box ...
1984 : ! **************************************************************************************************
1985 20 : SUBROUTINE collocate_gauge_new(current_env, qs_env, rs_grid_x, rs_grid_y, rs_grid_z, box)
1986 : !
1987 : TYPE(current_env_type) :: current_env
1988 : TYPE(qs_environment_type), POINTER :: qs_env
1989 : TYPE(realspace_grid_type), INTENT(IN) :: rs_grid_x, rs_grid_y, rs_grid_z
1990 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
1991 :
1992 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_gauge_new'
1993 :
1994 : INTEGER :: delta_lb(3), handle, i, iatom, ib, ibe, ibox, ibs, ie, is, j, jatom, jb, jbe, &
1995 : jbox, jbs, je, js, k, kb, kbe, kbox, kbs, ke, ks, lb(3), lb_local(3), natms, &
1996 : natms_local0, natms_local1, ng(3)
1997 20 : REAL(dp), DIMENSION(:, :), POINTER :: r_ptr
1998 : REAL(KIND=dp) :: ab, box_center(3), buf_tmp, dist, dr(3), &
1999 : gauge_atom_radius, offset(3), pa, pb, &
2000 : point(3), pra(3), r(3), res(3), summe, &
2001 : tmp, x, y, z
2002 20 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: buf, nrm_atms_pnt
2003 20 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: atms_pnt, ratom
2004 20 : REAL(kind=dp), DIMENSION(:, :, :), POINTER :: grid_x, grid_y, grid_z
2005 : TYPE(cell_type), POINTER :: cell
2006 20 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2007 20 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2008 :
2009 20 : CALL timeset(routineN, handle)
2010 :
2011 : !
2012 : CALL get_current_env(current_env=current_env, &
2013 20 : gauge_atom_radius=gauge_atom_radius)
2014 : !
2015 : CALL get_qs_env(qs_env=qs_env, &
2016 : qs_kind_set=qs_kind_set, &
2017 : cell=cell, &
2018 20 : particle_set=particle_set)
2019 : !
2020 20 : natms = SIZE(particle_set, 1)
2021 20 : dr(1) = rs_grid_x%desc%dh(1, 1)
2022 20 : dr(2) = rs_grid_x%desc%dh(2, 2)
2023 20 : dr(3) = rs_grid_x%desc%dh(3, 3)
2024 80 : lb(:) = rs_grid_x%desc%lb(:)
2025 80 : lb_local(:) = rs_grid_x%lb_local(:)
2026 20 : grid_x => rs_grid_x%r(:, :, :)
2027 20 : grid_y => rs_grid_y%r(:, :, :)
2028 20 : grid_z => rs_grid_z%r(:, :, :)
2029 80 : ng(:) = UBOUND(grid_x)
2030 80 : delta_lb(:) = lb_local(:) - lb(:)
2031 20 : offset(1) = cell%hmat(1, 1)*0.5_dp
2032 20 : offset(2) = cell%hmat(2, 2)*0.5_dp
2033 20 : offset(3) = cell%hmat(3, 3)*0.5_dp
2034 140 : ALLOCATE (buf(natms), ratom(3, natms), atms_pnt(3, natms), nrm_atms_pnt(natms))
2035 : !
2036 : ! find the boxes that match the grid
2037 20 : ibs = FLOOR(REAL(delta_lb(1), dp)*dr(1)/dbox(1))
2038 20 : ibe = FLOOR(REAL(ng(1) - 1 + delta_lb(1), dp)*dr(1)/dbox(1))
2039 20 : jbs = FLOOR(REAL(delta_lb(2), dp)*dr(2)/dbox(2))
2040 20 : jbe = FLOOR(REAL(ng(2) - 1 + delta_lb(2), dp)*dr(2)/dbox(2))
2041 20 : kbs = FLOOR(REAL(delta_lb(3), dp)*dr(3)/dbox(3))
2042 20 : kbe = FLOOR(REAL(ng(3) - 1 + delta_lb(3), dp)*dr(3)/dbox(3))
2043 : !
2044 : ! go over the box-list
2045 314 : DO kb = kbs, kbe
2046 5148 : DO jb = jbs, jbe
2047 89870 : DO ib = ibs, ibe
2048 84742 : ibox = MODULO(ib, nbox(1))
2049 84742 : jbox = MODULO(jb, nbox(2))
2050 84742 : kbox = MODULO(kb, nbox(3))
2051 : !
2052 84742 : is = MAX(CEILING(REAL(ib, dp)*dbox(1)/dr(1)), delta_lb(1)) - delta_lb(1) + 1
2053 84742 : ie = MIN(FLOOR(REAL(ib + 1, dp)*dbox(1)/dr(1)), ng(1) - 1 + delta_lb(1)) - delta_lb(1) + 1
2054 84742 : js = MAX(CEILING(REAL(jb, dp)*dbox(2)/dr(2)), delta_lb(2)) - delta_lb(2) + 1
2055 84742 : je = MIN(FLOOR(REAL(jb + 1, dp)*dbox(2)/dr(2)), ng(2) - 1 + delta_lb(2)) - delta_lb(2) + 1
2056 84742 : ks = MAX(CEILING(REAL(kb, dp)*dbox(3)/dr(3)), delta_lb(3)) - delta_lb(3) + 1
2057 84742 : ke = MIN(FLOOR(REAL(kb + 1, dp)*dbox(3)/dr(3)), ng(3) - 1 + delta_lb(3)) - delta_lb(3) + 1
2058 : !
2059 : ! sanity checks
2060 : IF (.TRUE.) THEN
2061 84742 : IF (REAL(ks - 1 + delta_lb(3), dp)*dr(3) .LT. REAL(kb, dp)*dbox(3) .OR. &
2062 : REAL(ke - 1 + delta_lb(3), dp)*dr(3) .GT. REAL(kb + 1, dp)*dbox(3)) THEN
2063 0 : WRITE (*, *) 'box_k', REAL(kb, dp)*dbox(3), REAL(kb + 1, dp)*dbox(3)
2064 0 : WRITE (*, *) 'point_k', REAL(ks - 1 + delta_lb(3), dp)*dr(3), REAL(ke - 1 + delta_lb(3), dp)*dr(3)
2065 0 : WRITE (*, *) 'ibox', ibox, 'jbox', jbox, 'kbox', kbox
2066 0 : WRITE (*, *) 'is,ie', is, ie, ' js,je', js, je, ' ks,ke', ks, ke
2067 0 : WRITE (*, *) 'ibs,ibe', ibs, ibe, ' jbs,jbe', jbs, jbe, ' kbs,kbe', kbs, kbe
2068 0 : CPABORT("we stop_k")
2069 : END IF
2070 84742 : IF (REAL(js - 1 + delta_lb(2), dp)*dr(2) .LT. REAL(jb, dp)*dbox(2) .OR. &
2071 : REAL(je - 1 + delta_lb(2), dp)*dr(2) .GT. REAL(jb + 1, dp)*dbox(2)) THEN
2072 0 : WRITE (*, *) 'box_j', REAL(jb, dp)*dbox(2), REAL(jb + 1, dp)*dbox(2)
2073 0 : WRITE (*, *) 'point_j', REAL(js - 1 + delta_lb(2), dp)*dr(2), REAL(je - 1 + delta_lb(2), dp)*dr(2)
2074 0 : WRITE (*, *) 'is,ie', is, ie, ' js,je', js, je, ' ks,ke', ks, ke
2075 0 : WRITE (*, *) 'ibs,ibe', ibs, ibe, ' jbs,jbe', jbs, jbe, ' kbs,kbe', kbs, kbe
2076 0 : CPABORT("we stop_j")
2077 : END IF
2078 84742 : IF (REAL(is - 1 + delta_lb(1), dp)*dr(1) .LT. REAL(ib, dp)*dbox(1) .OR. &
2079 : REAL(ie - 1 + delta_lb(1), dp)*dr(1) .GT. REAL(ib + 1, dp)*dbox(1)) THEN
2080 0 : WRITE (*, *) 'box_i', REAL(ib, dp)*dbox(1), REAL(ib + 1, dp)*dbox(1)
2081 0 : WRITE (*, *) 'point_i', REAL(is - 1 + delta_lb(1), dp)*dr(1), REAL(ie - 1 + delta_lb(1), dp)*dr(1)
2082 0 : WRITE (*, *) 'is,ie', is, ie, ' js,je', js, je, ' ks,ke', ks, ke
2083 0 : WRITE (*, *) 'ibs,ibe', ibs, ibe, ' jbs,jbe', jbs, jbe, ' kbs,kbe', kbs, kbe
2084 0 : CPABORT("we stop_i")
2085 : END IF
2086 : END IF
2087 : !
2088 : ! the center of the box
2089 84742 : box_center(1) = (REAL(ibox, dp) + 0.5_dp)*dbox(1)
2090 84742 : box_center(2) = (REAL(jbox, dp) + 0.5_dp)*dbox(2)
2091 84742 : box_center(3) = (REAL(kbox, dp) + 0.5_dp)*dbox(3)
2092 : !
2093 : ! find the atoms that are in the overlaping boxes
2094 84742 : natms_local0 = box(ibox, jbox, kbox)%n
2095 84742 : r_ptr => box(ibox, jbox, kbox)%r
2096 : !
2097 : ! go over the grid inside the box
2098 89576 : IF (natms_local0 .GT. 0) THEN
2099 : !
2100 : ! here there are some atoms...
2101 9850 : DO k = ks, ke
2102 31112 : DO j = js, je
2103 154636 : DO i = is, ie
2104 127152 : point(1) = REAL(i - 1 + delta_lb(1), dp)*dr(1)
2105 127152 : point(2) = REAL(j - 1 + delta_lb(2), dp)*dr(2)
2106 127152 : point(3) = REAL(k - 1 + delta_lb(3), dp)*dr(3)
2107 508608 : point = pbc(point, cell)
2108 : !
2109 : ! compute atom-point distances
2110 127152 : natms_local1 = 0
2111 364868 : DO iatom = 1, natms_local0
2112 1901728 : r(:) = pbc(r_ptr(:, iatom) - point(:), cell) + point(:) !needed?
2113 237716 : dist = (r(1) - point(1))**2 + (r(2) - point(2))**2 + (r(3) - point(3))**2
2114 364868 : IF (dist .LT. gauge_atom_radius**2) THEN
2115 154124 : natms_local1 = natms_local1 + 1
2116 616496 : ratom(:, natms_local1) = r(:)
2117 : !
2118 : ! compute the distance atoms-point
2119 154124 : x = point(1) - r(1)
2120 154124 : y = point(2) - r(2)
2121 154124 : z = point(3) - r(3)
2122 154124 : atms_pnt(1, natms_local1) = x
2123 154124 : atms_pnt(2, natms_local1) = y
2124 154124 : atms_pnt(3, natms_local1) = z
2125 154124 : nrm_atms_pnt(natms_local1) = SQRT(x*x + y*y + z*z)
2126 : END IF
2127 : END DO
2128 : !
2129 : !
2130 148414 : IF (natms_local1 .GT. 0) THEN
2131 : !
2132 : ! build the step
2133 238616 : DO iatom = 1, natms_local1
2134 154124 : buf_tmp = 1.0_dp
2135 154124 : pra(1) = atms_pnt(1, iatom)
2136 154124 : pra(2) = atms_pnt(2, iatom)
2137 154124 : pra(3) = atms_pnt(3, iatom)
2138 154124 : pa = nrm_atms_pnt(iatom)
2139 447512 : DO jatom = 1, natms_local1
2140 293388 : IF (iatom .EQ. jatom) CYCLE
2141 139264 : pb = nrm_atms_pnt(jatom)
2142 139264 : x = pra(1) - atms_pnt(1, jatom)
2143 139264 : y = pra(2) - atms_pnt(2, jatom)
2144 139264 : z = pra(3) - atms_pnt(3, jatom)
2145 139264 : ab = SQRT(x*x + y*y + z*z)
2146 : !
2147 139264 : tmp = (pa - pb)/ab
2148 139264 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2149 139264 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2150 139264 : tmp = 0.5_dp*(3.0_dp - tmp*tmp)*tmp
2151 447512 : buf_tmp = buf_tmp*0.5_dp*(1.0_dp - tmp)
2152 : END DO
2153 238616 : buf(iatom) = buf_tmp
2154 : END DO
2155 : res(1) = 0.0_dp
2156 : res(2) = 0.0_dp
2157 : res(3) = 0.0_dp
2158 : summe = 0.0_dp
2159 238616 : DO iatom = 1, natms_local1
2160 154124 : res(1) = res(1) + ratom(1, iatom)*buf(iatom)
2161 154124 : res(2) = res(2) + ratom(2, iatom)*buf(iatom)
2162 154124 : res(3) = res(3) + ratom(3, iatom)*buf(iatom)
2163 238616 : summe = summe + buf(iatom)
2164 : END DO
2165 84492 : res(1) = res(1)/summe
2166 84492 : res(2) = res(2)/summe
2167 84492 : res(3) = res(3)/summe
2168 84492 : grid_x(i, j, k) = point(1) - res(1)
2169 84492 : grid_y(i, j, k) = point(2) - res(2)
2170 84492 : grid_z(i, j, k) = point(3) - res(3)
2171 : ELSE
2172 42660 : grid_x(i, j, k) = 0.0_dp
2173 42660 : grid_y(i, j, k) = 0.0_dp
2174 42660 : grid_z(i, j, k) = 0.0_dp
2175 : END IF
2176 : END DO ! grid
2177 : END DO
2178 : END DO
2179 : !
2180 : ELSE
2181 : !
2182 : ! here there is no atom
2183 187142 : DO k = ks, ke
2184 376598 : DO j = js, je
2185 717810 : DO i = is, ie
2186 422326 : grid_x(i, j, k) = 0.0_dp
2187 422326 : grid_y(i, j, k) = 0.0_dp
2188 611782 : grid_z(i, j, k) = 0.0_dp
2189 : END DO ! grid
2190 : END DO
2191 : END DO
2192 : !
2193 : END IF
2194 : !
2195 : END DO ! list
2196 : END DO
2197 : END DO
2198 :
2199 20 : DEALLOCATE (buf, ratom, atms_pnt, nrm_atms_pnt)
2200 :
2201 20 : CALL timestop(handle)
2202 :
2203 20 : END SUBROUTINE collocate_gauge_new
2204 :
2205 : ! **************************************************************************************************
2206 : !> \brief ...
2207 : !> \param box ...
2208 : ! **************************************************************************************************
2209 28 : SUBROUTINE deallocate_box(box)
2210 : TYPE(box_type), DIMENSION(:, :, :), POINTER :: box
2211 :
2212 : INTEGER :: i, j, k
2213 :
2214 28 : IF (ASSOCIATED(box)) THEN
2215 100 : DO k = LBOUND(box, 3), UBOUND(box, 3)
2216 1594 : DO j = LBOUND(box, 2), UBOUND(box, 2)
2217 28390 : DO i = LBOUND(box, 1), UBOUND(box, 1)
2218 25624 : IF (ASSOCIATED(box(i, j, k)%r)) THEN
2219 1280 : DEALLOCATE (box(i, j, k)%r)
2220 : END IF
2221 : END DO
2222 : END DO
2223 : END DO
2224 6 : DEALLOCATE (box)
2225 : END IF
2226 28 : END SUBROUTINE deallocate_box
2227 : END SUBROUTINE current_set_gauge
2228 :
2229 : ! **************************************************************************************************
2230 : !> \brief ...
2231 : !> \param current_env ...
2232 : !> \param qs_env ...
2233 : !> \param iB ...
2234 : ! **************************************************************************************************
2235 522 : SUBROUTINE current_build_chi(current_env, qs_env, iB)
2236 : !
2237 : TYPE(current_env_type) :: current_env
2238 : TYPE(qs_environment_type), POINTER :: qs_env
2239 : INTEGER, INTENT(IN) :: iB
2240 :
2241 522 : IF (current_env%full) THEN
2242 426 : CALL current_build_chi_many_centers(current_env, qs_env, iB)
2243 96 : ELSE IF (current_env%nbr_center(1) > 1) THEN
2244 0 : CALL current_build_chi_many_centers(current_env, qs_env, iB)
2245 : ELSE
2246 96 : CALL current_build_chi_one_center(current_env, qs_env, iB)
2247 : END IF
2248 :
2249 522 : END SUBROUTINE current_build_chi
2250 :
2251 : ! **************************************************************************************************
2252 : !> \brief ...
2253 : !> \param current_env ...
2254 : !> \param qs_env ...
2255 : !> \param iB ...
2256 : ! **************************************************************************************************
2257 426 : SUBROUTINE current_build_chi_many_centers(current_env, qs_env, iB)
2258 : !
2259 : TYPE(current_env_type) :: current_env
2260 : TYPE(qs_environment_type), POINTER :: qs_env
2261 : INTEGER, INTENT(IN) :: iB
2262 :
2263 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_chi_many_centers'
2264 :
2265 : INTEGER :: handle, icenter, idir, idir2, ii, iiB, iii, iiiB, ispin, istate, j, jstate, &
2266 : max_states, nao, natom, nbr_center(2), nmo, nspins, nstate_loc, nstates(2), output_unit
2267 426 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
2268 426 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
2269 : LOGICAL :: chi_pbc, gapw
2270 : REAL(dp) :: chi(3), chi_tmp, contrib, contrib2, &
2271 : dk(3), int_current(3), &
2272 : int_current_tmp, maxocc
2273 : TYPE(cell_type), POINTER :: cell
2274 426 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
2275 426 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
2276 : TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct
2277 : TYPE(cp_fm_type) :: psi0, psi_D, psi_p1, psi_p2, psi_rxp
2278 4686 : TYPE(cp_fm_type), DIMENSION(3) :: p_rxp, r_p1, r_p2
2279 38766 : TYPE(cp_fm_type), DIMENSION(9, 3) :: rr_p1, rr_p2, rr_rxp
2280 426 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
2281 426 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_D, psi1_p, psi1_rxp
2282 : TYPE(cp_fm_type), POINTER :: mo_coeff
2283 : TYPE(cp_logger_type), POINTER :: logger
2284 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
2285 426 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_mom_ao, op_p_ao
2286 426 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_mom_der_ao
2287 : TYPE(dft_control_type), POINTER :: dft_control
2288 426 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2289 : TYPE(mp_para_env_type), POINTER :: para_env
2290 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2291 426 : POINTER :: sab_all, sab_orb
2292 426 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2293 426 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2294 :
2295 : !
2296 :
2297 426 : CALL timeset(routineN, handle)
2298 : !
2299 426 : NULLIFY (dft_control, mos, para_env, mo_coeff, op_mom_ao, &
2300 426 : op_mom_der_ao, center_list, centers_set, &
2301 426 : op_p_ao, psi1_p, psi1_rxp, psi1_D, &
2302 426 : cell, particle_set, qs_kind_set)
2303 :
2304 426 : logger => cp_get_default_logger()
2305 426 : output_unit = cp_logger_get_default_io_unit(logger)
2306 :
2307 : CALL get_qs_env(qs_env=qs_env, &
2308 : dft_control=dft_control, &
2309 : mos=mos, &
2310 : para_env=para_env, &
2311 : cell=cell, &
2312 : dbcsr_dist=dbcsr_dist, &
2313 : particle_set=particle_set, &
2314 : qs_kind_set=qs_kind_set, &
2315 : sab_all=sab_all, &
2316 426 : sab_orb=sab_orb)
2317 :
2318 426 : nspins = dft_control%nspins
2319 426 : gapw = dft_control%qs_control%gapw
2320 :
2321 : CALL get_current_env(current_env=current_env, &
2322 : chi_pbc=chi_pbc, &
2323 : nao=nao, &
2324 : nbr_center=nbr_center, &
2325 : center_list=center_list, &
2326 : centers_set=centers_set, &
2327 : psi1_p=psi1_p, &
2328 : psi1_rxp=psi1_rxp, &
2329 : psi1_D=psi1_D, &
2330 : nstates=nstates, &
2331 426 : psi0_order=psi0_order)
2332 : !
2333 : ! get max nbr of states per center
2334 426 : max_states = 0
2335 1044 : DO ispin = 1, nspins
2336 4662 : DO icenter = 1, nbr_center(ispin)
2337 : max_states = MAX(max_states, center_list(ispin)%array(1, icenter + 1)&
2338 4236 : & - center_list(ispin)%array(1, icenter))
2339 : END DO
2340 : END DO
2341 : !
2342 : ! Allocate sparse matrices for dipole, quadrupole and their derivatives => 9x3
2343 : ! Remember the derivatives are antisymmetric
2344 426 : CALL dbcsr_allocate_matrix_set(op_mom_ao, 9)
2345 426 : CALL dbcsr_allocate_matrix_set(op_mom_der_ao, 9, 3)
2346 : !
2347 : ! prepare for allocation
2348 426 : natom = SIZE(particle_set, 1)
2349 1278 : ALLOCATE (first_sgf(natom))
2350 852 : ALLOCATE (last_sgf(natom))
2351 : CALL get_particle_set(particle_set, qs_kind_set, &
2352 : first_sgf=first_sgf, &
2353 426 : last_sgf=last_sgf)
2354 852 : ALLOCATE (row_blk_sizes(natom))
2355 426 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
2356 426 : DEALLOCATE (first_sgf)
2357 426 : DEALLOCATE (last_sgf)
2358 : !
2359 : !
2360 426 : ALLOCATE (op_mom_ao(1)%matrix)
2361 : CALL dbcsr_create(matrix=op_mom_ao(1)%matrix, &
2362 : name="op_mom", &
2363 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
2364 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2365 426 : mutable_work=.TRUE.)
2366 426 : CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix, sab_all)
2367 :
2368 1704 : DO idir2 = 1, 3
2369 1278 : ALLOCATE (op_mom_der_ao(1, idir2)%matrix)
2370 : CALL dbcsr_copy(op_mom_der_ao(1, idir2)%matrix, op_mom_ao(1)%matrix, &
2371 1704 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir2))))
2372 : END DO
2373 :
2374 3834 : DO idir = 2, SIZE(op_mom_ao, 1)
2375 3408 : ALLOCATE (op_mom_ao(idir)%matrix)
2376 : CALL dbcsr_copy(op_mom_ao(idir)%matrix, op_mom_ao(1)%matrix, &
2377 3408 : "op_mom_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2378 14058 : DO idir2 = 1, 3
2379 10224 : ALLOCATE (op_mom_der_ao(idir, idir2)%matrix)
2380 : CALL dbcsr_copy(op_mom_der_ao(idir, idir2)%matrix, op_mom_ao(1)%matrix, &
2381 13632 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir*idir2))))
2382 : END DO
2383 : END DO
2384 : !
2385 426 : CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
2386 426 : ALLOCATE (op_p_ao(1)%matrix)
2387 : CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
2388 : name="op_p_ao", &
2389 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
2390 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2391 426 : mutable_work=.TRUE.)
2392 426 : CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
2393 :
2394 1278 : DO idir = 2, 3
2395 852 : ALLOCATE (op_p_ao(idir)%matrix)
2396 : CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
2397 1278 : "op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2398 : END DO
2399 : !
2400 : !
2401 426 : DEALLOCATE (row_blk_sizes)
2402 : !
2403 : !
2404 : ! Allocate full matrices for only one vector
2405 426 : mo_coeff => psi0_order(1)
2406 426 : NULLIFY (tmp_fm_struct)
2407 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
2408 : ncol_global=max_states, para_env=para_env, &
2409 426 : context=mo_coeff%matrix_struct%context)
2410 426 : CALL cp_fm_create(psi0, tmp_fm_struct)
2411 426 : CALL cp_fm_create(psi_D, tmp_fm_struct)
2412 426 : CALL cp_fm_create(psi_rxp, tmp_fm_struct)
2413 426 : CALL cp_fm_create(psi_p1, tmp_fm_struct)
2414 426 : CALL cp_fm_create(psi_p2, tmp_fm_struct)
2415 426 : CALL cp_fm_struct_release(tmp_fm_struct)
2416 : !
2417 : CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nao, &
2418 : ncol_global=max_states, para_env=para_env, &
2419 426 : context=mo_coeff%matrix_struct%context)
2420 1704 : DO idir = 1, 3
2421 1278 : CALL cp_fm_create(p_rxp(idir), tmp_fm_struct)
2422 1278 : CALL cp_fm_create(r_p1(idir), tmp_fm_struct)
2423 1278 : CALL cp_fm_create(r_p2(idir), tmp_fm_struct)
2424 13206 : DO idir2 = 1, 9
2425 11502 : CALL cp_fm_create(rr_rxp(idir2, idir), tmp_fm_struct)
2426 11502 : CALL cp_fm_create(rr_p1(idir2, idir), tmp_fm_struct)
2427 12780 : CALL cp_fm_create(rr_p2(idir2, idir), tmp_fm_struct)
2428 : END DO
2429 : END DO
2430 426 : CALL cp_fm_struct_release(tmp_fm_struct)
2431 : !
2432 : !
2433 : !
2434 : ! recompute the linear momentum matrices
2435 426 : CALL build_lin_mom_matrix(qs_env, op_p_ao)
2436 : !CALL p_xyz_ao(op_p_ao,qs_env,minimum_image=.FALSE.)
2437 : !
2438 : !
2439 : ! get iiB and iiiB
2440 426 : CALL set_vecp(iB, iiB, iiiB)
2441 1044 : DO ispin = 1, nspins
2442 : !
2443 : ! get ground state MOS
2444 618 : nmo = nstates(ispin)
2445 618 : mo_coeff => psi0_order(ispin)
2446 618 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
2447 : !
2448 : ! Initialize the temporary vector chi
2449 618 : chi = 0.0_dp
2450 : int_current = 0.0_dp
2451 : !
2452 : ! Start loop over the occupied states
2453 4236 : DO icenter = 1, nbr_center(ispin)
2454 : !
2455 : ! Get the Wannier center of the istate-th ground state orbital
2456 14472 : dk(1:3) = centers_set(ispin)%array(1:3, icenter)
2457 : !
2458 : ! Compute the multipole integrals for the state istate,
2459 : ! using as reference center the corresponding Wannier center
2460 36180 : DO idir = 1, 9
2461 32562 : CALL dbcsr_set(op_mom_ao(idir)%matrix, 0.0_dp)
2462 133866 : DO idir2 = 1, 3
2463 130248 : CALL dbcsr_set(op_mom_der_ao(idir, idir2)%matrix, 0.0_dp)
2464 : END DO
2465 : END DO
2466 : CALL rRc_xyz_der_ao(op_mom_ao, op_mom_der_ao, qs_env, dk, order=2, &
2467 3618 : minimum_image=.FALSE., soft=gapw)
2468 : !
2469 : ! collecte the states that belong to a given center
2470 3618 : CALL cp_fm_set_all(psi0, 0.0_dp)
2471 3618 : CALL cp_fm_set_all(psi_rxp, 0.0_dp)
2472 3618 : CALL cp_fm_set_all(psi_D, 0.0_dp)
2473 3618 : CALL cp_fm_set_all(psi_p1, 0.0_dp)
2474 3618 : CALL cp_fm_set_all(psi_p2, 0.0_dp)
2475 3618 : nstate_loc = center_list(ispin)%array(1, icenter + 1) - center_list(ispin)%array(1, icenter)
2476 3618 : jstate = 1
2477 7614 : DO j = center_list(ispin)%array(1, icenter), center_list(ispin)%array(1, icenter + 1) - 1
2478 3996 : istate = center_list(ispin)%array(2, j)
2479 : !
2480 : ! block the states that belong to this center
2481 3996 : CALL cp_fm_to_fm(mo_coeff, psi0, 1, istate, jstate)
2482 : !
2483 3996 : CALL cp_fm_to_fm(psi1_rxp(ispin, iB), psi_rxp, 1, istate, jstate)
2484 3996 : IF (current_env%full) CALL cp_fm_to_fm(psi1_D(ispin, iB), psi_D, 1, istate, jstate)
2485 : !
2486 : ! psi1_p_iiB_istate and psi1_p_iiiB_istate
2487 3996 : CALL cp_fm_to_fm(psi1_p(ispin, iiB), psi_p1, 1, istate, jstate)
2488 3996 : CALL cp_fm_to_fm(psi1_p(ispin, iiiB), psi_p2, 1, istate, jstate)
2489 : !
2490 7614 : jstate = jstate + 1
2491 : END DO ! istate
2492 : !
2493 : ! scale the ordered mos
2494 3618 : IF (current_env%full) CALL cp_fm_scale_and_add(1.0_dp, psi_rxp, -1.0_dp, psi_D)
2495 : !
2496 14472 : DO idir = 1, 3
2497 10854 : CALL set_vecp(idir, ii, iii)
2498 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix, psi_rxp, &
2499 10854 : p_rxp(idir), ncol=nstate_loc, alpha=1.e0_dp)
2500 10854 : IF (iiiB .EQ. iii .OR. iiiB .EQ. ii) THEN
2501 : CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, psi_p1, &
2502 7236 : r_p1(idir), ncol=nstate_loc, alpha=1.e0_dp)
2503 : END IF
2504 10854 : IF (iiB .EQ. iii .OR. iiB .EQ. ii) THEN
2505 : CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, psi_p2, &
2506 7236 : r_p2(idir), ncol=nstate_loc, alpha=1.e0_dp)
2507 : END IF
2508 123012 : DO idir2 = 1, 9
2509 97686 : IF (idir2 .EQ. ii .OR. idir2 .EQ. iii) THEN
2510 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_rxp, &
2511 21708 : rr_rxp(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2512 : END IF
2513 : !
2514 97686 : IF (idir2 .EQ. ind_m2(ii, iiiB) .OR. idir2 .EQ. ind_m2(iii, iiiB) .OR. idir2 .EQ. iiiB) THEN
2515 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_p1, &
2516 32562 : rr_p1(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2517 : END IF
2518 : !
2519 108540 : IF (idir2 .EQ. ind_m2(ii, iiB) .OR. idir2 .EQ. ind_m2(iii, iiB) .OR. idir2 .EQ. iiB) THEN
2520 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(idir2, idir)%matrix, psi_p2, &
2521 32562 : rr_p2(idir2, idir), ncol=nstate_loc, alpha=1.e0_dp)
2522 : END IF
2523 : END DO
2524 : END DO
2525 : !
2526 : ! Multuply left and right by the appropriate coefficients and sum into the
2527 : ! correct component of the chi tensor using the appropriate multiplicative factor
2528 : ! (don't forget the occupation number)
2529 : ! Loop over the cartesian components of the tensor
2530 : ! The loop over the components of the external field is external, thereby
2531 : ! only one column of the chi tensor is computed here
2532 15090 : DO idir = 1, 3
2533 10854 : chi_tmp = 0.0_dp
2534 10854 : int_current_tmp = 0.0_dp
2535 : !
2536 : ! get ii and iii
2537 10854 : CALL set_vecp(idir, ii, iii)
2538 : !
2539 : ! term: 2[C0| (r-dk)_ii |d_iii(C1(rxp-D))]-2[C0| (r-dk)_iii |d_ii(C1(rxp-D))]
2540 : ! the factor 2 should be already included in the matrix elements
2541 : contrib = 0.0_dp
2542 10854 : CALL cp_fm_trace(psi0, rr_rxp(ii, iii), contrib)
2543 10854 : chi_tmp = chi_tmp + 2.0_dp*contrib
2544 : !
2545 : contrib = 0.0_dp
2546 10854 : CALL cp_fm_trace(psi0, rr_rxp(iii, ii), contrib)
2547 10854 : chi_tmp = chi_tmp - 2.0_dp*contrib
2548 : !
2549 : ! correction: dk_ii*2[C0| d_iii(C1(rxp-D))] - dk_iii*2[C0| d_ii(C1(rxp-D))]
2550 : ! factor 2 not included in the matrix elements
2551 : contrib = 0.0_dp
2552 10854 : CALL cp_fm_trace(psi0, p_rxp(iii), contrib)
2553 10854 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(ii)*contrib
2554 10854 : int_current_tmp = int_current_tmp + 2.0_dp*contrib
2555 : !
2556 : contrib2 = 0.0_dp
2557 10854 : CALL cp_fm_trace(psi0, p_rxp(ii), contrib2)
2558 10854 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(iii)*contrib2
2559 : !
2560 : ! term: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] \
2561 : ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))]
2562 : ! the factor 2 should be already included in the matrix elements
2563 : contrib = 0.0_dp
2564 10854 : idir2 = ind_m2(ii, iiB)
2565 10854 : CALL cp_fm_trace(psi0, rr_p2(idir2, iii), contrib)
2566 10854 : chi_tmp = chi_tmp - 2.0_dp*contrib
2567 : contrib2 = 0.0_dp
2568 10854 : IF (iiB == iii) THEN
2569 3618 : CALL cp_fm_trace(psi0, r_p2(ii), contrib2)
2570 3618 : chi_tmp = chi_tmp - contrib2
2571 : END IF
2572 : !
2573 : contrib = 0.0_dp
2574 10854 : idir2 = ind_m2(iii, iiB)
2575 10854 : CALL cp_fm_trace(psi0, rr_p2(idir2, ii), contrib)
2576 10854 : chi_tmp = chi_tmp + 2.0_dp*contrib
2577 : contrib2 = 0.0_dp
2578 10854 : IF (iiB == ii) THEN
2579 3618 : CALL cp_fm_trace(psi0, r_p2(iii), contrib2)
2580 3618 : chi_tmp = chi_tmp + contrib2
2581 : END IF
2582 : !
2583 : ! correction: -dk_ii * 2[C0|(r-dk)_iiB | d_iii(C1(piiiB))] \
2584 : ! +dk_iii * 2[C0|(r-dk)_iiB | d_ii(C1(piiiB))]
2585 : ! the factor 2 should be already included in the matrix elements
2586 : ! no additional correction terms because of the orthogonality between C0 and C1
2587 : contrib = 0.0_dp
2588 10854 : CALL cp_fm_trace(psi0, rr_p2(iiB, iii), contrib)
2589 10854 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(ii)*contrib
2590 10854 : int_current_tmp = int_current_tmp - 2.0_dp*contrib
2591 : !
2592 : contrib2 = 0.0_dp
2593 10854 : CALL cp_fm_trace(psi0, rr_p2(iiB, ii), contrib2)
2594 10854 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(iii)*contrib2
2595 : !
2596 : ! term: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] \
2597 : ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))]
2598 : ! the factor 2 should be already included in the matrix elements
2599 : contrib = 0.0_dp
2600 10854 : idir2 = ind_m2(ii, iiiB)
2601 10854 : CALL cp_fm_trace(psi0, rr_p1(idir2, iii), contrib)
2602 10854 : chi_tmp = chi_tmp + 2.0_dp*contrib
2603 : contrib2 = 0.0_dp
2604 10854 : IF (iiiB == iii) THEN
2605 3618 : CALL cp_fm_trace(psi0, r_p1(ii), contrib2)
2606 3618 : chi_tmp = chi_tmp + contrib2
2607 : END IF
2608 : !
2609 : contrib = 0.0_dp
2610 10854 : idir2 = ind_m2(iii, iiiB)
2611 10854 : CALL cp_fm_trace(psi0, rr_p1(idir2, ii), contrib)
2612 10854 : chi_tmp = chi_tmp - 2.0_dp*contrib
2613 : contrib2 = 0.0_dp
2614 10854 : IF (iiiB == ii) THEN
2615 3618 : CALL cp_fm_trace(psi0, r_p1(iii), contrib2)
2616 3618 : chi_tmp = chi_tmp - contrib2
2617 : END IF
2618 : !
2619 : ! correction: +dk_ii * 2[C0|(r-dk)_iiiB | d_iii(C1(piiB))] +\
2620 : ! -dk_iii * 2[C0|(r-dk)_iiiB | d_ii(C1(piiB))]
2621 : ! the factor 2 should be already included in the matrix elements
2622 : contrib = 0.0_dp
2623 10854 : CALL cp_fm_trace(psi0, rr_p1(iiiB, iii), contrib)
2624 10854 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp + 2.0_dp*dk(ii)*contrib
2625 10854 : int_current_tmp = int_current_tmp + 2.0_dp*contrib
2626 : !
2627 : contrib2 = 0.0_dp
2628 10854 : CALL cp_fm_trace(psi0, rr_p1(iiiB, ii), contrib2)
2629 10854 : IF (.NOT. chi_pbc) chi_tmp = chi_tmp - 2.0_dp*dk(iii)*contrib2
2630 : !
2631 : ! accumulate
2632 10854 : chi(idir) = chi(idir) + maxocc*chi_tmp
2633 112158 : int_current(iii) = int_current(iii) + int_current_tmp
2634 : END DO ! idir
2635 :
2636 : END DO ! icenter
2637 : !
2638 3516 : DO idir = 1, 3
2639 : current_env%chi_tensor(idir, iB, ispin) = current_env%chi_tensor(idir, iB, ispin) + &
2640 1854 : chi(idir)
2641 618 : IF (output_unit > 0) THEN
2642 : !WRITE(output_unit,'(A,E12.6)') ' chi_'//ACHAR(119+idir)//ACHAR(119+iB)//&
2643 : ! & ' = ',chi(idir)
2644 : !WRITE(output_unit,'(A,E12.6)') ' analytic \int j_'//ACHAR(119+idir)//ACHAR(119+iB)//&
2645 : ! & '(r) d^3r = ',int_current(idir)
2646 : END IF
2647 : END DO
2648 : !
2649 : END DO ! ispin
2650 : !
2651 : ! deallocate the sparse matrices
2652 426 : CALL dbcsr_deallocate_matrix_set(op_mom_ao)
2653 426 : CALL dbcsr_deallocate_matrix_set(op_mom_der_ao)
2654 426 : CALL dbcsr_deallocate_matrix_set(op_p_ao)
2655 :
2656 426 : CALL cp_fm_release(psi0)
2657 426 : CALL cp_fm_release(psi_rxp)
2658 426 : CALL cp_fm_release(psi_D)
2659 426 : CALL cp_fm_release(psi_p1)
2660 426 : CALL cp_fm_release(psi_p2)
2661 1704 : DO idir = 1, 3
2662 1278 : CALL cp_fm_release(p_rxp(idir))
2663 1278 : CALL cp_fm_release(r_p1(idir))
2664 1278 : CALL cp_fm_release(r_p2(idir))
2665 13206 : DO idir2 = 1, 9
2666 11502 : CALL cp_fm_release(rr_rxp(idir2, idir))
2667 11502 : CALL cp_fm_release(rr_p1(idir2, idir))
2668 12780 : CALL cp_fm_release(rr_p2(idir2, idir))
2669 : END DO
2670 : END DO
2671 :
2672 426 : CALL timestop(handle)
2673 :
2674 2556 : END SUBROUTINE current_build_chi_many_centers
2675 :
2676 : ! **************************************************************************************************
2677 : !> \brief ...
2678 : !> \param current_env ...
2679 : !> \param qs_env ...
2680 : !> \param iB ...
2681 : ! **************************************************************************************************
2682 96 : SUBROUTINE current_build_chi_one_center(current_env, qs_env, iB)
2683 : !
2684 : TYPE(current_env_type) :: current_env
2685 : TYPE(qs_environment_type), POINTER :: qs_env
2686 : INTEGER, INTENT(IN) :: iB
2687 :
2688 : CHARACTER(LEN=*), PARAMETER :: routineN = 'current_build_chi_one_center'
2689 :
2690 : INTEGER :: handle, idir, idir2, iiB, iiiB, ispin, jdir, jjdir, kdir, max_states, nao, natom, &
2691 : nbr_center(2), nmo, nspins, nstates(2), output_unit
2692 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf
2693 96 : INTEGER, DIMENSION(:), POINTER :: row_blk_sizes
2694 : LOGICAL :: chi_pbc, gapw
2695 : REAL(dp) :: chi(3), contrib, dk(3), int_current(3), &
2696 : maxocc
2697 : TYPE(cell_type), POINTER :: cell
2698 96 : TYPE(cp_2d_i_p_type), DIMENSION(:), POINTER :: center_list
2699 96 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: centers_set
2700 : TYPE(cp_fm_type) :: buf
2701 96 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: psi0_order
2702 96 : TYPE(cp_fm_type), DIMENSION(:, :), POINTER :: psi1_p, psi1_rxp
2703 : TYPE(cp_fm_type), POINTER :: mo_coeff
2704 : TYPE(cp_logger_type), POINTER :: logger
2705 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
2706 96 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: op_mom_ao, op_p_ao
2707 96 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_mom_der_ao
2708 : TYPE(dft_control_type), POINTER :: dft_control
2709 96 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
2710 : TYPE(mp_para_env_type), POINTER :: para_env
2711 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
2712 96 : POINTER :: sab_all, sab_orb
2713 96 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2714 96 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2715 :
2716 : !
2717 :
2718 96 : CALL timeset(routineN, handle)
2719 : !
2720 96 : NULLIFY (dft_control, mos, para_env, mo_coeff, op_mom_ao, &
2721 96 : op_mom_der_ao, center_list, centers_set, &
2722 96 : op_p_ao, psi1_p, psi1_rxp, cell, psi0_order)
2723 :
2724 96 : logger => cp_get_default_logger()
2725 96 : output_unit = cp_logger_get_default_io_unit(logger)
2726 :
2727 : CALL get_qs_env(qs_env=qs_env, &
2728 : dft_control=dft_control, &
2729 : mos=mos, &
2730 : para_env=para_env, &
2731 : cell=cell, &
2732 : particle_set=particle_set, &
2733 : qs_kind_set=qs_kind_set, &
2734 : sab_all=sab_all, &
2735 : sab_orb=sab_orb, &
2736 96 : dbcsr_dist=dbcsr_dist)
2737 :
2738 96 : nspins = dft_control%nspins
2739 96 : gapw = dft_control%qs_control%gapw
2740 :
2741 : CALL get_current_env(current_env=current_env, &
2742 : chi_pbc=chi_pbc, &
2743 : nao=nao, &
2744 : nbr_center=nbr_center, &
2745 : center_list=center_list, &
2746 : centers_set=centers_set, &
2747 : psi1_p=psi1_p, &
2748 : psi1_rxp=psi1_rxp, &
2749 : nstates=nstates, &
2750 96 : psi0_order=psi0_order)
2751 : !
2752 228 : max_states = MAXVAL(nstates(1:nspins))
2753 : !
2754 : ! Allocate sparse matrices for dipole, quadrupole and their derivatives => 9x3
2755 : ! Remember the derivatives are antisymmetric
2756 96 : CALL dbcsr_allocate_matrix_set(op_mom_ao, 9)
2757 96 : CALL dbcsr_allocate_matrix_set(op_mom_der_ao, 9, 3)
2758 : !
2759 : ! prepare for allocation
2760 96 : natom = SIZE(particle_set, 1)
2761 288 : ALLOCATE (first_sgf(natom))
2762 192 : ALLOCATE (last_sgf(natom))
2763 : CALL get_particle_set(particle_set, qs_kind_set, &
2764 : first_sgf=first_sgf, &
2765 96 : last_sgf=last_sgf)
2766 192 : ALLOCATE (row_blk_sizes(natom))
2767 96 : CALL dbcsr_convert_offsets_to_sizes(first_sgf, row_blk_sizes, last_sgf)
2768 96 : DEALLOCATE (first_sgf)
2769 96 : DEALLOCATE (last_sgf)
2770 : !
2771 : !
2772 96 : ALLOCATE (op_mom_ao(1)%matrix)
2773 : CALL dbcsr_create(matrix=op_mom_ao(1)%matrix, &
2774 : name="op_mom", &
2775 : dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
2776 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2777 96 : mutable_work=.TRUE.)
2778 96 : CALL cp_dbcsr_alloc_block_from_nbl(op_mom_ao(1)%matrix, sab_all)
2779 :
2780 384 : DO idir2 = 1, 3
2781 288 : ALLOCATE (op_mom_der_ao(1, idir2)%matrix)
2782 : CALL dbcsr_copy(op_mom_der_ao(1, idir2)%matrix, op_mom_ao(1)%matrix, &
2783 384 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir2))))
2784 : END DO
2785 :
2786 864 : DO idir = 2, SIZE(op_mom_ao, 1)
2787 768 : ALLOCATE (op_mom_ao(idir)%matrix)
2788 : CALL dbcsr_copy(op_mom_ao(idir)%matrix, op_mom_ao(1)%matrix, &
2789 768 : "op_mom_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2790 3168 : DO idir2 = 1, 3
2791 2304 : ALLOCATE (op_mom_der_ao(idir, idir2)%matrix)
2792 : CALL dbcsr_copy(op_mom_der_ao(idir, idir2)%matrix, op_mom_ao(1)%matrix, &
2793 3072 : "op_mom_der_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir*idir2))))
2794 : END DO
2795 : END DO
2796 : !
2797 96 : CALL dbcsr_allocate_matrix_set(op_p_ao, 3)
2798 96 : ALLOCATE (op_p_ao(1)%matrix)
2799 : CALL dbcsr_create(matrix=op_p_ao(1)%matrix, &
2800 : name="op_p_ao", &
2801 : dist=dbcsr_dist, matrix_type=dbcsr_type_antisymmetric, &
2802 : row_blk_size=row_blk_sizes, col_blk_size=row_blk_sizes, &
2803 96 : mutable_work=.TRUE.)
2804 96 : CALL cp_dbcsr_alloc_block_from_nbl(op_p_ao(1)%matrix, sab_orb)
2805 :
2806 288 : DO idir = 2, 3
2807 192 : ALLOCATE (op_p_ao(idir)%matrix)
2808 : CALL dbcsr_copy(op_p_ao(idir)%matrix, op_p_ao(1)%matrix, &
2809 288 : "op_p_ao"//"-"//TRIM(ADJUSTL(cp_to_string(idir))))
2810 : END DO
2811 : !
2812 : !
2813 96 : DEALLOCATE (row_blk_sizes)
2814 : !
2815 : ! recompute the linear momentum matrices
2816 96 : CALL build_lin_mom_matrix(qs_env, op_p_ao)
2817 : !CALL p_xyz_ao(op_p_ao,qs_env,minimum_image=.FALSE.)
2818 : !
2819 : !
2820 : ! get iiB and iiiB
2821 96 : CALL set_vecp(iB, iiB, iiiB)
2822 228 : DO ispin = 1, nspins
2823 : !
2824 132 : CPASSERT(nbr_center(ispin) == 1)
2825 : !
2826 : ! get ground state MOS
2827 132 : nmo = nstates(ispin)
2828 132 : mo_coeff => psi0_order(ispin)
2829 132 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc)
2830 : !
2831 : ! Create buffer matrix
2832 132 : CALL cp_fm_create(buf, mo_coeff%matrix_struct)
2833 : !
2834 : ! Initialize the temporary vector chi
2835 132 : chi = 0.0_dp
2836 : int_current = 0.0_dp
2837 : !
2838 : !
2839 : ! Get the Wannier center of the istate-th ground state orbital
2840 528 : dk(1:3) = centers_set(ispin)%array(1:3, 1)
2841 : !
2842 : ! Compute the multipole integrals for the state istate,
2843 : ! using as reference center the corresponding Wannier center
2844 1320 : DO idir = 1, 9
2845 1188 : CALL dbcsr_set(op_mom_ao(idir)%matrix, 0.0_dp)
2846 4884 : DO idir2 = 1, 3
2847 4752 : CALL dbcsr_set(op_mom_der_ao(idir, idir2)%matrix, 0.0_dp)
2848 : END DO
2849 : END DO
2850 : CALL rRc_xyz_der_ao(op_mom_ao, op_mom_der_ao, qs_env, dk, order=2, &
2851 132 : minimum_image=.FALSE., soft=gapw)
2852 : !
2853 : !
2854 : ! Multuply left and right by the appropriate coefficients and sum into the
2855 : ! correct component of the chi tensor using the appropriate multiplicative factor
2856 : ! (don't forget the occupation number)
2857 : ! Loop over the cartesian components of the tensor
2858 : ! The loop over the components of the external field is external, thereby
2859 : ! only one column of the chi tensor is computed here
2860 528 : DO idir = 1, 3
2861 : !
2862 : !
2863 : !
2864 : ! term: dk_ii*2[C0| d_iii(C1(rxp-D))] - dk_iii*2[C0| d_ii(C1(rxp-D))]
2865 396 : IF (.NOT. chi_pbc) THEN
2866 : CALL cp_dbcsr_sm_fm_multiply(op_p_ao(idir)%matrix, mo_coeff, &
2867 396 : buf, ncol=nmo, alpha=1.e0_dp)
2868 1584 : DO jdir = 1, 3
2869 5148 : DO kdir = 1, 3
2870 3564 : IF (Levi_Civita(kdir, jdir, idir) .EQ. 0.0_dp) CYCLE
2871 792 : CALL cp_fm_trace(buf, psi1_rxp(ispin, iB), contrib)
2872 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jdir, idir)*2.0_dp*dk(jdir)*contrib
2873 : END DO
2874 : END DO
2875 : END IF
2876 : !
2877 : !
2878 : !
2879 : ! term: 2[C0| (r-dk)_ii |d_iii(C1(rxp-D))]-2[C0| (r-dk)_iii |d_ii(C1(rxp-D))]
2880 : ! and
2881 : ! term: -dk_ii * 2[C0|(r-dk)_iiB | d_iii(C1(piiiB))] +
2882 : ! +dk_iii * 2[C0|(r-dk)_iiB | d_ii(C1(piiiB))]
2883 : ! and
2884 : ! term: +dk_ii * 2[C0|(r-dk)_iiiB | d_iii(C1(piiB))] +
2885 : ! -dk_iii * 2[C0|(r-dk)_iiiB | d_ii(C1(piiB))]
2886 1584 : DO jdir = 1, 3
2887 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(jdir, idir)%matrix, mo_coeff, &
2888 1188 : buf, ncol=nmo, alpha=1.e0_dp)
2889 4752 : DO kdir = 1, 3
2890 3564 : IF (Levi_Civita(kdir, jdir, idir) .EQ. 0.0_dp) CYCLE
2891 792 : CALL cp_fm_trace(buf, psi1_rxp(ispin, iB), contrib)
2892 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jdir, idir)*2.0_dp*contrib
2893 : END DO
2894 : !
2895 1584 : IF (.NOT. chi_pbc) THEN
2896 1188 : IF (jdir .EQ. iiB) THEN
2897 1584 : DO jjdir = 1, 3
2898 5148 : DO kdir = 1, 3
2899 3564 : IF (Levi_Civita(kdir, jjdir, idir) .EQ. 0.0_dp) CYCLE
2900 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiiB), contrib)
2901 4752 : chi(kdir) = chi(kdir) + Levi_Civita(kdir, jjdir, idir)*2.0_dp*dk(jjdir)*contrib
2902 : END DO
2903 : END DO
2904 : END IF
2905 : !
2906 1188 : IF (jdir .EQ. iiiB) THEN
2907 1584 : DO jjdir = 1, 3
2908 5148 : DO kdir = 1, 3
2909 3564 : IF (Levi_Civita(kdir, jjdir, idir) .EQ. 0.0_dp) CYCLE
2910 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiB), contrib)
2911 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jjdir, idir)*2.0_dp*dk(jjdir)*contrib
2912 : END DO
2913 : END DO
2914 : END IF
2915 : END IF
2916 : END DO
2917 : !
2918 : !
2919 : !
2920 : ! term1: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] +
2921 : ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))]
2922 : ! and
2923 : ! term1: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] +
2924 : ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))]
2925 : ! HERE THERE IS ONE EXTRA MULTIPLY
2926 1584 : DO jdir = 1, 3
2927 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir, iiB), idir)%matrix, mo_coeff, &
2928 1188 : buf, ncol=nmo, alpha=1.e0_dp)
2929 4752 : DO kdir = 1, 3
2930 3564 : IF (Levi_Civita(kdir, jdir, idir) .EQ. 0.0_dp) CYCLE
2931 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiiB), contrib)
2932 4752 : chi(kdir) = chi(kdir) + Levi_Civita(kdir, jdir, idir)*2.0_dp*contrib
2933 : END DO
2934 : !
2935 : CALL cp_dbcsr_sm_fm_multiply(op_mom_der_ao(ind_m2(jdir, iiiB), idir)%matrix, mo_coeff, &
2936 1188 : buf, ncol=nmo, alpha=1.e0_dp)
2937 5148 : DO kdir = 1, 3
2938 3564 : IF (Levi_Civita(kdir, jdir, idir) .EQ. 0.0_dp) CYCLE
2939 792 : CALL cp_fm_trace(buf, psi1_p(ispin, iiB), contrib)
2940 4752 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, jdir, idir)*2.0_dp*contrib
2941 : END DO
2942 : END DO
2943 : !
2944 : !
2945 : !
2946 : ! term2: -2[C0| (r-dk)_ii (r-dk)_iiB | d_iii(C1(piiiB))] +
2947 : ! +2[C0| (r-dk)_iii (r-dk)_iiB | d_ii(C1(piiiB))]
2948 : ! and
2949 : ! term2: +2[C0| (r-dk)_ii (r-dk)_iiiB | d_iii(C1(piiB))] +
2950 : ! -2[C0| (r-dk)_iii (r-dk)_iiiB | d_ii(C1(piiB))]
2951 : CALL cp_dbcsr_sm_fm_multiply(op_mom_ao(idir)%matrix, mo_coeff, &
2952 396 : buf, ncol=nmo, alpha=1.e0_dp)
2953 1584 : DO jdir = 1, 3
2954 5148 : DO kdir = 1, 3
2955 3564 : IF (Levi_Civita(kdir, idir, jdir) .EQ. 0.0_dp) CYCLE
2956 1980 : IF (iiB == jdir) THEN
2957 264 : CALL cp_fm_trace(buf, psi1_p(ispin, iiiB), contrib)
2958 264 : chi(kdir) = chi(kdir) + Levi_Civita(kdir, idir, jdir)*contrib
2959 : END IF
2960 : END DO
2961 : END DO
2962 : !
2963 1716 : DO jdir = 1, 3
2964 5148 : DO kdir = 1, 3
2965 3564 : IF (Levi_Civita(kdir, idir, jdir) .EQ. 0.0_dp) CYCLE
2966 1980 : IF (iiiB == jdir) THEN
2967 264 : CALL cp_fm_trace(buf, psi1_p(ispin, iiB), contrib)
2968 264 : chi(kdir) = chi(kdir) - Levi_Civita(kdir, idir, jdir)*contrib
2969 : END IF
2970 : !
2971 : END DO
2972 : END DO
2973 : !
2974 : !
2975 : !
2976 : !
2977 : END DO ! idir
2978 : !
2979 528 : DO idir = 1, 3
2980 : current_env%chi_tensor(idir, iB, ispin) = current_env%chi_tensor(idir, iB, ispin) + &
2981 396 : maxocc*chi(idir)
2982 132 : IF (output_unit > 0) THEN
2983 : !WRITE(output_unit,'(A,E12.6)') ' chi_'//ACHAR(119+idir)//ACHAR(119+iB)//&
2984 : ! & ' = ',maxocc * chi(idir)
2985 : END IF
2986 : END DO
2987 : !
2988 360 : CALL cp_fm_release(buf)
2989 : END DO ! ispin
2990 : !
2991 : ! deallocate the sparse matrices
2992 96 : CALL dbcsr_deallocate_matrix_set(op_mom_ao)
2993 96 : CALL dbcsr_deallocate_matrix_set(op_mom_der_ao)
2994 96 : CALL dbcsr_deallocate_matrix_set(op_p_ao)
2995 :
2996 96 : CALL timestop(handle)
2997 :
2998 192 : END SUBROUTINE current_build_chi_one_center
2999 :
3000 0 : END MODULE qs_linres_current
|
