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