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