comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-29 06:35:28 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

GW: performance improvements

Browse source
This commit is contained in:
Max Graml 2026-01-28 17:47:31 +01:00 • committed by GitHub
parent e60cd36aed
commit 35340893dd
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 605 additions and 154 deletions
Download patch file Download diff file Expand all files Collapse all files

20
src/dbm/dbm_multiply_comm.c
View file

@ -67,6 +67,16 @@ typedef struct {
int col_size;
} plan_t;
/*******************************************************************************
* \brief Private routine for calculating tick indices in pack plans.
* \author Maximilian Graml
******************************************************************************/
static inline unsigned long long calculate_tick_index(int sum_index,
int nticks) {
// 1021 is used as a random prime to scramble the index
return ((unsigned long long)sum_index * 1021ULL) % (unsigned long long)nticks;
}
/*******************************************************************************
* \brief Private routine for planing packs.
* \author Ole Schuett
@ -94,8 +104,8 @@ static void create_pack_plans(const bool trans_matrix, const bool trans_dist,
for (int iblock = 0; iblock < shard->nblocks; iblock++) {
const dbm_block_t *blk = &shard->blocks[iblock];
const int sum_index = (trans_matrix) ? blk->row : blk->col;
const int itick = (1021 * sum_index) % nticks; // 1021 = a random prime
const int ipack = itick / dist_ticks->nranks;
unsigned long long itick64 = calculate_tick_index(sum_index, nticks);
const int ipack = itick64 / dist_ticks->nranks;
nblks_mythread[ipack]++;
}
}
@ -124,11 +134,11 @@ static void create_pack_plans(const bool trans_matrix, const bool trans_dist,
const dbm_block_t *blk = &shard->blocks[iblock];
const int free_index = (trans_matrix) ? blk->col : blk->row;
const int sum_index = (trans_matrix) ? blk->row : blk->col;
const int itick = (1021 * sum_index) % nticks; // Same mapping as above.
const int ipack = itick / dist_ticks->nranks;
unsigned long long itick64 = calculate_tick_index(sum_index, nticks);
const int ipack = itick64 / dist_ticks->nranks;
// Compute rank to which this block should be sent.
const int coord_free_idx = dist_indices->index2coord[free_index];
const int coord_sum_idx = itick % dist_ticks->nranks;
const int coord_sum_idx = itick64 % dist_ticks->nranks;
const int coords[2] = {(trans_dist) ? coord_sum_idx : coord_free_idx,
(trans_dist) ? coord_free_idx : coord_sum_idx};
const int rank = cp_mpi_cart_rank(comm, coords);

1
src/dbt/tas/dbt_tas_mm.F
View file

@ -1637,7 +1637,6 @@ CONTAINS
IF (move_data) CALL dbm_clear(matrix_in)
END IF
END IF
CALL timestop(handle)
END SUBROUTINE convert_to_new_pgrid

449
src/gw_large_cell_gamma.F
View file

@ -55,6 +55,7 @@ MODULE gw_large_cell_gamma
dbt_copy,&
dbt_create,&
dbt_destroy,&
dbt_filter,&
dbt_type
USE gw_communication, ONLY: fm_to_local_tensor,&
local_dbt_to_global_mat
@ -158,7 +159,10 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_mat_chi_Gamma_tau'
INTEGER :: handle, i_intval_idx, i_t, inner_loop_atoms_interval_index, ispin, j_intval_idx
INTEGER, DIMENSION(2) :: i_atoms, IL_atoms, j_atoms
INTEGER(KIND=int_8) :: flop
INTEGER, DIMENSION(2) :: bounds_P, bounds_Q, i_atoms, IL_atoms, &
j_atoms
INTEGER, DIMENSION(2, 2) :: bounds_comb
LOGICAL :: dist_too_long_i, dist_too_long_j
REAL(KIND=dp) :: t1, tau
TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_3c_for_Gocc, &
@ -180,10 +184,10 @@ CONTAINS
keep_sparsity=.FALSE.)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
'Read χ(iτ,k=0) from file for time point ', i_t, ' /', &
bs_env%num_time_freq_points, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
CYCLE
@ -219,6 +223,15 @@ CONTAINS
i_atoms = bs_env%i_atom_intervals(1:2, i_intval_idx)
j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
IF (bs_env%skip_chi(i_intval_idx, j_intval_idx)) THEN
! Do that only after first timestep to avoid skips due to vanishing G
! caused by gaps
IF (i_t == 2) THEN
bs_env%n_skip_chi = bs_env%n_skip_chi + 1
END IF
CYCLE
END IF
DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
@ -227,33 +240,63 @@ CONTAINS
CALL check_dist(j_atoms, IL_atoms, qs_env, bs_env, dist_too_long_j)
IF (dist_too_long_i .OR. dist_too_long_j) CYCLE
! 2. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gocc, i_atoms, IL_atoms)
! 2. compute 3-center integrals (Pν|µ) ("|": truncated Coulomb operator)
CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gocc, &
atoms_AO_1=i_atoms, atoms_AO_2=IL_atoms)
! 3. tensor operation M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
! 3. tensor operation M_Pνλ(iτ) = sum_µ (Pν|µ) G^occ_λµ(i|τ|,k=0)
CALL G_times_3c(t_3c_for_Gocc, t_2c_Gocc, t_3c_x_Gocc, bs_env, &
j_atoms, i_atoms, IL_atoms)
! 4. compute 3-center integrals (σλ|Q) ("|": truncated Coulomb operator)
CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gvir, j_atoms, IL_atoms)
! 4. compute 3-center integrals (Qλ|σ) ("|": truncated Coulomb operator)
CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_Gvir, &
atoms_AO_1=j_atoms, atoms_AO_2=IL_atoms)
! 5. tensor operation N_νλQ(iτ) = sum_σ (σλ|Q) G^vir_σν(i|τ|,k=0)
! 5. tensor operation N_Qλν(iτ) = sum_σ (Qλ|σ) G^vir_νσ(i|τ|,k=0)
CALL G_times_3c(t_3c_for_Gvir, t_2c_Gvir, t_3c_x_Gvir, bs_env, &
i_atoms, j_atoms, IL_atoms)
END DO ! IL_atoms
! 6. reorder tensors
! 6. reorder tensors: M_Pνλ -> M_Pλν
CALL dbt_copy(t_3c_x_Gocc, t_3c_x_Gocc_2, move_data=.TRUE., order=[1, 3, 2])
CALL dbt_copy(t_3c_x_Gvir, t_3c_x_Gvir_2, move_data=.TRUE.)
! 7. tensor operation χ_PQ(iτ,k=0) = sum_λν M_λνP(iτ) N_νλQ(iτ),
CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
tensor_1=t_3c_x_Gocc_2, tensor_2=t_3c_x_Gvir_2, &
beta=1.0_dp, tensor_3=bs_env%t_chi, &
contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
filter_eps=bs_env%eps_filter, move_data=.TRUE.)
! 7. tensor operation χ_PQ(iτ,k=0) = sum_λν M_Pλν(iτ) N_Qλν(iτ),
! Bounds:
! "comb" (combined index)
! -> λ bounds from j_atoms
! -> ν bounds from i_atoms
! P -> sparse in ν (see 3.)
! Q -> sparse in λ (see 5.)
bounds_comb(1:2, 1) = [bs_env%i_ao_start_from_atom(j_atoms(1)), &
bs_env%i_ao_end_from_atom(j_atoms(2))]
bounds_comb(1:2, 2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
bs_env%i_ao_end_from_atom(i_atoms(2))]
CALL get_bounds_from_atoms(bounds_P, i_atoms, [1, bs_env%n_atom], &
bs_env%min_RI_idx_from_AO_AO_atom, &
bs_env%max_RI_idx_from_AO_AO_atom)
CALL get_bounds_from_atoms(bounds_Q, [1, bs_env%n_atom], j_atoms, &
bs_env%min_RI_idx_from_AO_AO_atom, &
bs_env%max_RI_idx_from_AO_AO_atom)
IF (bounds_Q(1) > bounds_Q(2) .OR. bounds_P(1) > bounds_P(2)) THEN
flop = 0_int_8
ELSE
CALL dbt_contract(alpha=bs_env%spin_degeneracy, &
tensor_1=t_3c_x_Gocc_2, tensor_2=t_3c_x_Gvir_2, &
beta=1.0_dp, tensor_3=bs_env%t_chi, &
contract_1=[2, 3], notcontract_1=[1], map_1=[1], &
contract_2=[2, 3], notcontract_2=[1], map_2=[2], &
bounds_1=bounds_comb, &
bounds_2=bounds_P, &
bounds_3=bounds_Q, &
filter_eps=bs_env%eps_filter, move_data=.FALSE., flop=flop, &
unit_nr=bs_env%unit_nr_contract, &
log_verbose=bs_env%print_contract_verbose)
END IF
IF (flop == 0_int_8) bs_env%skip_chi(i_intval_idx, j_intval_idx) = .TRUE.
END DO ! j_atoms
END DO ! i_atoms
@ -272,9 +315,9 @@ CONTAINS
t_3c_x_Gocc, t_3c_x_Gvir, t_3c_x_Gocc_2, t_3c_x_Gvir_2)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I13,A,I3,A,F10.1,A)') &
'Computed χ(iτ,k=0) for time point', i_t, ' /', bs_env%num_time_freq_points, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
END DO ! i_t
@ -356,14 +399,14 @@ CONTAINS
CALL timeset(routineN, handle)
CALL dbt_create(bs_env%t_G, t_2c_Gocc)
CALL dbt_create(bs_env%t_G, t_2c_Gvir)
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gocc)
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gvir)
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gocc)
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gvir)
CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gocc_2)
CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gvir_2)
CALL dbt_create(bs_env%t_G, t_2c_Gocc, name="Gocc 2c (AO|AO)")
CALL dbt_create(bs_env%t_G, t_2c_Gvir, name="Gvir 2c (AO|AO)")
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gocc, name="Gocc 3c (RI AO|AO)")
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_for_Gvir, name="Gvir 3c (RI AO|AO)")
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gocc, name="xGocc 3c (RI AO|AO)")
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_Gvir, name="xGvir 3c (RI AO|AO)")
CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gocc_2, name="x2Gocc 3c (RI AO|AO)")
CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_Gvir_2, name="x2Gvir 3c (RI AO|AO)")
CALL timestop(handle)
@ -614,6 +657,8 @@ CONTAINS
bounds_k=atoms_AO_2, &
desymmetrize=.FALSE.)
CALL dbt_filter(t_3c_array(1, 1), bs_env%eps_filter)
CALL dbt_copy(t_3c_array(1, 1), t_3c, move_data=.TRUE.)
CALL dbt_destroy(t_3c_array(1, 1))
@ -641,34 +686,63 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'G_times_3c'
INTEGER :: handle
INTEGER, DIMENSION(2) :: bounds_IL, bounds_l
INTEGER, DIMENSION(2, 2) :: bounds_k
INTEGER(KIND=int_8) :: flop
INTEGER, DIMENSION(2) :: bounds_ao_1, bounds_IL
INTEGER, DIMENSION(2, 2) :: bounds_comb
CALL timeset(routineN, handle)
! JW bounds_IL and bounds_k do not safe any operations, but maybe communication
! maybe remove "bounds_1=bounds_IL, &" and "bounds_2=bounds_k, &" later and
! check whether performance improves
! Bounds reduce needed memory and therefore scaling behavior
! Operations are of the form, e.g, M_Pνλ = sum_µ (Pν|µ) G_λµ
! "comb" (combined index)
! -> P sparse in ν and µ
! -> λ bounds from j_atoms (via atoms_AO_1)
! µ bounds from inner loop "IL" indices and sparse in P and ν
! ν bounds from i_atoms (via atoms_AO_2) and sparse in P and µ
bounds_IL(1:2) = [bs_env%i_ao_start_from_atom(atoms_IL(1)), &
bs_env%i_ao_end_from_atom(atoms_IL(2))]
bounds_k(1:2, 1) = [1, bs_env%n_RI]
bounds_k(1:2, 2) = [bs_env%i_ao_start_from_atom(atoms_AO_2(1)), &
bs_env%i_ao_end_from_atom(atoms_AO_2(2))]
bounds_l(1:2) = [bs_env%i_ao_start_from_atom(atoms_AO_1(1)), &
bs_env%i_ao_end_from_atom(atoms_AO_1(2))]
! µ index
CALL get_bounds_from_atoms(bounds_IL, [1, bs_env%n_atom], atoms_AO_2, &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom, &
atoms_3=atoms_IL, &
indices_3_start=bs_env%i_ao_start_from_atom, &
indices_3_end=bs_env%i_ao_end_from_atom)
CALL dbt_contract(alpha=1.0_dp, &
tensor_1=t_3c_for_G, &
tensor_2=t_G, &
beta=1.0_dp, &
tensor_3=t_M, &
contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
contract_2=[2], notcontract_2=[1], map_2=[3], &
bounds_1=bounds_IL, &
bounds_2=bounds_k, &
bounds_3=bounds_l, &
filter_eps=bs_env%eps_filter)
! P index
CALL get_bounds_from_atoms(bounds_comb(:, 1), atoms_IL, atoms_AO_2, &
bs_env%min_RI_idx_from_AO_AO_atom, &
bs_env%max_RI_idx_from_AO_AO_atom)
! ν index
CALL get_bounds_from_atoms(bounds_comb(:, 2), [1, bs_env%n_atom], atoms_IL, &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom, &
atoms_3=atoms_AO_2, &
indices_3_start=bs_env%i_ao_start_from_atom, &
indices_3_end=bs_env%i_ao_end_from_atom)
! λ index
bounds_ao_1(1:2) = [bs_env%i_ao_start_from_atom(atoms_AO_1(1)), &
bs_env%i_ao_end_from_atom(atoms_AO_1(2))]
IF (bounds_IL(1) > bounds_IL(2) .OR. bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
flop = 0_int_8
ELSE
CALL dbt_contract(alpha=1.0_dp, &
tensor_1=t_3c_for_G, &
tensor_2=t_G, &
beta=1.0_dp, &
tensor_3=t_M, &
contract_1=[3], notcontract_1=[1, 2], map_1=[1, 2], &
contract_2=[2], notcontract_2=[1], map_2=[3], &
bounds_1=bounds_IL, &
bounds_2=bounds_comb, &
bounds_3=bounds_ao_1, &
flop=flop, &
filter_eps=bs_env%eps_filter, &
unit_nr=bs_env%unit_nr_contract, &
log_verbose=bs_env%print_contract_verbose)
END IF
CALL dbt_clear(t_3c_for_G)
@ -705,13 +779,11 @@ CONTAINS
min_dist_AO_atoms = 1.0E5_dp
DO atom_1 = atoms_1(1), atoms_1(2)
DO atom_2 = atoms_2(1), atoms_2(2)
rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
min_dist_AO_atoms = MIN(min_dist_AO_atoms, abs_rab)
END DO
END DO
@ -956,9 +1028,9 @@ CONTAINS
CALL fm_read(fm_W_MIC_time(i_t), bs_env, bs_env%W_time_name, i_t)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I5,A,I3,A,F10.1,A)') &
'Read W^MIC(iτ) from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
END DO
@ -972,9 +1044,9 @@ CONTAINS
t1 = m_walltime()
CALL fm_read(bs_env%fm_W_MIC_freq_zero, bs_env, "W_freq_rtp", 0)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I3,A,I3,A,F10.1,A)') &
'Read W^MIC(f=0) from file for freq. point ', 1, ' /', 1, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
END IF
@ -1047,10 +1119,10 @@ CONTAINS
DEALLOCATE (fm_V_kp)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I12,A,I3,A,F10.1,A)') &
'Computed W(iτ,k) for k-point batch', &
ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
END DO ! ikp_batch
@ -1086,10 +1158,10 @@ CONTAINS
CALL fm_write(bs_env%fm_W_MIC_freq_zero, 0, "W_freq_rtp", qs_env)
! Report calculation
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I11,A,I3,A,F10.1,A)') &
'Computed W(f=0,k) for k-point batch', &
1, ' /', 1, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
END IF
@ -1701,13 +1773,13 @@ CONTAINS
j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
! 4. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
! 5. M_νσQ(iτ) = sum_P (νσ|P) (M^-1(k=0)*V^tr(k=0)*M^-1(k=0))_PQ(iτ)
! 5. M_Qνσ(iτ) = sum_P (νσ|P) (M^-1(k=0)*V^tr(k=0)*M^-1(k=0))_QP(iτ)
CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_V, t_2c_V)
! 6. tensor operations with D and computation of Σ^x
! Σ^x_λσ(k=0) = sum_νQ M_νσQ(iτ) sum_µ (λµ|Q) D_µν
! Σ^x_λσ(k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) D_νµ
CALL contract_to_Sigma(t_2c_D, t_3c_x_V, t_2c_Sigma_x, i_atoms, j_atoms, &
qs_env, bs_env, occ=.TRUE., vir=.FALSE., clear_W=.TRUE.)
qs_env, bs_env, occ=.TRUE., vir=.FALSE.)
END DO ! j_atoms
END DO ! i_atoms
@ -1723,7 +1795,7 @@ CONTAINS
END DO ! ispin
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,T58,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
'Computed Σ^x(k=0),', ' Execution time', m_walltime() - t1, ' s'
WRITE (bs_env%unit_nr, '(A)') ' '
END IF
@ -1786,9 +1858,9 @@ CONTAINS
CALL copy_fm_to_dbcsr(bs_env%fm_work_mo(1), mat_Sigma_neg_tau(i_t, ispin)%matrix, &
keep_sparsity=.FALSE.)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F7.1,A)') 'Read Σ^c(iτ,k=0) ', &
WRITE (bs_env%unit_nr, '(T2,2A,I3,A,I3,A,F10.1,A)') 'Read Σ^c(iτ,k=0) ', &
'from file for time point ', i_t, ' /', bs_env%num_time_freq_points, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
CYCLE
@ -1819,21 +1891,28 @@ CONTAINS
j_atoms = bs_env%j_atom_intervals(1:2, j_intval_idx)
IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) CYCLE
bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
! Do that only after first timestep to avoid skips due to vanishing G
! caused by gaps
IF (i_t == 2) THEN
bs_env%n_skip_sigma = bs_env%n_skip_sigma + 1
END IF
CYCLE
END IF
! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
! 2. tensor operation M_νσQ(iτ) = sum_P (νσ|P) W^MIC_PQ(iτ)
! 2. tensor operation M_Qνσ(iτ) = sum_P (νσ|P) W^MIC_QP(iτ)
CALL compute_3c_and_contract_W(qs_env, bs_env, i_atoms, j_atoms, t_3c_x_W, t_2c_W)
! 3. Σ_λσ(iτ,k=0) = sum_νQ M_νσQ(iτ) sum_µ (λµ|Q) G^occ_µν(i|τ|) for τ < 0
! (recall M_νσQ(iτ) = M_νσQ(-iτ) because W^MIC_PQ(iτ) = W^MIC_PQ(-iτ) )
! 3. Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^occ_νµ(i|τ|) for τ < 0
! (recall M_Qνσ(iτ) = M_Qνσ(-iτ) because W^MIC_PQ(iτ) = W^MIC_PQ(-iτ) )
CALL contract_to_Sigma(t_2c_Gocc, t_3c_x_W, t_2c_Sigma_neg_tau, i_atoms, j_atoms, &
qs_env, bs_env, occ=.TRUE., vir=.FALSE., clear_W=.FALSE., &
qs_env, bs_env, occ=.TRUE., vir=.FALSE., &
can_skip=bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx))
! Σ_λσ(iτ,k=0) = sum_νQ M_νσQ(iτ) sum_µ (λµ|Q) G^vir_µν(i|τ|) for τ > 0
! Σ_λσ(iτ,k=0) = sum_νQ M_Qνσ(iτ) sum_µ (Qλ|µ) G^vir_νµ(i|τ|) for τ > 0
CALL contract_to_Sigma(t_2c_Gvir, t_3c_x_W, t_2c_Sigma_pos_tau, i_atoms, j_atoms, &
qs_env, bs_env, occ=.FALSE., vir=.TRUE., clear_W=.TRUE., &
qs_env, bs_env, occ=.FALSE., vir=.TRUE., &
can_skip=bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx))
END DO ! j_atoms
@ -1852,9 +1931,9 @@ CONTAINS
bs_env%Sigma_n_name, bs_env%fm_work_mo(1), qs_env)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F7.1,A)') &
WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F10.1,A)') &
'Computed Σ^c(iτ,k=0) for time point ', i_t, ' /', bs_env%num_time_freq_points, &
', Execution time', m_walltime() - t1, ' s'
', Execution time', m_walltime() - t1, ' s'
END IF
END DO ! ispin
@ -1949,7 +2028,9 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_3c_and_contract_W'
INTEGER :: handle, RI_intval_idx
INTEGER, DIMENSION(2) :: bounds_j, RI_atoms
INTEGER(KIND=int_8) :: flop
INTEGER, DIMENSION(2) :: bounds_P, bounds_Q, RI_atoms
INTEGER, DIMENSION(2, 2) :: bounds_ao
TYPE(dbt_type) :: t_3c_for_W, t_3c_x_W_tmp
CALL timeset(routineN, handle)
@ -1957,17 +2038,48 @@ CONTAINS
CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_x_W_tmp)
CALL dbt_create(bs_env%t_RI__AO_AO, t_3c_for_W)
bounds_j(1:2) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
! final layout will be: M_Qνσ(iτ) = sum_P (P|νσ) W^MIC_QP(iτ)
! Bounds:
! "AO"
! -> ν (AO_1 in compute_3c_integrals) bounds from i_atoms and sparse in σ and P
! -> σ (AO_2 in compute_3c_integrals) sparse in ν and P
! Q bounds from j_atoms
! P bounds from inner loop indices and sparse in ν and σ
bounds_Q(1:2) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
bs_env%i_RI_end_from_atom(j_atoms(2))]
DO RI_intval_idx = 1, bs_env%n_intervals_inner_loop_atoms
RI_atoms = bs_env%inner_loop_atom_intervals(1:2, RI_intval_idx)
! 1. compute 3-center integrals (µν|P) ("|": truncated Coulomb operator)
CALL get_bounds_from_atoms(bounds_P, i_atoms, [1, bs_env%n_atom], &
bs_env%min_RI_idx_from_AO_AO_atom, &
bs_env%max_RI_idx_from_AO_AO_atom, &
atoms_3=RI_atoms, &
indices_3_start=bs_env%i_RI_start_from_atom, &
indices_3_end=bs_env%i_RI_end_from_atom)
! σ
CALL get_bounds_from_atoms(bounds_ao(:, 2), RI_atoms, i_atoms, &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom)
! ν
CALL get_bounds_from_atoms(bounds_ao(:, 1), RI_atoms, [1, bs_env%n_atom], &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom, &
atoms_3=i_atoms, &
indices_3_start=bs_env%i_ao_start_from_atom, &
indices_3_end=bs_env%i_ao_end_from_atom)
IF (bounds_P(1) > bounds_P(2) .OR. bounds_ao(1, 2) > bounds_ao(2, 2)) THEN
CYCLE
END IF
! 1. compute 3-center integrals (P|µν) ("|": truncated Coulomb operator)
CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_W, &
atoms_AO_1=i_atoms, atoms_RI=RI_atoms)
! 2. tensor operation M_νσQ(iτ) = sum_P (νσ|P) W^MIC_PQ(iτ)
! 2. tensor operation M_Qνσ(iτ) = sum_P W^MIC_QP(iτ) (P|νσ)
CALL dbt_contract(alpha=1.0_dp, &
tensor_1=t_2c_W, &
tensor_2=t_3c_for_W, &
@ -1975,8 +2087,14 @@ CONTAINS
tensor_3=t_3c_x_W_tmp, &
contract_1=[2], notcontract_1=[1], map_1=[1], &
contract_2=[1], notcontract_2=[2, 3], map_2=[2, 3], &
bounds_2=bounds_j, &
filter_eps=bs_env%eps_filter)
bounds_1=bounds_P, &
bounds_2=bounds_Q, &
bounds_3=bounds_ao, &
flop=flop, &
move_data=.FALSE., &
filter_eps=bs_env%eps_filter, &
unit_nr=bs_env%unit_nr_contract, &
log_verbose=bs_env%print_contract_verbose)
END DO ! RI_atoms
@ -2001,23 +2119,24 @@ CONTAINS
!> \param bs_env ...
!> \param occ ...
!> \param vir ...
!> \param clear_W ...
!> \param can_skip ...
! **************************************************************************************************
SUBROUTINE contract_to_Sigma(t_2c_G, t_3c_x_W, t_2c_Sigma, i_atoms, j_atoms, qs_env, bs_env, &
occ, vir, clear_W, can_skip)
occ, vir, can_skip)
TYPE(dbt_type) :: t_2c_G, t_3c_x_W, t_2c_Sigma
INTEGER, DIMENSION(2) :: i_atoms, j_atoms
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
LOGICAL :: occ, vir, clear_W
LOGICAL :: occ, vir
LOGICAL, OPTIONAL :: can_skip
CHARACTER(LEN=*), PARAMETER :: routineN = 'contract_to_Sigma'
INTEGER :: handle, inner_loop_atoms_interval_index
INTEGER(KIND=int_8) :: flop
INTEGER, DIMENSION(2) :: bounds_i, IL_atoms
INTEGER, DIMENSION(2) :: bounds_lambda, bounds_mu, bounds_nu, &
bounds_sigma, IL_atoms
INTEGER, DIMENSION(2, 2) :: bounds_comb
REAL(KIND=dp) :: sign_Sigma
TYPE(dbt_type) :: t_3c_for_G, t_3c_x_G, t_3c_x_G_2
@ -2031,12 +2150,49 @@ CONTAINS
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G)
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_x_G_2)
bounds_i(1:2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
bs_env%i_ao_end_from_atom(i_atoms(2))]
! Here, in the first step e.g., is computed: N_Qλν = sum_µ (Qλ|µ) G_νµ
! Afterwards e.g., is computed: Σ_λσ = sum_νQ M_Qνσ N_Qνλ (after reordering)
! Bounds:
! "comb" (combined index)
! -> Q bounds from j_atoms and sparse in λ
! -> λ (AO_1 in compute_3c_integrals) sparse in Q and µ
! µ (AO_2 in compute_3c_integrals) bounds from inner loop "IL" indices and sparse in Q and λ
! ν bounds from i_atoms
! σ sparse in ν
! ν
bounds_nu(1:2) = [bs_env%i_ao_start_from_atom(i_atoms(1)), &
bs_env%i_ao_end_from_atom(i_atoms(2))]
DO inner_loop_atoms_interval_index = 1, bs_env%n_intervals_inner_loop_atoms
IL_atoms = bs_env%inner_loop_atom_intervals(1:2, inner_loop_atoms_interval_index)
! µ
CALL get_bounds_from_atoms(bounds_mu, j_atoms, [1, bs_env%n_atom], &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom, &
atoms_3=IL_atoms, &
indices_3_start=bs_env%i_ao_start_from_atom, &
indices_3_end=bs_env%i_ao_end_from_atom)
! Q
CALL get_bounds_from_atoms(bounds_comb(:, 1), IL_atoms, [1, bs_env%n_atom], &
bs_env%min_RI_idx_from_AO_AO_atom, &
bs_env%max_RI_idx_from_AO_AO_atom, &
atoms_3=j_atoms, &
indices_3_start=bs_env%i_RI_start_from_atom, &
indices_3_end=bs_env%i_RI_end_from_atom)
! λ
CALL get_bounds_from_atoms(bounds_comb(:, 2), j_atoms, IL_atoms, &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom)
IF (bounds_mu(1) > bounds_mu(2) .OR. bounds_comb(1, 1) > bounds_comb(2, 1) .OR. &
bounds_comb(1, 2) > bounds_comb(2, 2)) THEN
CYCLE
END IF
CALL compute_3c_integrals(qs_env, bs_env, t_3c_for_G, &
atoms_RI=j_atoms, atoms_AO_2=IL_atoms)
@ -2047,21 +2203,49 @@ CONTAINS
tensor_3=t_3c_x_G, &
contract_1=[2], notcontract_1=[1], map_1=[3], &
contract_2=[3], notcontract_2=[1, 2], map_2=[1, 2], &
bounds_2=bounds_i, &
filter_eps=bs_env%eps_filter)
bounds_1=bounds_mu, &
bounds_2=bounds_nu, &
bounds_3=bounds_comb, &
flop=flop, &
move_data=.FALSE., &
filter_eps=bs_env%eps_filter, &
unit_nr=bs_env%unit_nr_contract, &
log_verbose=bs_env%print_contract_verbose)
END DO ! IL_atoms
! Reordering: N_Qλν -> N_Qνλ
CALL dbt_copy(t_3c_x_G, t_3c_x_G_2, order=[1, 3, 2], move_data=.TRUE.)
CALL dbt_contract(alpha=sign_Sigma, &
tensor_1=t_3c_x_W, &
tensor_2=t_3c_x_G_2, &
beta=1.0_dp, &
tensor_3=t_2c_Sigma, &
contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
filter_eps=bs_env%eps_filter, move_data=clear_W, flop=flop)
! Here, the last contraction is done, e.g., Σ_λσ = sum_νQ M_Qνσ N_Qνλ
! Bounds as above, new "comb" with upper ingredients
bounds_comb(1:2, 1) = [bs_env%i_RI_start_from_atom(j_atoms(1)), &
bs_env%i_RI_end_from_atom(j_atoms(2))]
bounds_comb(1:2, 2) = bounds_nu(1:2)
CALL get_bounds_from_atoms(bounds_lambda, j_atoms, [1, bs_env%n_atom], &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom)
CALL get_bounds_from_atoms(bounds_sigma, [1, bs_env%n_atom], i_atoms, &
bs_env%min_AO_idx_from_RI_AO_atom, &
bs_env%max_AO_idx_from_RI_AO_atom)
IF (bounds_sigma(1) > bounds_sigma(2) .OR. bounds_lambda(1) > bounds_lambda(2)) THEN
flop = 0_int_8
ELSE
CALL dbt_contract(alpha=sign_Sigma, &
tensor_1=t_3c_x_W, &
tensor_2=t_3c_x_G_2, &
beta=1.0_dp, &
tensor_3=t_2c_Sigma, &
contract_1=[1, 2], notcontract_1=[3], map_1=[1], &
contract_2=[1, 2], notcontract_2=[3], map_2=[2], &
bounds_1=bounds_comb, &
bounds_2=bounds_sigma, &
bounds_3=bounds_lambda, &
filter_eps=bs_env%eps_filter, move_data=.FALSE., flop=flop, &
unit_nr=bs_env%unit_nr_contract, &
log_verbose=bs_env%print_contract_verbose)
END IF
IF (PRESENT(can_skip)) THEN
IF (flop == 0_int_8) can_skip = .TRUE.
@ -2123,26 +2307,23 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'print_skipping'
INTEGER :: handle, i_intval_idx, j_intval_idx, &
n_skip
INTEGER :: handle, n_pairs
CALL timeset(routineN, handle)
n_skip = 0
n_pairs = bs_env%n_intervals_i*bs_env%n_intervals_j
DO i_intval_idx = 1, bs_env%n_intervals_i
DO j_intval_idx = 1, bs_env%n_intervals_j
IF (bs_env%skip_Sigma_occ(i_intval_idx, j_intval_idx) .AND. &
bs_env%skip_Sigma_vir(i_intval_idx, j_intval_idx)) THEN
n_skip = n_skip + 1
END IF
END DO
END DO
CALL bs_env%para_env_tensor%sum(bs_env%n_skip_sigma)
CALL bs_env%para_env_tensor%sum(bs_env%n_skip_chi)
CALL bs_env%para_env_tensor%sum(n_pairs)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
'Sparsity of Σ^c(iτ,k=0): Percentage of skipped atom pairs:', &
REAL(100*n_skip, KIND=dp)/REAL(i_intval_idx*j_intval_idx, KIND=dp), ' %'
REAL(100*bs_env%n_skip_sigma, KIND=dp)/REAL(n_pairs, KIND=dp), ' %'
WRITE (bs_env%unit_nr, '(T2,A,T74,F7.1,A)') &
'Sparsity of χ(iτ,k=0): Percentage of skipped atom pairs:', &
REAL(100*bs_env%n_skip_chi, KIND=dp)/REAL(n_pairs, KIND=dp), ' %'
END IF
CALL timestop(handle)
@ -2427,4 +2608,52 @@ CONTAINS
END SUBROUTINE fm_Gamma_ao_to_cfm_ikp_mo
! **************************************************************************************************
!> \brief Computes bounds (AO or RI) for given atom intervals atoms_1 and atoms_2 from indices_min
!> and indices_max and returns them in bounds_out.
!> In case, atoms_3 and indices_3 are given, the bounds are computed as the intersection
!> \param bounds_out Bounds to be computed
!> \param atoms_1 First atom interval
!> \param atoms_2 Second atom interval
!> \param indices_min Minimum indices for each atom pair (typically from bs_env,
!> computed in get_i_j_atom_ranges in gw_utils.F, e.g. bs_env%min_RI_idx_from_AO_AO_atom)
!> \param indices_max Maximum indices for each atom pair (typically from bs_env,
!> computed in get_i_j_atom_ranges in gw_utils.F)
!> \param atoms_3 (Optional) Third atom interval for intersection
!> \param indices_3_start (Optional) Indices for third atom interval for intersection
!> \param indices_3_end (Optional) Indices for third atom interval for intersection
! **************************************************************************************************
SUBROUTINE get_bounds_from_atoms(bounds_out, atoms_1, atoms_2, indices_min, indices_max, &
atoms_3, indices_3_start, indices_3_end)
INTEGER, DIMENSION(2), INTENT(OUT) :: bounds_out
INTEGER, DIMENSION(2), INTENT(IN) :: atoms_1, atoms_2
INTEGER, DIMENSION(:, :) :: indices_min, indices_max
INTEGER, DIMENSION(2), INTENT(IN), OPTIONAL :: atoms_3
INTEGER, DIMENSION(:), OPTIONAL :: indices_3_start, indices_3_end
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_bounds_from_atoms'
INTEGER :: handle, i_at, j_at
CALL timeset(routineN, handle)
bounds_out(1) = HUGE(0)
bounds_out(2) = -1
!Loop over all atoms in the two intervals and find min/max indices
DO i_at = atoms_1(1), atoms_1(2)
DO j_at = atoms_2(1), atoms_2(2)
bounds_out(1) = MIN(bounds_out(1), indices_min(i_at, j_at))
bounds_out(2) = MAX(bounds_out(2), indices_max(i_at, j_at))
END DO
END DO
IF (PRESENT(atoms_3) .AND. PRESENT(indices_3_start) .AND. PRESENT(indices_3_end)) THEN
bounds_out(1) = MAX(bounds_out(1), indices_3_start(atoms_3(1)))
bounds_out(2) = MIN(bounds_out(2), indices_3_end(atoms_3(2)))
END IF
CALL timestop(handle)
END SUBROUTINE get_bounds_from_atoms
END MODULE gw_large_cell_gamma

126
src/gw_utils.F
View file

@ -289,7 +289,14 @@ CONTAINS
CALL section_vals_val_get(gw_sec, "KPOINTS_W", i_vals=bs_env%nkp_grid_chi_eps_W_input)
CALL section_vals_val_get(gw_sec, "HEDIN_SHIFT", l_val=bs_env%do_hedin_shift)
CALL section_vals_val_get(gw_sec, "FREQ_MAX_FIT", r_val=bs_env%freq_max_fit)
CALL section_vals_val_get(gw_sec, "PRINT%PRINT_DBT_CONTRACT", l_val=bs_env%print_contract)
CALL section_vals_val_get(gw_sec, "PRINT%PRINT_DBT_CONTRACT_VERBOSE", l_val=bs_env%print_contract_verbose)
IF (bs_env%print_contract) THEN
bs_env%unit_nr_contract = bs_env%unit_nr
ELSE
bs_env%unit_nr_contract = 0
END IF
CALL timestop(handle)
END SUBROUTINE read_gw_input_parameters
@ -786,26 +793,33 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'check_sparsity_3c'
INTEGER :: handle, n_atom_step, RI_atom
INTEGER(int_8) :: mem, non_zero_elements_sum, nze
INTEGER(int_8) :: non_zero_elements_sum, nze
REAL(dp) :: max_dist_AO_atoms, occ, occupation_sum
REAL(KIND=dp) :: t1, t2
TYPE(dbt_type) :: t_3c_global
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_global_array
TYPE(neighbor_list_3c_type) :: nl_3c_global
!TYPE(dbt_type) :: t_3c_global
!TYPE(neighbor_list_3c_type) :: nl_3c_global
CALL timeset(routineN, handle)
! check the sparsity of 3c integral tensor (µν|P); calculate maximum distance between
! AO atoms µ, ν where at least a single integral (µν|P) is larger than the filter threshold
CALL create_3c_t(t_3c_global, bs_env%para_env, "(RI AO | AO)", [1, 2], [3], &
bs_env%sizes_RI, bs_env%sizes_AO, &
create_nl_3c=.TRUE., nl_3c=nl_3c_global, qs_env=qs_env)
CALL m_memory(mem)
CALL bs_env%para_env%max(mem)
ALLOCATE (t_3c_global_array(1, 1))
CALL dbt_create(t_3c_global, t_3c_global_array(1, 1))
CALL dbt_create(bs_env%t_RI_AO__AO, t_3c_global_array(1, 1))
! Allocate arrays to store min/max indices for overlap with other AO/RI functions on each atom
! (Filled during loop)
ALLOCATE (bs_env%min_RI_idx_from_AO_AO_atom(bs_env%n_atom, bs_env%n_atom))
ALLOCATE (bs_env%max_RI_idx_from_AO_AO_atom(bs_env%n_atom, bs_env%n_atom))
ALLOCATE (bs_env%min_AO_idx_from_RI_AO_atom(bs_env%n_atom, bs_env%n_atom))
ALLOCATE (bs_env%max_AO_idx_from_RI_AO_atom(bs_env%n_atom, bs_env%n_atom))
bs_env%min_RI_idx_from_AO_AO_atom(:, :) = bs_env%n_RI
bs_env%max_RI_idx_from_AO_AO_atom(:, :) = 1
bs_env%min_AO_idx_from_RI_AO_atom(:, :) = bs_env%n_AO
bs_env%max_AO_idx_from_RI_AO_atom(:, :) = 1
CALL bs_env%para_env%sync()
t1 = m_walltime()
@ -820,7 +834,7 @@ CONTAINS
CALL build_3c_integrals(t_3c_global_array, &
bs_env%eps_filter, &
qs_env, &
nl_3c_global, &
bs_env%nl_3c, &
int_eps=bs_env%eps_filter, &
basis_i=bs_env%basis_set_RI, &
basis_j=bs_env%basis_set_AO, &
@ -839,21 +853,26 @@ CONTAINS
CALL get_max_dist_AO_atoms(t_3c_global_array(1, 1), max_dist_AO_atoms, qs_env)
! Extract indices per block
CALL get_i_j_atom_ranges(t_3c_global_array(1, 1), bs_env)
CALL dbt_clear(t_3c_global_array(1, 1))
END DO
t2 = m_walltime()
CALL bs_env%para_env%min(bs_env%min_RI_idx_from_AO_AO_atom)
CALL bs_env%para_env%max(bs_env%max_RI_idx_from_AO_AO_atom)
CALL bs_env%para_env%min(bs_env%min_AO_idx_from_RI_AO_atom)
CALL bs_env%para_env%max(bs_env%max_AO_idx_from_RI_AO_atom)
bs_env%occupation_3c_int = occupation_sum
bs_env%max_dist_AO_atoms = max_dist_AO_atoms
CALL dbt_destroy(t_3c_global)
CALL dbt_destroy(t_3c_global_array(1, 1))
DEALLOCATE (t_3c_global_array)
CALL neighbor_list_3c_destroy(nl_3c_global)
IF (bs_env%unit_nr > 0) THEN
WRITE (bs_env%unit_nr, '(T2,A)') ''
WRITE (bs_env%unit_nr, '(T2,A,F27.1,A)') &
@ -870,6 +889,64 @@ CONTAINS
END SUBROUTINE check_sparsity_3c
! **************************************************************************************************
!> \brief ...
!> \param t_3c ...
!> \param bs_env ...
! **************************************************************************************************
SUBROUTINE get_i_j_atom_ranges(t_3c, bs_env)
TYPE(dbt_type) :: t_3c
TYPE(post_scf_bandstructure_type), POINTER :: bs_env
CHARACTER(LEN=*), PARAMETER :: routineN = 'get_i_j_atom_ranges'
INTEGER :: handle, idx_AO_end, idx_AO_start, &
idx_RI_end, idx_RI_start
INTEGER, DIMENSION(3) :: atom_ind
TYPE(dbt_iterator_type) :: iter
CALL timeset(routineN, handle)
! Loop over blocks in 3c, for given min_atom: RI_min/max index from min_atom
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(t_3c, bs_env) &
!$OMP PRIVATE(iter, atom_ind, &
!$OMP idx_RI_start, idx_RI_end, idx_AO_start, idx_AO_end)
CALL dbt_iterator_start(iter, t_3c)
DO WHILE (dbt_iterator_blocks_left(iter))
CALL dbt_iterator_next_block(iter, atom_ind)
! Pre-fetch indices to avoid referencing 'bs_env' twice inside the ATOMIC blocks
idx_RI_start = bs_env%i_RI_start_from_atom(atom_ind(1))
idx_RI_end = bs_env%i_RI_end_from_atom(atom_ind(1))
idx_AO_start = bs_env%i_ao_start_from_atom(atom_ind(2))
idx_AO_end = bs_env%i_ao_end_from_atom(atom_ind(2))
! Update values safely inside ATOMIC blocks, otherwise race conditions occur
!$OMP ATOMIC UPDATE
bs_env%min_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)) = &
MIN(bs_env%min_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)), idx_RI_start)
!$OMP ATOMIC UPDATE
bs_env%max_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)) = &
MAX(bs_env%max_RI_idx_from_AO_AO_atom(atom_ind(2), atom_ind(3)), idx_RI_end)
!$OMP ATOMIC UPDATE
bs_env%min_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)) = &
MIN(bs_env%min_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)), idx_AO_start)
!$OMP ATOMIC UPDATE
bs_env%max_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)) = &
MAX(bs_env%max_AO_idx_from_RI_AO_atom(atom_ind(1), atom_ind(3)), idx_AO_end)
END DO
CALL dbt_iterator_stop(iter)
!$OMP END PARALLEL
CALL timestop(handle)
END SUBROUTINE get_i_j_atom_ranges
! **************************************************************************************************
!> \brief ...
!> \param bs_env ...
@ -1036,20 +1113,22 @@ CONTAINS
NULLIFY (cell, particle_set, para_env)
CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, para_env=para_env)
! IMPORTANT: Use thread-local copy for max_dist_AO_atoms via REDUCTION to avoid race conditions
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(t_3c_int, max_dist_AO_atoms, num_cells, index_to_cell, particle_set, cell) &
!$OMP PRIVATE(iter,atom_ind,rab, abs_rab, atom_1, atom_2)
!$OMP SHARED(t_3c_int, num_cells, index_to_cell, particle_set, cell) &
!$OMP PRIVATE(iter, atom_ind, rab, abs_rab, atom_1, atom_2) &
!$OMP REDUCTION(MAX:max_dist_AO_atoms)
CALL dbt_iterator_start(iter, t_3c_int)
DO WHILE (dbt_iterator_blocks_left(iter))
CALL dbt_iterator_next_block(iter, atom_ind)
atom_1 = atom_ind(2)
atom_2 = atom_ind(3)
rab = pbc(particle_set(atom_1)%r(1:3), particle_set(atom_2)%r(1:3), cell)
abs_rab = SQRT(rab(1)**2 + rab(2)**2 + rab(3)**2)
! Reduction takes care of using a thread-local copy
max_dist_AO_atoms = MAX(max_dist_AO_atoms, abs_rab)
END DO
@ -1116,6 +1195,11 @@ CONTAINS
ALLOCATE (bs_env%skip_Sigma_vir(n_intervals_i, n_intervals_j))
bs_env%skip_Sigma_occ(:, :) = .FALSE.
bs_env%skip_Sigma_vir(:, :) = .FALSE.
bs_env%n_skip_chi = 0
ALLOCATE (bs_env%skip_chi(n_intervals_i, n_intervals_j))
bs_env%skip_chi(:, :) = .FALSE.
bs_env%n_skip_sigma = 0
! choose atomic range for µ and σ ("inner loop (IL) atom") in
! M_λνP(iτ) = sum_µ (µν|P) G^occ_µλ(i|τ|,k=0)
@ -1288,7 +1372,6 @@ CONTAINS
INTEGER :: color_sub, dummy_1, dummy_2, handle, &
num_pe, num_t_groups, u
INTEGER(KIND=int_8) :: mem
TYPE(mp_para_env_type), POINTER :: para_env
CALL timeset(routineN, handle)
@ -1327,9 +1410,6 @@ CONTAINS
dummy_1, dummy_2, color_sub, bs_env)
END DO
CALL m_memory(mem)
CALL bs_env%para_env%max(mem)
u = bs_env%unit_nr
IF (u > 0) THEN
WRITE (u, '(T2,A,I47)') 'Group size for tensor operations', bs_env%group_size_tensor
@ -1593,6 +1673,8 @@ CONTAINS
WRITE (u, '(T2,A,ES27.1)') 'Input: Filter threshold for sparse tensor operations', &
bs_env%eps_filter
WRITE (u, "(T2,A,L55)") 'Input: Apply Hedin shift', bs_env%do_hedin_shift
WRITE (u, '(T2,A,F37.1,A)') 'Input: Available memory per MPI process', &
bs_env%input_memory_per_proc_GB, ' GB'
END IF
CALL timestop(handle)

18
src/input_cp2k_properties_dft.F
View file

@ -2412,7 +2412,7 @@ CONTAINS
NULLIFY (subsection, print_key)
CALL section_create(subsection, __LOCATION__, name="PRINT", &
description="Printing of GW restarts.", &
n_keywords=0, n_subsections=1, repeats=.FALSE.)
n_keywords=2, n_subsections=1, repeats=.FALSE.)
CALL cp_print_key_section_create(print_key, __LOCATION__, "RESTART", &
description="Controls the printing of restart files "// &
"for χ, W, Σ.", &
@ -2421,6 +2421,22 @@ CONTAINS
CALL section_add_subsection(subsection, print_key)
CALL section_release(print_key)
CALL keyword_create(keyword, __LOCATION__, name="PRINT_DBT_CONTRACT", &
description="Prints information of contraction routines.", &
usage="PRINT_DBT_CONTRACT", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, name="PRINT_DBT_CONTRACT_VERBOSE", &
description="Prints verbose information of contraction routines.", &
usage="PRINT_DBT_CONTRACT_VERBOSE", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(subsection, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, subsection)
CALL section_release(subsection)

38
src/mpiwrap/message_passing.F
View file

@ -281,28 +281,42 @@ MODULE message_passing
mp_sum_partial_cm, mp_sum_partial_zm
PROCEDURE, PRIVATE, PASS(comm), NON_OVERRIDABLE :: mp_max_i, mp_max_iv, &
mp_max_l, mp_max_lv, mp_max_r, mp_max_rv, &
mp_max_d, mp_max_dv, mp_max_c, mp_max_cv, &
mp_max_z, mp_max_zv, mp_max_root_i, mp_max_root_l, &
mp_max_im, &
mp_max_l, mp_max_lv, mp_max_lm, &
mp_max_r, mp_max_rv, mp_max_rm, &
mp_max_d, mp_max_dv, mp_max_dm, &
mp_max_c, mp_max_cv, mp_max_cm, &
mp_max_z, mp_max_zv, mp_max_zm, &
mp_max_root_i, mp_max_root_l, &
mp_max_root_r, mp_max_root_d, mp_max_root_c, mp_max_root_z, &
mp_max_root_im, mp_max_root_lm, mp_max_root_rm, mp_max_root_dm, &
mp_max_root_cm, mp_max_root_zm
GENERIC, PUBLIC :: max => mp_max_i, mp_max_iv, &
mp_max_l, mp_max_lv, mp_max_r, mp_max_rv, &
mp_max_d, mp_max_dv, mp_max_c, mp_max_cv, &
mp_max_z, mp_max_zv, mp_max_root_i, mp_max_root_l, &
mp_max_im, &
mp_max_l, mp_max_lv, mp_max_lm, &
mp_max_r, mp_max_rv, mp_max_rm, &
mp_max_d, mp_max_dv, mp_max_dm, &
mp_max_c, mp_max_cv, mp_max_cm, &
mp_max_z, mp_max_zv, mp_max_zm, &
mp_max_root_i, mp_max_root_l, &
mp_max_root_r, mp_max_root_d, mp_max_root_c, mp_max_root_z, &
mp_max_root_im, mp_max_root_lm, mp_max_root_rm, mp_max_root_dm, &
mp_max_root_cm, mp_max_root_zm
PROCEDURE, PRIVATE, PASS(comm), NON_OVERRIDABLE :: mp_min_i, mp_min_iv, &
mp_min_l, mp_min_lv, mp_min_r, mp_min_rv, &
mp_min_d, mp_min_dv, mp_min_c, mp_min_cv, &
mp_min_z, mp_min_zv
mp_min_im, &
mp_min_l, mp_min_lv, mp_min_lm, &
mp_min_r, mp_min_rv, mp_min_rm, &
mp_min_d, mp_min_dv, mp_min_dm, &
mp_min_c, mp_min_cv, mp_min_cm, &
mp_min_z, mp_min_zv, mp_min_zm
GENERIC, PUBLIC :: min => mp_min_i, mp_min_iv, &
mp_min_l, mp_min_lv, mp_min_r, mp_min_rv, &
mp_min_d, mp_min_dv, mp_min_c, mp_min_cv, &
mp_min_z, mp_min_zv
mp_min_im, &
mp_min_l, mp_min_lv, mp_min_lm, &
mp_min_r, mp_min_rv, mp_min_rm, &
mp_min_d, mp_min_dv, mp_min_dm, &
mp_min_c, mp_min_cv, mp_min_cm, &
mp_min_z, mp_min_zv, mp_min_zm
PROCEDURE, PUBLIC, PASS(comm), NON_OVERRIDABLE :: &
mp_sum_scatter_iv, mp_sum_scatter_lv, mp_sum_scatter_rv, &

84
src/mpiwrap/message_passing.fypp
View file

@ -1710,6 +1710,48 @@
CALL mp_timestop(handle)
END SUBROUTINE mp_max_${nametype1}$v
! **************************************************************************************************
!> \brief Finds the element-wise maximum of a rank2-array with the result left on
!> all processes.
!> \param[in] msg Matrix - Find maximum among these data (input) and
!> maximum (output)
!> \param comm ...
!> \note see mp_max_${nametype1}$
! **************************************************************************************************
SUBROUTINE mp_max_${nametype1}$m(msg, comm)
${type1}$, CONTIGUOUS, INTENT(INOUT) :: msg(:, :)
CLASS(mp_comm_type), INTENT(IN) :: comm
CHARACTER(len=*), PARAMETER :: routineN = 'mp_max_${nametype1}$m'
INTEGER :: handle
#if defined(__parallel)
INTEGER, PARAMETER :: max_msg = 2**25
INTEGER :: ierr, m1, msglen, step, msglensum
#endif
CALL mp_timeset(routineN, handle)
#if defined(__parallel)
! chunk up the call so that message sizes are limited, to avoid overflows in mpich triggered in large rpa calcs
step = MAX(1, SIZE(msg, 2)/MAX(1, SIZE(msg)/max_msg))
msglensum = 0
DO m1 = LBOUND(msg, 2), UBOUND(msg, 2), step
msglen = SIZE(msg, 1)*(MIN(UBOUND(msg, 2), m1 + step - 1) - m1 + 1)
msglensum = msglensum + msglen
IF (msglen > 0) THEN
CALL mpi_allreduce(MPI_IN_PLACE, msg(LBOUND(msg, 1), m1), msglen, ${mpi_type1}$, MPI_MAX, comm%handle, ierr)
IF (ierr /= 0) CALL mp_stop(ierr, "mpi_allreduce @ "//routineN)
END IF
END DO
CALL add_perf(perf_id=3, count=1, msg_size=msglensum*${bytes1}$)
#else
MARK_USED(msg)
MARK_USED(comm)
#endif
CALL mp_timestop(handle)
END SUBROUTINE mp_max_${nametype1}$m
! **************************************************************************************************
!> \brief Finds the element-wise maximum of a vector with the result left on
!> all processes.
@ -1815,6 +1857,48 @@
CALL mp_timestop(handle)
END SUBROUTINE mp_min_${nametype1}$v
! **************************************************************************************************
!> \brief Finds the element-wise minimum of a rank2-array with the result left on
!> all processes.
!> \param[in] msg Matrix - Find maximum among these data (input) and
!> minimum (output)
!> \param comm ...
!> \note see mp_min_${nametype1}$
! **************************************************************************************************
SUBROUTINE mp_min_${nametype1}$m(msg, comm)
${type1}$, CONTIGUOUS, INTENT(INOUT) :: msg(:, :)
CLASS(mp_comm_type), INTENT(IN) :: comm
CHARACTER(len=*), PARAMETER :: routineN = 'mp_min_${nametype1}$m'
INTEGER :: handle
#if defined(__parallel)
INTEGER, PARAMETER :: max_msg = 2**25
INTEGER :: ierr, m1, msglen, step, msglensum
#endif
CALL mp_timeset(routineN, handle)
#if defined(__parallel)
! chunk up the call so that message sizes are limited, to avoid overflows in mpich triggered in large rpa calcs
step = MAX(1, SIZE(msg, 2)/MAX(1, SIZE(msg)/max_msg))
msglensum = 0
DO m1 = LBOUND(msg, 2), UBOUND(msg, 2), step
msglen = SIZE(msg, 1)*(MIN(UBOUND(msg, 2), m1 + step - 1) - m1 + 1)
msglensum = msglensum + msglen
IF (msglen > 0) THEN
CALL mpi_allreduce(MPI_IN_PLACE, msg(LBOUND(msg, 1), m1), msglen, ${mpi_type1}$, MPI_MIN, comm%handle, ierr)
IF (ierr /= 0) CALL mp_stop(ierr, "mpi_allreduce @ "//routineN)
END IF
END DO
CALL add_perf(perf_id=3, count=1, msg_size=msglensum*${bytes1}$)
#else
MARK_USED(msg)
MARK_USED(comm)
#endif
CALL mp_timestop(handle)
END SUBROUTINE mp_min_${nametype1}$m
! **************************************************************************************************
!> \brief Multiplies a set of numbers scattered across a number of processes,
!> then replicates the result.

23
src/post_scf_bandstructure_types.F
View file

@ -89,6 +89,10 @@ MODULE post_scf_bandstructure_types
i_ao_end_from_atom, &
i_RI_start_from_atom, &
i_RI_end_from_atom
INTEGER, DIMENSION(:, :), ALLOCATABLE :: min_RI_idx_from_AO_AO_atom, &
max_RI_idx_from_AO_AO_atom, &
min_AO_idx_from_RI_AO_atom, &
max_AO_idx_from_RI_AO_atom
INTEGER, DIMENSION(2) :: n_occ = -1, &
n_vir = -1
REAL(KIND=dp) :: spin_degeneracy = -1.0_dp
@ -215,14 +219,17 @@ MODULE post_scf_bandstructure_types
n_intervals_j = -1, &
n_atom_per_interval_ij = -1, &
n_intervals_inner_loop_atoms = -1, &
n_atom_per_IL_interval = -1
n_atom_per_IL_interval = -1, &
n_skip_sigma = -1, &
n_skip_chi = -1
INTEGER, DIMENSION(:, :), ALLOCATABLE :: i_atom_intervals, &
j_atom_intervals, &
inner_loop_atom_intervals, &
atoms_i_t_group, &
atoms_j_t_group
LOGICAL, DIMENSION(:, :), ALLOCATABLE :: skip_Sigma_occ, &
skip_Sigma_vir
skip_Sigma_vir, &
skip_chi
! Marek : rtbse_method
INTEGER :: rtp_method = rtp_method_bse
@ -238,7 +245,8 @@ MODULE post_scf_bandstructure_types
CHARACTER(LEN=13) :: Sigma_p_name = "Sigma_pos_tau", &
Sigma_n_name = "Sigma_neg_tau"
CHARACTER(LEN=default_string_length) :: prefix = ""
INTEGER :: unit_nr = -1
INTEGER :: unit_nr = -1, &
unit_nr_contract = -1
! parameters and data for basis sets
TYPE(gto_basis_set_p_type), &
@ -311,6 +319,10 @@ MODULE post_scf_bandstructure_types
REAL(KIND=dp), DIMENSION(:, :, :), ALLOCATABLE :: v_xc_n
TYPE(dbt_type), ALLOCATABLE, DIMENSION(:, :) :: t_3c_int
!MG: Print options
LOGICAL :: print_contract = .FALSE., &
print_contract_verbose = .FALSE.
END TYPE post_scf_bandstructure_type
CONTAINS
@ -353,12 +365,17 @@ CONTAINS
IF (ALLOCATED(bs_env%i_ao_end_from_atom)) DEALLOCATE (bs_env%i_ao_end_from_atom)
IF (ALLOCATED(bs_env%i_RI_start_from_atom)) DEALLOCATE (bs_env%i_RI_start_from_atom)
IF (ALLOCATED(bs_env%i_RI_end_from_atom)) DEALLOCATE (bs_env%i_RI_end_from_atom)
IF (ALLOCATED(bs_env%min_RI_idx_from_AO_AO_atom)) DEALLOCATE (bs_env%min_RI_idx_from_AO_AO_atom)
IF (ALLOCATED(bs_env%max_RI_idx_from_AO_AO_atom)) DEALLOCATE (bs_env%max_RI_idx_from_AO_AO_atom)
IF (ALLOCATED(bs_env%min_AO_idx_from_RI_AO_atom)) DEALLOCATE (bs_env%min_AO_idx_from_RI_AO_atom)
IF (ALLOCATED(bs_env%max_AO_idx_from_RI_AO_atom)) DEALLOCATE (bs_env%max_AO_idx_from_RI_AO_atom)
IF (ALLOCATED(bs_env%i_atom_intervals)) DEALLOCATE (bs_env%i_atom_intervals)
IF (ALLOCATED(bs_env%j_atom_intervals)) DEALLOCATE (bs_env%j_atom_intervals)
IF (ALLOCATED(bs_env%atoms_i_t_group)) DEALLOCATE (bs_env%atoms_i_t_group)
IF (ALLOCATED(bs_env%atoms_j_t_group)) DEALLOCATE (bs_env%atoms_j_t_group)
IF (ALLOCATED(bs_env%skip_Sigma_occ)) DEALLOCATE (bs_env%skip_Sigma_occ)
IF (ALLOCATED(bs_env%skip_Sigma_vir)) DEALLOCATE (bs_env%skip_Sigma_vir)
IF (ALLOCATED(bs_env%skip_chi)) DEALLOCATE (bs_env%skip_chi)
IF (ALLOCATED(bs_env%read_chi)) DEALLOCATE (bs_env%read_chi)
IF (ALLOCATED(bs_env%calc_chi)) DEALLOCATE (bs_env%calc_chi)
IF (ALLOCATED(bs_env%Sigma_c_exists)) DEALLOCATE (bs_env%Sigma_c_exists)