34 dbt_batched_contract_finalize, dbt_batched_contract_init, dbt_contract, dbt_copy, &
35 dbt_copy_matrix_to_tensor, dbt_copy_tensor_to_matrix, dbt_create, dbt_destroy, dbt_filter, &
36 dbt_get_info, dbt_nblks_total, dbt_nd_mp_comm, dbt_pgrid_destroy, dbt_pgrid_type, dbt_type
61#include "./base/base_uses.f90"
67 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'rpa_im_time'
122 fm_scaled_dm_virt_tau, fm_mo_coeff_occ, fm_mo_coeff_virt, &
123 fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled, &
127 t_3c_M, t_3c_O, t_3c_O_compressed, t_3c_O_ind, &
128 starts_array_mc, ends_array_mc, &
129 starts_array_mc_block, ends_array_mc_block, &
130 weights_cos_tf_t_to_w, &
131 tj, tau_tj, e_fermi, eps_filter, &
132 alpha, eps_filter_im_time, Eigenval, nmo, &
133 num_integ_points, cut_memory, unit_nr, &
135 qs_env, do_kpoints_from_Gamma, &
136 index_to_cell_3c, cell_to_index_3c, &
137 has_mat_P_blocks, do_ri_sos_laplace_mp2, &
138 dbcsr_time, dbcsr_nflop)
139 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
INTENT(IN) :: mat_p_omega
140 TYPE(
cp_fm_type),
INTENT(IN) :: fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, &
141 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, fm_mo_coeff_virt_scaled
144 INTEGER,
INTENT(IN) :: ispin
145 TYPE(dbt_type),
INTENT(INOUT) :: t_3c_m
146 TYPE(dbt_type),
DIMENSION(:, :),
INTENT(INOUT) :: t_3c_o
148 INTENT(INOUT) :: t_3c_o_compressed
150 INTENT(INOUT) :: t_3c_o_ind
151 INTEGER,
DIMENSION(:),
INTENT(IN) :: starts_array_mc, ends_array_mc, &
152 starts_array_mc_block, &
154 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :), &
155 INTENT(IN) :: weights_cos_tf_t_to_w
156 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:), &
158 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: tau_tj
159 REAL(kind=
dp),
INTENT(IN) :: e_fermi, eps_filter, alpha, &
161 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
162 INTEGER,
INTENT(IN) :: nmo, num_integ_points, cut_memory, &
167 LOGICAL,
INTENT(IN) :: do_kpoints_from_gamma
168 INTEGER,
ALLOCATABLE,
DIMENSION(:, :),
INTENT(IN) :: index_to_cell_3c
169 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :), &
170 INTENT(IN) :: cell_to_index_3c
171 LOGICAL,
DIMENSION(:, :, :, :, :),
INTENT(INOUT) :: has_mat_p_blocks
172 LOGICAL,
INTENT(IN) :: do_ri_sos_laplace_mp2
173 REAL(
dp),
INTENT(INOUT) :: dbcsr_time
174 INTEGER(int_8),
INTENT(INOUT) :: dbcsr_nflop
176 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_mat_P_omega'
178 INTEGER :: comm_2d_handle, handle, handle2, handle3, i, i_cell, i_cell_r_1, &
179 i_cell_r_1_minus_s, i_cell_r_1_minus_t, i_cell_r_2, i_cell_r_2_minus_s_minus_t, i_cell_s, &
180 i_cell_t, i_mem, iquad, j, j_mem, jquad, num_3c_repl, num_cells_dm, unit_nr_dbcsr
181 INTEGER(int_8) :: nze, nze_dm_occ, nze_dm_virt, nze_m_occ, &
183 INTEGER(KIND=int_8) :: flops_1_occ, flops_1_virt, flops_2
184 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: dist_1, dist_2, mc_ranges, size_dm, &
186 INTEGER,
DIMENSION(2) :: pdims_2d
187 INTEGER,
DIMENSION(2, 1) :: ibounds_2, jbounds_2
188 INTEGER,
DIMENSION(2, 2) :: ibounds_1, jbounds_1
189 INTEGER,
DIMENSION(3) :: bounds_3c
190 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell_dm
191 LOGICAL :: do_gamma_rpa, do_kpoints_cubic_rpa, first_cycle_im_time, first_cycle_omega_loop, &
192 memory_info, r_1_minus_s_needed, r_1_minus_t_needed, r_2_minus_s_minus_t_needed
193 REAL(
dp) :: occ, occ_dm_occ, occ_dm_virt, occ_m_occ, &
194 occ_m_virt, occ_o, t1_flop
195 REAL(kind=
dp) :: omega, omega_old, t1, t2, tau, weight, &
198 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_occ_global, mat_dm_virt_global
199 TYPE(dbt_pgrid_type) :: pgrid_2d
200 TYPE(dbt_type) :: t_3c_m_occ, t_3c_m_occ_tmp, t_3c_m_virt, &
201 t_3c_m_virt_tmp, t_dm, t_dm_tmp, t_p, &
203 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:) :: t_dm_occ, t_dm_virt
204 TYPE(dbt_type),
ALLOCATABLE,
DIMENSION(:, :) :: t_3c_o_occ, t_3c_o_virt
207 CALL timeset(routinen, handle)
209 memory_info = mp2_env%ri_rpa_im_time%memory_info
210 IF (memory_info)
THEN
211 unit_nr_dbcsr = unit_nr
216 do_kpoints_cubic_rpa = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
217 do_gamma_rpa = .NOT. do_kpoints_cubic_rpa
218 num_3c_repl = maxval(cell_to_index_3c)
220 first_cycle_im_time = .true.
221 ALLOCATE (t_3c_o_occ(
SIZE(t_3c_o, 1),
SIZE(t_3c_o, 2)), t_3c_o_virt(
SIZE(t_3c_o, 1),
SIZE(t_3c_o, 2)))
222 DO i = 1,
SIZE(t_3c_o, 1)
223 DO j = 1,
SIZE(t_3c_o, 2)
224 CALL dbt_create(t_3c_o(i, j), t_3c_o_occ(i, j))
225 CALL dbt_create(t_3c_o(i, j), t_3c_o_virt(i, j))
229 CALL dbt_create(t_3c_m, t_3c_m_occ, name=
"M occ (RI | AO AO)")
230 CALL dbt_create(t_3c_m, t_3c_m_virt, name=
"M virt (RI | AO AO)")
232 ALLOCATE (mc_ranges(cut_memory + 1))
233 mc_ranges(:cut_memory) = starts_array_mc_block(:)
234 mc_ranges(cut_memory + 1) = ends_array_mc_block(cut_memory) + 1
236 DO jquad = 1, num_integ_points
241 CALL compute_mat_dm_global(fm_scaled_dm_occ_tau, fm_scaled_dm_virt_tau, tau_tj, num_integ_points, nmo, &
242 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
243 fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, &
245 eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
246 jquad, do_kpoints_cubic_rpa, do_kpoints_from_gamma, qs_env, &
247 num_cells_dm, index_to_cell_dm, para_env)
249 ALLOCATE (t_dm_virt(num_cells_dm))
250 ALLOCATE (t_dm_occ(num_cells_dm))
253 CALL comm_2d%set_handle(comm_2d_handle)
255 pgrid_2d = dbt_nd_mp_comm(comm_2d, [1], [2], pdims_2d=pdims_2d)
256 ALLOCATE (size_p(dbt_nblks_total(t_3c_m, 1)))
257 CALL dbt_get_info(t_3c_m, blk_size_1=size_p)
259 ALLOCATE (size_dm(dbt_nblks_total(t_3c_o(1, 1), 3)))
260 CALL dbt_get_info(t_3c_o(1, 1), blk_size_3=size_dm)
261 CALL create_2c_tensor(t_dm, dist_1, dist_2, pgrid_2d, size_dm, size_dm, name=
"D (AO | AO)")
263 DEALLOCATE (dist_1, dist_2)
264 CALL create_2c_tensor(t_p, dist_1, dist_2, pgrid_2d, size_p, size_p, name=
"P (RI | RI)")
266 DEALLOCATE (dist_1, dist_2)
267 CALL dbt_pgrid_destroy(pgrid_2d)
269 DO i_cell = 1, num_cells_dm
270 CALL dbt_create(t_dm, t_dm_virt(i_cell), name=
"D virt (AO | AO)")
271 CALL dbt_create(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
272 CALL dbt_copy_matrix_to_tensor(mat_dm_virt_global(jquad, i_cell)%matrix, t_dm_tmp)
273 CALL dbt_copy(t_dm_tmp, t_dm_virt(i_cell), move_data=.true.)
274 CALL dbcsr_clear(mat_dm_virt_global(jquad, i_cell)%matrix)
276 CALL dbt_create(t_dm, t_dm_occ(i_cell), name=
"D occ (AO | AO)")
277 CALL dbt_copy_matrix_to_tensor(mat_dm_occ_global(jquad, i_cell)%matrix, t_dm_tmp)
278 CALL dbt_copy(t_dm_tmp, t_dm_occ(i_cell), move_data=.true.)
279 CALL dbt_destroy(t_dm_tmp)
280 CALL dbcsr_clear(mat_dm_occ_global(jquad, i_cell)%matrix)
286 CALL dbt_destroy(t_dm)
288 CALL dbt_create(t_3c_o_occ(1, 1), t_3c_m_occ_tmp, name=
"M (RI AO | AO)")
289 CALL dbt_create(t_3c_o_virt(1, 1), t_3c_m_virt_tmp, name=
"M (RI AO | AO)")
291 CALL timeset(routinen//
"_contract", handle2)
296 DO i = 1,
SIZE(t_3c_o_occ, 1)
297 DO j = 1,
SIZE(t_3c_o_occ, 2)
298 CALL dbt_batched_contract_init(t_3c_o_occ(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
301 DO i = 1,
SIZE(t_3c_o_virt, 1)
302 DO j = 1,
SIZE(t_3c_o_virt, 2)
303 CALL dbt_batched_contract_init(t_3c_o_virt(i, j), batch_range_2=mc_ranges, batch_range_3=mc_ranges)
306 CALL dbt_batched_contract_init(t_3c_m_occ_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
307 CALL dbt_batched_contract_init(t_3c_m_virt_tmp, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
308 CALL dbt_batched_contract_init(t_3c_m_occ, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
309 CALL dbt_batched_contract_init(t_3c_m_virt, batch_range_2=mc_ranges, batch_range_3=mc_ranges)
311 DO i_cell_t = 1, num_cells_dm/2 + 1
313 IF (.NOT. any(has_mat_p_blocks(i_cell_t, :, :, :, :))) cycle
315 CALL dbt_batched_contract_init(t_p)
317 IF (do_gamma_rpa)
THEN
326 DO j_mem = 1, cut_memory
328 CALL dbt_get_info(t_3c_o_occ(1, 1), nfull_total=bounds_3c)
330 jbounds_1(:, 1) = [1, bounds_3c(1)]
331 jbounds_1(:, 2) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
333 jbounds_2(:, 1) = [starts_array_mc(j_mem), ends_array_mc(j_mem)]
335 IF (do_gamma_rpa)
CALL dbt_batched_contract_init(t_dm_virt(1))
337 DO i_mem = 1, cut_memory
339 IF (.NOT. any(has_mat_p_blocks(i_cell_t, i_mem, j_mem, :, :))) cycle
341 ibounds_1(:, 1) = [1, bounds_3c(1)]
342 ibounds_1(:, 2) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
344 ibounds_2(:, 1) = [starts_array_mc(i_mem), ends_array_mc(i_mem)]
346 IF (unit_nr_dbcsr > 0)
WRITE (unit=unit_nr_dbcsr, fmt=
"(T3,A,I3,1X,I3)") &
347 "RPA_LOW_SCALING_INFO| Memory Cut iteration", i_mem, j_mem
349 DO i_cell_r_1 = 1, num_3c_repl
351 DO i_cell_r_2 = 1, num_3c_repl
353 IF (.NOT. has_mat_p_blocks(i_cell_t, i_mem, j_mem, i_cell_r_1, i_cell_r_2)) cycle
355 CALL get_diff_index_3c(i_cell_r_1, i_cell_t, i_cell_r_1_minus_t, &
356 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
357 r_1_minus_t_needed, do_kpoints_cubic_rpa)
359 IF (do_gamma_rpa)
CALL dbt_batched_contract_init(t_dm_occ(1))
360 DO i_cell_s = 1, num_cells_dm
361 CALL get_diff_index_3c(i_cell_r_1, i_cell_s, i_cell_r_1_minus_s, index_to_cell_3c, &
362 cell_to_index_3c, index_to_cell_dm, r_1_minus_s_needed, &
363 do_kpoints_cubic_rpa)
364 IF (r_1_minus_s_needed)
THEN
366 CALL timeset(routinen//
"_calc_M_occ_t", handle3)
368 t_3c_o_ind(i_cell_r_1_minus_s, i_cell_r_2, j_mem)%ind, &
369 t_3c_o_compressed(i_cell_r_1_minus_s, i_cell_r_2, j_mem), &
370 qs_env%mp2_env%ri_rpa_im_time%eps_compress)
372 IF (do_gamma_rpa .AND. i_mem == 1)
THEN
378 CALL dbt_copy(t_3c_o(i_cell_r_1_minus_s, i_cell_r_2), &
379 t_3c_o_occ(i_cell_r_1_minus_s, i_cell_r_2), move_data=.true.)
381 CALL dbt_contract(alpha=1.0_dp, &
382 tensor_1=t_3c_o_occ(i_cell_r_1_minus_s, i_cell_r_2), &
383 tensor_2=t_dm_occ(i_cell_s), &
385 tensor_3=t_3c_m_occ_tmp, &
386 contract_1=[3], notcontract_1=[1, 2], &
387 contract_2=[2], notcontract_2=[1], &
388 map_1=[1, 2], map_2=[3], &
389 bounds_2=jbounds_1, bounds_3=ibounds_2, &
390 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
392 CALL timestop(handle3)
394 dbcsr_nflop = dbcsr_nflop + flops_1_occ
399 IF (do_gamma_rpa)
CALL dbt_batched_contract_finalize(t_dm_occ(1))
404 CALL timeset(routinen//
"_copy_M_occ_t", handle3)
405 CALL dbt_copy(t_3c_m_occ_tmp, t_3c_m_occ, order=[1, 3, 2], move_data=.true.)
406 CALL dbt_filter(t_3c_m_occ, eps_filter)
407 CALL timestop(handle3)
409 IF (do_gamma_rpa)
THEN
411 nze_m_occ = nze_m_occ + nze
412 occ_m_occ = occ_m_occ + occ
415 DO i_cell_s = 1, num_cells_dm
416 CALL get_diff_diff_index_3c(i_cell_r_2, i_cell_s, i_cell_t, i_cell_r_2_minus_s_minus_t, &
417 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
418 r_2_minus_s_minus_t_needed, do_kpoints_cubic_rpa)
420 IF (r_1_minus_t_needed .AND. r_2_minus_s_minus_t_needed)
THEN
422 t_3c_o_ind(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t, i_mem)%ind, &
423 t_3c_o_compressed(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t, i_mem), &
424 qs_env%mp2_env%ri_rpa_im_time%eps_compress)
426 CALL dbt_copy(t_3c_o(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), &
427 t_3c_o_virt(i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), move_data=.true.)
429 CALL timeset(routinen//
"_calc_M_virt_t", handle3)
430 CALL dbt_contract(alpha=alpha/2.0_dp, &
431 tensor_1=t_3c_o_virt( &
432 i_cell_r_2_minus_s_minus_t, i_cell_r_1_minus_t), &
433 tensor_2=t_dm_virt(i_cell_s), &
435 tensor_3=t_3c_m_virt_tmp, &
436 contract_1=[3], notcontract_1=[1, 2], &
437 contract_2=[2], notcontract_2=[1], &
438 map_1=[1, 2], map_2=[3], &
439 bounds_2=ibounds_1, bounds_3=jbounds_2, &
440 filter_eps=eps_filter, unit_nr=unit_nr_dbcsr, &
442 CALL timestop(handle3)
444 dbcsr_nflop = dbcsr_nflop + flops_1_virt
449 CALL timeset(routinen//
"_copy_M_virt_t", handle3)
450 CALL dbt_copy(t_3c_m_virt_tmp, t_3c_m_virt, move_data=.true.)
451 CALL dbt_filter(t_3c_m_virt, eps_filter)
452 CALL timestop(handle3)
454 IF (do_gamma_rpa)
THEN
456 nze_m_virt = nze_m_virt + nze
457 occ_m_virt = occ_m_virt + occ
462 CALL timeset(routinen//
"_calc_P_t", handle3)
464 CALL dbt_contract(alpha=1.0_dp, tensor_1=t_3c_m_occ, &
465 tensor_2=t_3c_m_virt, &
468 contract_1=[2, 3], notcontract_1=[1], &
469 contract_2=[2, 3], notcontract_2=[1], &
470 map_1=[1], map_2=[2], &
471 filter_eps=eps_filter_im_time/real(cut_memory**2, kind=
dp), &
474 unit_nr=unit_nr_dbcsr)
476 CALL timestop(handle3)
478 first_cycle_im_time = .false.
480 IF (jquad == 1 .AND. flops_2 == 0) &
481 has_mat_p_blocks(i_cell_t, i_mem, j_mem, i_cell_r_1, i_cell_r_2) = .false.
486 IF (do_gamma_rpa)
CALL dbt_batched_contract_finalize(t_dm_virt(1))
489 CALL dbt_batched_contract_finalize(t_p, unit_nr=unit_nr_dbcsr)
491 CALL dbt_create(mat_p_global%matrix, t_p_tmp)
492 CALL dbt_copy(t_p, t_p_tmp, move_data=.true.)
493 CALL dbt_copy_tensor_to_matrix(t_p_tmp, mat_p_global%matrix)
494 CALL dbt_destroy(t_p_tmp)
496 IF (do_ri_sos_laplace_mp2)
THEN
500 CALL dbcsr_add(mat_p_omega(jquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
502 CALL timeset(routinen//
"_Fourier_transform", handle3)
505 first_cycle_omega_loop = .true.
509 DO iquad = 1, num_integ_points
512 weight = weights_cos_tf_t_to_w(iquad, jquad)
514 IF (first_cycle_omega_loop)
THEN
517 CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)*weight)
519 CALL dbcsr_scale(mat_p_global%matrix, cos(omega*tau)/cos(omega_old*tau)*weight/weight_old)
522 CALL dbcsr_add(mat_p_omega(iquad, i_cell_t)%matrix, mat_p_global%matrix, 1.0_dp, 1.0_dp)
524 first_cycle_omega_loop = .false.
531 CALL timestop(handle3)
536 CALL timestop(handle2)
538 CALL dbt_batched_contract_finalize(t_3c_m_occ_tmp)
539 CALL dbt_batched_contract_finalize(t_3c_m_virt_tmp)
540 CALL dbt_batched_contract_finalize(t_3c_m_occ)
541 CALL dbt_batched_contract_finalize(t_3c_m_virt)
543 DO i = 1,
SIZE(t_3c_o_occ, 1)
544 DO j = 1,
SIZE(t_3c_o_occ, 2)
545 CALL dbt_batched_contract_finalize(t_3c_o_occ(i, j))
549 DO i = 1,
SIZE(t_3c_o_virt, 1)
550 DO j = 1,
SIZE(t_3c_o_virt, 2)
551 CALL dbt_batched_contract_finalize(t_3c_o_virt(i, j))
555 CALL dbt_destroy(t_p)
556 DO i_cell = 1, num_cells_dm
557 CALL dbt_destroy(t_dm_virt(i_cell))
558 CALL dbt_destroy(t_dm_occ(i_cell))
561 CALL dbt_destroy(t_3c_m_occ_tmp)
562 CALL dbt_destroy(t_3c_m_virt_tmp)
563 DEALLOCATE (t_dm_virt)
564 DEALLOCATE (t_dm_occ)
569 dbcsr_time = dbcsr_time + t2 - t1_flop
571 IF (unit_nr > 0)
THEN
572 WRITE (unit_nr,
'(/T3,A,1X,I3)') &
573 'RPA_LOW_SCALING_INFO| Info for time point', jquad
574 WRITE (unit_nr,
'(T6,A,T56,F25.1)') &
575 'Execution time (s):', t2 - t1
576 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
577 'Occupancy of D occ:', real(nze_dm_occ,
dp),
'/', occ_dm_occ*100,
'%'
578 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
579 'Occupancy of D virt:', real(nze_dm_virt,
dp),
'/', occ_dm_virt*100,
'%'
580 IF (do_gamma_rpa)
THEN
581 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
582 'Occupancy of 3c ints:', real(nze_o,
dp),
'/', occ_o*100,
'%'
583 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
584 'Occupancy of M occ:', real(nze_m_occ,
dp),
'/', occ_m_occ*100,
'%'
585 WRITE (unit_nr,
'(T6,A,T63,ES7.1,1X,A1,1X,F7.3,A1)') &
586 'Occupancy of M virt:', real(nze_m_virt,
dp),
'/', occ_m_virt*100,
'%'
594 CALL dbt_destroy(t_3c_m_occ)
595 CALL dbt_destroy(t_3c_m_virt)
597 DO i = 1,
SIZE(t_3c_o, 1)
598 DO j = 1,
SIZE(t_3c_o, 2)
599 CALL dbt_destroy(t_3c_o_occ(i, j))
600 CALL dbt_destroy(t_3c_o_virt(i, j))
604 CALL clean_up(mat_dm_occ_global, mat_dm_virt_global)
606 CALL timestop(handle)
615 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
INTENT(IN) :: mat_p_omega
617 INTEGER :: i_kp, jquad
619 DO jquad = 1,
SIZE(mat_p_omega, 1)
620 DO i_kp = 1,
SIZE(mat_p_omega, 2)
622 CALL dbcsr_set(mat_p_omega(jquad, i_kp)%matrix, 0.0_dp)
658 fm_mo_coeff_occ, fm_mo_coeff_virt, fm_mo_coeff_occ_scaled, &
659 fm_mo_coeff_virt_scaled, mat_dm_occ_global, mat_dm_virt_global, &
661 Eigenval, e_fermi, eps_filter, memory_info, unit_nr, &
662 jquad, do_kpoints_cubic_RPA, do_kpoints_from_Gamma, qs_env, &
663 num_cells_dm, index_to_cell_dm, para_env)
665 TYPE(
cp_fm_type),
INTENT(IN) :: fm_scaled_dm_occ_tau, &
666 fm_scaled_dm_virt_tau
667 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: tau_tj
668 INTEGER,
INTENT(IN) :: num_integ_points, nmo
669 TYPE(
cp_fm_type),
INTENT(IN) :: fm_mo_coeff_occ, fm_mo_coeff_virt, &
670 fm_mo_coeff_occ_scaled, &
671 fm_mo_coeff_virt_scaled
672 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_occ_global, mat_dm_virt_global
673 TYPE(
dbcsr_p_type),
DIMENSION(:),
INTENT(IN) :: matrix_s
674 INTEGER,
INTENT(IN) :: ispin
675 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: eigenval
676 REAL(kind=
dp),
INTENT(IN) :: e_fermi, eps_filter
677 LOGICAL,
INTENT(IN) :: memory_info
678 INTEGER,
INTENT(IN) :: unit_nr, jquad
679 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa, &
680 do_kpoints_from_gamma
682 INTEGER,
INTENT(OUT) :: num_cells_dm
683 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell_dm
686 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_mat_dm_global'
687 REAL(kind=
dp),
PARAMETER :: stabilize_exp = 70.0_dp
689 INTEGER :: handle, i_global, iib, iquad, jjb, &
690 ncol_local, nrow_local, size_dm_occ, &
692 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
695 CALL timeset(routinen, handle)
697 IF (memory_info .AND. unit_nr > 0)
WRITE (unit=unit_nr, fmt=
"(T3,A,T75,i6)") &
698 "RPA_LOW_SCALING_INFO| Started with time point: ", jquad
702 IF (do_kpoints_cubic_rpa)
THEN
705 ispin, num_integ_points, jquad, e_fermi, tau, &
706 eps_filter, num_cells_dm, index_to_cell_dm, &
707 remove_occ=.false., remove_virt=.true., first_jquad=1)
710 ispin, num_integ_points, jquad, e_fermi, tau, &
711 eps_filter, num_cells_dm, index_to_cell_dm, &
712 remove_occ=.true., remove_virt=.false., first_jquad=1)
714 ELSE IF (do_kpoints_from_gamma)
THEN
717 ispin, num_integ_points, jquad, e_fermi, tau, &
718 remove_occ=.false., remove_virt=.true., &
719 alloc_dm=(jquad == 1))
722 ispin, num_integ_points, jquad, e_fermi, tau, &
723 remove_occ=.true., remove_virt=.false., &
724 alloc_dm=(jquad == 1))
736 nrow_local=nrow_local, &
737 ncol_local=ncol_local, &
738 row_indices=row_indices, &
739 col_indices=col_indices)
746 DO jjb = 1, nrow_local
747 DO iib = 1, ncol_local
748 i_global = col_indices(iib)
752 IF (abs(tau*0.5_dp*(eigenval(i_global) - e_fermi)) < stabilize_exp)
THEN
753 fm_mo_coeff_occ_scaled%local_data(jjb, iib) = &
754 fm_mo_coeff_occ%local_data(jjb, iib)*exp(tau*0.5_dp*(eigenval(i_global) - e_fermi))
756 fm_mo_coeff_occ_scaled%local_data(jjb, iib) = 0.0_dp
764 nrow_local=nrow_local, &
765 ncol_local=ncol_local, &
766 row_indices=row_indices, &
767 col_indices=col_indices)
770 DO jjb = 1, nrow_local
771 DO iib = 1, ncol_local
772 i_global = col_indices(iib)
774 IF (abs(tau*0.5_dp*(eigenval(i_global) - e_fermi)) < stabilize_exp)
THEN
775 fm_mo_coeff_virt_scaled%local_data(jjb, iib) = &
776 fm_mo_coeff_virt%local_data(jjb, iib)*exp(-tau*0.5_dp*(eigenval(i_global) - e_fermi))
778 fm_mo_coeff_virt_scaled%local_data(jjb, iib) = 0.0_dp
789 CALL parallel_gemm(transa=
"N", transb=
"T", m=size_dm_occ, n=size_dm_occ, k=nmo, alpha=1.0_dp, &
790 matrix_a=fm_mo_coeff_occ_scaled, matrix_b=fm_mo_coeff_occ_scaled, beta=0.0_dp, &
791 matrix_c=fm_scaled_dm_occ_tau)
793 CALL parallel_gemm(transa=
"N", transb=
"T", m=size_dm_virt, n=size_dm_virt, k=nmo, alpha=1.0_dp, &
794 matrix_a=fm_mo_coeff_virt_scaled, matrix_b=fm_mo_coeff_virt_scaled, beta=0.0_dp, &
795 matrix_c=fm_scaled_dm_virt_tau)
800 NULLIFY (mat_dm_occ_global)
803 DO iquad = 1, num_integ_points
805 ALLOCATE (mat_dm_occ_global(iquad, 1)%matrix)
806 CALL dbcsr_create(matrix=mat_dm_occ_global(iquad, 1)%matrix, &
807 template=matrix_s(1)%matrix, &
808 matrix_type=dbcsr_type_no_symmetry)
815 mat_dm_occ_global(jquad, 1)%matrix, &
816 keep_sparsity=.false.)
818 CALL dbcsr_filter(mat_dm_occ_global(jquad, 1)%matrix, eps_filter)
822 NULLIFY (mat_dm_virt_global)
827 ALLOCATE (mat_dm_virt_global(jquad, 1)%matrix)
828 CALL dbcsr_create(matrix=mat_dm_virt_global(jquad, 1)%matrix, &
829 template=matrix_s(1)%matrix, &
830 matrix_type=dbcsr_type_no_symmetry)
832 mat_dm_virt_global(jquad, 1)%matrix, &
833 keep_sparsity=.false.)
835 CALL dbcsr_filter(mat_dm_virt_global(jquad, 1)%matrix, eps_filter)
839 CALL dbcsr_set(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
840 CALL dbcsr_set(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp)
841 CALL dbcsr_filter(mat_dm_occ_global(jquad - 1, 1)%matrix, 0.0_dp)
842 CALL dbcsr_filter(mat_dm_virt_global(jquad - 1, 1)%matrix, 0.0_dp)
847 CALL timestop(handle)
856 SUBROUTINE clean_up(mat_dm_occ_global, mat_dm_virt_global)
857 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_occ_global, mat_dm_virt_global
862 END SUBROUTINE clean_up
872 SUBROUTINE kpoint_density_matrices_rpa(kpoint, tau, e_fermi, remove_occ, remove_virt)
875 REAL(kind=
dp),
INTENT(IN) :: tau, e_fermi
876 LOGICAL,
INTENT(IN) :: remove_occ, remove_virt
878 CHARACTER(LEN=*),
PARAMETER :: routinen =
'kpoint_density_matrices_rpa'
879 REAL(kind=
dp),
PARAMETER :: stabilize_exp = 70.0_dp
881 INTEGER :: handle, i_mo, ikpgr, ispin, kplocal, &
883 INTEGER,
DIMENSION(2) :: kp_range
884 REAL(kind=
dp),
DIMENSION(:),
POINTER :: eigenvalues, exp_scaling, occupation
891 CALL timeset(routinen, handle)
894 cpassert(kpoint%use_real_wfn .EQV. .false.)
897 mo_set => kpoint%kp_env(1)%kpoint_env%mos(1, 1)
904 ALLOCATE (exp_scaling(nmo))
910 kplocal = kp_range(2) - kp_range(1) + 1
912 DO ikpgr = 1, kplocal
913 kp => kpoint%kp_env(ikpgr)%kpoint_env
914 nspin =
SIZE(kp%mos, 2)
916 mo_set => kp%mos(1, ispin)
917 CALL get_mo_set(mo_set, eigenvalues=eigenvalues)
918 rpmat => kp%wmat(1, ispin)
919 cpmat => kp%wmat(2, ispin)
920 CALL get_mo_set(mo_set, occupation_numbers=occupation)
923 IF (remove_virt)
THEN
937 IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp)
THEN
938 exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
940 exp_scaling(i_mo) = 0.0_dp
947 CALL parallel_gemm(
"N",
"T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, rpmat)
948 mo_set => kp%mos(2, ispin)
950 CALL parallel_gemm(
"N",
"T", nao, nao, nmo, -1.0_dp, mo_set%mo_coeff, fwork, 0.0_dp, cpmat)
952 CALL parallel_gemm(
"N",
"T", nao, nao, nmo, 1.0_dp, fwork, mo_set%mo_coeff, 1.0_dp, cpmat)
956 IF (remove_virt)
THEN
969 IF (abs(tau*0.5_dp*(eigenvalues(i_mo) - e_fermi)) < stabilize_exp)
THEN
970 exp_scaling(i_mo) = exp(-abs(tau*(eigenvalues(i_mo) - e_fermi)))
972 exp_scaling(i_mo) = 0.0_dp
978 CALL parallel_gemm(
"N",
"T", nao, nao, nmo, 1.0_dp, mo_set%mo_coeff, fwork, 1.0_dp, rpmat)
985 DEALLOCATE (exp_scaling)
987 CALL timestop(handle)
989 END SUBROUTINE kpoint_density_matrices_rpa
1008 eps_filter, num_cells_dm, index_to_cell_dm, remove_occ, remove_virt, &
1010 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_global
1012 INTEGER,
INTENT(IN) :: ispin, num_integ_points, jquad
1013 REAL(kind=
dp),
INTENT(IN) :: e_fermi, tau, eps_filter
1014 INTEGER,
INTENT(OUT) :: num_cells_dm
1015 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell_dm
1016 LOGICAL,
INTENT(IN) :: remove_occ, remove_virt
1017 INTEGER,
INTENT(IN) :: first_jquad
1019 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_transl_dm'
1021 INTEGER :: handle, i_dim, i_img, iquad, jspin, nspin
1022 INTEGER,
DIMENSION(3) :: cell_grid_dm
1024 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_global_work, matrix_s_kp
1028 CALL timeset(routinen, handle)
1031 matrix_s_kp=matrix_s_kp, &
1041 cell_grid_dm(i_dim) = (kpoints%nkp_grid(i_dim)/2)*2 - 1
1044 num_cells_dm = cell_grid_dm(1)*cell_grid_dm(2)*cell_grid_dm(3)
1046 NULLIFY (mat_dm_global_work)
1051 DO i_img = 1, num_cells_dm
1053 ALLOCATE (mat_dm_global_work(jspin, i_img)%matrix)
1054 CALL dbcsr_create(matrix=mat_dm_global_work(jspin, i_img)%matrix, &
1055 template=matrix_s_kp(1, 1)%matrix, &
1057 matrix_type=dbcsr_type_no_symmetry)
1061 CALL dbcsr_set(mat_dm_global_work(ispin, i_img)%matrix, 0.0_dp)
1068 CALL kpoint_density_matrices_rpa(kpoints, tau, e_fermi, &
1069 remove_occ=remove_occ, remove_virt=remove_virt)
1072 CALL init_cell_index_rpa(cell_grid_dm, kpoints%cell_to_index, kpoints%index_to_cell, cell)
1076 CALL density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, kpoints%index_to_cell)
1079 index_to_cell_dm => kpoints%index_to_cell
1082 IF (jquad == first_jquad)
THEN
1084 NULLIFY (mat_dm_global)
1085 ALLOCATE (mat_dm_global(first_jquad:num_integ_points, num_cells_dm))
1087 DO iquad = first_jquad, num_integ_points
1088 DO i_img = 1, num_cells_dm
1089 NULLIFY (mat_dm_global(iquad, i_img)%matrix)
1090 ALLOCATE (mat_dm_global(iquad, i_img)%matrix)
1091 CALL dbcsr_create(matrix=mat_dm_global(iquad, i_img)%matrix, &
1092 template=matrix_s_kp(1, 1)%matrix, &
1093 matrix_type=dbcsr_type_no_symmetry)
1100 DO i_img = 1, num_cells_dm
1103 CALL dbcsr_filter(mat_dm_global_work(ispin, i_img)%matrix, eps_filter)
1105 CALL dbcsr_copy(mat_dm_global(jquad, i_img)%matrix, &
1106 mat_dm_global_work(ispin, i_img)%matrix)
1112 CALL timestop(handle)
1130 remove_occ, remove_virt, alloc_dm)
1131 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_global
1133 INTEGER,
INTENT(IN) :: ispin, num_integ_points, jquad
1134 REAL(kind=
dp),
INTENT(IN) :: e_fermi, tau
1135 LOGICAL,
INTENT(IN) :: remove_occ, remove_virt, alloc_dm
1137 CHARACTER(LEN=*),
PARAMETER :: routinen =
'compute_periodic_dm'
1139 INTEGER :: handle, iquad, jspin, nspin, num_cells_dm
1140 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_global_work, matrix_s_kp
1144 CALL timeset(routinen, handle)
1146 NULLIFY (matrix_s_kp, mos)
1149 matrix_s_kp=matrix_s_kp, &
1152 kpoints_g => qs_env%mp2_env%ri_rpa_im_time%kpoints_G
1158 NULLIFY (mat_dm_global_work)
1164 NULLIFY (mat_dm_global)
1165 ALLOCATE (mat_dm_global(1:num_integ_points, num_cells_dm))
1167 DO iquad = 1, num_integ_points
1168 NULLIFY (mat_dm_global(iquad, 1)%matrix)
1169 ALLOCATE (mat_dm_global(iquad, 1)%matrix)
1170 CALL dbcsr_create(matrix=mat_dm_global(iquad, 1)%matrix, &
1171 template=matrix_s_kp(1, 1)%matrix, &
1172 matrix_type=dbcsr_type_no_symmetry)
1180 ALLOCATE (mat_dm_global_work(jspin, 1)%matrix)
1181 CALL dbcsr_create(matrix=mat_dm_global_work(jspin, 1)%matrix, &
1182 template=matrix_s_kp(1, 1)%matrix, &
1183 matrix_type=dbcsr_type_no_symmetry)
1187 CALL dbcsr_set(mat_dm_global_work(jspin, 1)%matrix, 0.0_dp)
1192 CALL kpoint_density_matrices_rpa(kpoints_g, tau, e_fermi, &
1193 remove_occ=remove_occ, remove_virt=remove_virt)
1195 CALL density_matrix_from_kp_to_mic(kpoints_g, mat_dm_global_work, qs_env)
1197 CALL dbcsr_copy(mat_dm_global(jquad, 1)%matrix, &
1198 mat_dm_global_work(ispin, 1)%matrix)
1202 CALL timestop(handle)
1212 SUBROUTINE density_matrix_from_kp_to_mic(kpoints_G, mat_dm_global_work, qs_env)
1215 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
POINTER :: mat_dm_global_work
1218 CHARACTER(LEN=*),
PARAMETER :: routinen =
'density_matrix_from_kp_to_mic'
1220 INTEGER :: handle, iatom, iatom_old, ik, irow, &
1221 ispin, jatom, jatom_old, jcol, nao, &
1222 ncol_local, nkp, nrow_local, nspin, &
1224 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_from_ao_index
1225 INTEGER,
DIMENSION(:),
POINTER :: col_indices, row_indices
1226 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1227 REAL(kind=
dp) :: contribution, weight_im, weight_re
1228 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat
1229 REAL(kind=
dp),
DIMENSION(:),
POINTER :: wkp
1230 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: xkp
1238 CALL timeset(routinen, handle)
1242 CALL cp_fm_create(fm_mat_work, kpoints_g%kp_env(1)%kpoint_env%wmat(1, 1)%matrix_struct)
1246 index_to_cell => kpoints_g%index_to_cell
1247 num_cells =
SIZE(index_to_cell, 2)
1249 nspin =
SIZE(mat_dm_global_work, 1)
1251 mo_set => kpoints_g%kp_env(1)%kpoint_env%mos(1, 1)
1254 ALLOCATE (atom_from_ao_index(nao))
1258 CALL cp_fm_get_info(matrix=kpoints_g%kp_env(1)%kpoint_env%wmat(1, 1), &
1259 nrow_local=nrow_local, &
1260 ncol_local=ncol_local, &
1261 row_indices=row_indices, &
1262 col_indices=col_indices)
1264 NULLIFY (cell, particle_set)
1265 CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set)
1273 CALL dbcsr_set(mat_dm_global_work(ispin, 1)%matrix, 0.0_dp)
1277 kp => kpoints_g%kp_env(ik)%kpoint_env
1278 rpmat => kp%wmat(1, ispin)
1279 cpmat => kp%wmat(2, ispin)
1281 DO irow = 1, nrow_local
1282 DO jcol = 1, ncol_local
1284 iatom = atom_from_ao_index(row_indices(irow))
1285 jatom = atom_from_ao_index(col_indices(jcol))
1287 IF (iatom .NE. iatom_old .OR. jatom .NE. jatom_old)
THEN
1290 num_cells, iatom, jatom, xkp(1:3, ik), wkp(ik), &
1291 cell, index_to_cell, hmat, particle_set)
1299 contribution = weight_re*rpmat%local_data(irow, jcol) - &
1300 weight_im*cpmat%local_data(irow, jcol)
1302 fm_mat_work%local_data(irow, jcol) = fm_mat_work%local_data(irow, jcol) + contribution
1309 CALL copy_fm_to_dbcsr(fm_mat_work, mat_dm_global_work(ispin, 1)%matrix, keep_sparsity=.false.)
1315 DEALLOCATE (atom_from_ao_index)
1317 CALL timestop(handle)
1319 END SUBROUTINE density_matrix_from_kp_to_mic
1327 SUBROUTINE density_matrix_from_kp_to_transl(kpoints, mat_dm_global_work, index_to_cell)
1330 TYPE(
dbcsr_p_type),
DIMENSION(:, :),
INTENT(IN) :: mat_dm_global_work
1331 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: index_to_cell
1333 CHARACTER(LEN=*),
PARAMETER :: routinen =
'density_matrix_from_kp_to_transl'
1335 INTEGER :: handle, icell, ik, ispin, nkp, nspin, &
1337 REAL(kind=
dp) :: arg, coskl, sinkl
1338 REAL(kind=
dp),
DIMENSION(:),
POINTER :: wkp
1339 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: xkp
1341 TYPE(
dbcsr_type),
POINTER :: mat_work_im, mat_work_re
1344 CALL timeset(routinen, handle)
1348 NULLIFY (mat_work_re)
1351 template=mat_dm_global_work(1, 1)%matrix, &
1352 matrix_type=dbcsr_type_no_symmetry)
1354 NULLIFY (mat_work_im)
1357 template=mat_dm_global_work(1, 1)%matrix, &
1358 matrix_type=dbcsr_type_no_symmetry)
1362 nspin =
SIZE(mat_dm_global_work, 1)
1364 cpassert(
SIZE(mat_dm_global_work, 2) ==
SIZE(index_to_cell, 2))
1368 DO icell = 1,
SIZE(mat_dm_global_work, 2)
1370 CALL dbcsr_set(mat_dm_global_work(ispin, icell)%matrix, 0.0_dp)
1380 kp => kpoints%kp_env(ik)%kpoint_env
1381 rpmat => kp%wmat(1, ispin)
1382 cpmat => kp%wmat(2, ispin)
1387 DO icell = 1,
SIZE(mat_dm_global_work, 2)
1389 xcell = index_to_cell(1, icell)
1390 ycell = index_to_cell(2, icell)
1391 zcell = index_to_cell(3, icell)
1393 arg = real(xcell,
dp)*xkp(1, ik) + real(ycell,
dp)*xkp(2, ik) + real(zcell,
dp)*xkp(3, ik)
1394 coskl = wkp(ik)*cos(
twopi*arg)
1395 sinkl = wkp(ik)*sin(
twopi*arg)
1397 CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_re, 1.0_dp, coskl)
1398 CALL dbcsr_add(mat_dm_global_work(ispin, icell)%matrix, mat_work_im, 1.0_dp, sinkl)
1408 CALL timestop(handle)
1410 END SUBROUTINE density_matrix_from_kp_to_transl
1420 INTEGER,
DIMENSION(3),
INTENT(IN) :: cell_grid
1421 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
1422 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell
1423 TYPE(
cell_type),
INTENT(IN),
POINTER :: cell
1425 CHARACTER(LEN=*),
PARAMETER :: routinen =
'init_cell_index_rpa'
1427 INTEGER :: cell_counter, handle, i_cell, &
1428 index_min_dist, num_cells, xcell, &
1430 INTEGER,
DIMENSION(3) :: itm
1431 INTEGER,
DIMENSION(:, :),
POINTER :: index_to_cell_unsorted
1432 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index_unsorted
1433 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: abs_cell_vectors
1434 REAL(kind=
dp),
DIMENSION(3) :: cell_vector
1435 REAL(kind=
dp),
DIMENSION(3, 3) :: hmat
1437 CALL timeset(routinen, handle)
1441 num_cells = cell_grid(1)*cell_grid(2)*cell_grid(3)
1442 itm(:) = cell_grid(:)/2
1446 cpassert(cell_grid(1) .NE. itm(1)*2)
1447 cpassert(cell_grid(2) .NE. itm(2)*2)
1448 cpassert(cell_grid(3) .NE. itm(3)*2)
1450 IF (
ASSOCIATED(cell_to_index))
DEALLOCATE (cell_to_index)
1451 IF (
ASSOCIATED(index_to_cell))
DEALLOCATE (index_to_cell)
1453 ALLOCATE (cell_to_index_unsorted(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1454 cell_to_index_unsorted(:, :, :) = 0
1456 ALLOCATE (index_to_cell_unsorted(3, num_cells))
1457 index_to_cell_unsorted(:, :) = 0
1459 ALLOCATE (cell_to_index(-itm(1):itm(1), -itm(2):itm(2), -itm(3):itm(3)))
1460 cell_to_index(:, :, :) = 0
1462 ALLOCATE (index_to_cell(3, num_cells))
1463 index_to_cell(:, :) = 0
1465 ALLOCATE (abs_cell_vectors(1:num_cells))
1469 DO xcell = -itm(1), itm(1)
1470 DO ycell = -itm(2), itm(2)
1471 DO zcell = -itm(3), itm(3)
1473 cell_counter = cell_counter + 1
1474 cell_to_index_unsorted(xcell, ycell, zcell) = cell_counter
1476 index_to_cell_unsorted(1, cell_counter) = xcell
1477 index_to_cell_unsorted(2, cell_counter) = ycell
1478 index_to_cell_unsorted(3, cell_counter) = zcell
1480 cell_vector(1:3) = matmul(hmat, real(index_to_cell_unsorted(1:3, cell_counter),
dp))
1482 abs_cell_vectors(cell_counter) = sqrt(cell_vector(1)**2 + cell_vector(2)**2 + cell_vector(3)**2)
1490 DO i_cell = 1, num_cells/2 + 1
1492 index_min_dist = minloc(abs_cell_vectors(1:num_cells/2 + 1), dim=1)
1494 xcell = index_to_cell_unsorted(1, index_min_dist)
1495 ycell = index_to_cell_unsorted(2, index_min_dist)
1496 zcell = index_to_cell_unsorted(3, index_min_dist)
1498 index_to_cell(1, i_cell) = xcell
1499 index_to_cell(2, i_cell) = ycell
1500 index_to_cell(3, i_cell) = zcell
1502 cell_to_index(xcell, ycell, zcell) = i_cell
1504 abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1509 DO i_cell = num_cells/2 + 2, num_cells
1511 index_min_dist = minloc(abs_cell_vectors(1:num_cells), dim=1)
1513 xcell = index_to_cell_unsorted(1, index_min_dist)
1514 ycell = index_to_cell_unsorted(2, index_min_dist)
1515 zcell = index_to_cell_unsorted(3, index_min_dist)
1517 index_to_cell(1, i_cell) = xcell
1518 index_to_cell(2, i_cell) = ycell
1519 index_to_cell(3, i_cell) = zcell
1521 cell_to_index(xcell, ycell, zcell) = i_cell
1523 abs_cell_vectors(index_min_dist) = 10000000000.0_dp
1527 DEALLOCATE (index_to_cell_unsorted, cell_to_index_unsorted, abs_cell_vectors)
1529 CALL timestop(handle)
1544 SUBROUTINE get_diff_index_3c(i_cell_R, i_cell_S, i_cell_R_minus_S, index_to_cell_3c, &
1545 cell_to_index_3c, index_to_cell_dm, R_minus_S_needed, &
1546 do_kpoints_cubic_RPA)
1548 INTEGER,
INTENT(IN) :: i_cell_r, i_cell_s
1549 INTEGER,
INTENT(OUT) :: i_cell_r_minus_s
1550 INTEGER,
ALLOCATABLE,
DIMENSION(:, :),
INTENT(IN) :: index_to_cell_3c
1551 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :), &
1552 INTENT(IN) :: cell_to_index_3c
1553 INTEGER,
DIMENSION(:, :),
INTENT(IN),
POINTER :: index_to_cell_dm
1554 LOGICAL,
INTENT(OUT) :: r_minus_s_needed
1555 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa
1557 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_diff_index_3c'
1559 INTEGER :: handle, x_cell_r, x_cell_r_minus_s, x_cell_s, y_cell_r, y_cell_r_minus_s, &
1560 y_cell_s, z_cell_r, z_cell_r_minus_s, z_cell_s
1562 CALL timeset(routinen, handle)
1564 IF (do_kpoints_cubic_rpa)
THEN
1566 x_cell_r = index_to_cell_3c(1, i_cell_r)
1567 y_cell_r = index_to_cell_3c(2, i_cell_r)
1568 z_cell_r = index_to_cell_3c(3, i_cell_r)
1570 x_cell_s = index_to_cell_dm(1, i_cell_s)
1571 y_cell_s = index_to_cell_dm(2, i_cell_s)
1572 z_cell_s = index_to_cell_dm(3, i_cell_s)
1574 x_cell_r_minus_s = x_cell_r - x_cell_s
1575 y_cell_r_minus_s = y_cell_r - y_cell_s
1576 z_cell_r_minus_s = z_cell_r - z_cell_s
1578 IF (x_cell_r_minus_s .GE. lbound(cell_to_index_3c, 1) .AND. &
1579 x_cell_r_minus_s .LE. ubound(cell_to_index_3c, 1) .AND. &
1580 y_cell_r_minus_s .GE. lbound(cell_to_index_3c, 2) .AND. &
1581 y_cell_r_minus_s .LE. ubound(cell_to_index_3c, 2) .AND. &
1582 z_cell_r_minus_s .GE. lbound(cell_to_index_3c, 3) .AND. &
1583 z_cell_r_minus_s .LE. ubound(cell_to_index_3c, 3))
THEN
1585 i_cell_r_minus_s = cell_to_index_3c(x_cell_r_minus_s, y_cell_r_minus_s, z_cell_r_minus_s)
1588 IF (i_cell_r_minus_s == 0)
THEN
1590 r_minus_s_needed = .false.
1591 i_cell_r_minus_s = 0
1595 r_minus_s_needed = .true.
1601 i_cell_r_minus_s = 0
1602 r_minus_s_needed = .false.
1608 r_minus_s_needed = .true.
1609 i_cell_r_minus_s = 1
1613 CALL timestop(handle)
1615 END SUBROUTINE get_diff_index_3c
1629 SUBROUTINE get_diff_diff_index_3c(i_cell_R, i_cell_S, i_cell_T, i_cell_R_minus_S_minus_T, &
1630 index_to_cell_3c, cell_to_index_3c, index_to_cell_dm, &
1631 R_minus_S_minus_T_needed, &
1632 do_kpoints_cubic_RPA)
1634 INTEGER,
INTENT(IN) :: i_cell_r, i_cell_s, i_cell_t
1635 INTEGER,
INTENT(OUT) :: i_cell_r_minus_s_minus_t
1636 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: index_to_cell_3c
1637 INTEGER,
ALLOCATABLE,
DIMENSION(:, :, :), &
1638 INTENT(IN) :: cell_to_index_3c
1639 INTEGER,
DIMENSION(:, :),
INTENT(IN) :: index_to_cell_dm
1640 LOGICAL,
INTENT(OUT) :: r_minus_s_minus_t_needed
1641 LOGICAL,
INTENT(IN) :: do_kpoints_cubic_rpa
1643 CHARACTER(LEN=*),
PARAMETER :: routinen =
'get_diff_diff_index_3c'
1645 INTEGER :: handle, x_cell_r, x_cell_r_minus_s_minus_t, x_cell_s, x_cell_t, y_cell_r, &
1646 y_cell_r_minus_s_minus_t, y_cell_s, y_cell_t, z_cell_r, z_cell_r_minus_s_minus_t, &
1649 CALL timeset(routinen, handle)
1651 IF (do_kpoints_cubic_rpa)
THEN
1653 x_cell_r = index_to_cell_3c(1, i_cell_r)
1654 y_cell_r = index_to_cell_3c(2, i_cell_r)
1655 z_cell_r = index_to_cell_3c(3, i_cell_r)
1657 x_cell_s = index_to_cell_dm(1, i_cell_s)
1658 y_cell_s = index_to_cell_dm(2, i_cell_s)
1659 z_cell_s = index_to_cell_dm(3, i_cell_s)
1661 x_cell_t = index_to_cell_dm(1, i_cell_t)
1662 y_cell_t = index_to_cell_dm(2, i_cell_t)
1663 z_cell_t = index_to_cell_dm(3, i_cell_t)
1665 x_cell_r_minus_s_minus_t = x_cell_r - x_cell_s - x_cell_t
1666 y_cell_r_minus_s_minus_t = y_cell_r - y_cell_s - y_cell_t
1667 z_cell_r_minus_s_minus_t = z_cell_r - z_cell_s - z_cell_t
1669 IF (x_cell_r_minus_s_minus_t .GE. lbound(cell_to_index_3c, 1) .AND. &
1670 x_cell_r_minus_s_minus_t .LE. ubound(cell_to_index_3c, 1) .AND. &
1671 y_cell_r_minus_s_minus_t .GE. lbound(cell_to_index_3c, 2) .AND. &
1672 y_cell_r_minus_s_minus_t .LE. ubound(cell_to_index_3c, 2) .AND. &
1673 z_cell_r_minus_s_minus_t .GE. lbound(cell_to_index_3c, 3) .AND. &
1674 z_cell_r_minus_s_minus_t .LE. ubound(cell_to_index_3c, 3))
THEN
1676 i_cell_r_minus_s_minus_t = cell_to_index_3c(x_cell_r_minus_s_minus_t, &
1677 y_cell_r_minus_s_minus_t, &
1678 z_cell_r_minus_s_minus_t)
1681 IF (i_cell_r_minus_s_minus_t == 0)
THEN
1683 r_minus_s_minus_t_needed = .false.
1687 r_minus_s_minus_t_needed = .true.
1693 i_cell_r_minus_s_minus_t = 0
1694 r_minus_s_minus_t_needed = .false.
1701 r_minus_s_minus_t_needed = .true.
1702 i_cell_r_minus_s_minus_t = 1
1706 CALL timestop(handle)
1708 END SUBROUTINE get_diff_diff_index_3c
Handles all functions related to the CELL.
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_clear(matrix)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
subroutine, public dbcsr_reserve_all_blocks(matrix)
Reserves all blocks.
DBCSR operations in CP2K.
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_column_scale(matrixa, scaling)
scales column i of matrix a with scaling(i)
subroutine, public cp_fm_scale(alpha, matrix_a)
scales a matrix matrix_a = alpha * matrix_b
represent the structure of a full matrix
represent a full matrix distributed on many processors
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
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....
Types and set/get functions for HFX.
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
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 twopi
Interface to the message passing library MPI.
Types needed for MP2 calculations.
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
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.
Utility methods to build 3-center integral tensors of various types.
subroutine, public create_2c_tensor(t2c, dist_1, dist_2, pgrid, sizes_1, sizes_2, order, name)
...
Utility methods to build 3-center integral tensors of various types.
subroutine, public decompress_tensor(tensor, blk_indices, compressed, eps)
...
subroutine, public get_tensor_occupancy(tensor, nze, occ)
...
Utility routines for GW with imaginary time.
subroutine, public compute_weight_re_im(weight_re, weight_im, num_cells, iatom, jatom, xkp, wkp_w, cell, index_to_cell, hmat, particle_set)
...
subroutine, public get_atom_index_from_basis_function_index(qs_env, atom_from_basis_index, basis_size, basis_type, first_bf_from_atom)
...
Routines for low-scaling RPA/GW with imaginary time.
subroutine, public zero_mat_p_omega(mat_p_omega)
...
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)
...
subroutine, public compute_mat_p_omega(mat_p_omega, 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, mat_p_global, matrix_s, ispin, t_3c_m, t_3c_o, t_3c_o_compressed, t_3c_o_ind, starts_array_mc, ends_array_mc, starts_array_mc_block, ends_array_mc_block, weights_cos_tf_t_to_w, tj, tau_tj, e_fermi, eps_filter, alpha, eps_filter_im_time, eigenval, nmo, num_integ_points, cut_memory, unit_nr, mp2_env, para_env, qs_env, do_kpoints_from_gamma, index_to_cell_3c, cell_to_index_3c, has_mat_p_blocks, do_ri_sos_laplace_mp2, dbcsr_time, dbcsr_nflop)
...
subroutine, public compute_periodic_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, remove_occ, remove_virt, alloc_dm)
...
subroutine, public init_cell_index_rpa(cell_grid, cell_to_index, index_to_cell, cell)
...
subroutine, public compute_transl_dm(mat_dm_global, qs_env, ispin, num_integ_points, jquad, e_fermi, tau, eps_filter, num_cells_dm, index_to_cell_dm, remove_occ, remove_virt, first_jquad)
...
Type defining parameters related to the simulation cell.
keeps the information about the structure of a full matrix
Keeps information about a specific k-point.
Contains information about kpoints.
stores all the informations relevant to an mpi environment