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