48 USE dbcsr_api,
ONLY: &
49 dbcsr_add, dbcsr_add_on_diag, dbcsr_clear, dbcsr_complete_redistribute, dbcsr_copy, &
50 dbcsr_create, dbcsr_distribution_new, dbcsr_distribution_release, dbcsr_distribution_type, &
51 dbcsr_frobenius_norm, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
52 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
53 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
54 dbcsr_type_antisymmetric, dbcsr_type_no_symmetry, dbcsr_type_symmetric
56 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_clear, dbt_contract, &
57 dbt_copy, dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, &
58 dbt_filter, dbt_get_info, dbt_mp_environ_pgrid, dbt_pgrid_create, dbt_pgrid_destroy, &
59 dbt_pgrid_type, dbt_scale, dbt_type
73 hfx_compression_type,&
115 USE pw_types,
ONLY: pw_c1d_gs_type,&
122 qs_environment_type,&
127 integrate_v_core_rspace,&
155 distribution_3d_type,&
156 neighbor_list_3c_type
165 #include "./base/base_uses.f90"
171 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rpa_im_time_force_methods'
189 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
190 TYPE(cp_fm_type),
INTENT(IN) :: fm_matrix_pq
191 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_m
192 INTEGER,
INTENT(IN) :: unit_nr
193 TYPE(mp2_type) :: mp2_env
194 TYPE(qs_environment_type),
POINTER :: qs_env
196 CHARACTER(LEN=*),
PARAMETER :: routinen =
'init_im_time_forces'
198 INTEGER :: handle, i_mem, i_xyz, ibasis, ispin, &
199 n_dependent, n_mem, n_rep, natom, &
201 INTEGER(int_8) :: nze, nze_tot
202 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist1, dist2, dist_ao_1, dist_ao_2, &
203 dist_ri, dummy_end, dummy_start, &
204 end_blocks, sizes_ao, sizes_ri, &
206 INTEGER,
DIMENSION(2) :: pdims_t2c
207 INTEGER,
DIMENSION(3) :: nblks_total, pcoord, pdims, pdims_t3c
208 INTEGER,
DIMENSION(:),
POINTER :: col_bsize, row_bsize
209 LOGICAL :: do_periodic, use_virial
210 REAL(
dp) :: compression_factor, eps_pgf_orb, &
211 eps_pgf_orb_old, memory, occ
212 TYPE(cell_type),
POINTER :: cell
213 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
214 TYPE(dbcsr_distribution_type) :: dbcsr_dist
215 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_s, rho_ao
216 TYPE(dbcsr_type) :: dbcsr_work, dbcsr_work2, dbcsr_work3
217 TYPE(dbcsr_type),
DIMENSION(1) :: t_2c_int_tmp
218 TYPE(dbcsr_type),
DIMENSION(1, 3) :: t_2c_der_tmp
219 TYPE(dbt_pgrid_type) :: pgrid_t2c, pgrid_t3c
220 TYPE(dbt_type) :: t_2c_template, t_2c_tmp, t_3c_template
221 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :, :) :: t_3c_der_ao_prv, t_3c_der_ri_prv
222 TYPE(dft_control_type),
POINTER :: dft_control
223 TYPE(distribution_2d_type),
POINTER :: dist_2d
224 TYPE(distribution_3d_type) :: dist_3d, dist_vir
225 TYPE(gto_basis_set_p_type),
ALLOCATABLE, &
226 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri_aux
227 TYPE(gto_basis_set_type),
POINTER :: orb_basis, ri_basis
228 TYPE(libint_potential_type) :: identity_pot
229 TYPE(mp_cart_type) :: mp_comm_t3c, mp_comm_vir
230 TYPE(mp_para_env_type),
POINTER :: para_env
231 TYPE(neighbor_list_3c_type) :: nl_3c
232 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
234 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
235 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
236 TYPE(qs_rho_type),
POINTER :: rho
237 TYPE(section_vals_type),
POINTER :: qs_section
238 TYPE(virial_type),
POINTER :: virial
240 NULLIFY (dft_control, para_env, particle_set, qs_kind_set, dist_2d, nl_2c, blacs_env, matrix_s, &
241 rho, rho_ao, cell, qs_section, orb_basis, ri_basis, virial)
245 CALL timeset(routinen, handle)
247 CALL get_qs_env(qs_env, natom=natom, nkind=nkind, dft_control=dft_control, para_env=para_env, &
248 particle_set=particle_set, qs_kind_set=qs_kind_set, cell=cell, virial=virial)
249 IF (dft_control%qs_control%gapw)
THEN
250 cpabort(
"Low-scaling RPA/SOS-MP2 forces only available with GPW")
253 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
255 do_periodic = .false.
256 IF (any(cell%perd == 1)) do_periodic = .true.
257 force_data%do_periodic = do_periodic
261 CALL dbt_pgrid_create(para_env, pdims_t3c, pgrid_t3c)
270 eps_pgf_orb = sqrt(eps_pgf_orb)
272 eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
274 ALLOCATE (sizes_ri(natom), sizes_ao(natom))
275 ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
277 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
279 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ao, basis=basis_set_ao)
281 DO ibasis = 1,
SIZE(basis_set_ao)
282 orb_basis => basis_set_ao(ibasis)%gto_basis_set
284 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
288 CALL create_3c_tensor(t_3c_template, dist_ri, dist_ao_1, dist_ao_2, pgrid_t3c, &
289 sizes_ri, sizes_ao, sizes_ao, map1=[1], map2=[2, 3], name=
"der (RI AO | AO)")
291 ALLOCATE (t_3c_der_ri_prv(1, 1, 3), t_3c_der_ao_prv(1, 1, 3))
293 CALL dbt_create(t_3c_template, t_3c_der_ri_prv(1, 1, i_xyz))
294 CALL dbt_create(t_3c_template, t_3c_der_ao_prv(1, 1, i_xyz))
298 ALLOCATE (force_data%t_3c_virial, force_data%t_3c_virial_split)
299 CALL dbt_create(t_3c_template, force_data%t_3c_virial)
300 CALL dbt_create(t_3c_m, force_data%t_3c_virial_split)
302 CALL dbt_destroy(t_3c_template)
304 CALL dbt_mp_environ_pgrid(pgrid_t3c, pdims, pcoord)
305 CALL mp_comm_t3c%create(pgrid_t3c%mp_comm_2d, 3, pdims)
307 nkind, particle_set, mp_comm_t3c, own_comm=.true.)
311 ALLOCATE (force_data%nl_3c)
312 CALL mp_comm_vir%create(pgrid_t3c%mp_comm_2d, 3, pdims)
314 nkind, particle_set, mp_comm_vir, own_comm=.true.)
316 dist_vir, mp2_env%ri_metric,
"RPA_3c_nl", qs_env, op_pos=1, &
317 sym_jk=.false., own_dist=.true.)
321 mp2_env%ri_metric,
"RPA_3c_nl", qs_env, op_pos=1, sym_jk=.true., &
323 DEALLOCATE (dist_ri, dist_ao_1, dist_ao_2)
326 CALL dbt_get_info(t_3c_m, nblks_total=nblks_total)
327 ALLOCATE (force_data%bsizes_RI_split(nblks_total(1)), force_data%bsizes_AO_split(nblks_total(2)))
328 CALL dbt_get_info(t_3c_m, blk_size_1=force_data%bsizes_RI_split, blk_size_2=force_data%bsizes_AO_split)
330 CALL dbt_create(t_3c_m, force_data%t_3c_der_RI(i_xyz))
331 CALL dbt_create(t_3c_m, force_data%t_3c_der_AO(i_xyz))
335 ALLOCATE (force_data%idx_to_at_RI(nblks_total(1)))
336 CALL get_idx_to_atom(force_data%idx_to_at_RI, force_data%bsizes_RI_split, sizes_ri)
338 ALLOCATE (force_data%idx_to_at_AO(nblks_total(2)))
339 CALL get_idx_to_atom(force_data%idx_to_at_AO, force_data%bsizes_AO_split, sizes_ao)
341 n_mem = mp2_env%ri_rpa_im_time%cut_memory
343 DEALLOCATE (dummy_start, dummy_end)
345 ALLOCATE (force_data%t_3c_der_AO_comp(n_mem, 3), force_data%t_3c_der_RI_comp(n_mem, 3))
346 ALLOCATE (force_data%t_3c_der_AO_ind(n_mem, 3), force_data%t_3c_der_RI_ind(n_mem, 3))
351 CALL build_3c_derivatives(t_3c_der_ri_prv, t_3c_der_ao_prv, mp2_env%ri_rpa_im_time%eps_filter, &
352 qs_env, nl_3c, basis_set_ri_aux, basis_set_ao, basis_set_ao, &
353 mp2_env%ri_metric, der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1, &
354 bounds_i=[start_blocks(i_mem), end_blocks(i_mem)])
357 CALL dbt_copy(t_3c_der_ri_prv(1, 1, i_xyz), force_data%t_3c_der_RI(i_xyz), move_data=.true.)
358 CALL dbt_filter(force_data%t_3c_der_RI(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
360 nze_tot = nze_tot + nze
363 CALL compress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
364 force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress, memory)
365 CALL dbt_clear(force_data%t_3c_der_RI(i_xyz))
367 CALL dbt_copy(t_3c_der_ao_prv(1, 1, i_xyz), force_data%t_3c_der_AO(i_xyz), move_data=.true.)
368 CALL dbt_filter(force_data%t_3c_der_AO(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
370 nze_tot = nze_tot + nze
373 CALL compress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
374 force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress, memory)
375 CALL dbt_clear(force_data%t_3c_der_AO(i_xyz))
380 CALL dbt_destroy(t_3c_der_ri_prv(1, 1, i_xyz))
381 CALL dbt_destroy(t_3c_der_ao_prv(1, 1, i_xyz))
384 CALL para_env%sum(memory)
385 compression_factor = real(nze_tot,
dp)*1.0e-06*8.0_dp/memory
386 IF (unit_nr > 0)
THEN
387 WRITE (unit=unit_nr, fmt=
"((T3,A,T66,F11.2,A4))") &
388 "MEMORY_INFO| Memory for 3-center derivatives (compressed):", memory,
' MiB'
390 WRITE (unit=unit_nr, fmt=
"((T3,A,T60,F21.2))") &
391 "MEMORY_INFO| Compression factor: ", compression_factor
395 CALL get_qs_env(qs_env, distribution_2d=dist_2d, blacs_env=blacs_env, matrix_s=matrix_s)
397 ALLOCATE (row_bsize(
SIZE(sizes_ri)))
398 ALLOCATE (col_bsize(
SIZE(sizes_ri)))
399 row_bsize(:) = sizes_ri(:)
400 col_bsize(:) = sizes_ri(:)
403 CALL dbt_pgrid_create(para_env, pdims_t2c, pgrid_t2c)
404 CALL create_2c_tensor(t_2c_template, dist1, dist2, pgrid_t2c, force_data%bsizes_RI_split, &
405 force_data%bsizes_RI_split, name=
'(RI| RI)')
406 DEALLOCATE (dist1, dist2)
408 CALL dbcsr_create(t_2c_int_tmp(1),
"(P|Q) RPA", dbcsr_dist, dbcsr_type_symmetric, row_bsize, col_bsize)
410 CALL dbcsr_create(t_2c_der_tmp(1, i_xyz),
"(P|Q) RPA der", dbcsr_dist, &
411 dbcsr_type_antisymmetric, row_bsize, col_bsize)
415 ALLOCATE (force_data%RI_virial_pot, force_data%RI_virial_met)
416 CALL dbcsr_create(force_data%RI_virial_pot,
"RI_virial", dbcsr_dist, &
417 dbcsr_type_no_symmetry, row_bsize, col_bsize)
418 CALL dbcsr_create(force_data%RI_virial_met,
"RI_virial", dbcsr_dist, &
419 dbcsr_type_no_symmetry, row_bsize, col_bsize)
424 CALL dbcsr_create(dbcsr_work, template=t_2c_int_tmp(1))
428 CALL dbcsr_create(dbcsr_work2, template=t_2c_int_tmp(1))
429 CALL dbcsr_create(dbcsr_work3, template=t_2c_int_tmp(1))
430 CALL dbcsr_copy(dbcsr_work2, dbcsr_work)
431 CALL cp_dbcsr_power(dbcsr_work, -0.5_dp, 1.0e-7_dp, n_dependent, para_env, blacs_env)
434 CALL dbt_create(dbcsr_work, t_2c_tmp)
435 CALL dbt_copy_matrix_to_tensor(dbcsr_work, t_2c_tmp)
436 CALL dbt_create(t_2c_template, force_data%t_2c_pot_msqrt)
437 CALL dbt_copy(t_2c_tmp, force_data%t_2c_pot_msqrt, move_data=.true.)
438 CALL dbt_filter(force_data%t_2c_pot_msqrt, mp2_env%ri_rpa_im_time%eps_filter)
440 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, dbcsr_work2, dbcsr_work, 0.0_dp, dbcsr_work3)
441 CALL dbt_copy_matrix_to_tensor(dbcsr_work3, t_2c_tmp)
442 CALL dbt_create(t_2c_template, force_data%t_2c_pot_psqrt)
443 CALL dbt_copy(t_2c_tmp, force_data%t_2c_pot_psqrt, move_data=.true.)
444 CALL dbt_filter(force_data%t_2c_pot_psqrt, mp2_env%ri_rpa_im_time%eps_filter)
445 CALL dbt_destroy(t_2c_tmp)
446 CALL dbcsr_release(dbcsr_work2)
447 CALL dbcsr_release(dbcsr_work3)
448 CALL dbcsr_clear(dbcsr_work)
451 IF (.NOT. do_periodic)
THEN
453 "RPA_2c_nl_pot", qs_env, sym_ij=.true., dist_2d=dist_2d)
455 basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
459 CALL dbt_create(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
460 CALL dbt_copy_matrix_to_tensor(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
461 CALL dbt_create(t_2c_template, force_data%t_2c_der_pot(i_xyz))
462 CALL dbt_copy(t_2c_tmp, force_data%t_2c_der_pot(i_xyz), move_data=.true.)
463 CALL dbt_filter(force_data%t_2c_der_pot(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
464 CALL dbt_destroy(t_2c_tmp)
465 CALL dbcsr_clear(t_2c_der_tmp(1, i_xyz))
470 mp2_env%potential_parameter,
"RPA_2c_nl_pot", qs_env, &
471 sym_ij=.false., dist_2d=dist_2d)
475 CALL dbcsr_create(force_data%G_PQ,
"G_PQ", dbcsr_dist, dbcsr_type_no_symmetry, row_bsize, col_bsize)
479 "RPA_2c_nl_metric", qs_env, sym_ij=.true., dist_2d=dist_2d)
480 CALL build_2c_integrals(t_2c_int_tmp, mp2_env%ri_rpa_im_time%eps_filter, qs_env, nl_2c, &
481 basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
483 basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
488 mp2_env%ri_metric,
"RPA_2c_nl_metric", qs_env, sym_ij=.false., &
492 CALL dbcsr_copy(dbcsr_work, t_2c_int_tmp(1))
496 CALL dbt_create(dbcsr_work, t_2c_tmp)
497 CALL dbt_copy_matrix_to_tensor(dbcsr_work, t_2c_tmp)
498 CALL dbt_create(t_2c_template, force_data%t_2c_inv_metric)
499 CALL dbt_copy(t_2c_tmp, force_data%t_2c_inv_metric, move_data=.true.)
500 CALL dbt_filter(force_data%t_2c_inv_metric, mp2_env%ri_rpa_im_time%eps_filter)
501 CALL dbt_destroy(t_2c_tmp)
502 CALL dbcsr_clear(dbcsr_work)
503 CALL dbcsr_clear(t_2c_int_tmp(1))
506 CALL dbt_create(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
507 CALL dbt_copy_matrix_to_tensor(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
508 CALL dbt_create(t_2c_template, force_data%t_2c_der_metric(i_xyz))
509 CALL dbt_copy(t_2c_tmp, force_data%t_2c_der_metric(i_xyz), move_data=.true.)
510 CALL dbt_filter(force_data%t_2c_der_metric(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
511 CALL dbt_destroy(t_2c_tmp)
512 CALL dbcsr_clear(t_2c_der_tmp(1, i_xyz))
516 CALL dbt_create(t_2c_template, force_data%t_2c_K)
517 CALL dbt_contract(1.0_dp, force_data%t_2c_inv_metric, force_data%t_2c_pot_psqrt, &
518 0.0_dp, force_data%t_2c_K, &
519 contract_1=[2], notcontract_1=[1], &
520 contract_2=[1], notcontract_2=[2], &
521 map_1=[1], map_2=[2], filter_eps=mp2_env%ri_rpa_im_time%eps_filter)
524 CALL dbt_destroy(t_2c_template)
525 CALL dbcsr_release(dbcsr_work)
526 CALL dbcsr_release(t_2c_int_tmp(1))
528 CALL dbcsr_release(t_2c_der_tmp(1, i_xyz))
531 DEALLOCATE (row_bsize, col_bsize)
532 ALLOCATE (row_bsize(
SIZE(sizes_ao)))
533 ALLOCATE (col_bsize(
SIZE(sizes_ao)))
534 row_bsize(:) = sizes_ao(:)
535 col_bsize(:) = sizes_ao(:)
537 CALL create_2c_tensor(t_2c_template, dist1, dist2, pgrid_t2c, force_data%bsizes_AO_split, &
538 force_data%bsizes_AO_split, name=
'(AO| AO)')
539 DEALLOCATE (dist1, dist2)
542 CALL dbcsr_create(t_2c_der_tmp(1, i_xyz),
"(P|Q) RPA der", dbcsr_dist, &
543 dbcsr_type_antisymmetric, row_bsize, col_bsize)
548 "RPA_2c_nl_metric", qs_env, sym_ij=.true., dist_2d=dist_2d)
550 basis_set_ao, basis_set_ao, identity_pot)
555 "RPA_2c_nl_metric", qs_env, sym_ij=.false., dist_2d=dist_2d)
558 CALL dbcsr_create(force_data%inv_ovlp, template=matrix_s(1)%matrix)
559 CALL dbcsr_copy(force_data%inv_ovlp, matrix_s(1)%matrix)
564 CALL dbt_create(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
565 CALL dbt_copy_matrix_to_tensor(t_2c_der_tmp(1, i_xyz), t_2c_tmp)
566 CALL dbt_create(t_2c_template, force_data%t_2c_der_ovlp(i_xyz))
567 CALL dbt_copy(t_2c_tmp, force_data%t_2c_der_ovlp(i_xyz), move_data=.true.)
568 CALL dbt_filter(force_data%t_2c_der_ovlp(i_xyz), mp2_env%ri_rpa_im_time%eps_filter)
569 CALL dbt_destroy(t_2c_tmp)
570 CALL dbcsr_clear(t_2c_der_tmp(1, i_xyz))
574 nspins = dft_control%nspins
575 ALLOCATE (force_data%P_virt(nspins), force_data%P_occ(nspins))
576 ALLOCATE (force_data%sum_YP_tau(nspins), force_data%sum_O_tau(nspins))
578 ALLOCATE (force_data%P_virt(ispin)%matrix, force_data%P_occ(ispin)%matrix)
579 ALLOCATE (force_data%sum_YP_tau(ispin)%matrix, force_data%sum_O_tau(ispin)%matrix)
580 CALL dbcsr_create(force_data%P_virt(ispin)%matrix, template=matrix_s(1)%matrix)
581 CALL dbcsr_create(force_data%P_occ(ispin)%matrix, template=matrix_s(1)%matrix)
582 CALL dbcsr_create(force_data%sum_O_tau(ispin)%matrix, template=matrix_s(1)%matrix)
583 CALL dbcsr_create(force_data%sum_YP_tau(ispin)%matrix, template=matrix_s(1)%matrix)
585 CALL dbcsr_copy(force_data%sum_O_tau(ispin)%matrix, matrix_s(1)%matrix)
586 CALL dbcsr_copy(force_data%sum_YP_tau(ispin)%matrix, matrix_s(1)%matrix)
588 CALL dbcsr_set(force_data%sum_O_tau(ispin)%matrix, 0.0_dp)
589 CALL dbcsr_set(force_data%sum_YP_tau(ispin)%matrix, 0.0_dp)
595 CALL dbcsr_copy(force_data%P_occ(1)%matrix, rho_ao(1)%matrix)
596 IF (nspins == 1)
THEN
597 CALL dbcsr_scale(force_data%P_occ(1)%matrix, 0.5_dp)
599 CALL dbcsr_copy(force_data%P_occ(2)%matrix, rho_ao(2)%matrix)
602 CALL dbcsr_copy(force_data%P_virt(ispin)%matrix, force_data%inv_ovlp)
603 CALL dbcsr_add(force_data%P_virt(ispin)%matrix, force_data%P_occ(ispin)%matrix, 1.0_dp, -1.0_dp)
606 DO ibasis = 1,
SIZE(basis_set_ao)
607 orb_basis => basis_set_ao(ibasis)%gto_basis_set
609 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
613 CALL dbt_destroy(t_2c_template)
614 CALL dbcsr_release(dbcsr_work)
616 CALL dbcsr_release(t_2c_der_tmp(1, i_xyz))
618 DEALLOCATE (row_bsize, col_bsize)
619 CALL dbt_pgrid_destroy(pgrid_t3c)
620 CALL dbt_pgrid_destroy(pgrid_t2c)
621 CALL dbcsr_distribution_release(dbcsr_dist)
622 CALL timestop(handle)
661 t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
662 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
663 fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, &
664 starts_array_mc_block, ends_array_mc_block, num_integ_points, &
665 nmo, Eigenval, tau_tj, tau_wj, cut_memory, Pspin, Qspin, &
666 open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
668 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
669 TYPE(dbcsr_p_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_p_omega
670 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_m, t_3c_o
671 TYPE(hfx_compression_type),
DIMENSION(:) :: t_3c_o_compressed
672 TYPE(block_ind_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_o_ind
673 TYPE(cp_fm_type),
INTENT(IN) :: fm_scaled_dm_occ_tau, &
674 fm_scaled_dm_virt_tau
675 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt
676 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff_occ_scaled, &
677 fm_mo_coeff_virt_scaled
678 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc, &
679 starts_array_mc_block, &
681 INTEGER,
INTENT(IN) :: num_integ_points, nmo
682 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: eigenval
683 REAL(kind=
dp),
DIMENSION(0:num_integ_points), &
685 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
687 INTEGER,
INTENT(IN) :: cut_memory, pspin, qspin
688 LOGICAL,
INTENT(IN) :: open_shell
689 INTEGER,
INTENT(IN) :: unit_nr
690 REAL(
dp),
INTENT(INOUT) :: dbcsr_time
691 INTEGER(int_8),
INTENT(INOUT) :: dbcsr_nflop
692 TYPE(mp2_type) :: mp2_env
693 TYPE(qs_environment_type),
POINTER :: qs_env
695 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_laplace_loop_forces'
697 INTEGER :: dummy_int, handle, handle2, i_mem, i_xyz, ibasis, ispin, j_xyz, jquad, k_xyz, &
698 n_mem_ri, n_rep, natom, nkind, nspins, unit_nr_dbcsr
699 INTEGER(int_8) :: flop, nze, nze_ddint, nze_der_ao, &
701 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, batch_blk_end, &
702 batch_blk_start, batch_end_ri, &
703 batch_start_ri, kind_of, mc_ranges, &
705 INTEGER,
DIMENSION(:, :),
POINTER :: dummy_ptr
706 LOGICAL :: memory_info, use_virial
707 REAL(
dp) :: eps_filter, eps_pgf_orb, &
708 eps_pgf_orb_old,
fac, occ, occ_ddint, &
709 occ_der_ao, occ_der_ri, occ_kqk, &
710 omega, pref, t1, t2, tau
711 REAL(
dp),
DIMENSION(3, 3) :: work_virial, work_virial_ovlp
712 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
713 TYPE(cell_type),
POINTER :: cell
714 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_s
715 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_occ, mat_dm_virt
716 TYPE(dbcsr_type) :: dbcsr_work1, dbcsr_work2, dbcsr_work3, &
717 exp_occ, exp_virt, r_occ, r_virt, &
718 virial_ovlp, y_1, y_2
719 TYPE(dbt_type) :: t_2c_ao, t_2c_ri, t_2c_ri_2, t_2c_tmp, t_3c_0, t_3c_1, t_3c_3, t_3c_4, &
720 t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_sparse, &
721 t_3c_work, t_dm_occ, t_dm_virt, t_kqkt, t_m_occ, t_m_virt, t_q, t_r_occ, t_r_virt
722 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_p
723 TYPE(dft_control_type),
POINTER :: dft_control
724 TYPE(gto_basis_set_p_type),
ALLOCATABLE, &
725 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri_aux
726 TYPE(gto_basis_set_type),
POINTER :: orb_basis, ri_basis
727 TYPE(libint_potential_type) :: identity_pot
728 TYPE(mp_para_env_type),
POINTER :: para_env
729 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
730 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
731 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
732 TYPE(section_vals_type),
POINTER :: qs_section
733 TYPE(virial_type),
POINTER :: virial
735 NULLIFY (matrix_s, dummy_ptr, atomic_kind_set, force, matrix_s, matrix_ks, mat_dm_occ, mat_dm_virt)
736 NULLIFY (dft_control, virial, particle_set, cell, para_env, orb_basis, ri_basis, qs_section)
737 NULLIFY (qs_kind_set)
739 CALL timeset(routinen, handle)
741 CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=natom, atomic_kind_set=atomic_kind_set, &
742 force=force, matrix_ks=matrix_ks, dft_control=dft_control, virial=virial, &
743 particle_set=particle_set, cell=cell, para_env=para_env, nkind=nkind, &
744 qs_kind_set=qs_kind_set)
745 eps_filter = mp2_env%ri_rpa_im_time%eps_filter
746 nspins = dft_control%nspins
748 memory_info = mp2_env%ri_rpa_im_time%memory_info
749 IF (memory_info)
THEN
750 unit_nr_dbcsr = unit_nr
755 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
757 IF (use_virial) virial%pv_calculate = .true.
766 eps_pgf_orb = sqrt(eps_pgf_orb)
768 eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
770 ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
774 DO ibasis = 1,
SIZE(basis_set_ao)
775 orb_basis => basis_set_ao(ibasis)%gto_basis_set
777 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
783 ALLOCATE (t_p(nspins))
784 CALL dbt_create(force_data%t_2c_K, t_2c_ri)
785 CALL dbt_create(force_data%t_2c_K, t_2c_ri_2)
786 CALL dbt_create(force_data%t_2c_der_ovlp(1), t_2c_ao)
792 ALLOCATE (mc_ranges(cut_memory + 1))
793 mc_ranges(:cut_memory) = starts_array_mc_block(:)
794 mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
797 n_mem_ri = cut_memory
799 batch_blk_start, batch_blk_end)
800 ALLOCATE (mc_ranges_ri(n_mem_ri + 1))
801 mc_ranges_ri(1:n_mem_ri) = batch_blk_start(1:n_mem_ri)
802 mc_ranges_ri(n_mem_ri + 1) = batch_blk_end(n_mem_ri) + 1
803 DEALLOCATE (batch_blk_start, batch_blk_end)
807 CALL dbt_create(t_2c_ri, t_p(ispin))
809 CALL dbt_create(t_2c_ri, t_q)
810 CALL dbt_create(t_2c_ri, t_kqkt)
811 CALL dbt_create(t_2c_ao, t_dm_occ)
812 CALL dbt_create(t_2c_ao, t_dm_virt)
815 CALL dbt_create(t_3c_o, t_m_occ)
816 CALL dbt_create(t_3c_o, t_m_virt)
817 CALL dbt_create(t_3c_o, t_3c_0)
819 CALL dbt_create(t_3c_o, t_3c_1)
820 CALL dbt_create(t_3c_o, t_3c_3)
821 CALL dbt_create(t_3c_o, t_3c_4)
822 CALL dbt_create(t_3c_o, t_3c_5)
823 CALL dbt_create(t_3c_m, t_3c_6)
824 CALL dbt_create(t_3c_m, t_3c_7)
825 CALL dbt_create(t_3c_m, t_3c_8)
826 CALL dbt_create(t_3c_m, t_3c_sparse)
827 CALL dbt_create(t_3c_o, t_3c_help_1)
828 CALL dbt_create(t_3c_o, t_3c_help_2)
829 CALL dbt_create(t_2c_ao, t_r_occ)
830 CALL dbt_create(t_2c_ao, t_r_virt)
831 CALL dbt_create(t_3c_m, t_3c_ints)
832 CALL dbt_create(t_3c_m, t_3c_work)
835 occ_der_ao = 0; nze_der_ao = 0
836 occ_der_ri = 0; nze_der_ri = 0
838 DO i_mem = 1, cut_memory
839 CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
840 force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
842 occ_der_ri = occ_der_ri + occ
843 nze_der_ri = nze_der_ri + nze
844 CALL dbt_copy(force_data%t_3c_der_RI(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
846 CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
847 force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
849 occ_der_ao = occ_der_ao + occ
850 nze_der_ao = nze_der_ao + nze
851 CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, order=[1, 3, 2], summation=.true.)
852 CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
855 occ_der_ri = occ_der_ri/3.0_dp
856 occ_der_ao = occ_der_ao/3.0_dp
857 nze_der_ri = nze_der_ri/3
858 nze_der_ao = nze_der_ao/3
860 CALL dbcsr_create(r_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
861 CALL dbcsr_create(r_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
862 CALL dbcsr_create(dbcsr_work1, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
863 CALL dbcsr_create(dbcsr_work2, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
864 CALL dbcsr_create(dbcsr_work3, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
865 CALL dbcsr_create(exp_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
866 CALL dbcsr_create(exp_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
867 IF (use_virial)
CALL dbcsr_create(virial_ovlp, template=dbcsr_work1)
869 CALL dbt_batched_contract_init(t_3c_0, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
870 CALL dbt_batched_contract_init(t_3c_1, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
871 CALL dbt_batched_contract_init(t_3c_3, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
872 CALL dbt_batched_contract_init(t_m_occ, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
873 CALL dbt_batched_contract_init(t_m_virt, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
875 CALL dbt_batched_contract_init(t_3c_ints, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
876 CALL dbt_batched_contract_init(t_3c_work, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
878 CALL dbt_batched_contract_init(t_3c_4, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
879 batch_range_3=mc_ranges)
880 CALL dbt_batched_contract_init(t_3c_5, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
881 batch_range_3=mc_ranges)
882 CALL dbt_batched_contract_init(t_3c_6, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
883 batch_range_3=mc_ranges)
884 CALL dbt_batched_contract_init(t_3c_7, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
885 batch_range_3=mc_ranges)
886 CALL dbt_batched_contract_init(t_3c_8, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
887 batch_range_3=mc_ranges)
888 CALL dbt_batched_contract_init(t_3c_sparse, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
889 batch_range_3=mc_ranges)
892 work_virial_ovlp = 0.0_dp
893 DO jquad = 1, num_integ_points
895 omega = tau_wj(jquad)
896 fac = -2.0_dp*omega*mp2_env%scale_S
897 IF (open_shell)
fac = 0.5_dp*
fac
898 occ_ddint = 0; nze_ddint = 0
906 CALL dbt_create(mat_p_omega(jquad, ispin)%matrix, t_2c_tmp)
907 CALL dbt_copy_matrix_to_tensor(mat_p_omega(jquad, ispin)%matrix, t_2c_tmp)
908 CALL dbt_copy(t_2c_tmp, t_p(ispin), move_data=.true.)
909 CALL dbt_filter(t_p(ispin), eps_filter)
910 CALL dbt_destroy(t_2c_tmp)
914 CALL dbt_contract(1.0_dp, t_p(qspin), force_data%t_2c_K, 0.0_dp, t_2c_ri, &
915 contract_1=[2], notcontract_1=[1], &
916 contract_2=[1], notcontract_2=[2], &
917 map_1=[1], map_2=[2], filter_eps=eps_filter, &
918 flop=flop, unit_nr=unit_nr_dbcsr)
919 dbcsr_nflop = dbcsr_nflop + flop
920 CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_2c_ri, 0.0_dp, t_q, &
921 contract_1=[1], notcontract_1=[2], &
922 contract_2=[1], notcontract_2=[2], &
923 map_1=[1], map_2=[2], filter_eps=eps_filter, &
924 flop=flop, unit_nr=unit_nr_dbcsr)
925 dbcsr_nflop = dbcsr_nflop + flop
926 CALL dbt_clear(t_2c_ri)
928 CALL perform_2c_ops(force, t_kqkt, force_data,
fac, t_q, t_p(pspin), t_2c_ri, t_2c_ri_2, &
929 use_virial, atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
934 nmo, fm_mo_coeff_occ(pspin), fm_mo_coeff_virt(pspin), &
935 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ, mat_dm_virt, &
936 matrix_s, pspin, eigenval(:, pspin), 0.0_dp, eps_filter, &
937 mp2_env%ri_rpa_im_time%memory_info, unit_nr, &
938 jquad, .false., .false., qs_env, dummy_int, dummy_ptr, para_env)
940 CALL dbt_create(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
941 CALL dbt_copy_matrix_to_tensor(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
942 CALL dbt_copy(t_2c_tmp, t_dm_occ, move_data=.true.)
943 CALL dbt_filter(t_dm_occ, eps_filter)
944 CALL dbt_destroy(t_2c_tmp)
946 CALL dbt_create(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
947 CALL dbt_copy_matrix_to_tensor(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
948 CALL dbt_copy(t_2c_tmp, t_dm_virt, move_data=.true.)
949 CALL dbt_filter(t_dm_virt, eps_filter)
950 CALL dbt_destroy(t_2c_tmp)
953 CALL perform_3c_ops(force, t_r_occ, t_r_virt, force_data,
fac, cut_memory, n_mem_ri, &
954 t_kqkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, t_m_virt, t_3c_0, t_3c_1, &
955 t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, t_3c_help_1, t_3c_help_2, &
956 t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_ri, &
957 batch_end_ri, t_3c_o_compressed, t_3c_o_ind, use_virial, &
958 atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
959 unit_nr_dbcsr, mp2_env)
961 CALL timeset(routinen//
"_dbcsr", handle2)
964 CALL dbt_create(matrix_s(1)%matrix, t_2c_tmp)
965 CALL dbt_copy(t_r_occ, t_2c_tmp, move_data=.true.)
966 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_occ)
968 CALL dbt_copy(t_r_virt, t_2c_tmp, move_data=.true.)
969 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_virt)
972 CALL dbcsr_multiply(
'N',
'N', tau, force_data%P_occ(pspin)%matrix, &
973 matrix_ks(pspin)%matrix, 0.0_dp, dbcsr_work1)
974 CALL build_y_matrix(y_1, dbcsr_work1, force_data%P_occ(pspin)%matrix, r_virt, eps_filter)
977 CALL dbcsr_multiply(
'N',
'N', -tau, force_data%P_virt(pspin)%matrix, &
978 matrix_ks(pspin)%matrix, 0.0_dp, dbcsr_work1)
979 CALL build_y_matrix(y_2, dbcsr_work1, force_data%P_virt(pspin)%matrix, r_occ, eps_filter)
984 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, r_occ, force_data%inv_ovlp, 0.0_dp, dbcsr_work1)
985 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, exp_virt, dbcsr_work1, 0.0_dp, dbcsr_work3)
986 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, force_data%inv_ovlp, dbcsr_work3, 0.0_dp, dbcsr_work2)
988 CALL dbcsr_multiply(
'N',
'T', tau, force_data%inv_ovlp, y_2, 0.0_dp, dbcsr_work3)
989 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, dbcsr_work3, matrix_ks(pspin)%matrix, 0.0_dp, dbcsr_work1)
990 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, dbcsr_work1, force_data%inv_ovlp, 0.0_dp, dbcsr_work3)
992 CALL dbcsr_add(dbcsr_work2, dbcsr_work3, 1.0_dp, -1.0_dp)
994 CALL dbt_copy_matrix_to_tensor(dbcsr_work2, t_2c_tmp)
995 CALL dbt_copy(t_2c_tmp, t_2c_ao, move_data=.true.)
998 CALL get_2c_der_force(force, t_2c_ao, force_data%t_2c_der_ovlp, atom_of_kind, &
999 kind_of, force_data%idx_to_at_AO, pref, do_ovlp=.true.)
1001 IF (use_virial)
CALL dbcsr_add(virial_ovlp, dbcsr_work2, 1.0_dp, pref)
1004 CALL dbcsr_multiply(
'N',
'N', tau*
fac, y_1, force_data%P_occ(pspin)%matrix, 1.0_dp, &
1005 force_data%sum_YP_tau(pspin)%matrix, retain_sparsity=.true.)
1006 CALL dbcsr_multiply(
'N',
'N', -tau*
fac, y_2, force_data%P_virt(pspin)%matrix, 1.0_dp, &
1007 force_data%sum_YP_tau(pspin)%matrix, retain_sparsity=.true.)
1010 pref = -omega*mp2_env%scale_S
1011 CALL dbcsr_multiply(
'N',
'N', pref, r_virt, exp_occ, 1.0_dp, &
1012 force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1013 CALL dbcsr_multiply(
'N',
'N', -pref, r_occ, exp_virt, 1.0_dp, &
1014 force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1015 CALL dbcsr_multiply(
'N',
'N', pref*tau, matrix_ks(pspin)%matrix, y_1, 1.0_dp, &
1016 force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1017 CALL dbcsr_multiply(
'N',
'N', pref*tau, matrix_ks(pspin)%matrix, y_2, 1.0_dp, &
1018 force_data%sum_O_tau(pspin)%matrix, retain_sparsity=.true.)
1020 CALL timestop(handle2)
1023 CALL para_env%sync()
1025 dbcsr_time = dbcsr_time + t2 - t1
1027 IF (unit_nr > 0)
THEN
1028 WRITE (unit_nr,
'(/T3,A,1X,I3,A)') &
1029 'RPA_LOW_SCALING_INFO| Info for time point', jquad,
' (gradients)'
1030 WRITE (unit_nr,
'(T6,A,T56,F25.6)') &
1031 'Execution time (s):', t2 - t1
1032 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1033 'Occupancy of 3c AO derivs:', real(nze_der_ao,
dp),
'/', occ_der_ao*100,
'%'
1034 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1035 'Occupancy of 3c RI derivs:', real(nze_der_ri,
dp),
'/', occ_der_ri*100,
'%'
1036 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1037 'Occupancy of the Docc * Dvirt * 3c-int tensor', real(nze_ddint,
dp),
'/', occ_ddint*100,
'%'
1038 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1039 'Occupancy of KQK^T 2c-tensor:', real(nze_kqk,
dp),
'/', occ_kqk*100,
'%'
1044 CALL dbcsr_release(y_1)
1045 CALL dbcsr_release(y_2)
1046 CALL dbt_destroy(t_2c_tmp)
1049 CALL dbt_batched_contract_finalize(t_3c_0)
1050 CALL dbt_batched_contract_finalize(t_3c_1)
1051 CALL dbt_batched_contract_finalize(t_3c_3)
1052 CALL dbt_batched_contract_finalize(t_m_occ)
1053 CALL dbt_batched_contract_finalize(t_m_virt)
1055 CALL dbt_batched_contract_finalize(t_3c_ints)
1056 CALL dbt_batched_contract_finalize(t_3c_work)
1058 CALL dbt_batched_contract_finalize(t_3c_4)
1059 CALL dbt_batched_contract_finalize(t_3c_5)
1060 CALL dbt_batched_contract_finalize(t_3c_6)
1061 CALL dbt_batched_contract_finalize(t_3c_7)
1062 CALL dbt_batched_contract_finalize(t_3c_8)
1063 CALL dbt_batched_contract_finalize(t_3c_sparse)
1066 IF (use_virial)
THEN
1067 CALL dbt_copy(force_data%t_3c_virial_split, force_data%t_3c_virial, move_data=.true.)
1068 CALL calc_3c_virial(work_virial, force_data%t_3c_virial, 1.0_dp, qs_env, force_data%nl_3c, &
1069 basis_set_ri_aux, basis_set_ao, basis_set_ao, mp2_env%ri_metric, &
1070 der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1)
1072 CALL calc_2c_virial(work_virial, force_data%RI_virial_met, 1.0_dp, qs_env, force_data%nl_2c_met, &
1073 basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
1074 CALL dbcsr_clear(force_data%RI_virial_met)
1076 IF (.NOT. force_data%do_periodic)
THEN
1077 CALL calc_2c_virial(work_virial, force_data%RI_virial_pot, 1.0_dp, qs_env, force_data%nl_2c_pot, &
1078 basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
1079 CALL dbcsr_clear(force_data%RI_virial_pot)
1083 CALL calc_2c_virial(work_virial_ovlp, virial_ovlp, 1.0_dp, qs_env, force_data%nl_2c_ovlp, &
1084 basis_set_ao, basis_set_ao, identity_pot)
1085 CALL dbcsr_release(virial_ovlp)
1090 virial%pv_mp2(i_xyz, j_xyz) = virial%pv_mp2(i_xyz, j_xyz) &
1091 - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1092 virial%pv_overlap(i_xyz, j_xyz) = virial%pv_overlap(i_xyz, j_xyz) &
1093 - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1094 virial%pv_virial(i_xyz, j_xyz) = virial%pv_virial(i_xyz, j_xyz) &
1095 - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz) &
1096 - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1103 work_virial = 0.0_dp
1104 IF (force_data%do_periodic)
THEN
1106 CALL get_2c_gpw_forces(force_data%G_PQ, force, work_virial, use_virial, mp2_env, qs_env)
1107 ELSE IF (mp2_env%eri_method ==
do_eri_mme)
THEN
1108 CALL get_2c_mme_forces(force_data%G_PQ, force, mp2_env, qs_env)
1109 IF (use_virial) cpabort(
"Stress tensor not available with MME intrgrals")
1111 cpabort(
"Periodic case not possible with OS integrals")
1113 CALL dbcsr_clear(force_data%G_PQ)
1116 IF (use_virial)
THEN
1117 virial%pv_mp2 = virial%pv_mp2 + work_virial
1118 virial%pv_virial = virial%pv_virial + work_virial
1119 virial%pv_calculate = .false.
1121 DO ibasis = 1,
SIZE(basis_set_ao)
1122 orb_basis => basis_set_ao(ibasis)%gto_basis_set
1124 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1130 IF (
ASSOCIATED(dummy_ptr))
DEALLOCATE (dummy_ptr)
1131 DO ispin = 1, nspins
1132 CALL dbt_destroy(t_p(ispin))
1134 CALL dbt_destroy(t_3c_0)
1135 CALL dbt_destroy(t_3c_1)
1136 CALL dbt_destroy(t_3c_3)
1137 CALL dbt_destroy(t_3c_4)
1138 CALL dbt_destroy(t_3c_5)
1139 CALL dbt_destroy(t_3c_6)
1140 CALL dbt_destroy(t_3c_7)
1141 CALL dbt_destroy(t_3c_8)
1142 CALL dbt_destroy(t_3c_sparse)
1143 CALL dbt_destroy(t_3c_help_1)
1144 CALL dbt_destroy(t_3c_help_2)
1145 CALL dbt_destroy(t_3c_ints)
1146 CALL dbt_destroy(t_3c_work)
1147 CALL dbt_destroy(t_r_occ)
1148 CALL dbt_destroy(t_r_virt)
1149 CALL dbt_destroy(t_dm_occ)
1150 CALL dbt_destroy(t_dm_virt)
1151 CALL dbt_destroy(t_q)
1152 CALL dbt_destroy(t_kqkt)
1153 CALL dbt_destroy(t_m_occ)
1154 CALL dbt_destroy(t_m_virt)
1155 CALL dbcsr_release(r_occ)
1156 CALL dbcsr_release(r_virt)
1157 CALL dbcsr_release(dbcsr_work1)
1158 CALL dbcsr_release(dbcsr_work2)
1159 CALL dbcsr_release(dbcsr_work3)
1160 CALL dbcsr_release(exp_occ)
1161 CALL dbcsr_release(exp_virt)
1163 CALL dbt_destroy(t_2c_ri)
1164 CALL dbt_destroy(t_2c_ri_2)
1165 CALL dbt_destroy(t_2c_ao)
1169 CALL timestop(handle)
1211 t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
1212 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
1213 fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, &
1214 starts_array_mc_block, ends_array_mc_block, num_integ_points, &
1215 nmo, Eigenval, e_fermi, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, &
1216 tj, wj, tau_tj, cut_memory, ispin, open_shell, unit_nr, dbcsr_time, &
1217 dbcsr_nflop, mp2_env, qs_env)
1219 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
1220 TYPE(dbcsr_p_type),
DIMENSION(:, :),
INTENT(INOUT) :: mat_p_omega
1221 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_m, t_3c_o
1222 TYPE(hfx_compression_type),
DIMENSION(:) :: t_3c_o_compressed
1223 TYPE(block_ind_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_o_ind
1224 TYPE(cp_fm_type),
INTENT(IN) :: fm_scaled_dm_occ_tau, &
1225 fm_scaled_dm_virt_tau
1226 TYPE(cp_fm_type),
DIMENSION(:),
INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt
1227 TYPE(cp_fm_type),
INTENT(IN) :: fm_mo_coeff_occ_scaled, &
1228 fm_mo_coeff_virt_scaled
1229 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc, &
1230 starts_array_mc_block, &
1232 INTEGER,
INTENT(IN) :: num_integ_points, nmo
1233 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(IN) :: eigenval
1234 REAL(kind=
dp),
INTENT(IN) :: e_fermi
1235 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
1236 INTENT(IN) :: weights_cos_tf_t_to_w, &
1237 weights_cos_tf_w_to_t
1238 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
1239 INTENT(IN) :: tj, wj
1240 REAL(kind=
dp),
DIMENSION(0:num_integ_points), &
1241 INTENT(IN) :: tau_tj
1242 INTEGER,
INTENT(IN) :: cut_memory, ispin
1243 LOGICAL,
INTENT(IN) :: open_shell
1244 INTEGER,
INTENT(IN) :: unit_nr
1245 REAL(
dp),
INTENT(INOUT) :: dbcsr_time
1246 INTEGER(int_8),
INTENT(INOUT) :: dbcsr_nflop
1247 TYPE(mp2_type) :: mp2_env
1248 TYPE(qs_environment_type),
POINTER :: qs_env
1250 CHARACTER(LEN=*),
PARAMETER :: routinen =
'calc_rpa_loop_forces'
1252 INTEGER :: dummy_int, handle, handle2, i_mem, i_xyz, ibasis, iquad, j_xyz, jquad, k_xyz, &
1253 n_mem_ri, n_rep, natom, nkind, nspins, unit_nr_dbcsr
1254 INTEGER(int_8) :: flop, nze, nze_ddint, nze_der_ao, &
1256 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, batch_blk_end, &
1257 batch_blk_start, batch_end_ri, &
1258 batch_start_ri, kind_of, mc_ranges, &
1260 INTEGER,
DIMENSION(:, :),
POINTER :: dummy_ptr
1261 LOGICAL :: memory_info, use_virial
1262 REAL(
dp) :: eps_filter, eps_pgf_orb, eps_pgf_orb_old,
fac, occ, occ_ddint, occ_der_ao, &
1263 occ_der_ri, occ_kbk, omega, pref, spin_fac, t1, t2, tau, weight
1264 REAL(
dp),
DIMENSION(3, 3) :: work_virial, work_virial_ovlp
1265 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
1266 TYPE(cell_type),
POINTER :: cell
1267 TYPE(cp_blacs_env_type),
POINTER :: blacs_env
1268 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: mat_p_tau, matrix_ks, matrix_s
1269 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_occ, mat_dm_virt
1270 TYPE(dbcsr_type) :: dbcsr_work1, dbcsr_work2, dbcsr_work3, &
1271 dbcsr_work_symm, exp_occ, exp_virt, &
1272 r_occ, r_virt, virial_ovlp, y_1, y_2
1273 TYPE(dbt_type) :: t_2c_ao, t_2c_ri, t_2c_ri_2, t_2c_tmp, t_3c_0, t_3c_1, t_3c_3, t_3c_4, &
1274 t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_sparse, &
1275 t_3c_work, t_dm_occ, t_dm_virt, t_kbkt, t_m_occ, t_m_virt, t_p, t_r_occ, t_r_virt
1276 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_b
1277 TYPE(dft_control_type),
POINTER :: dft_control
1278 TYPE(gto_basis_set_p_type),
ALLOCATABLE, &
1279 DIMENSION(:),
TARGET :: basis_set_ao, basis_set_ri_aux
1280 TYPE(gto_basis_set_type),
POINTER :: orb_basis, ri_basis
1281 TYPE(libint_potential_type) :: identity_pot
1282 TYPE(mp_para_env_type),
POINTER :: para_env
1283 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
1284 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1285 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
1286 TYPE(section_vals_type),
POINTER :: qs_section
1287 TYPE(virial_type),
POINTER :: virial
1289 NULLIFY (matrix_s, dummy_ptr, atomic_kind_set, force, matrix_s, matrix_ks, mat_dm_occ, mat_dm_virt)
1290 NULLIFY (dft_control, virial, particle_set, cell, blacs_env, para_env, orb_basis, ri_basis)
1291 NULLIFY (qs_kind_set)
1293 CALL timeset(routinen, handle)
1295 CALL get_qs_env(qs_env, matrix_s=matrix_s, natom=natom, atomic_kind_set=atomic_kind_set, &
1296 force=force, matrix_ks=matrix_ks, dft_control=dft_control, virial=virial, &
1297 particle_set=particle_set, cell=cell, blacs_env=blacs_env, para_env=para_env, &
1298 qs_kind_set=qs_kind_set, nkind=nkind)
1299 eps_filter = mp2_env%ri_rpa_im_time%eps_filter
1300 nspins = dft_control%nspins
1302 memory_info = mp2_env%ri_rpa_im_time%memory_info
1303 IF (memory_info)
THEN
1304 unit_nr_dbcsr = unit_nr
1309 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1311 IF (use_virial) virial%pv_calculate = .true.
1313 IF (use_virial)
THEN
1316 IF (n_rep /= 0)
THEN
1320 eps_pgf_orb = sqrt(eps_pgf_orb)
1322 eps_pgf_orb_old = dft_control%qs_control%eps_pgf_orb
1324 ALLOCATE (basis_set_ri_aux(nkind), basis_set_ao(nkind))
1328 DO ibasis = 1,
SIZE(basis_set_ao)
1329 orb_basis => basis_set_ao(ibasis)%gto_basis_set
1331 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1337 CALL dbt_create(force_data%t_2c_K, t_2c_ri)
1338 CALL dbt_create(force_data%t_2c_K, t_2c_ri_2)
1339 CALL dbt_create(force_data%t_2c_der_ovlp(1), t_2c_ao)
1345 ALLOCATE (mc_ranges(cut_memory + 1))
1346 mc_ranges(:cut_memory) = starts_array_mc_block(:)
1347 mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
1350 n_mem_ri = cut_memory
1352 batch_blk_start, batch_blk_end)
1353 ALLOCATE (mc_ranges_ri(n_mem_ri + 1))
1354 mc_ranges_ri(1:n_mem_ri) = batch_blk_start(1:n_mem_ri)
1355 mc_ranges_ri(n_mem_ri + 1) = batch_blk_end(n_mem_ri) + 1
1356 DEALLOCATE (batch_blk_start, batch_blk_end)
1359 CALL dbt_create(t_2c_ri, t_p)
1360 CALL dbt_create(t_2c_ri, t_kbkt)
1361 CALL dbt_create(t_2c_ao, t_dm_occ)
1362 CALL dbt_create(t_2c_ao, t_dm_virt)
1365 CALL dbt_create(t_3c_o, t_m_occ)
1366 CALL dbt_create(t_3c_o, t_m_virt)
1367 CALL dbt_create(t_3c_o, t_3c_0)
1369 CALL dbt_create(t_3c_o, t_3c_1)
1370 CALL dbt_create(t_3c_o, t_3c_3)
1371 CALL dbt_create(t_3c_o, t_3c_4)
1372 CALL dbt_create(t_3c_o, t_3c_5)
1373 CALL dbt_create(t_3c_m, t_3c_6)
1374 CALL dbt_create(t_3c_m, t_3c_7)
1375 CALL dbt_create(t_3c_m, t_3c_8)
1376 CALL dbt_create(t_3c_m, t_3c_sparse)
1377 CALL dbt_create(t_3c_o, t_3c_help_1)
1378 CALL dbt_create(t_3c_o, t_3c_help_2)
1379 CALL dbt_create(t_2c_ao, t_r_occ)
1380 CALL dbt_create(t_2c_ao, t_r_virt)
1381 CALL dbt_create(t_3c_m, t_3c_ints)
1382 CALL dbt_create(t_3c_m, t_3c_work)
1386 ALLOCATE (t_b(num_integ_points))
1387 DO jquad = 1, num_integ_points
1388 CALL dbt_create(t_2c_ri, t_b(jquad))
1391 ALLOCATE (mat_p_tau(num_integ_points))
1392 DO jquad = 1, num_integ_points
1393 ALLOCATE (mat_p_tau(jquad)%matrix)
1394 CALL dbcsr_create(mat_p_tau(jquad)%matrix, template=mat_p_omega(jquad, ispin)%matrix)
1397 CALL dbcsr_create(dbcsr_work_symm, template=force_data%G_PQ, matrix_type=dbcsr_type_symmetric)
1398 CALL dbt_create(dbcsr_work_symm, t_2c_tmp)
1401 DO iquad = 1, num_integ_points
1406 CALL dbcsr_copy(dbcsr_work_symm, mat_p_omega(iquad, 1)%matrix)
1407 CALL dbt_copy_matrix_to_tensor(dbcsr_work_symm, t_2c_tmp)
1408 CALL dbt_copy(t_2c_tmp, t_2c_ri, move_data=.true.)
1410 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_2c_ri_2, &
1411 contract_1=[2], notcontract_1=[1], &
1412 contract_2=[1], notcontract_2=[2], &
1413 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1414 flop=flop, unit_nr=unit_nr_dbcsr)
1415 dbcsr_nflop = dbcsr_nflop + flop
1416 CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_2c_ri_2, 0.0_dp, t_2c_ri, &
1417 contract_1=[1], notcontract_1=[2], &
1418 contract_2=[1], notcontract_2=[2], &
1419 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1420 flop=flop, unit_nr=unit_nr_dbcsr)
1421 CALL dbt_copy(t_2c_ri, t_2c_tmp, move_data=.true.)
1422 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, dbcsr_work_symm)
1423 CALL dbcsr_add_on_diag(dbcsr_work_symm, 1.0_dp)
1428 CALL dbcsr_add_on_diag(dbcsr_work_symm, -1.0_dp)
1430 DO jquad = 1, num_integ_points
1434 weight = weights_cos_tf_w_to_t(jquad, iquad)*cos(tau*omega)
1435 IF (open_shell)
THEN
1436 IF (ispin == 1)
THEN
1438 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 1)%matrix, 1.0_dp, weight)
1439 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 2)%matrix, 1.0_dp, -weight)
1441 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 2)%matrix, 1.0_dp, weight)
1445 weight = 0.5_dp*weight
1446 CALL dbcsr_add(mat_p_tau(jquad)%matrix, mat_p_omega(iquad, 1)%matrix, 1.0_dp, weight)
1450 weight = weights_cos_tf_t_to_w(iquad, jquad)*cos(tau*omega)*wj(iquad)
1451 CALL dbt_copy_matrix_to_tensor(dbcsr_work_symm, t_2c_tmp)
1452 CALL dbt_scale(t_2c_tmp, weight)
1453 CALL dbt_copy(t_2c_tmp, t_b(jquad), summation=.true., move_data=.true.)
1456 CALL dbt_destroy(t_2c_tmp)
1457 CALL dbcsr_release(dbcsr_work_symm)
1458 CALL dbt_clear(t_2c_ri)
1459 CALL dbt_clear(t_2c_ri_2)
1462 occ_der_ao = 0; nze_der_ao = 0
1463 occ_der_ri = 0; nze_der_ri = 0
1465 DO i_mem = 1, cut_memory
1466 CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(i_mem, i_xyz)%ind, &
1467 force_data%t_3c_der_RI_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
1469 occ_der_ri = occ_der_ri + occ
1470 nze_der_ri = nze_der_ri + nze
1471 CALL dbt_copy(force_data%t_3c_der_RI(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
1473 CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(i_mem, i_xyz)%ind, &
1474 force_data%t_3c_der_AO_comp(i_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
1476 occ_der_ao = occ_der_ao + occ
1477 nze_der_ao = nze_der_ao + nze
1478 CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, order=[1, 3, 2], summation=.true.)
1479 CALL dbt_copy(force_data%t_3c_der_AO(i_xyz), t_3c_sparse, summation=.true., move_data=.true.)
1482 occ_der_ri = occ_der_ri/3.0_dp
1483 occ_der_ao = occ_der_ao/3.0_dp
1484 nze_der_ri = nze_der_ri/3
1485 nze_der_ao = nze_der_ao/3
1487 CALL dbcsr_create(r_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1488 CALL dbcsr_create(r_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1489 CALL dbcsr_create(dbcsr_work_symm, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
1490 CALL dbcsr_create(dbcsr_work1, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1491 CALL dbcsr_create(dbcsr_work2, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1492 CALL dbcsr_create(dbcsr_work3, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1493 CALL dbcsr_create(exp_occ, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1494 CALL dbcsr_create(exp_virt, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1495 IF (use_virial)
CALL dbcsr_create(virial_ovlp, template=dbcsr_work1)
1497 CALL dbt_batched_contract_init(t_3c_0, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
1498 CALL dbt_batched_contract_init(t_3c_1, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
1499 CALL dbt_batched_contract_init(t_3c_3, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1500 CALL dbt_batched_contract_init(t_m_occ, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1501 CALL dbt_batched_contract_init(t_m_virt, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1503 CALL dbt_batched_contract_init(t_3c_ints, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1504 CALL dbt_batched_contract_init(t_3c_work, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges)
1506 CALL dbt_batched_contract_init(t_3c_4, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1507 batch_range_3=mc_ranges)
1508 CALL dbt_batched_contract_init(t_3c_5, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1509 batch_range_3=mc_ranges)
1510 CALL dbt_batched_contract_init(t_3c_6, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1511 batch_range_3=mc_ranges)
1512 CALL dbt_batched_contract_init(t_3c_7, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1513 batch_range_3=mc_ranges)
1514 CALL dbt_batched_contract_init(t_3c_8, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1515 batch_range_3=mc_ranges)
1516 CALL dbt_batched_contract_init(t_3c_sparse, batch_range_1=mc_ranges_ri, batch_range_2=mc_ranges, &
1517 batch_range_3=mc_ranges)
1519 fac = 1.0_dp/
fourpi*mp2_env%ri_rpa%scale_rpa
1520 IF (open_shell)
fac = 0.5_dp*
fac
1522 work_virial = 0.0_dp
1523 work_virial_ovlp = 0.0_dp
1524 DO jquad = 1, num_integ_points
1526 occ_ddint = 0; nze_ddint = 0
1528 CALL para_env%sync()
1533 CALL dbt_create(mat_p_tau(jquad)%matrix, t_2c_tmp)
1534 CALL dbt_copy_matrix_to_tensor(mat_p_tau(jquad)%matrix, t_2c_tmp)
1535 CALL dbt_copy(t_2c_tmp, t_p, move_data=.true.)
1536 CALL dbt_filter(t_p, eps_filter)
1537 CALL dbt_destroy(t_2c_tmp)
1539 CALL perform_2c_ops(force, t_kbkt, force_data,
fac, t_b(jquad), t_p, t_2c_ri, t_2c_ri_2, &
1540 use_virial, atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
1545 nmo, fm_mo_coeff_occ(ispin), fm_mo_coeff_virt(ispin), &
1546 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ, mat_dm_virt, &
1547 matrix_s, ispin, eigenval(:, ispin), e_fermi, eps_filter, &
1548 mp2_env%ri_rpa_im_time%memory_info, unit_nr, &
1549 jquad, .false., .false., qs_env, dummy_int, dummy_ptr, para_env)
1551 CALL dbt_create(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
1552 CALL dbt_copy_matrix_to_tensor(mat_dm_occ(jquad, 1)%matrix, t_2c_tmp)
1553 CALL dbt_copy(t_2c_tmp, t_dm_occ, move_data=.true.)
1554 CALL dbt_filter(t_dm_occ, eps_filter)
1555 CALL dbt_destroy(t_2c_tmp)
1557 CALL dbt_create(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
1558 CALL dbt_copy_matrix_to_tensor(mat_dm_virt(jquad, 1)%matrix, t_2c_tmp)
1559 CALL dbt_copy(t_2c_tmp, t_dm_virt, move_data=.true.)
1560 CALL dbt_filter(t_dm_virt, eps_filter)
1561 CALL dbt_destroy(t_2c_tmp)
1564 CALL perform_3c_ops(force, t_r_occ, t_r_virt, force_data,
fac, cut_memory, n_mem_ri, &
1565 t_kbkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, t_m_virt, t_3c_0, t_3c_1, &
1566 t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, t_3c_help_1, t_3c_help_2, &
1567 t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_ri, &
1568 batch_end_ri, t_3c_o_compressed, t_3c_o_ind, use_virial, &
1569 atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
1570 unit_nr_dbcsr, mp2_env)
1572 CALL timeset(routinen//
"_dbcsr", handle2)
1575 CALL dbt_create(matrix_s(1)%matrix, t_2c_tmp)
1576 CALL dbt_copy(t_r_occ, t_2c_tmp, move_data=.true.)
1577 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_occ)
1579 CALL dbt_copy(t_r_virt, t_2c_tmp, move_data=.true.)
1580 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, r_virt)
1583 CALL dbcsr_copy(dbcsr_work_symm, matrix_ks(ispin)%matrix)
1584 CALL dbcsr_add(dbcsr_work_symm, matrix_s(1)%matrix, 1.0_dp, -e_fermi)
1585 CALL dbcsr_multiply(
'N',
'N', tau, force_data%P_occ(ispin)%matrix, &
1586 dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1587 CALL build_y_matrix(y_1, dbcsr_work1, force_data%P_occ(ispin)%matrix, r_virt, eps_filter)
1590 CALL dbcsr_multiply(
'N',
'N', -tau, force_data%P_virt(ispin)%matrix, &
1591 dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1592 CALL build_y_matrix(y_2, dbcsr_work1, force_data%P_virt(ispin)%matrix, r_occ, eps_filter)
1598 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, r_occ, force_data%inv_ovlp, 0.0_dp, dbcsr_work1)
1599 CALL dbcsr_multiply(
'T',
'N', 1.0_dp, exp_virt, dbcsr_work1, 0.0_dp, dbcsr_work3)
1600 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, force_data%inv_ovlp, dbcsr_work3, 0.0_dp, dbcsr_work2)
1602 CALL dbcsr_multiply(
'N',
'T', tau, force_data%inv_ovlp, y_2, 0.0_dp, dbcsr_work3)
1603 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, dbcsr_work3, dbcsr_work_symm, 0.0_dp, dbcsr_work1)
1604 CALL dbcsr_multiply(
'N',
'N', -1.0_dp, dbcsr_work1, force_data%inv_ovlp, 1.0_dp, dbcsr_work2)
1606 CALL dbcsr_multiply(
'N',
'T', tau*e_fermi, force_data%P_occ(ispin)%matrix, y_1, 1.0_dp, dbcsr_work2)
1607 CALL dbcsr_multiply(
'N',
'T', -tau*e_fermi, force_data%P_virt(ispin)%matrix, y_2, 1.0_dp, dbcsr_work2)
1609 CALL dbt_copy_matrix_to_tensor(dbcsr_work2, t_2c_tmp)
1610 CALL dbt_copy(t_2c_tmp, t_2c_ao, move_data=.true.)
1613 CALL get_2c_der_force(force, t_2c_ao, force_data%t_2c_der_ovlp, atom_of_kind, &
1614 kind_of, force_data%idx_to_at_AO, pref, do_ovlp=.true.)
1616 IF (use_virial)
CALL dbcsr_add(virial_ovlp, dbcsr_work2, 1.0_dp, pref)
1619 CALL dbcsr_multiply(
'N',
'N',
fac*tau, y_1, force_data%P_occ(ispin)%matrix, 1.0_dp, &
1620 force_data%sum_YP_tau(ispin)%matrix, retain_sparsity=.true.)
1621 CALL dbcsr_multiply(
'N',
'N', -
fac*tau, y_2, force_data%P_virt(ispin)%matrix, 1.0_dp, &
1622 force_data%sum_YP_tau(ispin)%matrix, retain_sparsity=.true.)
1624 spin_fac = 0.5_dp*
fac
1625 IF (open_shell) spin_fac = 2.0_dp*spin_fac
1627 CALL dbcsr_multiply(
'N',
'N', 1.0_dp*spin_fac, r_virt, exp_occ, 1.0_dp, &
1628 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1629 CALL dbcsr_multiply(
'N',
'N', -1.0_dp*spin_fac, r_occ, exp_virt, 1.0_dp, &
1630 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1631 CALL dbcsr_multiply(
'N',
'N', tau*spin_fac, dbcsr_work_symm, y_1, 1.0_dp, &
1632 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1633 CALL dbcsr_multiply(
'N',
'N', tau*spin_fac, dbcsr_work_symm, y_2, 1.0_dp, &
1634 force_data%sum_O_tau(ispin)%matrix, retain_sparsity=.true.)
1636 CALL timestop(handle2)
1639 CALL para_env%sync()
1641 dbcsr_time = dbcsr_time + t2 - t1
1643 IF (unit_nr > 0)
THEN
1644 WRITE (unit_nr,
'(/T3,A,1X,I3,A)') &
1645 'RPA_LOW_SCALING_INFO| Info for time point', jquad,
' (gradients)'
1646 WRITE (unit_nr,
'(T6,A,T56,F25.6)') &
1648 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1649 'Occupancy of 3c AO derivs:', real(nze_der_ao,
dp),
'/', occ_der_ao*100,
'%'
1650 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1651 'Occupancy of 3c RI derivs:', real(nze_der_ri,
dp),
'/', occ_der_ri*100,
'%'
1652 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1653 'Occupancy of the Docc * Dvirt * 3c-int tensor', real(nze_ddint,
dp),
'/', occ_ddint*100,
'%'
1654 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
1655 'Occupancy of KBK^T 2c-tensor:', real(nze_kbk,
dp),
'/', occ_kbk*100,
'%'
1660 CALL dbcsr_release(y_1)
1661 CALL dbcsr_release(y_2)
1662 CALL dbt_destroy(t_2c_tmp)
1666 CALL dbt_batched_contract_finalize(t_3c_0)
1667 CALL dbt_batched_contract_finalize(t_3c_1)
1668 CALL dbt_batched_contract_finalize(t_3c_3)
1669 CALL dbt_batched_contract_finalize(t_m_occ)
1670 CALL dbt_batched_contract_finalize(t_m_virt)
1672 CALL dbt_batched_contract_finalize(t_3c_ints)
1673 CALL dbt_batched_contract_finalize(t_3c_work)
1675 CALL dbt_batched_contract_finalize(t_3c_4)
1676 CALL dbt_batched_contract_finalize(t_3c_5)
1677 CALL dbt_batched_contract_finalize(t_3c_6)
1678 CALL dbt_batched_contract_finalize(t_3c_7)
1679 CALL dbt_batched_contract_finalize(t_3c_8)
1680 CALL dbt_batched_contract_finalize(t_3c_sparse)
1683 IF (use_virial)
THEN
1684 CALL dbt_copy(force_data%t_3c_virial_split, force_data%t_3c_virial, move_data=.true.)
1685 CALL calc_3c_virial(work_virial, force_data%t_3c_virial, 1.0_dp, qs_env, force_data%nl_3c, &
1686 basis_set_ri_aux, basis_set_ao, basis_set_ao, mp2_env%ri_metric, &
1687 der_eps=mp2_env%ri_rpa_im_time%eps_filter, op_pos=1)
1689 CALL calc_2c_virial(work_virial, force_data%RI_virial_met, 1.0_dp, qs_env, force_data%nl_2c_met, &
1690 basis_set_ri_aux, basis_set_ri_aux, mp2_env%ri_metric)
1691 CALL dbcsr_clear(force_data%RI_virial_met)
1693 IF (.NOT. force_data%do_periodic)
THEN
1694 CALL calc_2c_virial(work_virial, force_data%RI_virial_pot, 1.0_dp, qs_env, force_data%nl_2c_pot, &
1695 basis_set_ri_aux, basis_set_ri_aux, mp2_env%potential_parameter)
1696 CALL dbcsr_clear(force_data%RI_virial_pot)
1700 CALL calc_2c_virial(work_virial_ovlp, virial_ovlp, 1.0_dp, qs_env, force_data%nl_2c_ovlp, &
1701 basis_set_ao, basis_set_ao, identity_pot)
1702 CALL dbcsr_release(virial_ovlp)
1707 virial%pv_mp2(i_xyz, j_xyz) = virial%pv_mp2(i_xyz, j_xyz) &
1708 - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1709 virial%pv_overlap(i_xyz, j_xyz) = virial%pv_overlap(i_xyz, j_xyz) &
1710 - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1711 virial%pv_virial(i_xyz, j_xyz) = virial%pv_virial(i_xyz, j_xyz) &
1712 - work_virial(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz) &
1713 - work_virial_ovlp(i_xyz, k_xyz)*cell%hmat(j_xyz, k_xyz)
1720 work_virial = 0.0_dp
1721 IF (force_data%do_periodic)
THEN
1723 CALL get_2c_gpw_forces(force_data%G_PQ, force, work_virial, use_virial, mp2_env, qs_env)
1724 ELSE IF (mp2_env%eri_method ==
do_eri_mme)
THEN
1725 CALL get_2c_mme_forces(force_data%G_PQ, force, mp2_env, qs_env)
1726 IF (use_virial) cpabort(
"Stress tensor not available with MME intrgrals")
1728 cpabort(
"Periodic case not possible with OS integrals")
1730 CALL dbcsr_clear(force_data%G_PQ)
1733 IF (use_virial)
THEN
1734 virial%pv_mp2 = virial%pv_mp2 + work_virial
1735 virial%pv_virial = virial%pv_virial + work_virial
1736 virial%pv_calculate = .false.
1738 DO ibasis = 1,
SIZE(basis_set_ao)
1739 orb_basis => basis_set_ao(ibasis)%gto_basis_set
1741 ri_basis => basis_set_ri_aux(ibasis)%gto_basis_set
1747 IF (
ASSOCIATED(dummy_ptr))
DEALLOCATE (dummy_ptr)
1748 DO jquad = 1, num_integ_points
1749 CALL dbt_destroy(t_b(jquad))
1751 CALL dbt_destroy(t_p)
1752 CALL dbt_destroy(t_3c_0)
1753 CALL dbt_destroy(t_3c_1)
1754 CALL dbt_destroy(t_3c_3)
1755 CALL dbt_destroy(t_3c_4)
1756 CALL dbt_destroy(t_3c_5)
1757 CALL dbt_destroy(t_3c_6)
1758 CALL dbt_destroy(t_3c_7)
1759 CALL dbt_destroy(t_3c_8)
1760 CALL dbt_destroy(t_3c_sparse)
1761 CALL dbt_destroy(t_3c_help_1)
1762 CALL dbt_destroy(t_3c_help_2)
1763 CALL dbt_destroy(t_3c_ints)
1764 CALL dbt_destroy(t_3c_work)
1765 CALL dbt_destroy(t_r_occ)
1766 CALL dbt_destroy(t_r_virt)
1767 CALL dbt_destroy(t_dm_occ)
1768 CALL dbt_destroy(t_dm_virt)
1769 CALL dbt_destroy(t_kbkt)
1770 CALL dbt_destroy(t_m_occ)
1771 CALL dbt_destroy(t_m_virt)
1772 CALL dbcsr_release(r_occ)
1773 CALL dbcsr_release(r_virt)
1774 CALL dbcsr_release(dbcsr_work_symm)
1775 CALL dbcsr_release(dbcsr_work1)
1776 CALL dbcsr_release(dbcsr_work2)
1777 CALL dbcsr_release(dbcsr_work3)
1778 CALL dbcsr_release(exp_occ)
1779 CALL dbcsr_release(exp_virt)
1781 CALL dbt_destroy(t_2c_ri)
1782 CALL dbt_destroy(t_2c_ri_2)
1783 CALL dbt_destroy(t_2c_ao)
1788 CALL timestop(handle)
1810 SUBROUTINE perform_2c_ops(force, t_KBKT, force_data, fac, t_B, t_P, t_2c_RI, t_2c_RI_2, use_virial, &
1811 atom_of_kind, kind_of, eps_filter, dbcsr_nflop, unit_nr_dbcsr)
1813 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1814 TYPE(dbt_type),
INTENT(INOUT) :: t_kbkt
1815 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
1816 REAL(
dp),
INTENT(IN) ::
fac
1817 TYPE(dbt_type),
INTENT(INOUT) :: t_b, t_p, t_2c_ri, t_2c_ri_2
1818 LOGICAL,
INTENT(IN) :: use_virial
1819 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
1820 REAL(
dp),
INTENT(IN) :: eps_filter
1821 INTEGER(int_8),
INTENT(INOUT) :: dbcsr_nflop
1822 INTEGER,
INTENT(IN) :: unit_nr_dbcsr
1824 CHARACTER(LEN=*),
PARAMETER :: routinen =
'perform_2c_ops'
1827 INTEGER(int_8) :: flop
1829 TYPE(dbt_type) :: t_2c_tmp, t_2c_virial
1831 CALL timeset(routinen, handle)
1833 IF (use_virial)
CALL dbt_create(force_data%RI_virial_pot, t_2c_virial)
1836 CALL dbt_contract(1.0_dp, force_data%t_2c_K, t_b, 0.0_dp, t_2c_ri, &
1837 contract_1=[2], notcontract_1=[1], &
1838 contract_2=[1], notcontract_2=[2], &
1839 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1840 flop=flop, unit_nr=unit_nr_dbcsr)
1841 dbcsr_nflop = dbcsr_nflop + flop
1843 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_kbkt, &
1844 contract_1=[2], notcontract_1=[1], &
1845 contract_2=[2], notcontract_2=[1], &
1846 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1847 flop=flop, unit_nr=unit_nr_dbcsr)
1848 dbcsr_nflop = dbcsr_nflop + flop
1850 CALL dbt_contract(2.0_dp, t_p, t_2c_ri, 0.0_dp, t_2c_ri_2, &
1851 contract_1=[2], notcontract_1=[1], &
1852 contract_2=[1], notcontract_2=[2], &
1853 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1854 flop=flop, unit_nr=unit_nr_dbcsr)
1855 dbcsr_nflop = dbcsr_nflop + flop
1856 CALL dbt_clear(t_2c_ri)
1860 CALL dbt_contract(1.0_dp, force_data%t_2c_inv_metric, t_2c_ri_2, 0.0_dp, t_2c_ri, &
1861 contract_1=[2], notcontract_1=[1], &
1862 contract_2=[1], notcontract_2=[2], &
1863 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1864 flop=flop, unit_nr=unit_nr_dbcsr)
1865 dbcsr_nflop = dbcsr_nflop + flop
1867 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_K, 0.0_dp, t_2c_ri_2, &
1868 contract_1=[2], notcontract_1=[1], &
1869 contract_2=[2], notcontract_2=[1], &
1870 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1871 flop=flop, unit_nr=unit_nr_dbcsr)
1872 dbcsr_nflop = dbcsr_nflop + flop
1876 CALL get_2c_der_force(force, t_2c_ri_2, force_data%t_2c_der_metric, atom_of_kind, &
1877 kind_of, force_data%idx_to_at_RI, pref, do_mp2=.true.)
1878 IF (use_virial)
THEN
1879 CALL dbt_copy(t_2c_ri_2, t_2c_virial)
1880 CALL dbt_scale(t_2c_virial, pref)
1881 CALL dbt_copy_tensor_to_matrix(t_2c_virial, force_data%RI_virial_met, summation=.true.)
1882 CALL dbt_clear(t_2c_virial)
1887 CALL dbt_contract(1.0_dp, t_2c_ri, force_data%t_2c_pot_msqrt, 0.0_dp, t_2c_ri_2, &
1888 contract_1=[2], notcontract_1=[1], &
1889 contract_2=[1], notcontract_2=[2], &
1890 map_1=[1], map_2=[2], filter_eps=eps_filter, &
1891 flop=flop, unit_nr=unit_nr_dbcsr)
1892 dbcsr_nflop = dbcsr_nflop + flop
1896 IF (force_data%do_periodic)
THEN
1897 CALL dbt_scale(t_2c_ri_2, pref)
1898 CALL dbt_create(force_data%G_PQ, t_2c_tmp)
1899 CALL dbt_copy(t_2c_ri_2, t_2c_tmp, move_data=.true.)
1900 CALL dbt_copy_tensor_to_matrix(t_2c_tmp, force_data%G_PQ, summation=.true.)
1901 CALL dbt_destroy(t_2c_tmp)
1903 CALL get_2c_der_force(force, t_2c_ri_2, force_data%t_2c_der_pot, atom_of_kind, &
1904 kind_of, force_data%idx_to_at_RI, pref, do_mp2=.true.)
1906 IF (use_virial)
THEN
1907 CALL dbt_copy(t_2c_ri_2, t_2c_virial)
1908 CALL dbt_scale(t_2c_virial, pref)
1909 CALL dbt_copy_tensor_to_matrix(t_2c_virial, force_data%RI_virial_pot, summation=.true.)
1910 CALL dbt_clear(t_2c_virial)
1914 CALL dbt_clear(t_2c_ri)
1915 CALL dbt_clear(t_2c_ri_2)
1917 IF (use_virial)
CALL dbt_destroy(t_2c_virial)
1919 CALL timestop(handle)
1921 END SUBROUTINE perform_2c_ops
1969 SUBROUTINE perform_3c_ops(force, t_R_occ, t_R_virt, force_data, fac, cut_memory, n_mem_RI, &
1970 t_KBKT, t_dm_occ, t_dm_virt, t_3c_O, t_3c_M, t_M_occ, t_M_virt, t_3c_0, t_3c_1, &
1971 t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, t_3c_help_1, t_3c_help_2, &
1972 t_3c_ints, t_3c_work, starts_array_mc, ends_array_mc, batch_start_RI, &
1973 batch_end_RI, t_3c_O_compressed, t_3c_O_ind, use_virial, &
1974 atom_of_kind, kind_of, eps_filter, occ_ddint, nze_ddint, dbcsr_nflop, &
1975 unit_nr_dbcsr, mp2_env)
1977 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
1978 TYPE(dbt_type),
INTENT(INOUT) :: t_r_occ, t_r_virt
1979 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
1980 REAL(
dp),
INTENT(IN) ::
fac
1981 INTEGER,
INTENT(IN) :: cut_memory, n_mem_ri
1982 TYPE(dbt_type),
INTENT(INOUT) :: t_kbkt, t_dm_occ, t_dm_virt, t_3c_o, t_3c_m, t_m_occ, &
1983 t_m_virt, t_3c_0, t_3c_1, t_3c_3, t_3c_4, t_3c_5, t_3c_6, t_3c_7, t_3c_8, t_3c_sparse, &
1984 t_3c_help_1, t_3c_help_2, t_3c_ints, t_3c_work
1985 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc, &
1986 batch_start_ri, batch_end_ri
1987 TYPE(hfx_compression_type),
DIMENSION(:) :: t_3c_o_compressed
1988 TYPE(block_ind_type),
DIMENSION(:),
INTENT(INOUT) :: t_3c_o_ind
1989 LOGICAL,
INTENT(IN) :: use_virial
1990 INTEGER,
DIMENSION(:),
INTENT(IN) :: atom_of_kind, kind_of
1991 REAL(
dp),
INTENT(IN) :: eps_filter
1992 REAL(
dp),
INTENT(INOUT) :: occ_ddint
1993 INTEGER(int_8),
INTENT(INOUT) :: nze_ddint, dbcsr_nflop
1994 INTEGER,
INTENT(IN) :: unit_nr_dbcsr
1995 TYPE(mp2_type) :: mp2_env
1997 CHARACTER(LEN=*),
PARAMETER :: routinen =
'perform_3c_ops'
1999 INTEGER :: dummy_int, handle, handle2, i_mem, &
2001 INTEGER(int_8) :: flop, nze
2002 INTEGER,
DIMENSION(2, 1) :: ibounds, jbounds, kbounds
2003 INTEGER,
DIMENSION(2, 2) :: bounds_2c
2004 INTEGER,
DIMENSION(2, 3) :: bounds_cpy
2005 INTEGER,
DIMENSION(3) :: bounds_3c
2006 REAL(
dp) :: memory, occ, pref
2007 TYPE(block_ind_type),
ALLOCATABLE,
DIMENSION(:, :) :: blk_indices
2008 TYPE(hfx_compression_type),
ALLOCATABLE, &
2009 DIMENSION(:, :) :: store_3c
2011 CALL timeset(routinen, handle)
2013 CALL dbt_get_info(t_3c_m, nfull_total=bounds_3c)
2016 ALLOCATE (store_3c(n_mem_ri, cut_memory))
2017 ALLOCATE (blk_indices(n_mem_ri, cut_memory))
2019 CALL timeset(routinen//
"_pre_3c", handle2)
2021 CALL dbt_copy(t_3c_o, t_3c_0)
2022 DO i_mem = 1, cut_memory
2023 CALL decompress_tensor(t_3c_o, t_3c_o_ind(i_mem)%ind, t_3c_o_compressed(i_mem), &
2024 mp2_env%ri_rpa_im_time%eps_compress)
2025 CALL dbt_copy(t_3c_o, t_3c_ints)
2026 CALL dbt_copy(t_3c_o, t_3c_0, move_data=.true., summation=.true.)
2028 DO k_mem = 1, n_mem_ri
2029 kbounds(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2034 CALL dbt_batched_contract_init(t_kbkt)
2035 CALL dbt_contract(1.0_dp, t_kbkt, t_3c_ints, 0.0_dp, t_3c_work, &
2036 contract_1=[2], notcontract_1=[1], &
2037 contract_2=[1], notcontract_2=[2, 3], &
2038 map_1=[1], map_2=[2, 3], filter_eps=eps_filter, &
2039 bounds_2=kbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2040 CALL dbt_batched_contract_finalize(t_kbkt)
2041 dbcsr_nflop = dbcsr_nflop + flop
2043 CALL dbt_copy(t_3c_work, t_3c_m, move_data=.true.)
2044 CALL compress_tensor(t_3c_m, blk_indices(k_mem, i_mem)%ind, store_3c(k_mem, i_mem), &
2045 mp2_env%ri_rpa_im_time%eps_compress, memory)
2048 CALL dbt_clear(t_3c_m)
2049 CALL dbt_copy(t_3c_m, t_3c_ints)
2050 CALL timestop(handle2)
2052 CALL dbt_batched_contract_init(t_r_occ)
2053 CALL dbt_batched_contract_init(t_r_virt)
2054 DO i_mem = 1, cut_memory
2055 ibounds(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2058 CALL timeset(routinen//
"_3c_M", handle2)
2059 CALL dbt_batched_contract_init(t_dm_occ)
2060 CALL dbt_contract(1.0_dp, t_3c_0, t_dm_occ, 0.0_dp, t_3c_1, &
2061 contract_1=[3], notcontract_1=[1, 2], &
2062 contract_2=[1], notcontract_2=[2], &
2063 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2064 bounds_3=ibounds, flop=flop, unit_nr=unit_nr_dbcsr)
2065 dbcsr_nflop = dbcsr_nflop + flop
2066 CALL dbt_batched_contract_finalize(t_dm_occ)
2067 CALL dbt_copy(t_3c_1, t_m_occ, order=[1, 3, 2], move_data=.true.)
2069 CALL dbt_batched_contract_init(t_dm_virt)
2070 CALL dbt_contract(1.0_dp, t_3c_0, t_dm_virt, 0.0_dp, t_3c_1, &
2071 contract_1=[3], notcontract_1=[1, 2], &
2072 contract_2=[1], notcontract_2=[2], &
2073 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2074 bounds_3=ibounds, flop=flop, unit_nr=unit_nr_dbcsr)
2075 dbcsr_nflop = dbcsr_nflop + flop
2076 CALL dbt_batched_contract_finalize(t_dm_virt)
2077 CALL dbt_copy(t_3c_1, t_m_virt, order=[1, 3, 2], move_data=.true.)
2078 CALL timestop(handle2)
2081 CALL timeset(routinen//
"_3c_R", handle2)
2082 DO k_mem = 1, n_mem_ri
2083 CALL decompress_tensor(t_3c_m, blk_indices(k_mem, i_mem)%ind, store_3c(k_mem, i_mem), &
2084 mp2_env%ri_rpa_im_time%eps_compress)
2085 CALL dbt_copy(t_3c_m, t_3c_3, move_data=.true.)
2087 CALL dbt_contract(1.0_dp, t_m_occ, t_3c_3, 1.0_dp, t_r_occ, &
2088 contract_1=[1, 2], notcontract_1=[3], &
2089 contract_2=[1, 2], notcontract_2=[3], &
2090 map_1=[1], map_2=[2], filter_eps=eps_filter, &
2091 flop=flop, unit_nr=unit_nr_dbcsr)
2092 dbcsr_nflop = dbcsr_nflop + flop
2094 CALL dbt_contract(1.0_dp, t_m_virt, t_3c_3, 1.0_dp, t_r_virt, &
2095 contract_1=[1, 2], notcontract_1=[3], &
2096 contract_2=[1, 2], notcontract_2=[3], &
2097 map_1=[1], map_2=[2], filter_eps=eps_filter, &
2098 flop=flop, unit_nr=unit_nr_dbcsr)
2099 dbcsr_nflop = dbcsr_nflop + flop
2101 CALL dbt_copy(t_3c_m, t_3c_3)
2102 CALL dbt_copy(t_3c_m, t_m_virt)
2103 CALL timestop(handle2)
2105 CALL dbt_copy(t_m_occ, t_3c_4, move_data=.true.)
2107 DO j_mem = 1, cut_memory
2108 jbounds(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2110 bounds_cpy(:, 1) = [1, bounds_3c(1)]
2111 bounds_cpy(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2112 bounds_cpy(:, 3) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
2113 CALL dbt_copy(t_3c_sparse, t_3c_7, bounds=bounds_cpy)
2115 CALL dbt_batched_contract_init(t_dm_virt)
2116 DO k_mem = 1, n_mem_ri
2117 bounds_2c(:, 1) = [batch_start_ri(k_mem), batch_end_ri(k_mem)]
2118 bounds_2c(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
2120 CALL timeset(routinen//
"_3c_dm", handle2)
2124 CALL dbt_contract(2.0_dp, t_3c_4, t_dm_virt, 0.0_dp, t_3c_5, &
2125 contract_1=[3], notcontract_1=[1, 2], &
2126 contract_2=[1], notcontract_2=[2], &
2127 map_1=[1, 2], map_2=[3], filter_eps=eps_filter, &
2128 bounds_2=bounds_2c, bounds_3=jbounds, flop=flop, unit_nr=unit_nr_dbcsr)
2129 dbcsr_nflop = dbcsr_nflop + flop
2132 nze_ddint = nze_ddint + nze
2133 occ_ddint = occ_ddint + occ
2135 CALL dbt_copy(t_3c_5, t_3c_6, move_data=.true.)
2136 CALL timestop(handle2)
2139 CALL timeset(routinen//
"_3c_KBK", handle2)
2140 CALL dbt_batched_contract_init(t_kbkt)
2141 CALL dbt_contract(1.0_dp, t_kbkt, t_3c_6, 0.0_dp, t_3c_7, &
2142 contract_1=[2], notcontract_1=[1], &
2143 contract_2=[1], notcontract_2=[2, 3], &
2144 map_1=[1], map_2=[2, 3], &
2145 retain_sparsity=.true., flop=flop, unit_nr=unit_nr_dbcsr)
2146 dbcsr_nflop = dbcsr_nflop + flop
2147 CALL dbt_batched_contract_finalize(t_kbkt)
2148 CALL timestop(handle2)
2149 CALL dbt_copy(t_3c_7, t_3c_8, summation=.true.)
2152 CALL dbt_batched_contract_finalize(t_dm_virt)
2155 CALL dbt_copy(t_3c_8, t_3c_help_1, move_data=.true.)
2158 DO k_mem = 1, cut_memory
2160 CALL dbt_clear(force_data%t_3c_der_RI(i_xyz))
2161 CALL decompress_tensor(force_data%t_3c_der_RI(i_xyz), force_data%t_3c_der_RI_ind(k_mem, i_xyz)%ind, &
2162 force_data%t_3c_der_RI_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2165 force_data%idx_to_at_RI, pref, do_mp2=.true., deriv_dim=1)
2168 IF (use_virial)
THEN
2169 CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2170 CALL dbt_scale(t_3c_help_2, pref)
2171 CALL dbt_copy(t_3c_help_2, force_data%t_3c_virial_split, summation=.true., move_data=.true.)
2174 CALL dbt_copy(t_3c_help_1, t_3c_help_2)
2175 CALL dbt_copy(t_3c_help_1, t_3c_help_2, order=[1, 3, 2], move_data=.true., summation=.true.)
2176 DO k_mem = 1, cut_memory
2178 CALL dbt_clear(force_data%t_3c_der_AO(i_xyz))
2179 CALL decompress_tensor(force_data%t_3c_der_AO(i_xyz), force_data%t_3c_der_AO_ind(k_mem, i_xyz)%ind, &
2180 force_data%t_3c_der_AO_comp(k_mem, i_xyz), mp2_env%ri_rpa_im_time%eps_compress)
2183 force_data%idx_to_at_AO, pref, do_mp2=.true., deriv_dim=3)
2186 CALL dbt_clear(t_3c_help_2)
2188 CALL dbt_batched_contract_finalize(t_r_occ)
2189 CALL dbt_batched_contract_finalize(t_r_virt)
2191 DO k_mem = 1, n_mem_ri
2192 DO i_mem = 1, cut_memory
2196 DEALLOCATE (store_3c, blk_indices)
2198 CALL timestop(handle)
2200 END SUBROUTINE perform_3c_ops
2212 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
2213 INTEGER,
INTENT(IN) :: unit_nr
2214 TYPE(qs_environment_type),
POINTER :: qs_env
2216 CHARACTER(len=*),
PARAMETER :: routinen =
'calc_post_loop_forces'
2218 INTEGER :: handle, ispin, nao, nao_aux, nocc, nspins
2221 TYPE(admm_type),
POINTER :: admm_env
2222 TYPE(cp_fm_struct_type),
POINTER :: fm_struct
2223 TYPE(cp_fm_type),
ALLOCATABLE,
DIMENSION(:) :: cpmos, mo_occ
2224 TYPE(cp_fm_type),
POINTER :: mo_coeff
2225 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: dbcsr_p_work, matrix_p_mp2, &
2226 matrix_p_mp2_admm, matrix_s, &
2227 matrix_s_aux, work_admm, yp_admm
2228 TYPE(dft_control_type),
POINTER :: dft_control
2229 TYPE(linres_control_type),
POINTER :: linres_control
2230 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
2231 TYPE(qs_p_env_type),
POINTER :: p_env
2232 TYPE(section_vals_type),
POINTER :: hfx_section, lr_section
2234 NULLIFY (linres_control, p_env, dft_control, matrix_s, mos, mo_coeff, fm_struct, lr_section, &
2235 dbcsr_p_work, yp_admm, matrix_p_mp2, admm_env, work_admm, matrix_s_aux, matrix_p_mp2_admm)
2237 CALL timeset(routinen, handle)
2239 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_s=matrix_s, mos=mos)
2240 nspins = dft_control%nspins
2247 ALLOCATE (linres_control)
2250 CALL section_vals_val_get(lr_section,
"PRECONDITIONER", i_val=linres_control%preconditioner_type)
2253 linres_control%do_kernel = .true.
2254 linres_control%lr_triplet = .false.
2255 linres_control%converged = .false.
2256 linres_control%eps_filter = qs_env%mp2_env%ri_rpa_im_time%eps_filter
2258 CALL set_qs_env(qs_env, linres_control=linres_control)
2260 IF (unit_nr > 0)
THEN
2262 WRITE (unit_nr,
'(T3,A)')
'MP2_CPHF| Iterative solution of Z-Vector equations'
2263 WRITE (unit_nr,
'(T3,A,T45,ES8.1)')
'MP2_CPHF| Convergence threshold:', linres_control%eps
2264 WRITE (unit_nr,
'(T3,A,T45,I8)')
'MP2_CPHF| Maximum number of iterations: ', linres_control%max_iter
2268 CALL p_env_create(p_env, qs_env, orthogonal_orbitals=.true., linres_control=linres_control)
2275 DO ispin = 1, nspins
2276 ALLOCATE (p_env%p1(ispin)%matrix, p_env%w1(ispin)%matrix, dbcsr_p_work(ispin)%matrix)
2277 CALL dbcsr_create(matrix=p_env%p1(ispin)%matrix, template=matrix_s(1)%matrix)
2278 CALL dbcsr_create(matrix=p_env%w1(ispin)%matrix, template=matrix_s(1)%matrix)
2279 CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2280 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_s(1)%matrix)
2281 CALL dbcsr_copy(p_env%w1(ispin)%matrix, matrix_s(1)%matrix)
2282 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2283 CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
2284 CALL dbcsr_set(p_env%w1(ispin)%matrix, 0.0_dp)
2285 CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2288 IF (dft_control%do_admm)
THEN
2289 CALL get_admm_env(qs_env%admm_env, matrix_s_aux_fit=matrix_s_aux)
2292 DO ispin = 1, nspins
2293 ALLOCATE (p_env%p1_admm(ispin)%matrix, work_admm(ispin)%matrix)
2294 CALL dbcsr_create(p_env%p1_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2295 CALL dbcsr_copy(p_env%p1_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2296 CALL dbcsr_set(p_env%p1_admm(ispin)%matrix, 0.0_dp)
2297 CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2298 CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2299 CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2304 CALL prepare_for_response(force_data, qs_env)
2305 ALLOCATE (cpmos(nspins), mo_occ(nspins))
2306 DO ispin = 1, nspins
2307 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, homo=nocc)
2310 template_fmstruct=mo_coeff%matrix_struct)
2314 CALL cp_fm_to_fm(mo_coeff, mo_occ(ispin), nocc)
2329 ext_hfx_section=hfx_section, &
2330 x_data=qs_env%mp2_env%ri_rpa%x_data, &
2331 recalc_integrals=.false., &
2332 do_admm=qs_env%mp2_env%ri_rpa%do_admm, &
2334 reuse_hfx=qs_env%mp2_env%ri_rpa%reuse_hfx)
2338 IF (nspins == 1) focc = 4.0_dp
2339 DO ispin = 1, nspins
2340 CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, homo=nocc)
2342 cpmos(ispin), nocc, &
2343 alpha=focc, beta=0.0_dp)
2350 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2)
2352 DO ispin = 1, nspins
2353 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, p_env%p1(ispin)%matrix)
2354 CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, force_data%sum_YP_tau(ispin)%matrix, 1.0_dp, 1.0_dp)
2356 CALL set_ks_env(qs_env%ks_env, matrix_p_mp2=dbcsr_p_work)
2358 IF (dft_control%do_admm)
THEN
2360 CALL get_qs_env(qs_env, matrix_p_mp2_admm=matrix_p_mp2_admm, admm_env=admm_env)
2361 nao = admm_env%nao_orb
2362 nao_aux = admm_env%nao_aux_fit
2364 DO ispin = 1, nspins
2367 CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2368 CALL parallel_gemm(
'N',
'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2369 0.0_dp, admm_env%work_aux_orb)
2370 CALL parallel_gemm(
'N',
'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2371 0.0_dp, admm_env%work_aux_aux)
2372 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, work_admm(ispin)%matrix, keep_sparsity=.true.)
2375 ALLOCATE (yp_admm(ispin)%matrix)
2376 CALL dbcsr_create(yp_admm(ispin)%matrix, template=work_admm(ispin)%matrix)
2377 CALL dbcsr_copy(yp_admm(ispin)%matrix, work_admm(ispin)%matrix)
2379 CALL dbcsr_add(work_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
2382 CALL set_ks_env(qs_env%ks_env, matrix_p_mp2_admm=work_admm)
2386 CALL update_im_time_forces(p_env, force_data%sum_O_tau, force_data%sum_YP_tau, yp_admm, qs_env)
2391 CALL cp_fm_release(cpmos)
2392 CALL cp_fm_release(mo_occ)
2396 CALL timestop(handle)
2406 SUBROUTINE prepare_for_response(force_data, qs_env)
2408 TYPE(im_time_force_type),
INTENT(INOUT) :: force_data
2409 TYPE(qs_environment_type),
POINTER :: qs_env
2411 CHARACTER(LEN=*),
PARAMETER :: routinen =
'prepare_for_response'
2413 INTEGER :: handle, ispin, nao, nao_aux, nspins
2414 LOGICAL :: do_hfx, do_tau, do_tau_admm
2415 REAL(
dp) :: ehartree
2416 TYPE(admm_type),
POINTER :: admm_env
2417 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: dbcsr_p_work, ker_tau_admm, matrix_s, &
2418 matrix_s_aux, work_admm
2419 TYPE(dbcsr_type) :: dbcsr_work
2420 TYPE(dft_control_type),
POINTER :: dft_control
2421 TYPE(pw_c1d_gs_type) :: rhoz_tot_gspace, zv_hartree_gspace
2422 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rhoz_g
2423 TYPE(pw_env_type),
POINTER :: pw_env
2424 TYPE(pw_poisson_type),
POINTER :: poisson_env
2425 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
2426 TYPE(pw_r3d_rs_type) :: zv_hartree_rspace
2427 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau
2428 TYPE(qs_rho_type),
POINTER :: rho, rho_aux_fit
2429 TYPE(section_vals_type),
POINTER :: hfx_section, xc_section
2430 TYPE(task_list_type),
POINTER :: task_list_aux_fit
2432 NULLIFY (pw_env, rhoz_r, rhoz_g, tauz_r, v_xc, v_xc_tau, &
2433 poisson_env, auxbas_pw_pool, dft_control, admm_env, xc_section, rho, rho_aux_fit, &
2434 task_list_aux_fit, ker_tau_admm, work_admm, dbcsr_p_work, matrix_s, hfx_section)
2436 CALL timeset(routinen, handle)
2438 CALL get_qs_env(qs_env, dft_control=dft_control, pw_env=pw_env, rho=rho, matrix_s=matrix_s)
2439 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
2440 nspins = dft_control%nspins
2443 DO ispin = 1, nspins
2444 ALLOCATE (dbcsr_p_work(ispin)%matrix)
2445 CALL dbcsr_create(matrix=dbcsr_p_work(ispin)%matrix, template=matrix_s(1)%matrix)
2446 CALL dbcsr_copy(dbcsr_p_work(ispin)%matrix, matrix_s(1)%matrix)
2447 CALL dbcsr_set(dbcsr_p_work(ispin)%matrix, 0.0_dp)
2451 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
2452 DO ispin = 1, nspins
2453 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
2454 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
2456 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
2457 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
2458 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
2460 CALL pw_zero(rhoz_tot_gspace)
2461 DO ispin = 1, nspins
2462 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2463 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
2464 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
2467 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
2470 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
2471 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
2476 TYPE(pw_c1d_gs_type) :: tauz_g
2477 ALLOCATE (tauz_r(nspins))
2478 CALL auxbas_pw_pool%create_pw(tauz_g)
2479 DO ispin = 1, nspins
2480 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
2482 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=force_data%sum_YP_tau(ispin)%matrix, &
2483 rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
2485 CALL auxbas_pw_pool%give_back_pw(tauz_g)
2489 IF (dft_control%do_admm)
THEN
2491 xc_section => admm_env%xc_section_primary
2497 CALL create_kernel(qs_env, v_xc, v_xc_tau, rho, rhoz_r, rhoz_g, tauz_r, xc_section)
2499 DO ispin = 1, nspins
2500 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2501 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
2502 CALL integrate_v_rspace(qs_env=qs_env, &
2503 v_rspace=v_xc(ispin), &
2504 hmat=dbcsr_p_work(ispin), &
2505 calculate_forces=.false.)
2507 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2509 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
2510 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
2511 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
2515 DO ispin = 1, nspins
2516 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2517 CALL integrate_v_rspace(qs_env=qs_env, &
2518 v_rspace=v_xc_tau(ispin), &
2519 hmat=dbcsr_p_work(ispin), &
2520 compute_tau=.true., &
2521 calculate_forces=.false.)
2522 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2524 DEALLOCATE (v_xc_tau)
2528 IF (dft_control%do_admm)
THEN
2530 CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux, &
2531 task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit)
2533 CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
2537 DO ispin = 1, nspins
2538 ALLOCATE (work_admm(ispin)%matrix, ker_tau_admm(ispin)%matrix)
2539 CALL dbcsr_create(work_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2540 CALL dbcsr_copy(work_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2541 CALL dbcsr_set(work_admm(ispin)%matrix, 0.0_dp)
2542 CALL dbcsr_create(ker_tau_admm(ispin)%matrix, template=matrix_s_aux(1)%matrix)
2543 CALL dbcsr_copy(ker_tau_admm(ispin)%matrix, matrix_s_aux(1)%matrix)
2544 CALL dbcsr_set(ker_tau_admm(ispin)%matrix, 0.0_dp)
2548 cpassert(
ASSOCIATED(admm_env%work_orb_orb))
2549 cpassert(
ASSOCIATED(admm_env%work_aux_orb))
2550 cpassert(
ASSOCIATED(admm_env%work_aux_aux))
2551 nao = admm_env%nao_orb
2552 nao_aux = admm_env%nao_aux_fit
2553 DO ispin = 1, nspins
2554 CALL copy_dbcsr_to_fm(force_data%sum_YP_tau(ispin)%matrix, admm_env%work_orb_orb)
2555 CALL parallel_gemm(
'N',
'N', nao_aux, nao, nao, 1.0_dp, admm_env%A, admm_env%work_orb_orb, &
2556 0.0_dp, admm_env%work_aux_orb)
2557 CALL parallel_gemm(
'N',
'T', nao_aux, nao_aux, nao, 1.0_dp, admm_env%work_aux_orb, admm_env%A, &
2558 0.0_dp, admm_env%work_aux_aux)
2559 CALL copy_fm_to_dbcsr(admm_env%work_aux_aux, ker_tau_admm(ispin)%matrix, keep_sparsity=.true.)
2563 DO ispin = 1, nspins
2564 CALL pw_zero(rhoz_r(ispin))
2565 CALL pw_zero(rhoz_g(ispin))
2567 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
2568 basis_type=
"AUX_FIT", task_list_external=task_list_aux_fit)
2571 IF (do_tau_admm)
THEN
2573 TYPE(pw_c1d_gs_type) :: tauz_g
2574 CALL auxbas_pw_pool%create_pw(tauz_g)
2575 DO ispin = 1, nspins
2576 CALL pw_zero(tauz_r(ispin))
2578 rho=tauz_r(ispin), rho_gspace=tauz_g, &
2579 basis_type=
"AUX_FIT", task_list_external=task_list_aux_fit, &
2582 CALL auxbas_pw_pool%give_back_pw(tauz_g)
2586 xc_section => admm_env%xc_section_aux
2587 CALL create_kernel(qs_env, v_xc, v_xc_tau, rho_aux_fit, rhoz_r, rhoz_g, tauz_r, xc_section)
2589 DO ispin = 1, nspins
2590 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
2591 CALL integrate_v_rspace(qs_env=qs_env, &
2592 v_rspace=v_xc(ispin), &
2593 hmat=work_admm(ispin), &
2594 calculate_forces=.false., &
2595 basis_type=
"AUX_FIT", &
2596 task_list_external=task_list_aux_fit)
2597 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
2601 IF (do_tau_admm)
THEN
2602 DO ispin = 1, nspins
2603 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
2604 CALL integrate_v_rspace(qs_env=qs_env, &
2605 v_rspace=v_xc_tau(ispin), &
2606 hmat=work_admm(ispin), &
2607 calculate_forces=.false., &
2608 basis_type=
"AUX_FIT", &
2609 task_list_external=task_list_aux_fit, &
2611 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
2613 DEALLOCATE (v_xc_tau)
2618 DO ispin = 1, nspins
2619 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
2620 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
2622 DEALLOCATE (rhoz_r, rhoz_g)
2625 DO ispin = 1, nspins
2626 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
2632 hfx_section => section_vals_get_subs_vals(qs_env%input,
"DFT%XC%HF")
2633 CALL section_vals_get(hfx_section, explicit=do_hfx)
2635 IF (dft_control%do_admm)
THEN
2636 CALL tddft_hfx_matrix(work_admm, ker_tau_admm, qs_env, .false., .false.)
2639 CALL dbcsr_create(dbcsr_work, template=dbcsr_p_work(1)%matrix)
2640 CALL dbcsr_copy(dbcsr_work, dbcsr_p_work(1)%matrix)
2641 CALL dbcsr_set(dbcsr_work, 0.0_dp)
2642 DO ispin = 1, nspins
2643 CALL copy_dbcsr_to_fm(work_admm(ispin)%matrix, admm_env%work_aux_aux)
2644 CALL parallel_gemm(
'N',
'N', nao_aux, nao, nao_aux, 1.0_dp, admm_env%work_aux_aux, admm_env%A, &
2645 0.0_dp, admm_env%work_aux_orb)
2646 CALL parallel_gemm(
'T',
'N', nao, nao, nao_aux, 1.0_dp, admm_env%A, admm_env%work_aux_orb, &
2647 0.0_dp, admm_env%work_orb_orb)
2648 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
2649 CALL dbcsr_add(dbcsr_p_work(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
2651 CALL dbcsr_release(dbcsr_work)
2652 CALL dbcsr_deallocate_matrix_set(ker_tau_admm)
2654 CALL tddft_hfx_matrix(dbcsr_p_work, force_data%sum_YP_tau, qs_env, .false., .false.)
2658 DO ispin = 1, nspins
2659 CALL dbcsr_add(force_data%sum_O_tau(ispin)%matrix, dbcsr_p_work(ispin)%matrix, 1.0_dp, 1.0_dp)
2662 CALL dbcsr_deallocate_matrix_set(dbcsr_p_work)
2663 CALL dbcsr_deallocate_matrix_set(work_admm)
2665 CALL timestop(handle)
2667 END SUBROUTINE prepare_for_response
2678 SUBROUTINE get_2c_gpw_forces(G_PQ, force, h_stress, use_virial, mp2_env, qs_env)
2680 TYPE(dbcsr_type),
INTENT(INOUT) :: g_pq
2681 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
2682 REAL(dp),
DIMENSION(3, 3),
INTENT(INOUT) :: h_stress
2683 LOGICAL,
INTENT(IN) :: use_virial
2684 TYPE(mp2_type),
INTENT(INOUT) :: mp2_env
2685 TYPE(qs_environment_type),
POINTER :: qs_env
2687 CHARACTER(len=*),
PARAMETER :: routinen =
'get_2c_gpw_forces'
2689 INTEGER :: atom_a, color, handle, i, i_ri, i_xyz, iatom, igrid_level, ikind, ipgf, iset, j, &
2690 j_ri, jatom, lb_ri, n_ri, natom, ncoa, ncoms, nkind, nproc, nseta, o1, offset, pdims(2), &
2692 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, iproc_map, kind_of, &
2694 INTEGER,
DIMENSION(:),
POINTER :: col_dist, la_max, la_min, npgfa, nsgfa, &
2696 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, pgrid
2697 LOGICAL :: found, one_proc_group
2698 REAL(dp) :: cutoff_old, radius, relative_cutoff_old
2699 REAL(dp),
ALLOCATABLE,
DIMENSION(:) :: e_cutoff_old, wf_vector
2700 REAL(dp),
DIMENSION(3) :: force_a, force_b, ra
2701 REAL(dp),
DIMENSION(3, 3) :: my_virial_a, my_virial_b
2702 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: h_tmp, i_ab, pab, pblock, sphi_a, zeta
2703 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
2704 TYPE(cell_type),
POINTER :: cell
2705 TYPE(dbcsr_distribution_type) :: dbcsr_dist
2706 TYPE(dbcsr_type) :: tmp_g_pq
2707 TYPE(dft_control_type),
POINTER :: dft_control
2708 TYPE(gto_basis_set_p_type),
ALLOCATABLE, &
2709 DIMENSION(:),
TARGET :: basis_set_ri_aux
2710 TYPE(gto_basis_set_type),
POINTER :: basis_set_a
2711 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
2712 TYPE(mp_para_env_type),
POINTER :: para_env, para_env_ext
2713 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
2715 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
2716 TYPE(pw_c1d_gs_type) :: dvg(3), pot_g, rho_g, rho_g_copy
2717 TYPE(pw_env_type),
POINTER :: pw_env_ext
2718 TYPE(pw_poisson_type),
POINTER :: poisson_env
2719 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
2720 TYPE(pw_r3d_rs_type) :: psi_l, rho_r
2721 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
2722 TYPE(realspace_grid_type),
DIMENSION(:),
POINTER :: rs_v
2723 TYPE(task_list_type),
POINTER :: task_list_ext
2725 NULLIFY (sab_orb, task_list_ext, particle_set, qs_kind_set, dft_control, pw_env_ext, auxbas_pw_pool, &
2726 poisson_env, atomic_kind_set, para_env, cell, rs_v, mos, basis_set_a)
2728 CALL timeset(routinen, handle)
2730 CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env, sab_orb=sab_orb, &
2731 natom=natom, nkind=nkind, qs_kind_set=qs_kind_set, particle_set=particle_set, &
2732 mos=mos, cell=cell, atomic_kind_set=atomic_kind_set)
2740 IF (para_env%num_pe <= natom)
THEN
2742 pdims(2) = para_env%num_pe
2745 IF (
modulo(para_env%num_pe, i) == 0)
THEN
2746 pdims(1) = para_env%num_pe/i
2753 ALLOCATE (row_dist(natom), col_dist(natom))
2755 row_dist(iatom) =
modulo(iatom, pdims(1))
2758 col_dist(jatom) =
modulo(jatom, pdims(2))
2761 ALLOCATE (pgrid(0:pdims(1) - 1, 0:pdims(2) - 1))
2763 DO i = 0, pdims(1) - 1
2764 DO j = 0, pdims(2) - 1
2770 CALL dbcsr_distribution_new(dbcsr_dist, group=para_env%get_handle(), pgrid=pgrid, row_dist=row_dist, col_dist=col_dist)
2773 CALL dbcsr_create(tmp_g_pq, template=g_pq, matrix_type=dbcsr_type_no_symmetry, dist=dbcsr_dist)
2774 CALL dbcsr_complete_redistribute(g_pq, tmp_g_pq)
2776 ALLOCATE (basis_set_ri_aux(nkind), sizes_ri(natom))
2777 CALL basis_set_list_setup(basis_set_ri_aux,
"RI_AUX", qs_kind_set)
2778 CALL get_particle_set(particle_set, qs_kind_set, nsgf=sizes_ri, basis=basis_set_ri_aux)
2779 n_ri = sum(sizes_ri)
2781 one_proc_group = mp2_env%mp2_num_proc == 1
2782 ALLOCATE (para_env_ext)
2783 IF (one_proc_group)
THEN
2785 CALL para_env_ext%from_split(para_env, para_env%mepos)
2788 ncoms = min(pdims(2), para_env%num_pe/mp2_env%mp2_num_proc)
2789 DO i = 0, pdims(1) - 1
2790 DO j = 0, pdims(2) - 1
2791 IF (pgrid(i, j) == para_env%mepos) color =
modulo(j + 1, ncoms)
2794 CALL para_env_ext%from_split(para_env, color)
2798 CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2799 auxbas_pw_pool, poisson_env, task_list_ext, rho_r, rho_g, pot_g, psi_l, sab_orb)
2801 IF (use_virial)
THEN
2802 CALL auxbas_pw_pool%create_pw(rho_g_copy)
2804 CALL auxbas_pw_pool%create_pw(dvg(i_xyz))
2808 ALLOCATE (wf_vector(n_ri))
2810 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
2812 ALLOCATE (iproc_map(natom))
2818 IF (one_proc_group)
THEN
2821 IF (pgrid(row_dist(iatom), col_dist(jatom)) == para_env%mepos) iproc_map(iatom) = 1
2823 IF (.NOT. any(iproc_map == 1)) cycle
2825 IF (.NOT.
modulo(col_dist(jatom) + 1, ncoms) == color) cycle
2828 lb_ri = sum(sizes_ri(1:jatom - 1))
2829 ub_ri = lb_ri + sizes_ri(jatom)
2830 DO j_ri = lb_ri + 1, ub_ri
2833 wf_vector(j_ri) = 1.0_dp
2835 CALL calculate_wavefunction(mos(1)%mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
2836 qs_kind_set, cell, dft_control, particle_set, pw_env_ext, &
2837 basis_type=
"RI_AUX", external_vector=wf_vector)
2839 IF (use_virial)
THEN
2840 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, dvg)
2845 IF (one_proc_group)
THEN
2846 IF (.NOT. iproc_map(iatom) == 1) cycle
2849 CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2850 IF (.NOT. found) cycle
2852 i_ri = sum(sizes_ri(1:iatom - 1))
2853 wf_vector(i_ri + 1:i_ri + sizes_ri(iatom)) = pblock(:, j_ri - lb_ri)
2856 CALL pw_copy(rho_g, rho_g_copy)
2857 CALL calculate_wavefunction(mos(1)%mo_coeff, 1, psi_l, rho_g, atomic_kind_set, &
2858 qs_kind_set, cell, dft_control, particle_set, pw_env_ext, &
2859 basis_type=
"RI_AUX", external_vector=wf_vector)
2861 CALL calc_potential_gpw(psi_l, rho_g, poisson_env, pot_g, mp2_env%potential_parameter, &
2863 CALL virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, &
2864 mp2_env%potential_parameter, para_env_ext)
2866 CALL calc_potential_gpw(rho_r, rho_g, poisson_env, pot_g, mp2_env%potential_parameter)
2870 CALL pw_env_get(pw_env_ext, rs_grids=rs_v)
2871 CALL potential_pw2rs(rs_v, rho_r, pw_env_ext)
2876 IF (one_proc_group)
THEN
2877 IF (.NOT. iproc_map(iatom) == 1) cycle
2882 IF (use_virial)
THEN
2883 my_virial_a = 0.0_dp
2884 my_virial_b = 0.0_dp
2887 ikind = kind_of(iatom)
2888 atom_a = atom_of_kind(iatom)
2890 basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
2891 first_sgfa => basis_set_a%first_sgf
2892 la_max => basis_set_a%lmax
2893 la_min => basis_set_a%lmin
2894 nseta = basis_set_a%nset
2895 nsgfa => basis_set_a%nsgf_set
2896 sphi_a => basis_set_a%sphi
2897 zeta => basis_set_a%zet
2898 npgfa => basis_set_a%npgf
2900 ra(:) =
pbc(particle_set(iatom)%r, cell)
2902 CALL dbcsr_get_block_p(tmp_g_pq, iatom, jatom, pblock, found)
2903 IF (.NOT. found) cycle
2907 ncoa = npgfa(iset)*
ncoset(la_max(iset))
2908 sgfa = first_sgfa(1, iset)
2910 ALLOCATE (h_tmp(ncoa, 1)); h_tmp = 0.0_dp
2911 ALLOCATE (i_ab(nsgfa(iset), 1)); i_ab = 0.0_dp
2912 ALLOCATE (pab(ncoa, 1)); pab = 0.0_dp
2914 i_ab(1:nsgfa(iset), 1) = 2.0_dp*pblock(offset + 1:offset + nsgfa(iset), j_ri - lb_ri)
2915 CALL dgemm(
"N",
"N", ncoa, 1, nsgfa(iset), 1.0_dp, sphi_a(1, sgfa),
SIZE(sphi_a, 1), &
2916 i_ab(1, 1), nsgfa(iset), 0.0_dp, pab(1, 1), ncoa)
2918 igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, minval(zeta(:, iset)))
2922 IF (map_gaussian_here(rs_v(igrid_level), cell%h_inv, ra, 0, 1, 0))
THEN
2923 DO ipgf = 1, npgfa(iset)
2924 o1 = (ipgf - 1)*
ncoset(la_max(iset))
2925 igrid_level = gaussian_gridlevel(pw_env_ext%gridlevel_info, zeta(ipgf, iset))
2927 radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
2928 lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
2929 zetp=zeta(ipgf, iset), &
2930 eps=dft_control%qs_control%eps_gvg_rspace, &
2931 prefactor=1.0_dp, cutoff=1.0_dp)
2933 CALL integrate_pgf_product( &
2934 la_max=la_max(iset), zeta=zeta(ipgf, iset), la_min=la_min(iset), &
2935 lb_max=0, zetb=0.0_dp, lb_min=0, &
2936 ra=ra, rab=(/0.0_dp, 0.0_dp, 0.0_dp/), &
2937 rsgrid=rs_v(igrid_level), &
2938 hab=h_tmp, pab=pab, &
2942 calculate_forces=.true., &
2943 force_a=force_a, force_b=force_b, &
2944 use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
2950 offset = offset + nsgfa(iset)
2951 DEALLOCATE (pab, h_tmp, i_ab)
2954 force(ikind)%mp2_non_sep(:, atom_a) = force(ikind)%mp2_non_sep(:, atom_a) + force_a + force_b
2955 IF (use_virial) h_stress = h_stress + my_virial_a + my_virial_b
2961 IF (use_virial)
THEN
2962 CALL auxbas_pw_pool%give_back_pw(rho_g_copy)
2964 CALL auxbas_pw_pool%give_back_pw(dvg(i_xyz))
2968 CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_ext, pw_env_ext, &
2969 task_list_ext, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
2971 CALL dbcsr_release(tmp_g_pq)
2972 CALL dbcsr_distribution_release(dbcsr_dist)
2973 DEALLOCATE (col_dist, row_dist, pgrid)
2975 CALL mp_para_env_release(para_env_ext)
2977 CALL timestop(handle)
2979 END SUBROUTINE get_2c_gpw_forces
2988 SUBROUTINE get_2c_mme_forces(G_PQ, force, mp2_env, qs_env)
2990 TYPE(dbcsr_type),
INTENT(INOUT) :: g_pq
2991 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
2992 TYPE(mp2_type),
INTENT(INOUT) :: mp2_env
2993 TYPE(qs_environment_type),
POINTER :: qs_env
2995 CHARACTER(len=*),
PARAMETER :: routinen =
'get_2c_mme_forces'
2997 INTEGER :: atom_a, atom_b, blk, g_count, handle, i_xyz, iatom, ikind, iset, jatom, jkind, &
2998 jset, natom, nkind, nseta, nsetb, offset_hab_a, offset_hab_b, r_count, sgfa, sgfb
2999 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
3000 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
3002 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
3004 REAL(dp) :: new_force, pref
3005 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :) :: hab
3006 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: hdab
3007 REAL(dp),
DIMENSION(:, :),
POINTER :: pblock
3008 REAL(kind=dp),
DIMENSION(3) :: ra, rb
3009 REAL(kind=dp),
DIMENSION(:, :),
POINTER :: sphi_a, sphi_b, zeta, zetb
3010 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
3011 TYPE(cell_type),
POINTER :: cell
3012 TYPE(dbcsr_iterator_type) :: iter
3013 TYPE(gto_basis_set_p_type),
ALLOCATABLE, &
3014 DIMENSION(:),
TARGET :: basis_set_ri_aux
3015 TYPE(gto_basis_set_type),
POINTER :: basis_set_a, basis_set_b
3016 TYPE(mp_para_env_type),
POINTER :: para_env
3017 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
3018 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3020 NULLIFY (qs_kind_set, basis_set_a, basis_set_b, pblock, particle_set, &
3021 cell, la_max, la_min, lb_min, npgfa, lb_max, npgfb, nsgfa, &
3022 nsgfb, first_sgfa, first_sgfb, sphi_a, sphi_b, zeta, zetb, &
3023 atomic_kind_set, para_env)
3025 CALL timeset(routinen, handle)
3027 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, nkind=nkind, particle_set=particle_set, &
3028 cell=cell, atomic_kind_set=atomic_kind_set, natom=natom, para_env=para_env)
3030 CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
3032 ALLOCATE (basis_set_ri_aux(nkind))
3033 CALL basis_set_list_setup(basis_set_ri_aux,
"RI_AUX", qs_kind_set)
3035 g_count = 0; r_count = 0
3037 CALL dbcsr_iterator_start(iter, g_pq)
3038 DO WHILE (dbcsr_iterator_blocks_left(iter))
3040 CALL dbcsr_iterator_next_block(iter, row=iatom, column=jatom, blk=blk)
3041 CALL dbcsr_get_block_p(g_pq, iatom, jatom, pblock, found)
3042 IF (.NOT. found) cycle
3043 IF (iatom > jatom) cycle
3045 IF (iatom == jatom) pref = 1.0_dp
3047 ikind = kind_of(iatom)
3048 jkind = kind_of(jatom)
3050 atom_a = atom_of_kind(iatom)
3051 atom_b = atom_of_kind(jatom)
3053 basis_set_a => basis_set_ri_aux(ikind)%gto_basis_set
3054 first_sgfa => basis_set_a%first_sgf
3055 la_max => basis_set_a%lmax
3056 la_min => basis_set_a%lmin
3057 nseta = basis_set_a%nset
3058 nsgfa => basis_set_a%nsgf_set
3059 sphi_a => basis_set_a%sphi
3060 zeta => basis_set_a%zet
3061 npgfa => basis_set_a%npgf
3063 basis_set_b => basis_set_ri_aux(jkind)%gto_basis_set
3064 first_sgfb => basis_set_b%first_sgf
3065 lb_max => basis_set_b%lmax
3066 lb_min => basis_set_b%lmin
3067 nsetb = basis_set_b%nset
3068 nsgfb => basis_set_b%nsgf_set
3069 sphi_b => basis_set_b%sphi
3070 zetb => basis_set_b%zet
3071 npgfb => basis_set_b%npgf
3073 ra(:) =
pbc(particle_set(iatom)%r, cell)
3074 rb(:) =
pbc(particle_set(jatom)%r, cell)
3076 ALLOCATE (hab(basis_set_a%nsgf, basis_set_b%nsgf))
3077 ALLOCATE (hdab(3, basis_set_a%nsgf, basis_set_b%nsgf))
3079 hdab(:, :, :) = 0.0_dp
3083 sgfa = first_sgfa(1, iset)
3087 sgfb = first_sgfb(1, jset)
3089 CALL integrate_set_2c(mp2_env%eri_mme_param%par, mp2_env%potential_parameter, la_min(iset), &
3090 la_max(iset), lb_min(jset), lb_max(jset), npgfa(iset), npgfb(jset), &
3091 zeta(:, iset), zetb(:, jset), ra, rb, hab, nsgfa(iset), nsgfb(jset), &
3092 offset_hab_a, offset_hab_b, 0, 0, sphi_a, sphi_b, sgfa, sgfb, &
3093 nsgfa(iset), nsgfb(jset), do_eri_mme, hdab=hdab, &
3094 g_count=g_count, r_count=r_count)
3096 offset_hab_b = offset_hab_b + nsgfb(jset)
3098 offset_hab_a = offset_hab_a + nsgfa(iset)
3102 new_force = pref*sum(pblock(:, :)*hdab(i_xyz, :, :))
3103 force(ikind)%mp2_non_sep(i_xyz, atom_a) = force(ikind)%mp2_non_sep(i_xyz, atom_a) + new_force
3104 force(jkind)%mp2_non_sep(i_xyz, atom_b) = force(jkind)%mp2_non_sep(i_xyz, atom_b) - new_force
3107 DEALLOCATE (hab, hdab)
3109 CALL dbcsr_iterator_stop(iter)
3111 CALL cp_eri_mme_update_local_counts(mp2_env%eri_mme_param, para_env, g_count_2c=g_count, r_count_2c=r_count)
3113 CALL timestop(handle)
3115 END SUBROUTINE get_2c_mme_forces
3127 SUBROUTINE update_im_time_forces(p_env, matrix_hz, matrix_p_F, matrix_p_F_admm, qs_env)
3129 TYPE(qs_p_env_type),
POINTER :: p_env
3130 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_hz, matrix_p_f, matrix_p_f_admm
3131 TYPE(qs_environment_type),
POINTER :: qs_env
3133 CHARACTER(len=*),
PARAMETER :: routinen =
'update_im_time_forces'
3135 INTEGER :: handle, i, idens, ispin, n_rep_hf, nao, &
3136 nao_aux, nder, nimages, nocc, nspins
3137 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
3138 LOGICAL :: do_exx, do_hfx, do_tau, do_tau_admm, &
3140 REAL(dp) :: dummy_real1, dummy_real2, ehartree, &
3142 REAL(dp),
DIMENSION(3, 3) :: h_stress, pv_loc
3143 TYPE(admm_type),
POINTER :: admm_env
3144 TYPE(atomic_kind_type),
DIMENSION(:),
POINTER :: atomic_kind_set
3145 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: current_density, current_density_admm, &
3146 current_mat_h, matrix_p_mp2, matrix_p_mp2_admm, matrix_s, matrix_s_aux_fit, matrix_w, &
3147 rho_ao, rho_ao_aux, scrm, scrm_admm
3148 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: dbcsr_work_h, dbcsr_work_p
3149 TYPE(dbcsr_type) :: dbcsr_work
3150 TYPE(dft_control_type),
POINTER :: dft_control
3151 TYPE(hfx_type),
DIMENSION(:, :),
POINTER :: x_data
3152 TYPE(mo_set_type),
DIMENSION(:),
POINTER :: mos
3153 TYPE(mp_para_env_type),
POINTER :: para_env
3154 TYPE(neighbor_list_set_p_type),
DIMENSION(:), &
3155 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
3156 TYPE(particle_type),
DIMENSION(:),
POINTER :: particle_set
3157 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, rhoz_tot_gspace, &
3159 TYPE(pw_c1d_gs_type),
DIMENSION(:),
POINTER :: rhoz_g
3160 TYPE(pw_env_type),
POINTER :: pw_env
3161 TYPE(pw_poisson_type),
POINTER :: poisson_env
3162 TYPE(pw_pool_type),
POINTER :: auxbas_pw_pool
3163 TYPE(pw_r3d_rs_type) :: vh_rspace, vhxc_rspace, zv_hartree_rspace
3164 TYPE(pw_r3d_rs_type),
DIMENSION(:),
POINTER :: rhoz_r, tauz_r, v_xc, v_xc_tau, &
3165 vadmm_rspace, vtau_rspace, vxc_rspace
3166 TYPE(qs_force_type),
DIMENSION(:),
POINTER :: force
3167 TYPE(qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
3168 TYPE(qs_rho_type),
POINTER :: rho, rho_aux_fit
3169 TYPE(section_vals_type),
POINTER :: hfx_section, xc_section
3170 TYPE(task_list_type),
POINTER :: task_list_aux_fit
3171 TYPE(virial_type),
POINTER :: virial
3173 NULLIFY (scrm, rho, dft_control, matrix_p_mp2, matrix_s, matrix_p_mp2_admm, admm_env, sab_orb, &
3174 cell_to_index, dbcsr_work_p, dbcsr_work_h, sac_ae, sac_ppl, sap_ppnl, force, virial, &
3175 qs_kind_set, atomic_kind_set, particle_set, pw_env, poisson_env, auxbas_pw_pool, &
3176 task_list_aux_fit, matrix_s_aux_fit, scrm_admm, rho_aux_fit, rho_ao_aux, x_data, &
3177 hfx_section, xc_section, para_env, rhoz_g, rhoz_r, tauz_r, v_xc, v_xc_tau, &
3178 vxc_rspace, vtau_rspace, vadmm_rspace, rho_ao, matrix_w)
3180 CALL timeset(routinen, handle)
3182 CALL get_qs_env(qs_env, rho=rho, dft_control=dft_control, matrix_s=matrix_s, admm_env=admm_env, &
3183 sab_orb=sab_orb, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl, force=force, &
3184 virial=virial, particle_set=particle_set, qs_kind_set=qs_kind_set, &
3185 atomic_kind_set=atomic_kind_set, x_data=x_data, para_env=para_env)
3186 nspins = dft_control%nspins
3188 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
3189 IF (use_virial) virial%pv_calculate = .true.
3193 IF (qs_env%mp2_env%method == ri_rpa_method_gpw)
THEN
3194 hfx_section => section_vals_get_subs_vals(qs_env%input,
"DFT%XC%WF_CORRELATION%RI_RPA%HF")
3195 CALL section_vals_get(hfx_section, explicit=do_exx)
3199 CALL get_qs_env(qs_env, matrix_p_mp2=matrix_p_mp2, matrix_p_mp2_admm=matrix_p_mp2_admm)
3203 IF (nspins == 2)
CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, 1.0_dp)
3204 CALL build_kinetic_matrix(qs_env%ks_env, matrix_t=scrm, &
3205 matrix_name=
"KINETIC ENERGY MATRIX", &
3207 sab_nl=sab_orb, calculate_forces=.true., &
3208 matrix_p=matrix_p_mp2(1)%matrix)
3209 IF (nspins == 2)
CALL dbcsr_add(matrix_p_mp2(1)%matrix, matrix_p_mp2(2)%matrix, 1.0_dp, -1.0_dp)
3210 CALL dbcsr_deallocate_matrix_set(scrm)
3213 CALL dbcsr_allocate_matrix_set(scrm, nspins)
3214 DO ispin = 1, nspins
3215 ALLOCATE (scrm(ispin)%matrix)
3216 CALL dbcsr_create(scrm(ispin)%matrix, template=matrix_s(1)%matrix)
3217 CALL dbcsr_copy(scrm(ispin)%matrix, matrix_s(1)%matrix)
3218 CALL dbcsr_set(scrm(ispin)%matrix, 0.0_dp)
3223 NULLIFY (cell_to_index)
3224 ALLOCATE (dbcsr_work_p(nspins, 1), dbcsr_work_h(nspins, 1))
3225 DO ispin = 1, nspins
3226 dbcsr_work_p(ispin, 1)%matrix => matrix_p_mp2(ispin)%matrix
3227 dbcsr_work_h(ispin, 1)%matrix => scrm(ispin)%matrix
3230 IF (
ASSOCIATED(sac_ae))
THEN
3231 CALL build_core_ae(dbcsr_work_h, dbcsr_work_p, force, &
3232 virial, .true., use_virial, nder, &
3233 qs_kind_set, atomic_kind_set, particle_set, &
3234 sab_orb, sac_ae, nimages, cell_to_index)
3237 IF (
ASSOCIATED(sac_ppl))
THEN
3238 CALL build_core_ppl(dbcsr_work_h, dbcsr_work_p, force, &
3239 virial, .true., use_virial, nder, &
3240 qs_kind_set, atomic_kind_set, particle_set, &
3241 sab_orb, sac_ppl, nimages, cell_to_index,
"ORB")
3244 IF (
ASSOCIATED(sap_ppnl))
THEN
3245 eps_ppnl = dft_control%qs_control%eps_ppnl
3246 CALL build_core_ppnl(dbcsr_work_h, dbcsr_work_p, force, &
3247 virial, .true., use_virial, nder, &
3248 qs_kind_set, atomic_kind_set, particle_set, &
3249 sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index,
"ORB")
3251 DEALLOCATE (dbcsr_work_p, dbcsr_work_h)
3253 IF (use_virial)
THEN
3255 virial%pv_xc = 0.0_dp
3256 NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
3257 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, &
3258 dummy_real1, dummy_real2, h_stress)
3259 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
3260 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
3261 IF (.NOT. do_exx)
THEN
3263 virial%pv_exc = virial%pv_exc - virial%pv_xc
3264 virial%pv_virial = virial%pv_virial - virial%pv_xc
3267 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, dummy_real1, dummy_real2)
3269 do_tau =
ASSOCIATED(vtau_rspace)
3272 CALL integrate_v_core_rspace(vh_rspace, qs_env)
3276 CALL qs_rho_get(rho, rho_ao=rho_ao)
3277 DO ispin = 1, nspins
3278 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3281 CALL get_qs_env(qs_env, pw_env=pw_env)
3282 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
3283 poisson_env=poisson_env)
3284 CALL auxbas_pw_pool%create_pw(vhxc_rspace)
3286 IF (use_virial) pv_loc = virial%pv_virial
3290 DO ispin = 1, nspins
3292 CALL pw_transfer(vh_rspace, vhxc_rspace)
3293 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3294 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3295 qs_env=qs_env, calculate_forces=.true.)
3297 CALL pw_transfer(vxc_rspace(ispin), vhxc_rspace)
3298 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3299 hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3300 qs_env=qs_env, calculate_forces=.true.)
3302 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3303 hmat=scrm(ispin), pmat=matrix_p_mp2(ispin), &
3304 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3308 DO ispin = 1, nspins
3309 CALL pw_transfer(vh_rspace, vhxc_rspace)
3310 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
3311 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
3312 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3313 qs_env=qs_env, calculate_forces=.true.)
3315 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
3316 hmat=scrm(ispin), pmat=rho_ao(ispin), &
3317 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
3321 CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
3323 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3326 IF (dft_control%do_admm)
THEN
3327 CALL get_admm_env(admm_env, task_list_aux_fit=task_list_aux_fit, rho_aux_fit=rho_aux_fit, &
3328 matrix_s_aux_fit=matrix_s_aux_fit)
3329 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux)
3330 CALL dbcsr_allocate_matrix_set(scrm_admm, nspins)
3331 DO ispin = 1, nspins
3332 ALLOCATE (scrm_admm(ispin)%matrix)
3333 CALL dbcsr_create(scrm_admm(ispin)%matrix, template=matrix_s_aux_fit(1)%matrix)
3334 CALL dbcsr_copy(scrm_admm(ispin)%matrix, matrix_s_aux_fit(1)%matrix)
3335 CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3338 IF (use_virial) pv_loc = virial%pv_virial
3339 IF (.NOT. qs_env%admm_env%aux_exch_func == do_admm_aux_exch_func_none)
THEN
3340 DO ispin = 1, nspins
3342 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3343 hmat=scrm_admm(ispin), pmat=matrix_p_mp2_admm(ispin), &
3344 qs_env=qs_env, calculate_forces=.true., &
3345 basis_type=
"AUX_FIT", task_list_external=task_list_aux_fit)
3347 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3348 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), &
3349 hmat=scrm_admm(ispin), pmat=rho_ao_aux(ispin), &
3350 qs_env=qs_env, calculate_forces=.true., &
3351 basis_type=
"AUX_FIT", task_list_external=task_list_aux_fit)
3352 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3356 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3358 CALL tddft_hfx_matrix(scrm_admm, rho_ao_aux, qs_env, .false., .false.)
3361 CALL admm_projection_derivative(qs_env, scrm_admm, matrix_p_mp2)
3363 CALL admm_projection_derivative(qs_env, scrm_admm, rho_ao)
3368 xc_section => section_vals_get_subs_vals(qs_env%input,
"DFT%XC")
3369 hfx_section => section_vals_get_subs_vals(xc_section,
"HF")
3370 CALL section_vals_get(hfx_section, explicit=do_hfx)
3373 CALL section_vals_get(hfx_section, n_repetition=n_rep_hf)
3374 cpassert(n_rep_hf == 1)
3375 IF (use_virial) virial%pv_fock_4c = 0.0_dp
3379 IF (dft_control%do_admm)
THEN
3380 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3381 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3382 IF (x_data(1, 1)%do_hfx_ri)
THEN
3384 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3385 x_data(1, 1)%general_parameter%fraction, &
3386 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3387 use_virial=use_virial, resp_only=.true.)
3389 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3390 1, use_virial, resp_only=.true.)
3393 DO ispin = 1, nspins
3394 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3396 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3397 IF (x_data(1, 1)%do_hfx_ri)
THEN
3399 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3400 x_data(1, 1)%general_parameter%fraction, &
3401 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2, &
3402 use_virial=use_virial, resp_only=.true.)
3404 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3405 1, use_virial, resp_only=.true.)
3407 DO ispin = 1, nspins
3408 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
3413 IF (dft_control%do_admm)
THEN
3414 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3415 CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux, rho_ao_kp=dbcsr_work_p)
3416 DO ispin = 1, nspins
3417 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
3419 IF (x_data(1, 1)%do_hfx_ri)
THEN
3421 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3422 x_data(1, 1)%general_parameter%fraction, &
3423 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2_admm, &
3424 use_virial=use_virial, resp_only=.false.)
3426 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2_admm, hfx_section, para_env, &
3427 1, use_virial, resp_only=.false.)
3429 DO ispin = 1, nspins
3430 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, matrix_p_mp2_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
3433 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3434 IF (x_data(1, 1)%do_hfx_ri)
THEN
3436 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
3437 x_data(1, 1)%general_parameter%fraction, &
3438 rho_ao=dbcsr_work_p, rho_ao_resp=matrix_p_mp2, &
3439 use_virial=use_virial, resp_only=.false.)
3441 CALL derivatives_four_center(qs_env, dbcsr_work_p, matrix_p_mp2, hfx_section, para_env, &
3442 1, use_virial, resp_only=.false.)
3447 IF (use_virial)
THEN
3448 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
3449 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
3454 CALL qs_rho_get(rho, rho_ao=rho_ao)
3455 DO ispin = 1, nspins
3456 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
3464 ALLOCATE (current_density(nspins), current_mat_h(nspins), current_density_admm(nspins))
3466 DO ispin = 1, nspins
3467 IF (idens == 1)
THEN
3468 current_density(ispin)%matrix => matrix_p_f(ispin)%matrix
3469 current_mat_h(ispin)%matrix => scrm(ispin)%matrix
3470 IF (dft_control%do_admm) current_density_admm(ispin)%matrix => matrix_p_f_admm(ispin)%matrix
3472 current_density(ispin)%matrix => p_env%p1(ispin)%matrix
3473 current_mat_h(ispin)%matrix => matrix_hz(ispin)%matrix
3474 IF (dft_control%do_admm) current_density_admm(ispin)%matrix => p_env%p1_admm(ispin)%matrix
3479 ALLOCATE (rhoz_r(nspins), rhoz_g(nspins))
3480 DO ispin = 1, nspins
3481 CALL auxbas_pw_pool%create_pw(rhoz_r(ispin))
3482 CALL auxbas_pw_pool%create_pw(rhoz_g(ispin))
3484 CALL auxbas_pw_pool%create_pw(rhoz_tot_gspace)
3485 CALL auxbas_pw_pool%create_pw(zv_hartree_rspace)
3486 CALL auxbas_pw_pool%create_pw(zv_hartree_gspace)
3488 CALL pw_zero(rhoz_tot_gspace)
3489 DO ispin = 1, nspins
3490 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3491 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin))
3492 CALL pw_axpy(rhoz_g(ispin), rhoz_tot_gspace)
3495 IF (use_virial)
THEN
3497 CALL get_qs_env(qs_env, rho=rho)
3498 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
3500 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
3502 h_stress(:, :) = 0.0_dp
3503 CALL pw_poisson_solve(poisson_env, &
3504 density=rhoz_tot_gspace, &
3505 ehartree=ehartree, &
3506 vhartree=zv_hartree_gspace, &
3507 h_stress=h_stress, &
3508 aux_density=rho_tot_gspace)
3510 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
3513 virial%pv_ehartree = virial%pv_ehartree + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3514 virial%pv_virial = virial%pv_virial + 2.0_dp*h_stress/real(para_env%num_pe, dp)
3517 CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, &
3521 CALL pw_transfer(zv_hartree_gspace, zv_hartree_rspace)
3522 CALL pw_scale(zv_hartree_rspace, zv_hartree_rspace%pw_grid%dvol)
3523 CALL integrate_v_core_rspace(zv_hartree_rspace, qs_env)
3527 TYPE(pw_c1d_gs_type) :: tauz_g
3528 CALL auxbas_pw_pool%create_pw(tauz_g)
3529 ALLOCATE (tauz_r(nspins))
3530 DO ispin = 1, nspins
3531 CALL auxbas_pw_pool%create_pw(tauz_r(ispin))
3533 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3534 rho=tauz_r(ispin), rho_gspace=tauz_g, compute_tau=.true.)
3536 CALL auxbas_pw_pool%give_back_pw(tauz_g)
3541 IF (use_virial)
THEN
3544 DO ispin = 1, nspins
3545 exc = exc + pw_integral_ab(rhoz_r(ispin), vxc_rspace(ispin))/ &
3546 vxc_rspace(ispin)%pw_grid%dvol
3548 IF (
ASSOCIATED(vtau_rspace))
THEN
3549 DO ispin = 1, nspins
3550 exc = exc + pw_integral_ab(tauz_r(ispin), vtau_rspace(ispin))/ &
3551 vtau_rspace(ispin)%pw_grid%dvol
3555 virial%pv_ehartree(i, i) = virial%pv_ehartree(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp)
3556 virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3557 virial%pv_virial(i, i) = virial%pv_virial(i, i) - 4.0_dp*ehartree/real(para_env%num_pe, dp) &
3558 - exc/real(para_env%num_pe, dp)
3563 IF (dft_control%do_admm)
THEN
3564 CALL get_qs_env(qs_env, admm_env=admm_env)
3565 xc_section => admm_env%xc_section_primary
3567 xc_section => section_vals_get_subs_vals(qs_env%input,
"DFT%XC")
3570 IF (use_virial) virial%pv_xc = 0.0_dp
3572 CALL create_kernel(qs_env, &
3579 xc_section=xc_section, &
3580 compute_virial=use_virial, &
3581 virial_xc=virial%pv_xc)
3583 IF (use_virial)
THEN
3584 virial%pv_exc = virial%pv_exc + virial%pv_xc
3585 virial%pv_virial = virial%pv_virial + virial%pv_xc
3587 pv_loc = virial%pv_virial
3590 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3591 DO ispin = 1, nspins
3592 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3593 CALL pw_axpy(zv_hartree_rspace, v_xc(ispin))
3594 CALL integrate_v_rspace(qs_env=qs_env, &
3595 v_rspace=v_xc(ispin), &
3596 hmat=current_mat_h(ispin), &
3597 pmat=dbcsr_work_p(ispin, 1), &
3598 calculate_forces=.true.)
3599 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3601 CALL auxbas_pw_pool%give_back_pw(rhoz_tot_gspace)
3602 CALL auxbas_pw_pool%give_back_pw(zv_hartree_rspace)
3603 CALL auxbas_pw_pool%give_back_pw(zv_hartree_gspace)
3607 DO ispin = 1, nspins
3608 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3609 CALL integrate_v_rspace(qs_env=qs_env, &
3610 v_rspace=v_xc_tau(ispin), &
3611 hmat=current_mat_h(ispin), &
3612 pmat=dbcsr_work_p(ispin, 1), &
3613 compute_tau=.true., &
3614 calculate_forces=.true.)
3615 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3617 DEALLOCATE (v_xc_tau)
3620 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3623 IF (dft_control%do_admm)
THEN
3624 DO ispin = 1, nspins
3625 CALL dbcsr_set(scrm_admm(ispin)%matrix, 0.0_dp)
3627 CALL qs_rho_get(rho_aux_fit, tau_r_valid=do_tau_admm)
3629 IF (.NOT. admm_env%aux_exch_func == do_admm_aux_exch_func_none)
THEN
3630 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux_fit)
3631 DO ispin = 1, nspins
3632 CALL pw_zero(rhoz_r(ispin))
3633 CALL pw_zero(rhoz_g(ispin))
3634 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density_admm(ispin)%matrix, &
3635 rho=rhoz_r(ispin), rho_gspace=rhoz_g(ispin), &
3636 basis_type=
"AUX_FIT", task_list_external=task_list_aux_fit)
3639 IF (do_tau_admm)
THEN
3641 TYPE(pw_c1d_gs_type) :: tauz_g
3642 CALL auxbas_pw_pool%create_pw(tauz_g)
3643 DO ispin = 1, nspins
3644 CALL pw_zero(tauz_r(ispin))
3645 CALL calculate_rho_elec(ks_env=qs_env%ks_env, matrix_p=current_density(ispin)%matrix, &
3646 rho=tauz_r(ispin), rho_gspace=tauz_g, &
3647 basis_type=
"AUX_FIT", task_list_external=task_list_aux_fit, &
3650 CALL auxbas_pw_pool%give_back_pw(tauz_g)
3655 IF (use_virial)
THEN
3657 DO ispin = 1, nspins
3658 exc = exc + pw_integral_ab(rhoz_r(ispin), vadmm_rspace(ispin))/ &
3659 vadmm_rspace(ispin)%pw_grid%dvol
3662 virial%pv_exc(i, i) = virial%pv_exc(i, i) - exc/real(para_env%num_pe, dp)
3663 virial%pv_virial(i, i) = virial%pv_virial(i, i) - exc/real(para_env%num_pe, dp)
3666 virial%pv_xc = 0.0_dp
3669 xc_section => admm_env%xc_section_aux
3670 CALL create_kernel(qs_env, v_xc, v_xc_tau, rho_aux_fit, rhoz_r, rhoz_g, tauz_r, xc_section, &
3671 compute_virial=use_virial, virial_xc=virial%pv_xc)
3673 IF (use_virial)
THEN
3674 virial%pv_exc = virial%pv_exc + virial%pv_xc
3675 virial%pv_virial = virial%pv_virial + virial%pv_xc
3677 pv_loc = virial%pv_virial
3680 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=dbcsr_work_p)
3681 DO ispin = 1, nspins
3682 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
3683 CALL integrate_v_rspace(qs_env=qs_env, &
3684 v_rspace=v_xc(ispin), &
3685 hmat=scrm_admm(ispin), &
3686 pmat=dbcsr_work_p(ispin, 1), &
3687 calculate_forces=.true., &
3688 basis_type=
"AUX_FIT", &
3689 task_list_external=task_list_aux_fit)
3690 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
3694 IF (do_tau_admm)
THEN
3695 DO ispin = 1, nspins
3696 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
3697 CALL integrate_v_rspace(qs_env=qs_env, &
3698 v_rspace=v_xc_tau(ispin), &
3699 hmat=scrm_admm(ispin), &
3700 pmat=dbcsr_work_p(ispin, 1), &
3701 calculate_forces=.true., &
3702 basis_type=
"AUX_FIT", &
3703 task_list_external=task_list_aux_fit, &
3705 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
3707 DEALLOCATE (v_xc_tau)
3710 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
3713 CALL tddft_hfx_matrix(scrm_admm, current_density_admm, qs_env, .false., .false.)
3715 CALL qs_rho_get(rho, rho_ao_kp=dbcsr_work_p)
3716 CALL admm_projection_derivative(qs_env, scrm_admm, dbcsr_work_p(:, 1))
3719 CALL dbcsr_create(dbcsr_work, template=matrix_s(1)%matrix)
3720 IF (idens == 2)
THEN
3721 nao = admm_env%nao_orb
3722 nao_aux = admm_env%nao_aux_fit
3723 DO ispin = 1, nspins
3724 CALL dbcsr_copy(dbcsr_work, matrix_hz(ispin)%matrix)
3725 CALL dbcsr_set(dbcsr_work, 0.0_dp)
3727 CALL cp_dbcsr_sm_fm_multiply(scrm_admm(ispin)%matrix, admm_env%A, &
3728 admm_env%work_aux_orb, nao)
3729 CALL parallel_gemm(
'T',
'N', nao, nao, nao_aux, &
3730 1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
3731 admm_env%work_orb_orb)
3732 CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, dbcsr_work, keep_sparsity=.true.)
3733 CALL dbcsr_add(matrix_hz(ispin)%matrix, dbcsr_work, 1.0_dp, 1.0_dp)
3737 CALL dbcsr_release(dbcsr_work)
3741 IF (idens == 2)
THEN
3742 CALL tddft_hfx_matrix(matrix_hz, current_density, qs_env, .false., .false.)
3747 DO ispin = 1, nspins
3748 CALL auxbas_pw_pool%give_back_pw(rhoz_r(ispin))
3749 CALL auxbas_pw_pool%give_back_pw(rhoz_g(ispin))
3751 DEALLOCATE (rhoz_r, rhoz_g)
3754 DO ispin = 1, nspins
3755 CALL auxbas_pw_pool%give_back_pw(tauz_r(ispin))
3760 CALL dbcsr_deallocate_matrix_set(scrm_admm)
3762 DEALLOCATE (current_density, current_mat_h, current_density_admm)
3763 CALL dbcsr_deallocate_matrix_set(scrm)
3767 IF (nspins == 2) focc = 1.0_dp
3768 CALL get_qs_env(qs_env, mos=mos)
3769 DO ispin = 1, nspins
3770 CALL get_mo_set(mo_set=mos(ispin), homo=nocc)
3771 CALL calculate_whz_matrix(mos(ispin)%mo_coeff, matrix_hz(ispin)%matrix, &
3772 p_env%w1(ispin)%matrix, focc, nocc)
3774 IF (nspins == 2)
CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, 1.0_dp)
3777 IF (.NOT. do_exx)
THEN
3778 CALL qs_energies_compute_w(qs_env, calc_forces=.true.)
3779 CALL get_qs_env(qs_env, matrix_w=matrix_w)
3780 CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, 1.0_dp)
3781 IF (nspins == 2)
CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, 1.0_dp)
3785 CALL build_overlap_matrix(qs_env%ks_env, matrix_s=scrm, &
3786 matrix_name=
"OVERLAP MATRIX", &
3787 basis_type_a=
"ORB", basis_type_b=
"ORB", &
3788 sab_nl=sab_orb, calculate_forces=.true., &
3789 matrix_p=p_env%w1(1)%matrix)
3791 IF (.NOT. do_exx)
THEN
3792 CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(1)%matrix, 1.0_dp, -1.0_dp)
3793 IF (nspins == 2)
CALL dbcsr_add(p_env%w1(1)%matrix, matrix_w(2)%matrix, 1.0_dp, -1.0_dp)
3794 DO ispin = 1, nspins
3795 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp)
3799 IF (nspins == 2)
CALL dbcsr_add(p_env%w1(1)%matrix, p_env%w1(2)%matrix, 1.0_dp, -1.0_dp)
3800 CALL dbcsr_deallocate_matrix_set(scrm)
3802 IF (use_virial) virial%pv_calculate = .false.
3805 CALL auxbas_pw_pool%give_back_pw(vh_rspace)
3807 DO ispin = 1, nspins
3808 CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
3809 IF (
ASSOCIATED(vtau_rspace))
THEN
3810 CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
3812 IF (
ASSOCIATED(vadmm_rspace))
THEN
3813 CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
3816 DEALLOCATE (vxc_rspace)
3817 IF (
ASSOCIATED(vtau_rspace))
DEALLOCATE (vtau_rspace)
3818 IF (
ASSOCIATED(vadmm_rspace))
DEALLOCATE (vadmm_rspace)
3820 CALL timestop(handle)
3822 END SUBROUTINE update_im_time_forces
3835 SUBROUTINE build_y_matrix(Y, A, P, R, filter_eps)
3837 TYPE(dbcsr_type),
INTENT(OUT) :: y
3838 TYPE(dbcsr_type),
INTENT(INOUT) :: a, p, r
3839 REAL(dp),
INTENT(IN) :: filter_eps
3841 CHARACTER(LEN=*),
PARAMETER :: routinen =
'build_Y_matrix'
3843 INTEGER :: handle, k, n
3844 REAL(dp) :: norm_scalar, threshold
3845 TYPE(dbcsr_type) :: a2n, exp_a2n, prn, work, work2, yk
3847 CALL timeset(routinen, handle)
3849 threshold = 1.0e-16_dp
3852 norm_scalar = dbcsr_frobenius_norm(a)
3857 IF ((norm_scalar/2.0_dp**n) < 1.0_dp)
EXIT
3862 CALL dbcsr_create(prn, template=a, matrix_type=dbcsr_type_no_symmetry)
3863 CALL dbcsr_create(work, template=a, matrix_type=dbcsr_type_no_symmetry)
3864 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, p, r, 0.0_dp, work, filter_eps=filter_eps)
3865 CALL dbcsr_create(exp_a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3868 CALL matrix_exponential(exp_a2n, a, 1.0_dp, 0.5_dp**k, threshold)
3869 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, exp_a2n, work, 0.0_dp, prn, filter_eps=filter_eps)
3870 CALL dbcsr_multiply(
'N',
'N', 1.0_dp, work, exp_a2n, 1.0_dp, prn, filter_eps=filter_eps)
3871 CALL dbcsr_copy(work, prn)
3873 CALL dbcsr_release(exp_a2n)
3876 CALL dbcsr_create(a2n, template=a, matrix_type=dbcsr_type_no_symmetry)
3877 CALL dbcsr_copy(a2n, a)
3878 CALL dbcsr_scale(a2n, 0.5_dp**n)
3879 CALL dbcsr_create(y, template=a, matrix_type=dbcsr_type_no_symmetry)
3880 CALL dbcsr_create(yk, template=a, matrix_type=dbcsr_type_no_symmetry)
3881 CALL dbcsr_create(work2, template=a, matrix_type=dbcsr_type_no_symmetry)
3884 CALL dbcsr_scale(prn, 0.5_dp**n)
3885 CALL dbcsr_copy(work, prn)
3886 CALL dbcsr_copy(work2, prn)
3887 CALL dbcsr_add(y, prn, 1.0_dp, 1.0_dp)
3892 CALL dbcsr_multiply(
'N',
'N', 1.0_dp/real(k, dp), a2n, work, 0.0_dp, yk, filter_eps=filter_eps)
3893 CALL dbcsr_multiply(
'N',
'N', 1.0_dp/real(k, dp), work2, a2n, 0.0_dp, prn, filter_eps=filter_eps)
3895 CALL dbcsr_add(yk, prn, 1.0_dp, 1.0_dp)
3896 CALL dbcsr_add(y, yk, 1.0_dp, 1.0_dp)
3898 IF (dbcsr_frobenius_norm(yk) < threshold)
EXIT
3899 CALL dbcsr_copy(work, yk)
3900 CALL dbcsr_copy(work2, prn)
3903 CALL dbcsr_release(work)
3904 CALL dbcsr_release(work2)
3905 CALL dbcsr_release(prn)
3906 CALL dbcsr_release(a2n)
3907 CALL dbcsr_release(yk)
3909 CALL timestop(handle)
3911 END SUBROUTINE build_y_matrix
3927 SUBROUTINE keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, &
3928 do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
3930 REAL(dp),
ALLOCATABLE,
DIMENSION(:),
INTENT(INOUT) :: tj, wj, tau_tj, tau_wj
3931 REAL(dp),
ALLOCATABLE,
DIMENSION(:, :), &
3932 INTENT(INOUT) :: weights_cos_tf_t_to_w, &
3933 weights_cos_tf_w_to_t
3934 LOGICAL,
INTENT(IN) :: do_laplace, do_im_time
3935 INTEGER,
INTENT(IN) :: num_integ_points, unit_nr
3936 TYPE(qs_environment_type),
POINTER :: qs_env
3940 IF (do_laplace .OR. do_im_time)
THEN
3941 IF (.NOT.
ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tau_tj))
THEN
3942 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tau_tj(0:num_integ_points))
3943 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tau_wj(num_integ_points))
3944 qs_env%mp2_env%ri_rpa_im_time%tau_tj(:) = tau_tj(:)
3945 qs_env%mp2_env%ri_rpa_im_time%tau_wj(:) = tau_wj(:)
3948 tau_tj(:) = qs_env%mp2_env%ri_rpa_im_time%tau_tj(:)
3949 tau_wj(:) = qs_env%mp2_env%ri_rpa_im_time%tau_wj(:)
3952 IF (.NOT. do_laplace)
THEN
3953 IF (.NOT.
ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tj))
THEN
3954 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%tj(num_integ_points))
3955 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%wj(num_integ_points))
3956 qs_env%mp2_env%ri_rpa_im_time%tj(:) = tj(:)
3957 qs_env%mp2_env%ri_rpa_im_time%wj(:) = wj(:)
3958 IF (do_im_time)
THEN
3959 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(num_integ_points, num_integ_points))
3960 ALLOCATE (qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(num_integ_points, num_integ_points))
3961 qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(:, :) = weights_cos_tf_t_to_w(:, :)
3962 qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(:, :) = weights_cos_tf_w_to_t(:, :)
3965 tj(:) = qs_env%mp2_env%ri_rpa_im_time%tj(:)
3966 wj(:) = qs_env%mp2_env%ri_rpa_im_time%wj(:)
3967 IF (do_im_time)
THEN
3968 weights_cos_tf_t_to_w(:, :) = qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_t_to_w(:, :)
3969 weights_cos_tf_w_to_t(:, :) = qs_env%mp2_env%ri_rpa_im_time%weights_cos_tf_w_to_t(:, :)
3973 IF (unit_nr > 0)
THEN
3975 IF (
ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tj) .AND. (.NOT. do_laplace))
THEN
3976 WRITE (unit=unit_nr, fmt=
"(T3,A,T75,i6)") &
3977 "MINIMAX_INFO| Number of integration points:", num_integ_points
3978 WRITE (unit=unit_nr, fmt=
"(T3,A,T54,A,T72,A)") &
3979 "MINIMAX_INFO| Minimax params (freq grid, scaled):",
"Weights",
"Abscissas"
3980 DO jquad = 1, num_integ_points
3981 WRITE (unit=unit_nr, fmt=
"(T41,F20.10,F20.10)") wj(jquad), tj(jquad)
3983 CALL m_flush(unit_nr)
3985 IF (
ASSOCIATED(qs_env%mp2_env%ri_rpa_im_time%tau_tj))
THEN
3986 WRITE (unit=unit_nr, fmt=
"(T3,A,T75,i6)") &
3987 "MINIMAX_INFO| Number of integration points:", num_integ_points
3988 WRITE (unit=unit_nr, fmt=
"(T3,A,T54,A,T72,A)") &
3989 "MINIMAX_INFO| Minimax params (time grid, scaled):",
"Weights",
"Abscissas"
3990 DO jquad = 1, num_integ_points
3991 WRITE (unit=unit_nr, fmt=
"(T41,F20.10,F20.10)") tau_wj(jquad), tau_tj(jquad)
3993 CALL m_flush(unit_nr)
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
static GRID_HOST_DEVICE int ncoset(const int l)
Number of Cartesian orbitals up to given angular momentum quantum.
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Contains ADMM methods which require molecular orbitals.
subroutine, public admm_projection_derivative(qs_env, matrix_hz, matrix_pz, fval)
Calculate derivatives terms from overlap matrices.
Types and set/get functions for auxiliary density matrix methods.
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
All kind of helpful little routines.
real(kind=dp) function, public exp_radius_very_extended(la_min, la_max, lb_min, lb_max, pab, o1, o2, ra, rb, rp, zetp, eps, prefactor, cutoff, epsabs)
computes the radius of the Gaussian outside of which it is smaller than eps
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public bussy2023
Handles all functions related to the CELL.
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
...
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltaR)
...
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
methods related to the blacs parallel environment
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_cholesky_decompose(matrix, n, para_env, blacs_env)
used to replace a symmetric positive def. matrix M with its cholesky decomposition U: M = U^T * U,...
subroutine, public cp_dbcsr_cholesky_invert(matrix, n, para_env, blacs_env, upper_to_full)
used to replace the cholesky decomposition by the inverse
Interface to (sca)lapack for the Cholesky based procedures.
subroutine, public cp_dbcsr_power(matrix, exponent, threshold, n_dependent, para_env, blacs_env, verbose, eigenvectors, eigenvalues)
...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public cp_dbcsr_dist2d_to_dist(dist2d, dist)
Creates a DBCSR distribution from a distribution_2d.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
subroutine, public cp_eri_mme_update_local_counts(param, para_env, G_count_2c, R_count_2c, GG_count_3c, GR_count_3c, RR_count_3c)
Update local counters to gather statistics on different paths taken in MME algorithm (each Ewald sum ...
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
This is the start of a dbt_api, all publically needed functions are exported here....
stores a mapping of 2D info (e.g. matrix) on a 2D processor distribution (i.e. blacs grid) where cpus...
Routines used for Harris functional Kohn-Sham calculation.
subroutine, public create_kernel(qs_env, vxc, vxc_tau, rho, rho1_r, rho1_g, tau1_r, xc_section, compute_virial, virial_xc)
Creation of second derivative xc-potential.
integer function, public gaussian_gridlevel(gridlevel_info, exponent)
...
Utilities for hfx and admm methods.
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data)
Add the hfx contributions to the Hamiltonian.
Routines to calculate derivatives with respect to basis function origin.
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Routines to calculate EXX in RPA and energy correction methods.
subroutine, public add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian.
subroutine, public get_force_from_3c_trace(force, t_3c_contr, t_3c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, deriv_dim)
This routines calculates the force contribution from a trace over 3D tensors, i.e....
subroutine, public get_idx_to_atom(idx_to_at, bsizes_split, bsizes_orig)
a small utility function that returns the atom corresponding to a block of a split tensor
subroutine, public hfx_ri_update_forces(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, mos, use_virial, resp_only, rescale_factor)
the general routine that calls the relevant force code
subroutine, public get_2c_der_force(force, t_2c_contr, t_2c_der, atom_of_kind, kind_of, idx_to_at, pref, do_mp2, do_ovlp)
Update the forces due to the derivative of the a 2-center product d/dR (Q|R)
Types and set/get functions for HFX.
subroutine, public dealloc_containers(DATA, memory_usage)
...
subroutine, public alloc_containers(DATA, bin_size)
...
Routines useful for iterative matrix calculations.
subroutine, public matrix_exponential(matrix_exp, matrix, omega, alpha, threshold)
...
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Machine interface based on Fortran 2003 and POSIX.
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition of mathematical constants and functions.
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
Routines to calculate 2c- and 3c-integrals for RI with GPW.
subroutine, public cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
Cleanup GPW integration for RI-MP2/RI-RPA.
subroutine, public virial_gpw_potential(rho_g_copy, pot_g, rho_g, dvg, h_stress, potential_parameter, para_env_sub)
Calculates stress tensor contribution from the operator.
subroutine, public prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
Prepares GPW calculation for RI-MP2/RI-RPA.
subroutine, public calc_potential_gpw(pot_r, rho_g, poisson_env, pot_g, potential_parameter, dvg, no_transfer)
Calculates potential from a given density in g-space.
Interface to direct methods for electron repulsion integrals for MP2.
subroutine, public integrate_set_2c(param, potential_parameter, la_min, la_max, lb_min, lb_max, npgfa, npgfb, zeta, zetb, ra, rb, hab, n_hab_a, n_hab_b, offset_hab_a, offset_hab_b, offset_set_a, offset_set_b, sphi_a, sphi_b, sgfa, sgfb, nsgfa, nsgfb, eri_method, pab, force_a, force_b, hdab, hadb, G_count, R_count, do_reflection_a, do_reflection_b)
Integrate set pair and contract with sphi matrix.
Types needed for MP2 calculations.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
basic linear algebra operations for full matrixes
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Define the data structure for the particle information.
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
collects routines that calculate density matrices
subroutine, public calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbers....
Utility subroutine for qs energy calculation.
subroutine, public qs_energies_compute_w(qs_env, calc_forces)
Refactoring of qs_energies_scf. Moves computation of matrix_w into separate subroutine.
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Integrate single or product functions over a potential on a RS grid.
Calculate the interaction radii for the operator matrix calculation.
subroutine, public init_interaction_radii_orb_basis(orb_basis_set, eps_pgf_orb, eps_pgf_short)
...
Define the quickstep kind type and their sub types.
Calculation of kinetic energy matrix and forces.
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Calculate the KS reference potentials.
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Type definitiona for linear response calculations.
Definition and initialisation of the mo data type.
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Define the neighbor list data types and the corresponding functionality.
subroutine, public release_neighbor_list_sets(nlists)
releases an array of neighbor_list_sets
Calculation of overlap matrix, its derivatives and forces.
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Utility functions for the perturbation calculations.
subroutine, public p_env_psi0_changed(p_env, qs_env)
To be called after the value of psi0 has changed. Recalculates the quantities S_psi0 and m_epsilon.
subroutine, public p_env_create(p_env, qs_env, p1_option, p1_admm_option, orthogonal_orbitals, linres_control)
allocates and initializes the perturbation environment (no setup)
basis types for the calculation of the perturbation of density theory.
subroutine, public p_env_release(p_env)
relases the given p_env (see doc/ReferenceCounting.html)
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Utility methods to build 3-center integral tensors of various types.
subroutine, public distribution_3d_create(dist_3d, dist1, dist2, dist3, nkind, particle_set, mp_comm_3d, own_comm)
Create a 3d distribution.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
subroutine, public create_tensor_batches(sizes, nbatches, starts_array, ends_array, starts_array_block, ends_array_block)
...
subroutine, public create_3c_tensor(t3c, dist_1, dist_2, dist_3, pgrid, sizes_1, sizes_2, sizes_3, map1, map2, name)
...
Utility methods to build 3-center integral tensors of various types.
subroutine, public calc_3c_virial(work_virial, t3c_trace, pref, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos)
Calculates the 3c virial contributions on the fly.
subroutine, public build_2c_neighbor_lists(ij_list, basis_i, basis_j, potential_parameter, name, qs_env, sym_ij, molecular, dist_2d, pot_to_rad)
Build 2-center neighborlists adapted to different operators This mainly wraps build_neighbor_lists fo...
subroutine, public build_3c_derivatives(t3c_der_i, t3c_der_k, filter_eps, qs_env, nl_3c, basis_i, basis_j, basis_k, potential_parameter, der_eps, op_pos, do_kpoints, do_hfx_kpoints, bounds_i, bounds_j, bounds_k, RI_range, img_to_RI_cell)
Build 3-center derivative tensors.
subroutine, public compress_tensor(tensor, blk_indices, compressed, eps, memory)
...
subroutine, public neighbor_list_3c_destroy(ijk_list)
Destroy 3c neighborlist.
subroutine, public calc_2c_virial(work_virial, t2c_trace, pref, qs_env, nl_2c, basis_i, basis_j, potential_parameter)
Calculates the virial coming from 2c derivatives on the fly.
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
subroutine, public build_2c_derivatives(t2c_der, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints)
Calculates the derivatives of 2-center integrals, wrt to the first center.
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
subroutine, public build_3c_neighbor_lists(ijk_list, basis_i, basis_j, basis_k, dist_3d, potential_parameter, name, qs_env, sym_ij, sym_jk, sym_ik, molecular, op_pos, own_dist)
Build a 3-center neighbor list.
subroutine, public build_2c_integrals(t2c, filter_eps, qs_env, nl_2c, basis_i, basis_j, potential_parameter, do_kpoints, do_hfx_kpoints, ext_kpoints, regularization_RI)
...
pure logical function, public map_gaussian_here(rs_grid, h_inv, ra, offset, group_size, my_pos)
...
Calculate the CPKS equation and the resulting forces.
subroutine, public response_equation_new(qs_env, p_env, cpmos, iounit)
Initializes vectors for MO-coefficient based linear response solver and calculates response density,...
Routines needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
subroutine, public calc_rpa_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, Eigenval, e_fermi, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, tj, wj, tau_tj, cut_memory, ispin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling RPA contribution to the forces at each quadrature point....
subroutine, public calc_laplace_loop_forces(force_data, mat_P_omega, t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, num_integ_points, nmo, Eigenval, tau_tj, tau_wj, cut_memory, Pspin, Qspin, open_shell, unit_nr, dbcsr_time, dbcsr_nflop, mp2_env, qs_env)
Updates the cubic-scaling SOS-Laplace-MP2 contribution to the forces at each quadrature point.
subroutine, public calc_post_loop_forces(force_data, unit_nr, qs_env)
All the forces that can be calculated after the loop on the Laplace quaradture, using terms collected...
subroutine, public init_im_time_forces(force_data, fm_matrix_PQ, t_3c_M, unit_nr, mp2_env, qs_env)
Initializes and pre-calculates all needed tensors for the forces.
subroutine, public keep_initial_quad(tj, wj, tau_tj, tau_wj, weights_cos_tf_t_to_w, weights_cos_tf_w_to_t, do_laplace, do_im_time, num_integ_points, unit_nr, qs_env)
Overwrites the "optimal" Laplace quadrature with that of the first step.
Types needed for cubic-scaling RPA and SOS-Laplace-MP2 forces.
Routines for low-scaling RPA/GW with imaginary time.
subroutine, public compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, matrix_s, ispin, Eigenval, e_fermi, eps_filter, memory_info, unit_nr, jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, num_cells_dm, index_to_cell_dm, para_env)
...
Transfers densities from PW to RS grids and potentials from PW to RS.
subroutine, public potential_pw2rs(rs_v, v_rspace, pw_env)
transfers a potential from a pw_grid to a vector of realspace multigrids