(git:0de0cc2)
mp2_ri_grad.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Routines to calculate gradients of RI-GPW-MP2 energy using pw
10 !> \par History
11 !> 10.2013 created [Mauro Del Ben]
12 ! **************************************************************************************************
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
16  USE cell_types, ONLY: cell_type
17  USE cp_blacs_env, ONLY: cp_blacs_env_type
18  USE cp_control_types, ONLY: dft_control_type
21  USE cp_eri_mme_interface, ONLY: cp_eri_mme_param
24  cp_fm_struct_type
25  USE cp_fm_types, ONLY: &
28  USE dbcsr_api, ONLY: &
29  dbcsr_add, dbcsr_copy, dbcsr_copy_into_existing, dbcsr_create, dbcsr_multiply, &
30  dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_transposed, dbcsr_type, &
31  dbcsr_type_no_symmetry, dbcsr_type_symmetric
32  USE input_constants, ONLY: do_eri_gpw,&
33  do_eri_mme,&
37  USE kinds, ONLY: dp
40  mp_para_env_type,&
42  mp_request_type,&
43  mp_waitall
44  USE mp2_eri, ONLY: mp2_eri_2c_integrate,&
47  mp2_eri_force
48  USE mp2_eri_gpw, ONLY: cleanup_gpw,&
53  USE mp2_types, ONLY: integ_mat_buffer_type,&
54  integ_mat_buffer_type_2d,&
55  mp2_type
56  USE parallel_gemm_api, ONLY: parallel_gemm
57  USE particle_types, ONLY: particle_type
58  USE pw_env_types, ONLY: pw_env_type
59  USE pw_poisson_types, ONLY: pw_poisson_type
60  USE pw_pool_types, ONLY: pw_pool_type
61  USE pw_types, ONLY: pw_c1d_gs_type,&
62  pw_r3d_rs_type
63  USE qs_environment_types, ONLY: get_qs_env,&
64  qs_environment_type
66  qs_force_type,&
68  USE qs_kind_types, ONLY: qs_kind_type
69  USE qs_ks_types, ONLY: qs_ks_env_type
70  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
71  USE task_list_types, ONLY: task_list_type
72  USE util, ONLY: get_limit
73  USE virial_types, ONLY: virial_type
74 #include "./base/base_uses.f90"
75 
76  IMPLICIT NONE
77 
78  PRIVATE
79 
80  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_ri_grad'
81 
82  PUBLIC :: calc_ri_mp2_nonsep
83 
84 CONTAINS
85 
86 ! **************************************************************************************************
87 !> \brief Calculate the non-separable part of the gradients and update the
88 !> Lagrangian
89 !> \param qs_env ...
90 !> \param mp2_env ...
91 !> \param para_env ...
92 !> \param para_env_sub ...
93 !> \param cell ...
94 !> \param particle_set ...
95 !> \param atomic_kind_set ...
96 !> \param qs_kind_set ...
97 !> \param mo_coeff ...
98 !> \param nmo ...
99 !> \param homo ...
100 !> \param dimen_RI ...
101 !> \param Eigenval ...
102 !> \param my_group_L_start ...
103 !> \param my_group_L_end ...
104 !> \param my_group_L_size ...
105 !> \param sab_orb_sub ...
106 !> \param mat_munu ...
107 !> \param blacs_env_sub ...
108 !> \author Mauro Del Ben
109 ! **************************************************************************************************
110  SUBROUTINE calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, particle_set, &
111  atomic_kind_set, qs_kind_set, mo_coeff, nmo, homo, dimen_RI, Eigenval, &
112  my_group_L_start, my_group_L_end, my_group_L_size, sab_orb_sub, mat_munu, &
113  blacs_env_sub)
114  TYPE(qs_environment_type), POINTER :: qs_env
115  TYPE(mp2_type) :: mp2_env
116  TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
117  TYPE(cell_type), POINTER :: cell
118  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
119  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
120  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
121  TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
122  INTEGER, INTENT(IN) :: nmo
123  INTEGER, DIMENSION(:), INTENT(IN) :: homo
124  INTEGER, INTENT(IN) :: dimen_ri
125  REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: eigenval
126  INTEGER, INTENT(IN) :: my_group_l_start, my_group_l_end, &
127  my_group_l_size
128  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
129  POINTER :: sab_orb_sub
130  TYPE(dbcsr_p_type), INTENT(INOUT) :: mat_munu
131  TYPE(cp_blacs_env_type), POINTER :: blacs_env_sub
132 
133  CHARACTER(LEN=*), PARAMETER :: routinen = 'calc_ri_mp2_nonsep'
134 
135  INTEGER :: dimen, eri_method, handle, handle2, i, &
136  ikind, ispin, itmp(2), l_counter, lll, &
137  my_p_end, my_p_size, my_p_start, nspins
138  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, natom_of_kind, &
139  virtual
140  LOGICAL :: alpha_beta, use_virial
141  REAL(kind=dp) :: cutoff_old, eps_filter, factor, &
142  factor_2c, relative_cutoff_old
143  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: e_cutoff_old
144  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: g_pq_local, g_pq_local_2
145  REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_virial
146  REAL(kind=dp), DIMENSION(:, :), POINTER :: i_tmp2
147  TYPE(cp_eri_mme_param), POINTER :: eri_param
148  TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
149  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: l1_mu_i, l2_nu_a
150  TYPE(dbcsr_p_type) :: matrix_p_munu
151  TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff_o, mo_coeff_v
152  TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:, :) :: g_p_ia
153  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_munu_local, matrix_p_munu_local
154  TYPE(dbcsr_type) :: matrix_p_munu_nosym
155  TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: lag_mu_i_1, lag_nu_a_2, matrix_p_inu
156  TYPE(dft_control_type), POINTER :: dft_control
157  TYPE(mp2_eri_force), ALLOCATABLE, DIMENSION(:) :: force_2c, force_2c_ri, force_3c_aux, &
158  force_3c_orb_mu, force_3c_orb_nu
159  TYPE(pw_c1d_gs_type) :: dvg(3), pot_g, rho_g, rho_g_copy
160  TYPE(pw_env_type), POINTER :: pw_env_sub
161  TYPE(pw_poisson_type), POINTER :: poisson_env
162  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
163  TYPE(pw_r3d_rs_type) :: psi_l, rho_r
164  TYPE(qs_force_type), DIMENSION(:), POINTER :: force, mp2_force
165  TYPE(qs_ks_env_type), POINTER :: ks_env
166  TYPE(task_list_type), POINTER :: task_list_sub
167  TYPE(virial_type), POINTER :: virial
168 
169  CALL timeset(routinen, handle)
170 
171  eri_method = mp2_env%eri_method
172  eri_param => mp2_env%eri_mme_param
173 
174  ! Find out whether we have a closed or open shell
175  nspins = SIZE(homo)
176  alpha_beta = (nspins == 2)
177 
178  dimen = nmo
179  ALLOCATE (virtual(nspins))
180  virtual(:) = dimen - homo(:)
181  eps_filter = mp2_env%mp2_gpw%eps_filter
182  ALLOCATE (mo_coeff_o(nspins), mo_coeff_v(nspins), g_p_ia(nspins, my_group_l_size))
183  DO ispin = 1, nspins
184  mo_coeff_o(ispin)%matrix => mp2_env%ri_grad%mo_coeff_o(ispin)%matrix
185  mo_coeff_v(ispin)%matrix => mp2_env%ri_grad%mo_coeff_v(ispin)%matrix
186  DO lll = 1, my_group_l_size
187  g_p_ia(ispin, lll)%matrix => mp2_env%ri_grad%G_P_ia(lll, ispin)%matrix
188  END DO
189  END DO
190  DEALLOCATE (mp2_env%ri_grad%G_P_ia)
191 
192  itmp = get_limit(dimen_ri, para_env_sub%num_pe, para_env_sub%mepos)
193  my_p_start = itmp(1)
194  my_p_end = itmp(2)
195  my_p_size = itmp(2) - itmp(1) + 1
196 
197  ALLOCATE (g_pq_local(dimen_ri, my_group_l_size))
198  g_pq_local = 0.0_dp
199  g_pq_local(my_p_start:my_p_end, :) = mp2_env%ri_grad%Gamma_PQ
200  DEALLOCATE (mp2_env%ri_grad%Gamma_PQ)
201  g_pq_local(my_p_start:my_p_end, :) = g_pq_local(my_p_start:my_p_end, :)/real(nspins, dp)
202  CALL para_env_sub%sum(g_pq_local)
203  IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
204  ALLOCATE (g_pq_local_2(dimen_ri, my_group_l_size))
205  g_pq_local_2 = 0.0_dp
206  g_pq_local_2(my_p_start:my_p_end, :) = mp2_env%ri_grad%Gamma_PQ_2
207  DEALLOCATE (mp2_env%ri_grad%Gamma_PQ_2)
208  g_pq_local_2(my_p_start:my_p_end, :) = g_pq_local_2(my_p_start:my_p_end, :)/real(nspins, dp)
209  CALL para_env_sub%sum(g_pq_local_2)
210  END IF
211 
212  ! create matrix holding the back transformation (G_P_inu)
213  ALLOCATE (matrix_p_inu(nspins))
214  DO ispin = 1, nspins
215  CALL dbcsr_create(matrix_p_inu(ispin), template=mo_coeff_o(ispin)%matrix)
216  END DO
217 
218  ! non symmetric matrix
219  CALL dbcsr_create(matrix_p_munu_nosym, template=mat_munu%matrix, &
220  matrix_type=dbcsr_type_no_symmetry)
221 
222  ! create Lagrangian matrices in mixed AO/MO formalism
223  ALLOCATE (lag_mu_i_1(nspins))
224  DO ispin = 1, nspins
225  CALL dbcsr_create(lag_mu_i_1(ispin), template=mo_coeff_o(ispin)%matrix)
226  CALL dbcsr_set(lag_mu_i_1(ispin), 0.0_dp)
227  END DO
228 
229  ALLOCATE (lag_nu_a_2(nspins))
230  DO ispin = 1, nspins
231  CALL dbcsr_create(lag_nu_a_2(ispin), template=mo_coeff_v(ispin)%matrix)
232  CALL dbcsr_set(lag_nu_a_2(ispin), 0.0_dp)
233  END DO
234 
235  ! get forces
236  NULLIFY (force, virial)
237  CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
238 
239  ! check if we want to calculate the virial
240  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
241 
242  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, natom_of_kind=natom_of_kind)
243  NULLIFY (mp2_force)
244  CALL allocate_qs_force(mp2_force, natom_of_kind)
245  DEALLOCATE (natom_of_kind)
246  CALL zero_qs_force(mp2_force)
247  mp2_env%ri_grad%mp2_force => mp2_force
248 
249  factor_2c = -4.0_dp
250  IF (mp2_env%method == ri_rpa_method_gpw) THEN
251  factor_2c = -1.0_dp
252  IF (alpha_beta) factor_2c = -2.0_dp
253  END IF
254 
255  ! prepare integral derivatives with mme method
256  IF (eri_method .EQ. do_eri_mme) THEN
257  ALLOCATE (matrix_p_munu_local(my_group_l_size))
258  ALLOCATE (mat_munu_local(my_group_l_size))
259  l_counter = 0
260  DO lll = my_group_l_start, my_group_l_end
261  l_counter = l_counter + 1
262  ALLOCATE (mat_munu_local(l_counter)%matrix)
263  CALL dbcsr_create(mat_munu_local(l_counter)%matrix, template=mat_munu%matrix, &
264  matrix_type=dbcsr_type_symmetric)
265  CALL dbcsr_copy(mat_munu_local(l_counter)%matrix, mat_munu%matrix)
266  CALL dbcsr_set(mat_munu_local(l_counter)%matrix, 0.0_dp)
267 
268  CALL g_p_transform_mo_to_ao(matrix_p_munu_local(l_counter)%matrix, matrix_p_munu_nosym, mat_munu%matrix, &
269  g_p_ia(:, l_counter), matrix_p_inu, &
270  mo_coeff_v, mo_coeff_o, eps_filter)
271  END DO
272 
273  ALLOCATE (i_tmp2(dimen_ri, my_group_l_size))
274  i_tmp2(:, :) = 0.0_dp
275  CALL mp2_eri_2c_integrate(eri_param, mp2_env%potential_parameter, para_env_sub, qs_env, &
276  basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
277  hab=i_tmp2, first_b=my_group_l_start, last_b=my_group_l_end, &
278  eri_method=eri_method, pab=g_pq_local, force_a=force_2c)
279  IF (.NOT. compare_potential_types(mp2_env%potential_parameter, mp2_env%ri_metric)) THEN
280  i_tmp2(:, :) = 0.0_dp
281  CALL mp2_eri_2c_integrate(eri_param, mp2_env%ri_metric, para_env_sub, qs_env, &
282  basis_type_a="RI_AUX", basis_type_b="RI_AUX", &
283  hab=i_tmp2, first_b=my_group_l_start, last_b=my_group_l_end, &
284  eri_method=eri_method, pab=g_pq_local_2, force_a=force_2c_ri)
285  END IF
286  DEALLOCATE (i_tmp2)
287 
288  CALL mp2_eri_3c_integrate(eri_param, mp2_env%ri_metric, para_env_sub, qs_env, &
289  first_c=my_group_l_start, last_c=my_group_l_end, mat_ab=mat_munu_local, &
290  basis_type_a="ORB", basis_type_b="ORB", basis_type_c="RI_AUX", &
291  sab_nl=sab_orb_sub, eri_method=eri_method, &
292  pabc=matrix_p_munu_local, &
293  force_a=force_3c_orb_mu, force_b=force_3c_orb_nu, force_c=force_3c_aux)
294 
295  l_counter = 0
296  DO lll = my_group_l_start, my_group_l_end
297  l_counter = l_counter + 1
298  DO ispin = 1, nspins
299  CALL dbcsr_multiply("N", "T", 1.0_dp, mo_coeff_v(ispin)%matrix, g_p_ia(ispin, l_counter)%matrix, &
300  0.0_dp, matrix_p_inu(ispin), filter_eps=eps_filter)
301  END DO
302 
303  ! The matrices of G_P_ia are deallocated here
304  CALL update_lagrangian(mat_munu_local(l_counter)%matrix, matrix_p_inu, lag_mu_i_1, &
305  g_p_ia(:, l_counter), mo_coeff_o, lag_nu_a_2, &
306  eps_filter)
307  END DO
308 
309  DO ikind = 1, SIZE(force)
310  mp2_force(ikind)%mp2_non_sep(:, :) = factor_2c*force_2c(ikind)%forces(:, :) + &
311  force_3c_orb_mu(ikind)%forces(:, :) + &
312  force_3c_orb_nu(ikind)%forces(:, :) + &
313  force_3c_aux(ikind)%forces(:, :)
314 
315  IF (.NOT. compare_potential_types(mp2_env%potential_parameter, mp2_env%ri_metric)) THEN
316  mp2_force(ikind)%mp2_non_sep(:, :) = mp2_force(ikind)%mp2_non_sep(:, :) + factor_2c*force_2c_ri(ikind)%forces
317  END IF
318  END DO
319 
320  CALL mp2_eri_deallocate_forces(force_2c)
321  IF (.NOT. compare_potential_types(mp2_env%potential_parameter, mp2_env%ri_metric)) THEN
322  CALL mp2_eri_deallocate_forces(force_2c_ri)
323  END IF
324  CALL mp2_eri_deallocate_forces(force_3c_aux)
325  CALL mp2_eri_deallocate_forces(force_3c_orb_mu)
326  CALL mp2_eri_deallocate_forces(force_3c_orb_nu)
327  CALL dbcsr_deallocate_matrix_set(matrix_p_munu_local)
328  CALL dbcsr_deallocate_matrix_set(mat_munu_local)
329 
330  ELSEIF (eri_method == do_eri_gpw) THEN
331  CALL get_qs_env(qs_env, ks_env=ks_env)
332 
333  CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of, atom_of_kind=atom_of_kind)
334 
335  ! Supporting stuff for GPW
336  CALL prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
337  auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_l, sab_orb_sub)
338 
339  ! in case virial is required we need auxiliary pw
340  ! for calculate the MP2-volume contribution to the virial
341  ! (hartree potential derivatives)
342  IF (use_virial) THEN
343  CALL auxbas_pw_pool%create_pw(rho_g_copy)
344  DO i = 1, 3
345  CALL auxbas_pw_pool%create_pw(dvg(i))
346  END DO
347  END IF
348 
349  ! start main loop over auxiliary basis functions
350  CALL timeset(routinen//"_loop", handle2)
351 
352  IF (use_virial) h_stress = 0.0_dp
353 
354  l_counter = 0
355  DO lll = my_group_l_start, my_group_l_end
356  l_counter = l_counter + 1
357 
358  CALL g_p_transform_mo_to_ao(matrix_p_munu%matrix, matrix_p_munu_nosym, mat_munu%matrix, &
359  g_p_ia(:, l_counter), matrix_p_inu, &
360  mo_coeff_v, mo_coeff_o, eps_filter)
361 
362  CALL integrate_potential_forces_2c(rho_r, lll, mo_coeff(1), rho_g, atomic_kind_set, &
363  qs_kind_set, particle_set, cell, pw_env_sub, poisson_env, &
364  pot_g, mp2_env%potential_parameter, use_virial, &
365  rho_g_copy, dvg, kind_of, atom_of_kind, g_pq_local(:, l_counter), &
366  mp2_force, h_stress, para_env_sub, dft_control, psi_l, factor_2c)
367 
368  IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) THEN
369 
370  CALL integrate_potential_forces_2c(rho_r, lll, mo_coeff(1), rho_g, atomic_kind_set, &
371  qs_kind_set, particle_set, cell, pw_env_sub, poisson_env, &
372  pot_g, mp2_env%ri_metric, use_virial, &
373  rho_g_copy, dvg, kind_of, atom_of_kind, g_pq_local_2(:, l_counter), &
374  mp2_force, h_stress, para_env_sub, dft_control, psi_l, factor_2c)
375  END IF
376 
377  IF (use_virial) pv_virial = virial%pv_virial
378  CALL integrate_potential_forces_3c_1c(mat_munu, rho_r, matrix_p_munu, qs_env, pw_env_sub, &
379  task_list_sub)
380  IF (use_virial) THEN
381  h_stress = h_stress + (virial%pv_virial - pv_virial)
382  virial%pv_virial = pv_virial
383  END IF
384 
385  ! The matrices of G_P_ia are deallocated here
386  CALL update_lagrangian(mat_munu%matrix, matrix_p_inu, lag_mu_i_1, &
387  g_p_ia(:, l_counter), mo_coeff_o, lag_nu_a_2, &
388  eps_filter)
389 
390  CALL integrate_potential_forces_3c_2c(matrix_p_munu, rho_r, rho_g, task_list_sub, pw_env_sub, &
391  mp2_env%ri_metric, &
392  ks_env, poisson_env, pot_g, use_virial, rho_g_copy, dvg, &
393  h_stress, para_env_sub, kind_of, atom_of_kind, &
394  qs_kind_set, particle_set, cell, lll, mp2_force, dft_control)
395  END DO
396 
397  CALL timestop(handle2)
398 
399  DEALLOCATE (kind_of)
400  DEALLOCATE (atom_of_kind)
401 
402  IF (use_virial) THEN
403  CALL auxbas_pw_pool%give_back_pw(rho_g_copy)
404  DO i = 1, 3
405  CALL auxbas_pw_pool%give_back_pw(dvg(i))
406  END DO
407  END IF
408 
409  CALL cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, &
410  task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
411 
412  CALL dbcsr_release(matrix_p_munu%matrix)
413  DEALLOCATE (matrix_p_munu%matrix)
414 
415  END IF
416 
417  IF (use_virial) mp2_env%ri_grad%mp2_virial = h_stress
418 
419  DEALLOCATE (g_pq_local)
420  IF (.NOT. compare_potential_types(mp2_env%ri_metric, mp2_env%potential_parameter)) DEALLOCATE (g_pq_local_2)
421 
422  CALL dbcsr_release(matrix_p_munu_nosym)
423 
424  DO ispin = 1, nspins
425  CALL dbcsr_release(matrix_p_inu(ispin))
426  END DO
427  DEALLOCATE (matrix_p_inu, g_p_ia)
428 
429  ! move the forces in the correct place
430  IF (eri_method .EQ. do_eri_gpw) THEN
431  DO ikind = 1, SIZE(mp2_force)
432  mp2_force(ikind)%mp2_non_sep(:, :) = force(ikind)%rho_elec(:, :)
433  force(ikind)%rho_elec(:, :) = 0.0_dp
434  END DO
435  END IF
436 
437  ! Now we move from the local matrices to the global ones
438  ! defined over all MPI tasks
439  ! Start with moving from the DBCSR to FM for the lagrangians
440 
441  ALLOCATE (l1_mu_i(nspins), l2_nu_a(nspins))
442  DO ispin = 1, nspins
443  ! Now we move from the local matrices to the global ones
444  ! defined over all MPI tasks
445  ! Start with moving from the DBCSR to FM for the lagrangians
446  NULLIFY (fm_struct_tmp)
447  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env_sub, context=blacs_env_sub, &
448  nrow_global=dimen, ncol_global=homo(ispin))
449  CALL cp_fm_create(l1_mu_i(ispin), fm_struct_tmp, name="Lag_mu_i")
450  CALL cp_fm_struct_release(fm_struct_tmp)
451  CALL cp_fm_set_all(l1_mu_i(ispin), 0.0_dp)
452  CALL copy_dbcsr_to_fm(matrix=lag_mu_i_1(ispin), fm=l1_mu_i(ispin))
453 
454  ! release Lag_mu_i_1
455  CALL dbcsr_release(lag_mu_i_1(ispin))
456 
457  NULLIFY (fm_struct_tmp)
458  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env_sub, context=blacs_env_sub, &
459  nrow_global=dimen, ncol_global=virtual(ispin))
460  CALL cp_fm_create(l2_nu_a(ispin), fm_struct_tmp, name="Lag_nu_a")
461  CALL cp_fm_struct_release(fm_struct_tmp)
462  CALL cp_fm_set_all(l2_nu_a(ispin), 0.0_dp)
463  CALL copy_dbcsr_to_fm(matrix=lag_nu_a_2(ispin), fm=l2_nu_a(ispin))
464 
465  ! release Lag_nu_a_2
466  CALL dbcsr_release(lag_nu_a_2(ispin))
467  END DO
468  DEALLOCATE (lag_mu_i_1, lag_nu_a_2)
469 
470  ! Set the factor to multiply P_ij (depends on the open or closed shell)
471  factor = 1.0_dp
472  IF (alpha_beta) factor = 0.50_dp
473 
474  DO ispin = 1, nspins
475  CALL create_w_p(qs_env, mp2_env, mo_coeff(ispin), homo(ispin), virtual(ispin), dimen, para_env, &
476  para_env_sub, eigenval(:, ispin), l1_mu_i(ispin), l2_nu_a(ispin), &
477  factor, ispin)
478  END DO
479  DEALLOCATE (l1_mu_i, l2_nu_a)
480 
481  CALL timestop(handle)
482 
483  END SUBROUTINE calc_ri_mp2_nonsep
484 
485 ! **************************************************************************************************
486 !> \brief Transforms G_P_ia to G_P_munu
487 !> \param G_P_munu The container for G_P_munu, will be allocated and created if not allocated on entry
488 !> \param G_P_munu_nosym ...
489 !> \param mat_munu ...
490 !> \param G_P_ia ...
491 !> \param G_P_inu ...
492 !> \param mo_coeff_v ...
493 !> \param mo_coeff_o ...
494 !> \param eps_filter ...
495 ! **************************************************************************************************
496  SUBROUTINE g_p_transform_mo_to_ao(G_P_munu, G_P_munu_nosym, mat_munu, G_P_ia, G_P_inu, &
497  mo_coeff_v, mo_coeff_o, eps_filter)
498  TYPE(dbcsr_type), POINTER :: g_p_munu
499  TYPE(dbcsr_type), INTENT(INOUT) :: g_p_munu_nosym, mat_munu
500  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: g_p_ia
501  TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: g_p_inu
502  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mo_coeff_v, mo_coeff_o
503  REAL(kind=dp), INTENT(IN) :: eps_filter
504 
505  CHARACTER(LEN=*), PARAMETER :: routinen = 'G_P_transform_MO_to_AO'
506 
507  INTEGER :: handle
508 
509  IF (.NOT. ASSOCIATED(g_p_munu)) THEN
510  ALLOCATE (g_p_munu)
511  CALL dbcsr_create(g_p_munu, template=mat_munu, &
512  matrix_type=dbcsr_type_symmetric)
513  END IF
514 
515  CALL g_p_transform_alpha_beta(g_p_ia, g_p_inu, g_p_munu_nosym, mo_coeff_v, mo_coeff_o, eps_filter)
516 
517  ! symmetrize
518  CALL timeset(routinen//"_symmetrize", handle)
519  CALL dbcsr_set(g_p_munu, 0.0_dp)
520  CALL dbcsr_transposed(g_p_munu, g_p_munu_nosym)
521  CALL dbcsr_add(g_p_munu, g_p_munu_nosym, &
522  alpha_scalar=2.0_dp, beta_scalar=2.0_dp)
523  ! this is a trick to avoid that integrate_v_rspace starts to cry
524  CALL dbcsr_copy_into_existing(mat_munu, g_p_munu)
525  CALL dbcsr_copy(g_p_munu, mat_munu)
526 
527  CALL timestop(handle)
528 
529  END SUBROUTINE g_p_transform_mo_to_ao
530 
531 ! **************************************************************************************************
532 !> \brief ...
533 !> \param G_P_ia ...
534 !> \param G_P_inu ...
535 !> \param G_P_munu ...
536 !> \param mo_coeff_v ...
537 !> \param mo_coeff_o ...
538 !> \param eps_filter ...
539 ! **************************************************************************************************
540  SUBROUTINE g_p_transform_alpha_beta(G_P_ia, G_P_inu, G_P_munu, mo_coeff_v, mo_coeff_o, eps_filter)
541  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: g_p_ia
542  TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: g_p_inu
543  TYPE(dbcsr_type), INTENT(INOUT) :: g_p_munu
544  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mo_coeff_v, mo_coeff_o
545  REAL(kind=dp), INTENT(IN) :: eps_filter
546 
547  CHARACTER(LEN=*), PARAMETER :: routinen = 'G_P_transform_alpha_beta'
548 
549  INTEGER :: handle, ispin
550  REAL(kind=dp) :: factor
551 
552  CALL timeset(routinen, handle)
553 
554  factor = 1.0_dp/real(SIZE(g_p_ia), dp)
555 
556  CALL dbcsr_set(g_p_munu, 0.0_dp)
557 
558  DO ispin = 1, SIZE(g_p_ia)
559  ! first back-transformation a->nu
560  CALL dbcsr_multiply("N", "T", 1.0_dp, mo_coeff_v(ispin)%matrix, g_p_ia(ispin)%matrix, &
561  0.0_dp, g_p_inu(ispin), filter_eps=eps_filter)
562 
563  ! second back-transformation i->mu
564  CALL dbcsr_multiply("N", "T", factor, g_p_inu(ispin), mo_coeff_o(ispin)%matrix, &
565  1.0_dp, g_p_munu, filter_eps=eps_filter)
566  END DO
567 
568  CALL timestop(handle)
569 
570  END SUBROUTINE g_p_transform_alpha_beta
571 
572 ! **************************************************************************************************
573 !> \brief ...
574 !> \param mat_munu ...
575 !> \param matrix_P_inu ...
576 !> \param Lag_mu_i_1 ...
577 !> \param G_P_ia ...
578 !> \param mo_coeff_o ...
579 !> \param Lag_nu_a_2 ...
580 !> \param eps_filter ...
581 ! **************************************************************************************************
582  SUBROUTINE update_lagrangian(mat_munu, matrix_P_inu, Lag_mu_i_1, &
583  G_P_ia, mo_coeff_o, Lag_nu_a_2, &
584  eps_filter)
585  TYPE(dbcsr_type), INTENT(IN) :: mat_munu
586  TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: matrix_p_inu, lag_mu_i_1
587  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT) :: g_p_ia
588  TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN) :: mo_coeff_o
589  TYPE(dbcsr_type), DIMENSION(:), INTENT(INOUT) :: lag_nu_a_2
590  REAL(kind=dp), INTENT(IN) :: eps_filter
591 
592  CHARACTER(LEN=*), PARAMETER :: routinen = 'update_lagrangian'
593 
594  INTEGER :: handle, ispin
595 
596  ! update lagrangian
597  CALL timeset(routinen, handle)
598 
599  DO ispin = 1, SIZE(g_p_ia)
600  ! first contract mat_munu with the half back transformed Gamma_i_nu
601  ! in order to update Lag_mu_i_1
602  CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu, matrix_p_inu(ispin), &
603  1.0_dp, lag_mu_i_1(ispin), filter_eps=eps_filter)
604 
605  ! transform first index of mat_munu and store the result into matrix_P_inu
606  CALL dbcsr_set(matrix_p_inu(ispin), 0.0_dp)
607  CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu, mo_coeff_o(ispin)%matrix, &
608  0.0_dp, matrix_p_inu(ispin), filter_eps=eps_filter)
609 
610  ! contract the transformend matrix_P_inu with the untransformend Gamma_i_a
611  ! in order to update Lag_nu_a_2
612  CALL dbcsr_multiply("N", "N", -1.0_dp, matrix_p_inu(ispin), g_p_ia(ispin)%matrix, &
613  1.0_dp, lag_nu_a_2(ispin), filter_eps=eps_filter)
614 
615  ! release the actual gamma_P_ia
616  CALL dbcsr_release(g_p_ia(ispin)%matrix)
617  DEALLOCATE (g_p_ia(ispin)%matrix)
618  END DO
619 
620  CALL timestop(handle)
621 
622  END SUBROUTINE update_lagrangian
623 
624 ! **************************************************************************************************
625 !> \brief ...
626 !> \param qs_env ...
627 !> \param mp2_env ...
628 !> \param mo_coeff ...
629 !> \param homo ...
630 !> \param virtual ...
631 !> \param dimen ...
632 !> \param para_env ...
633 !> \param para_env_sub ...
634 !> \param Eigenval ...
635 !> \param L1_mu_i ...
636 !> \param L2_nu_a ...
637 !> \param factor ...
638 !> \param kspin ...
639 ! **************************************************************************************************
640  SUBROUTINE create_w_p(qs_env, mp2_env, mo_coeff, homo, virtual, dimen, para_env, para_env_sub, &
641  Eigenval, L1_mu_i, L2_nu_a, factor, kspin)
642  TYPE(qs_environment_type), POINTER :: qs_env
643  TYPE(mp2_type) :: mp2_env
644  TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
645  INTEGER, INTENT(IN) :: homo, virtual, dimen
646  TYPE(mp_para_env_type), POINTER :: para_env, para_env_sub
647  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: eigenval
648  TYPE(cp_fm_type), INTENT(INOUT) :: l1_mu_i, l2_nu_a
649  REAL(kind=dp), INTENT(IN) :: factor
650  INTEGER, INTENT(IN) :: kspin
651 
652  CHARACTER(LEN=*), PARAMETER :: routinen = 'create_W_P'
653 
654  INTEGER :: color_exchange, dummy_proc, handle, handle2, handle3, i_global, i_local, iib, &
655  iii, iproc, itmp(2), j_global, j_local, jjb, max_col_size, max_row_size, &
656  my_b_virtual_end, my_b_virtual_start, mypcol, myprow, ncol_block, ncol_block_1i, &
657  ncol_block_2a, ncol_local, ncol_local_1i, ncol_local_2a, npcol, npcol_1i, npcol_2a, &
658  nprow, nprow_1i, nprow_2a, nrow_block, nrow_block_1i, nrow_block_2a, nrow_local, &
659  nrow_local_1i, nrow_local_2a, number_of_rec, number_of_send, proc_receive, &
660  proc_receive_static, proc_send, proc_send_ex, proc_send_static, proc_send_sub, &
661  proc_shift, rec_col_size
662  INTEGER :: rec_counter, rec_row_size, send_col_size, send_counter, send_pcol, send_prow, &
663  send_row_size, size_rec_buffer, size_send_buffer
664  INTEGER, ALLOCATABLE, DIMENSION(:) :: iii_vet, map_rec_size, map_send_size, &
665  pos_info, pos_info_ex, proc_2_send_pos
666  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: grid_2_mepos, mepos_2_grid, my_col_indeces_info_1i, &
667  my_col_indeces_info_2a, my_row_indeces_info_1i, my_row_indeces_info_2a, sizes, sizes_1i, &
668  sizes_2a
669  INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: col_indeces_info_1i, &
670  col_indeces_info_2a, &
671  row_indeces_info_1i, &
672  row_indeces_info_2a
673  INTEGER, DIMENSION(:), POINTER :: col_indices, col_indices_1i, &
674  col_indices_2a, row_indices, &
675  row_indices_1i, row_indices_2a
676  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: ab_rec, ab_send, mat_rec, mat_send
677  TYPE(cp_blacs_env_type), POINTER :: blacs_env
678  TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
679  TYPE(cp_fm_type) :: fm_p_ij, l_mu_q
680  TYPE(integ_mat_buffer_type), ALLOCATABLE, &
681  DIMENSION(:) :: buffer_rec, buffer_send
682  TYPE(integ_mat_buffer_type_2d), ALLOCATABLE, &
683  DIMENSION(:) :: buffer_cyclic
684  TYPE(mp_para_env_type), POINTER :: para_env_exchange
685  TYPE(mp_request_type), ALLOCATABLE, DIMENSION(:) :: req_send
686 
687  CALL timeset(routinen, handle)
688 
689  ! create the globally distributed mixed lagrangian
690  NULLIFY (blacs_env)
691  CALL get_qs_env(qs_env, blacs_env=blacs_env)
692 
693  NULLIFY (fm_struct_tmp)
694  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
695  nrow_global=dimen, ncol_global=dimen)
696  CALL cp_fm_create(l_mu_q, fm_struct_tmp, name="Lag_mu_q")
697  CALL cp_fm_struct_release(fm_struct_tmp)
698  CALL cp_fm_set_all(l_mu_q, 0.0_dp)
699 
700  ! create all information array
701  ALLOCATE (pos_info(0:para_env%num_pe - 1))
702  CALL para_env%allgather(para_env_sub%mepos, pos_info)
703 
704  ! get matrix information for the global
705  CALL cp_fm_get_info(matrix=l_mu_q, &
706  nrow_local=nrow_local, &
707  ncol_local=ncol_local, &
708  row_indices=row_indices, &
709  col_indices=col_indices, &
710  nrow_block=nrow_block, &
711  ncol_block=ncol_block)
712  myprow = l_mu_q%matrix_struct%context%mepos(1)
713  mypcol = l_mu_q%matrix_struct%context%mepos(2)
714  nprow = l_mu_q%matrix_struct%context%num_pe(1)
715  npcol = l_mu_q%matrix_struct%context%num_pe(2)
716 
717  ALLOCATE (grid_2_mepos(0:nprow - 1, 0:npcol - 1))
718  grid_2_mepos = 0
719  grid_2_mepos(myprow, mypcol) = para_env%mepos
720  CALL para_env%sum(grid_2_mepos)
721 
722  ! get matrix information for L1_mu_i
723  CALL cp_fm_get_info(matrix=l1_mu_i, &
724  nrow_local=nrow_local_1i, &
725  ncol_local=ncol_local_1i, &
726  row_indices=row_indices_1i, &
727  col_indices=col_indices_1i, &
728  nrow_block=nrow_block_1i, &
729  ncol_block=ncol_block_1i)
730  nprow_1i = l1_mu_i%matrix_struct%context%num_pe(1)
731  npcol_1i = l1_mu_i%matrix_struct%context%num_pe(2)
732 
733  ALLOCATE (sizes_1i(2, 0:para_env_sub%num_pe - 1))
734  CALL para_env_sub%allgather([nrow_local_1i, ncol_local_1i], sizes_1i)
735 
736  ! get matrix information for L2_nu_a
737  CALL cp_fm_get_info(matrix=l2_nu_a, &
738  nrow_local=nrow_local_2a, &
739  ncol_local=ncol_local_2a, &
740  row_indices=row_indices_2a, &
741  col_indices=col_indices_2a, &
742  nrow_block=nrow_block_2a, &
743  ncol_block=ncol_block_2a)
744  nprow_2a = l2_nu_a%matrix_struct%context%num_pe(1)
745  npcol_2a = l2_nu_a%matrix_struct%context%num_pe(2)
746 
747  ALLOCATE (sizes_2a(2, 0:para_env_sub%num_pe - 1))
748  CALL para_env_sub%allgather([nrow_local_2a, ncol_local_2a], sizes_2a)
749 
750  ! Here we perform a ring communication scheme taking into account
751  ! for the sub-group distribution of the source matrices.
752  ! as a first step we need to redistribute the data within
753  ! the subgroup.
754  ! In order to do so we have to allocate the structure
755  ! that will hold the local data involved in communication, this
756  ! structure will be the same for processes in different subgroups
757  ! sharing the same position in the subgroup.
758  ! -1) create the exchange para_env
759  color_exchange = para_env_sub%mepos
760  ALLOCATE (para_env_exchange)
761  CALL para_env_exchange%from_split(para_env, color_exchange)
762  ALLOCATE (pos_info_ex(0:para_env%num_pe - 1))
763  CALL para_env%allgather(para_env_exchange%mepos, pos_info_ex)
764  ALLOCATE (sizes(2, 0:para_env_exchange%num_pe - 1))
765  CALL para_env_exchange%allgather([nrow_local, ncol_local], sizes)
766 
767  ! 0) store some info about indeces of the fm matrices (subgroup)
768  CALL timeset(routinen//"_inx", handle2)
769  ! matrix L1_mu_i
770  max_row_size = maxval(sizes_1i(1, :))
771  max_col_size = maxval(sizes_1i(2, :))
772  ALLOCATE (row_indeces_info_1i(2, max_row_size, 0:para_env_sub%num_pe - 1))
773  ALLOCATE (col_indeces_info_1i(2, max_col_size, 0:para_env_sub%num_pe - 1))
774  ALLOCATE (my_row_indeces_info_1i(2, max_row_size))
775  ALLOCATE (my_col_indeces_info_1i(2, max_col_size))
776  row_indeces_info_1i = 0
777  col_indeces_info_1i = 0
778  dummy_proc = 0
779  ! row
780  DO iib = 1, nrow_local_1i
781  i_global = row_indices_1i(iib)
782  send_prow = cp_fm_indxg2p(i_global, nrow_block, dummy_proc, &
783  l_mu_q%matrix_struct%first_p_pos(1), nprow)
784  i_local = cp_fm_indxg2l(i_global, nrow_block, dummy_proc, &
785  l_mu_q%matrix_struct%first_p_pos(1), nprow)
786  my_row_indeces_info_1i(1, iib) = send_prow
787  my_row_indeces_info_1i(2, iib) = i_local
788  END DO
789  ! col
790  DO jjb = 1, ncol_local_1i
791  j_global = col_indices_1i(jjb)
792  send_pcol = cp_fm_indxg2p(j_global, ncol_block, dummy_proc, &
793  l_mu_q%matrix_struct%first_p_pos(2), npcol)
794  j_local = cp_fm_indxg2l(j_global, ncol_block, dummy_proc, &
795  l_mu_q%matrix_struct%first_p_pos(2), npcol)
796  my_col_indeces_info_1i(1, jjb) = send_pcol
797  my_col_indeces_info_1i(2, jjb) = j_local
798  END DO
799  CALL para_env_sub%allgather(my_row_indeces_info_1i, row_indeces_info_1i)
800  CALL para_env_sub%allgather(my_col_indeces_info_1i, col_indeces_info_1i)
801  DEALLOCATE (my_row_indeces_info_1i, my_col_indeces_info_1i)
802 
803  ! matrix L2_nu_a
804  max_row_size = maxval(sizes_2a(1, :))
805  max_col_size = maxval(sizes_2a(2, :))
806  ALLOCATE (row_indeces_info_2a(2, max_row_size, 0:para_env_sub%num_pe - 1))
807  ALLOCATE (col_indeces_info_2a(2, max_col_size, 0:para_env_sub%num_pe - 1))
808  ALLOCATE (my_row_indeces_info_2a(2, max_row_size))
809  ALLOCATE (my_col_indeces_info_2a(2, max_col_size))
810  row_indeces_info_2a = 0
811  col_indeces_info_2a = 0
812  ! row
813  DO iib = 1, nrow_local_2a
814  i_global = row_indices_2a(iib)
815  send_prow = cp_fm_indxg2p(i_global, nrow_block, dummy_proc, &
816  l_mu_q%matrix_struct%first_p_pos(1), nprow)
817  i_local = cp_fm_indxg2l(i_global, nrow_block, dummy_proc, &
818  l_mu_q%matrix_struct%first_p_pos(1), nprow)
819  my_row_indeces_info_2a(1, iib) = send_prow
820  my_row_indeces_info_2a(2, iib) = i_local
821  END DO
822  ! col
823  DO jjb = 1, ncol_local_2a
824  j_global = col_indices_2a(jjb) + homo
825  send_pcol = cp_fm_indxg2p(j_global, ncol_block, dummy_proc, &
826  l_mu_q%matrix_struct%first_p_pos(2), npcol)
827  j_local = cp_fm_indxg2l(j_global, ncol_block, dummy_proc, &
828  l_mu_q%matrix_struct%first_p_pos(2), npcol)
829  my_col_indeces_info_2a(1, jjb) = send_pcol
830  my_col_indeces_info_2a(2, jjb) = j_local
831  END DO
832  CALL para_env_sub%allgather(my_row_indeces_info_2a, row_indeces_info_2a)
833  CALL para_env_sub%allgather(my_col_indeces_info_2a, col_indeces_info_2a)
834  DEALLOCATE (my_row_indeces_info_2a, my_col_indeces_info_2a)
835  CALL timestop(handle2)
836 
837  ! 1) define the map for sending data in the subgroup starting with L1_mu_i
838  CALL timeset(routinen//"_subinfo", handle2)
839  ALLOCATE (map_send_size(0:para_env_sub%num_pe - 1))
840  map_send_size = 0
841  DO jjb = 1, ncol_local_1i
842  send_pcol = col_indeces_info_1i(1, jjb, para_env_sub%mepos)
843  DO iib = 1, nrow_local_1i
844  send_prow = row_indeces_info_1i(1, iib, para_env_sub%mepos)
845  proc_send = grid_2_mepos(send_prow, send_pcol)
846  proc_send_sub = pos_info(proc_send)
847  map_send_size(proc_send_sub) = map_send_size(proc_send_sub) + 1
848  END DO
849  END DO
850  ! and the same for L2_nu_a
851  DO jjb = 1, ncol_local_2a
852  send_pcol = col_indeces_info_2a(1, jjb, para_env_sub%mepos)
853  DO iib = 1, nrow_local_2a
854  send_prow = row_indeces_info_2a(1, iib, para_env_sub%mepos)
855  proc_send = grid_2_mepos(send_prow, send_pcol)
856  proc_send_sub = pos_info(proc_send)
857  map_send_size(proc_send_sub) = map_send_size(proc_send_sub) + 1
858  END DO
859  END DO
860  ! and exchange data in order to create map_rec_size
861  ALLOCATE (map_rec_size(0:para_env_sub%num_pe - 1))
862  map_rec_size = 0
863  CALL para_env_sub%alltoall(map_send_size, map_rec_size, 1)
864  CALL timestop(handle2)
865 
866  ! 2) reorder data in sending buffer
867  CALL timeset(routinen//"_sub_Bsend", handle2)
868  ! count the number of messages (include myself)
869  number_of_send = 0
870  DO proc_shift = 0, para_env_sub%num_pe - 1
871  proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
872  IF (map_send_size(proc_send) > 0) THEN
873  number_of_send = number_of_send + 1
874  END IF
875  END DO
876  ! allocate the structure that will hold the messages to be sent
877  ALLOCATE (buffer_send(number_of_send))
878  send_counter = 0
879  ALLOCATE (proc_2_send_pos(0:para_env_sub%num_pe - 1))
880  proc_2_send_pos = 0
881  DO proc_shift = 0, para_env_sub%num_pe - 1
882  proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
883  size_send_buffer = map_send_size(proc_send)
884  IF (map_send_size(proc_send) > 0) THEN
885  send_counter = send_counter + 1
886  ! allocate the sending buffer (msg)
887  ALLOCATE (buffer_send(send_counter)%msg(size_send_buffer))
888  buffer_send(send_counter)%msg = 0.0_dp
889  buffer_send(send_counter)%proc = proc_send
890  proc_2_send_pos(proc_send) = send_counter
891  END IF
892  END DO
893  ! loop over the locally held data and fill the buffer_send
894  ! for doing that we need an array that keep track if the
895  ! sequential increase of the index for each message
896  ALLOCATE (iii_vet(number_of_send))
897  iii_vet = 0
898  DO jjb = 1, ncol_local_1i
899  send_pcol = col_indeces_info_1i(1, jjb, para_env_sub%mepos)
900  DO iib = 1, nrow_local_1i
901  send_prow = row_indeces_info_1i(1, iib, para_env_sub%mepos)
902  proc_send = grid_2_mepos(send_prow, send_pcol)
903  proc_send_sub = pos_info(proc_send)
904  send_counter = proc_2_send_pos(proc_send_sub)
905  iii_vet(send_counter) = iii_vet(send_counter) + 1
906  iii = iii_vet(send_counter)
907  buffer_send(send_counter)%msg(iii) = l1_mu_i%local_data(iib, jjb)
908  END DO
909  END DO
910  ! release the local data of L1_mu_i
911  DEALLOCATE (l1_mu_i%local_data)
912  ! and the same for L2_nu_a
913  DO jjb = 1, ncol_local_2a
914  send_pcol = col_indeces_info_2a(1, jjb, para_env_sub%mepos)
915  DO iib = 1, nrow_local_2a
916  send_prow = row_indeces_info_2a(1, iib, para_env_sub%mepos)
917  proc_send = grid_2_mepos(send_prow, send_pcol)
918  proc_send_sub = pos_info(proc_send)
919  send_counter = proc_2_send_pos(proc_send_sub)
920  iii_vet(send_counter) = iii_vet(send_counter) + 1
921  iii = iii_vet(send_counter)
922  buffer_send(send_counter)%msg(iii) = l2_nu_a%local_data(iib, jjb)
923  END DO
924  END DO
925  DEALLOCATE (l2_nu_a%local_data)
926  DEALLOCATE (proc_2_send_pos)
927  DEALLOCATE (iii_vet)
928  CALL timestop(handle2)
929 
930  ! 3) create the buffer for receive, post the message with irecv
931  ! and send the messages non-blocking
932  CALL timeset(routinen//"_sub_isendrecv", handle2)
933  ! count the number of messages to be received
934  number_of_rec = 0
935  DO proc_shift = 0, para_env_sub%num_pe - 1
936  proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
937  IF (map_rec_size(proc_receive) > 0) THEN
938  number_of_rec = number_of_rec + 1
939  END IF
940  END DO
941  ALLOCATE (buffer_rec(number_of_rec))
942  rec_counter = 0
943  DO proc_shift = 0, para_env_sub%num_pe - 1
944  proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
945  size_rec_buffer = map_rec_size(proc_receive)
946  IF (map_rec_size(proc_receive) > 0) THEN
947  rec_counter = rec_counter + 1
948  ! prepare the buffer for receive
949  ALLOCATE (buffer_rec(rec_counter)%msg(size_rec_buffer))
950  buffer_rec(rec_counter)%msg = 0.0_dp
951  buffer_rec(rec_counter)%proc = proc_receive
952  ! post the message to be received (not need to send to myself)
953  IF (proc_receive /= para_env_sub%mepos) THEN
954  CALL para_env_sub%irecv(buffer_rec(rec_counter)%msg, proc_receive, &
955  buffer_rec(rec_counter)%msg_req)
956  END IF
957  END IF
958  END DO
959  ! send messages
960  ALLOCATE (req_send(number_of_send))
961  req_send = mp_request_null
962  send_counter = 0
963  DO proc_shift = 0, para_env_sub%num_pe - 1
964  proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
965  IF (map_send_size(proc_send) > 0) THEN
966  send_counter = send_counter + 1
967  IF (proc_send == para_env_sub%mepos) THEN
968  buffer_rec(send_counter)%msg(:) = buffer_send(send_counter)%msg
969  ELSE
970  CALL para_env_sub%isend(buffer_send(send_counter)%msg, proc_send, &
971  buffer_send(send_counter)%msg_req)
972  req_send(send_counter) = buffer_send(send_counter)%msg_req
973  END IF
974  END IF
975  END DO
976  DEALLOCATE (map_send_size)
977  CALL timestop(handle2)
978 
979  ! 4) (if memory is a problem we should move this part after point 5)
980  ! Here we create the new buffer for cyclic(ring) communication and
981  ! we fill it with the data received from the other member of the
982  ! subgroup
983  CALL timeset(routinen//"_Bcyclic", handle2)
984  ! first allocata new structure
985  ALLOCATE (buffer_cyclic(0:para_env_exchange%num_pe - 1))
986  DO iproc = 0, para_env_exchange%num_pe - 1
987  rec_row_size = sizes(1, iproc)
988  rec_col_size = sizes(2, iproc)
989  ALLOCATE (buffer_cyclic(iproc)%msg(rec_row_size, rec_col_size))
990  buffer_cyclic(iproc)%msg = 0.0_dp
991  END DO
992  ! now collect data from other member of the subgroup and fill
993  ! buffer_cyclic
994  rec_counter = 0
995  DO proc_shift = 0, para_env_sub%num_pe - 1
996  proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
997  size_rec_buffer = map_rec_size(proc_receive)
998  IF (map_rec_size(proc_receive) > 0) THEN
999  rec_counter = rec_counter + 1
1000 
1001  ! wait for the message
1002  IF (proc_receive /= para_env_sub%mepos) CALL buffer_rec(rec_counter)%msg_req%wait()
1003 
1004  CALL timeset(routinen//"_fill", handle3)
1005  iii = 0
1006  DO jjb = 1, sizes_1i(2, proc_receive)
1007  send_pcol = col_indeces_info_1i(1, jjb, proc_receive)
1008  j_local = col_indeces_info_1i(2, jjb, proc_receive)
1009  DO iib = 1, sizes_1i(1, proc_receive)
1010  send_prow = row_indeces_info_1i(1, iib, proc_receive)
1011  proc_send = grid_2_mepos(send_prow, send_pcol)
1012  proc_send_sub = pos_info(proc_send)
1013  IF (proc_send_sub /= para_env_sub%mepos) cycle
1014  iii = iii + 1
1015  i_local = row_indeces_info_1i(2, iib, proc_receive)
1016  proc_send_ex = pos_info_ex(proc_send)
1017  buffer_cyclic(proc_send_ex)%msg(i_local, j_local) = buffer_rec(rec_counter)%msg(iii)
1018  END DO
1019  END DO
1020  ! and the same for L2_nu_a
1021  DO jjb = 1, sizes_2a(2, proc_receive)
1022  send_pcol = col_indeces_info_2a(1, jjb, proc_receive)
1023  j_local = col_indeces_info_2a(2, jjb, proc_receive)
1024  DO iib = 1, sizes_2a(1, proc_receive)
1025  send_prow = row_indeces_info_2a(1, iib, proc_receive)
1026  proc_send = grid_2_mepos(send_prow, send_pcol)
1027  proc_send_sub = pos_info(proc_send)
1028  IF (proc_send_sub /= para_env_sub%mepos) cycle
1029  iii = iii + 1
1030  i_local = row_indeces_info_2a(2, iib, proc_receive)
1031  proc_send_ex = pos_info_ex(proc_send)
1032  buffer_cyclic(proc_send_ex)%msg(i_local, j_local) = buffer_rec(rec_counter)%msg(iii)
1033  END DO
1034  END DO
1035  CALL timestop(handle3)
1036 
1037  ! deallocate the received message
1038  DEALLOCATE (buffer_rec(rec_counter)%msg)
1039  END IF
1040  END DO
1041  DEALLOCATE (row_indeces_info_1i)
1042  DEALLOCATE (col_indeces_info_1i)
1043  DEALLOCATE (row_indeces_info_2a)
1044  DEALLOCATE (col_indeces_info_2a)
1045  DEALLOCATE (buffer_rec)
1046  DEALLOCATE (map_rec_size)
1047  CALL timestop(handle2)
1048 
1049  ! 5) Wait for all messeges to be sent in the subgroup
1050  CALL timeset(routinen//"_sub_waitall", handle2)
1051  CALL mp_waitall(req_send(:))
1052  DO send_counter = 1, number_of_send
1053  DEALLOCATE (buffer_send(send_counter)%msg)
1054  END DO
1055  DEALLOCATE (buffer_send)
1056  DEALLOCATE (req_send)
1057  CALL timestop(handle2)
1058 
1059  ! 6) Start with ring communication
1060  CALL timeset(routinen//"_ring", handle2)
1061  proc_send_static = modulo(para_env_exchange%mepos + 1, para_env_exchange%num_pe)
1062  proc_receive_static = modulo(para_env_exchange%mepos - 1, para_env_exchange%num_pe)
1063  max_row_size = maxval(sizes(1, :))
1064  max_col_size = maxval(sizes(2, :))
1065  ALLOCATE (mat_send(max_row_size, max_col_size))
1066  ALLOCATE (mat_rec(max_row_size, max_col_size))
1067  mat_send = 0.0_dp
1068  mat_send(1:nrow_local, 1:ncol_local) = buffer_cyclic(para_env_exchange%mepos)%msg(:, :)
1069  DEALLOCATE (buffer_cyclic(para_env_exchange%mepos)%msg)
1070  DO proc_shift = 1, para_env_exchange%num_pe - 1
1071  proc_receive = modulo(para_env_exchange%mepos - proc_shift, para_env_exchange%num_pe)
1072 
1073  rec_row_size = sizes(1, proc_receive)
1074  rec_col_size = sizes(2, proc_receive)
1075 
1076  mat_rec = 0.0_dp
1077  CALL para_env_exchange%sendrecv(mat_send, proc_send_static, &
1078  mat_rec, proc_receive_static)
1079 
1080  mat_send = 0.0_dp
1081  mat_send(1:rec_row_size, 1:rec_col_size) = mat_rec(1:rec_row_size, 1:rec_col_size) + &
1082  buffer_cyclic(proc_receive)%msg(:, :)
1083 
1084  DEALLOCATE (buffer_cyclic(proc_receive)%msg)
1085  END DO
1086  ! and finally
1087  CALL para_env_exchange%sendrecv(mat_send, proc_send_static, &
1088  mat_rec, proc_receive_static)
1089  l_mu_q%local_data(1:nrow_local, 1:ncol_local) = mat_rec(1:nrow_local, 1:ncol_local)
1090  DEALLOCATE (buffer_cyclic)
1091  DEALLOCATE (mat_send)
1092  DEALLOCATE (mat_rec)
1093  CALL timestop(handle2)
1094 
1095  ! release para_env_exchange
1096  CALL mp_para_env_release(para_env_exchange)
1097 
1098  CALL cp_fm_release(l1_mu_i)
1099  CALL cp_fm_release(l2_nu_a)
1100  DEALLOCATE (pos_info_ex)
1101  DEALLOCATE (grid_2_mepos)
1102  DEALLOCATE (sizes)
1103  DEALLOCATE (sizes_1i)
1104  DEALLOCATE (sizes_2a)
1105 
1106  ! update the P_ij block of P_MP2 with the
1107  ! non-singular ij pairs
1108  CALL timeset(routinen//"_Pij", handle2)
1109  NULLIFY (fm_struct_tmp)
1110  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
1111  nrow_global=homo, ncol_global=homo)
1112  CALL cp_fm_create(fm_p_ij, fm_struct_tmp, name="fm_P_ij")
1113  CALL cp_fm_struct_release(fm_struct_tmp)
1114  CALL cp_fm_set_all(fm_p_ij, 0.0_dp)
1115 
1116  ! we have it, update P_ij local
1117  CALL cp_fm_get_info(matrix=fm_p_ij, &
1118  nrow_local=nrow_local, &
1119  ncol_local=ncol_local, &
1120  row_indices=row_indices, &
1121  col_indices=col_indices)
1122 
1123  IF (.NOT. mp2_env%method == ri_rpa_method_gpw) THEN
1124  CALL parallel_gemm('T', 'N', homo, homo, dimen, 1.0_dp, &
1125  mo_coeff, l_mu_q, 0.0_dp, fm_p_ij, &
1126  a_first_col=1, &
1127  a_first_row=1, &
1128  b_first_col=1, &
1129  b_first_row=1, &
1130  c_first_col=1, &
1131  c_first_row=1)
1132  CALL parallel_gemm('T', 'N', homo, homo, dimen, -2.0_dp, &
1133  l_mu_q, mo_coeff, 2.0_dp, fm_p_ij, &
1134  a_first_col=1, &
1135  a_first_row=1, &
1136  b_first_col=1, &
1137  b_first_row=1, &
1138  c_first_col=1, &
1139  c_first_row=1)
1140 
1141  DO jjb = 1, ncol_local
1142  j_global = col_indices(jjb)
1143  DO iib = 1, nrow_local
1144  i_global = row_indices(iib)
1145  ! diagonal elements and nearly degenerate ij pairs already updated
1146  IF (abs(eigenval(j_global) - eigenval(i_global)) < mp2_env%ri_grad%eps_canonical) THEN
1147  fm_p_ij%local_data(iib, jjb) = mp2_env%ri_grad%P_ij(kspin)%array(i_global, j_global)
1148  ELSE
1149  fm_p_ij%local_data(iib, jjb) = &
1150  factor*fm_p_ij%local_data(iib, jjb)/(eigenval(j_global) - eigenval(i_global))
1151  END IF
1152  END DO
1153  END DO
1154  ELSE
1155  DO jjb = 1, ncol_local
1156  j_global = col_indices(jjb)
1157  DO iib = 1, nrow_local
1158  i_global = row_indices(iib)
1159  fm_p_ij%local_data(iib, jjb) = mp2_env%ri_grad%P_ij(kspin)%array(i_global, j_global)
1160  END DO
1161  END DO
1162  END IF
1163  ! deallocate the local P_ij
1164  DEALLOCATE (mp2_env%ri_grad%P_ij(kspin)%array)
1165  CALL timestop(handle2)
1166 
1167  ! Now create and fill the P matrix (MO)
1168  ! FOR NOW WE ASSUME P_ab AND P_ij ARE REPLICATED OVER EACH MPI RANK
1169  IF (.NOT. ALLOCATED(mp2_env%ri_grad%P_mo)) THEN
1170  ALLOCATE (mp2_env%ri_grad%P_mo(SIZE(mp2_env%ri_grad%mo_coeff_o)))
1171  END IF
1172 
1173  CALL timeset(routinen//"_PMO", handle2)
1174  NULLIFY (fm_struct_tmp)
1175  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
1176  nrow_global=dimen, ncol_global=dimen)
1177  CALL cp_fm_create(mp2_env%ri_grad%P_mo(kspin), fm_struct_tmp, name="P_MP2_MO")
1178  CALL cp_fm_set_all(mp2_env%ri_grad%P_mo(kspin), 0.0_dp)
1179 
1180  ! start with the (easy) occ-occ block and locally held P_ab elements
1181  itmp = get_limit(virtual, para_env_sub%num_pe, para_env_sub%mepos)
1182  my_b_virtual_start = itmp(1)
1183  my_b_virtual_end = itmp(2)
1184 
1185  ! Fill occ-occ block
1186  CALL cp_fm_to_fm_submat(fm_p_ij, mp2_env%ri_grad%P_mo(kspin), homo, homo, 1, 1, 1, 1)
1187  CALL cp_fm_release(fm_p_ij)
1188 
1189  CALL cp_fm_get_info(mp2_env%ri_grad%P_mo(kspin), &
1190  nrow_local=nrow_local, &
1191  ncol_local=ncol_local, &
1192  row_indices=row_indices, &
1193  col_indices=col_indices, &
1194  nrow_block=nrow_block, &
1195  ncol_block=ncol_block)
1196  myprow = mp2_env%ri_grad%P_mo(kspin)%matrix_struct%context%mepos(1)
1197  mypcol = mp2_env%ri_grad%P_mo(kspin)%matrix_struct%context%mepos(2)
1198  nprow = mp2_env%ri_grad%P_mo(kspin)%matrix_struct%context%num_pe(1)
1199  npcol = mp2_env%ri_grad%P_mo(kspin)%matrix_struct%context%num_pe(2)
1200 
1201  IF (mp2_env%method == ri_mp2_laplace) THEN
1202  CALL parallel_gemm('T', 'N', virtual, virtual, dimen, 1.0_dp, &
1203  mo_coeff, l_mu_q, 0.0_dp, mp2_env%ri_grad%P_mo(kspin), &
1204  a_first_col=homo + 1, &
1205  a_first_row=1, &
1206  b_first_col=homo + 1, &
1207  b_first_row=1, &
1208  c_first_col=homo + 1, &
1209  c_first_row=homo + 1)
1210  CALL parallel_gemm('T', 'N', virtual, virtual, dimen, -2.0_dp, &
1211  l_mu_q, mo_coeff, 2.0_dp, mp2_env%ri_grad%P_mo(kspin), &
1212  a_first_col=homo + 1, &
1213  a_first_row=1, &
1214  b_first_col=homo + 1, &
1215  b_first_row=1, &
1216  c_first_col=homo + 1, &
1217  c_first_row=homo + 1)
1218  END IF
1219 
1220  IF (mp2_env%method == ri_mp2_method_gpw .OR. mp2_env%method == ri_rpa_method_gpw) THEN
1221  ! With MP2 and RPA, we have already calculated the density matrix elements
1222  DO jjb = 1, ncol_local
1223  j_global = col_indices(jjb)
1224  IF (j_global <= homo) cycle
1225  DO iib = 1, nrow_local
1226  i_global = row_indices(iib)
1227  IF (my_b_virtual_start <= i_global - homo .AND. i_global - homo <= my_b_virtual_end) THEN
1228  mp2_env%ri_grad%P_mo(kspin)%local_data(iib, jjb) = &
1229  mp2_env%ri_grad%P_ab(kspin)%array(i_global - homo - my_b_virtual_start + 1, j_global - homo)
1230  END IF
1231  END DO
1232  END DO
1233  ELSE IF (mp2_env%method == ri_mp2_laplace) THEN
1234  ! With Laplace-SOS-MP2, we still have to calculate the matrix elements of the non-degenerate pairs
1235  DO jjb = 1, ncol_local
1236  j_global = col_indices(jjb)
1237  IF (j_global <= homo) cycle
1238  DO iib = 1, nrow_local
1239  i_global = row_indices(iib)
1240  IF (abs(eigenval(i_global) - eigenval(j_global)) < mp2_env%ri_grad%eps_canonical) THEN
1241  IF (my_b_virtual_start <= i_global - homo .AND. i_global - homo <= my_b_virtual_end) THEN
1242  mp2_env%ri_grad%P_mo(kspin)%local_data(iib, jjb) = &
1243  mp2_env%ri_grad%P_ab(kspin)%array(i_global - homo - my_b_virtual_start + 1, j_global - homo)
1244  ELSE
1245  mp2_env%ri_grad%P_mo(kspin)%local_data(iib, jjb) = 0.0_dp
1246  END IF
1247  ELSE
1248  mp2_env%ri_grad%P_mo(kspin)%local_data(iib, jjb) = &
1249  factor*mp2_env%ri_grad%P_mo(kspin)%local_data(iib, jjb)/ &
1250  (eigenval(i_global) - eigenval(j_global))
1251  END IF
1252  END DO
1253  END DO
1254  ELSE
1255  cpabort("Calculation of virt-virt block of density matrix is dealt with elsewhere!")
1256  END IF
1257 
1258  ! send around the sub_group the local data and check if we
1259  ! have to update our block with external elements
1260  ALLOCATE (mepos_2_grid(2, 0:para_env_sub%num_pe - 1))
1261  CALL para_env_sub%allgather([myprow, mypcol], mepos_2_grid)
1262 
1263  ALLOCATE (sizes(2, 0:para_env_sub%num_pe - 1))
1264  CALL para_env_sub%allgather([nrow_local, ncol_local], sizes)
1265 
1266  ALLOCATE (ab_rec(nrow_local, ncol_local))
1267  DO proc_shift = 1, para_env_sub%num_pe - 1
1268  proc_send = modulo(para_env_sub%mepos + proc_shift, para_env_sub%num_pe)
1269  proc_receive = modulo(para_env_sub%mepos - proc_shift, para_env_sub%num_pe)
1270 
1271  send_prow = mepos_2_grid(1, proc_send)
1272  send_pcol = mepos_2_grid(2, proc_send)
1273 
1274  send_row_size = sizes(1, proc_send)
1275  send_col_size = sizes(2, proc_send)
1276 
1277  ALLOCATE (ab_send(send_row_size, send_col_size))
1278  ab_send = 0.0_dp
1279 
1280  ! first loop over row since in this way we can cycle
1281  DO iib = 1, send_row_size
1282  i_global = cp_fm_indxl2g(iib, nrow_block, send_prow, &
1283  mp2_env%ri_grad%P_mo(kspin)%matrix_struct%first_p_pos(1), nprow)
1284  IF (i_global <= homo) cycle
1285  i_global = i_global - homo
1286  IF (.NOT. (my_b_virtual_start <= i_global .AND. i_global <= my_b_virtual_end)) cycle
1287  DO jjb = 1, send_col_size
1288  j_global = cp_fm_indxl2g(jjb, ncol_block, send_pcol, &
1289  mp2_env%ri_grad%P_mo(kspin)%matrix_struct%first_p_pos(2), npcol)
1290  IF (j_global <= homo) cycle
1291  j_global = j_global - homo
1292  ab_send(iib, jjb) = mp2_env%ri_grad%P_ab(kspin)%array(i_global - my_b_virtual_start + 1, j_global)
1293  END DO
1294  END DO
1295 
1296  ab_rec = 0.0_dp
1297  CALL para_env_sub%sendrecv(ab_send, proc_send, &
1298  ab_rec, proc_receive)
1299  mp2_env%ri_grad%P_mo(kspin)%local_data(1:nrow_local, 1:ncol_local) = &
1300  mp2_env%ri_grad%P_mo(kspin)%local_data(1:nrow_local, 1:ncol_local) + &
1301  ab_rec(1:nrow_local, 1:ncol_local)
1302 
1303  DEALLOCATE (ab_send)
1304  END DO
1305  DEALLOCATE (ab_rec)
1306  DEALLOCATE (mepos_2_grid)
1307  DEALLOCATE (sizes)
1308 
1309  ! deallocate the local P_ab
1310  DEALLOCATE (mp2_env%ri_grad%P_ab(kspin)%array)
1311  CALL timestop(handle2)
1312 
1313  ! create also W_MP2_MO
1314  CALL timeset(routinen//"_WMO", handle2)
1315  IF (.NOT. ALLOCATED(mp2_env%ri_grad%W_mo)) THEN
1316  ALLOCATE (mp2_env%ri_grad%W_mo(SIZE(mp2_env%ri_grad%mo_coeff_o)))
1317  END IF
1318 
1319  CALL cp_fm_create(mp2_env%ri_grad%W_mo(kspin), fm_struct_tmp, name="W_MP2_MO")
1320  CALL cp_fm_struct_release(fm_struct_tmp)
1321 
1322  ! all block
1323  CALL parallel_gemm('T', 'N', dimen, dimen, dimen, 2.0_dp*factor, &
1324  l_mu_q, mo_coeff, 0.0_dp, mp2_env%ri_grad%W_mo(kspin), &
1325  a_first_col=1, &
1326  a_first_row=1, &
1327  b_first_col=1, &
1328  b_first_row=1, &
1329  c_first_col=1, &
1330  c_first_row=1)
1331 
1332  ! occ-occ block
1333  CALL parallel_gemm('T', 'N', homo, homo, dimen, -2.0_dp*factor, &
1334  l_mu_q, mo_coeff, 0.0_dp, mp2_env%ri_grad%W_mo(kspin), &
1335  a_first_col=1, &
1336  a_first_row=1, &
1337  b_first_col=1, &
1338  b_first_row=1, &
1339  c_first_col=1, &
1340  c_first_row=1)
1341 
1342  ! occ-virt block
1343  CALL parallel_gemm('T', 'N', homo, virtual, dimen, 2.0_dp*factor, &
1344  mo_coeff, l_mu_q, 0.0_dp, mp2_env%ri_grad%W_mo(kspin), &
1345  a_first_col=1, &
1346  a_first_row=1, &
1347  b_first_col=homo + 1, &
1348  b_first_row=1, &
1349  c_first_col=homo + 1, &
1350  c_first_row=1)
1351  CALL timestop(handle2)
1352 
1353  ! Calculate occ-virt block of the lagrangian in MO
1354  CALL timeset(routinen//"_Ljb", handle2)
1355  IF (.NOT. ALLOCATED(mp2_env%ri_grad%L_jb)) THEN
1356  ALLOCATE (mp2_env%ri_grad%L_jb(SIZE(mp2_env%ri_grad%mo_coeff_o)))
1357  END IF
1358 
1359  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
1360  nrow_global=homo, ncol_global=virtual)
1361  CALL cp_fm_create(mp2_env%ri_grad%L_jb(kspin), fm_struct_tmp, name="fm_L_jb")
1362  CALL cp_fm_struct_release(fm_struct_tmp)
1363 
1364  ! first Virtual
1365  CALL parallel_gemm('T', 'N', homo, virtual, dimen, 2.0_dp*factor, &
1366  l_mu_q, mo_coeff, 0.0_dp, mp2_env%ri_grad%L_jb(kspin), &
1367  a_first_col=1, &
1368  a_first_row=1, &
1369  b_first_col=homo + 1, &
1370  b_first_row=1, &
1371  c_first_col=1, &
1372  c_first_row=1)
1373  ! then occupied
1374  CALL parallel_gemm('T', 'N', homo, virtual, dimen, 2.0_dp*factor, &
1375  mo_coeff, l_mu_q, 1.0_dp, mp2_env%ri_grad%L_jb(kspin), &
1376  a_first_col=1, &
1377  a_first_row=1, &
1378  b_first_col=homo + 1, &
1379  b_first_row=1, &
1380  c_first_col=1, &
1381  c_first_row=1)
1382 
1383  ! finally release L_mu_q
1384  CALL cp_fm_release(l_mu_q)
1385  CALL timestop(handle2)
1386 
1387  ! here we should be done next CPHF
1388 
1389  DEALLOCATE (pos_info)
1390 
1391  CALL timestop(handle)
1392 
1393  END SUBROUTINE create_w_p
1394 
1395 END MODULE mp2_ri_grad
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition: grid_common.h:117
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Handles all functions related to the CELL.
Definition: cell_types.F:15
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
integer function, public cp_fm_indxl2g(INDXLOC, NB, IPROC, ISRCPROC, NPROCS)
wrapper to scalapack function INDXL2G that computes the global index of a distributed matrix entry po...
Definition: cp_fm_types.F:2585
integer function, public cp_fm_indxg2l(INDXGLOB, NB, IPROC, ISRCPROC, NPROCS)
wrapper to scalapack function INDXG2L that computes the local index of a distributed matrix entry poi...
Definition: cp_fm_types.F:2525
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
Definition: cp_fm_types.F:1016
integer function, public cp_fm_indxg2p(INDXGLOB, NB, IPROC, ISRCPROC, NPROCS)
wrapper to scalapack function INDXG2P that computes the process coordinate which possesses the entry ...
Definition: cp_fm_types.F:2466
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
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
Definition: cp_fm_types.F:535
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_eri_mme
integer, parameter, public ri_rpa_method_gpw
integer, parameter, public ri_mp2_method_gpw
integer, parameter, public ri_mp2_laplace
integer, parameter, public do_eri_gpw
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition: libint_2c_3c.F:14
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare libint_potential_types.
Definition: libint_2c_3c.F:963
Interface to the message passing library MPI.
subroutine, public mp_para_env_release(para_env)
releases the para object (to be called when you don't want anymore the shared copy of this object)
type(mp_request_type), parameter, public mp_request_null
Routines to calculate 2c- and 3c-integrals for RI with GPW.
Definition: mp2_eri_gpw.F:13
subroutine, public cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_L)
Cleanup GPW integration for RI-MP2/RI-RPA.
Definition: mp2_eri_gpw.F:1060
subroutine, public integrate_potential_forces_3c_1c(mat_munu, rho_r, matrix_P_munu, qs_env, pw_env_sub, task_list_sub)
Takes the precomputed potential of an RI wave-function and determines matrix element and gradients wi...
Definition: mp2_eri_gpw.F:472
subroutine, public integrate_potential_forces_2c(rho_r, LLL, matrix, rho_g, atomic_kind_set, qs_kind_set, particle_set, cell, pw_env_sub, poisson_env, pot_g, potential_parameter, use_virial, rho_g_copy, dvg, kind_of, atom_of_kind, G_PQ_local, force, h_stress, para_env_sub, dft_control, psi_L, factor)
Integrates the potential of a RI function obtaining the forces and stress tensor.
Definition: mp2_eri_gpw.F:383
subroutine, public integrate_potential_forces_3c_2c(matrix_P_munu, rho_r, rho_g, task_list_sub, pw_env_sub, potential_parameter, ks_env, poisson_env, pot_g, use_virial, rho_g_copy, dvg, h_stress, para_env_sub, kind_of, atom_of_kind, qs_kind_set, particle_set, cell, LLL, force, dft_control)
Integrates potential of two Gaussians to a potential.
Definition: mp2_eri_gpw.F:527
subroutine, public prepare_gpw(qs_env, dft_control, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, auxbas_pw_pool, poisson_env, task_list_sub, rho_r, rho_g, pot_g, psi_L, sab_orb_sub)
Prepares GPW calculation for RI-MP2/RI-RPA.
Definition: mp2_eri_gpw.F:967
Interface to direct methods for electron repulsion integrals for MP2.
Definition: mp2_eri.F:12
subroutine, public mp2_eri_2c_integrate(param, potential_parameter, para_env, qs_env, basis_type_a, basis_type_b, hab, first_b, last_b, eri_method, pab, force_a, force_b, hdab, hadb, reflection_z_a, reflection_z_b, do_reflection_a, do_reflection_b)
high-level integration routine for 2c integrals over CP2K basis sets. Contiguous column-wise distribu...
Definition: mp2_eri.F:115
subroutine, public mp2_eri_3c_integrate(param, potential_parameter, para_env, qs_env, first_c, last_c, mat_ab, basis_type_a, basis_type_b, basis_type_c, sab_nl, eri_method, pabc, force_a, force_b, force_c, mat_dabc, mat_adbc, mat_abdc)
high-level integration routine for 3c integrals (ab|c) over CP2K basis sets. For each local function ...
Definition: mp2_eri.F:680
pure subroutine, public mp2_eri_deallocate_forces(force)
...
Definition: mp2_eri.F:5842
Routines to calculate gradients of RI-GPW-MP2 energy using pw.
Definition: mp2_ri_grad.F:13
subroutine, public calc_ri_mp2_nonsep(qs_env, mp2_env, para_env, para_env_sub, cell, particle_set, atomic_kind_set, qs_kind_set, mo_coeff, nmo, homo, dimen_RI, Eigenval, my_group_L_start, my_group_L_end, my_group_L_size, sab_orb_sub, mat_munu, blacs_env_sub)
Calculate the non-separable part of the gradients and update the Lagrangian.
Definition: mp2_ri_grad.F:114
Types needed for MP2 calculations.
Definition: mp2_types.F:14
basic linear algebra operations for full matrixes
Define the data structure for the particle information.
container for various plainwaves related things
Definition: pw_env_types.F:14
functions related to the poisson solver on regular grids
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
subroutine, public allocate_qs_force(qs_force, natom_of_kind)
Allocate a Quickstep force data structure.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
Define the neighbor list data types and the corresponding functionality.
types for task lists
All kind of helpful little routines.
Definition: util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333