(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE mp2_ri_grad
14 USE atomic_kind_types, ONLY: atomic_kind_type,&
15 get_atomic_kind_set
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
19 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
20 dbcsr_deallocate_matrix_set
21 USE cp_eri_mme_interface, ONLY: cp_eri_mme_param
22 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
23 cp_fm_struct_release,&
24 cp_fm_struct_type
25 USE cp_fm_types, ONLY: &
26 cp_fm_create, cp_fm_get_info, cp_fm_indxg2l, cp_fm_indxg2p, cp_fm_indxl2g, cp_fm_release, &
27 cp_fm_set_all, cp_fm_to_fm_submat, cp_fm_type
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,&
34 ri_mp2_laplace,&
35 ri_mp2_method_gpw,&
36 ri_rpa_method_gpw
37 USE kinds, ONLY: dp
38 USE libint_2c_3c, ONLY: compare_potential_types
39 USE message_passing, ONLY: mp_para_env_release,&
40 mp_para_env_type,&
41 mp_request_null,&
42 mp_request_type,&
43 mp_waitall
44 USE mp2_eri, ONLY: mp2_eri_2c_integrate,&
45 mp2_eri_3c_integrate,&
46 mp2_eri_deallocate_forces,&
47 mp2_eri_force
48 USE mp2_eri_gpw, ONLY: cleanup_gpw,&
49 integrate_potential_forces_2c,&
50 integrate_potential_forces_3c_1c,&
51 integrate_potential_forces_3c_2c,&
52 prepare_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
65 USE qs_force_types, ONLY: allocate_qs_force,&
66 qs_force_type,&
67 zero_qs_force
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
84CONTAINS
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
1395END MODULE mp2_ri_grad
modulo
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
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:89
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
message_passing::mp_waitall
Definition message_passing.F:871
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
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.
Definition atomic_kind_types.F:223
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_blacs_env
methods related to the blacs parallel environment
Definition cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:208
cp_eri_mme_interface
Interface to Minimax-Ewald method for periodic ERI's to be used in CP2K.
Definition cp_eri_mme_interface.F:12
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
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
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_indxl2g
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
cp_fm_types::cp_fm_get_info
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
cp_fm_types::cp_fm_indxg2l
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
cp_fm_types::cp_fm_to_fm_submat
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
cp_fm_types::cp_fm_set_all
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
cp_fm_types::cp_fm_indxg2p
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
cp_fm_types::cp_fm_create
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
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_eri_mme
integer, parameter, public do_eri_mme
Definition input_constants.F:1226
input_constants::ri_rpa_method_gpw
integer, parameter, public ri_rpa_method_gpw
Definition input_constants.F:1046
input_constants::ri_mp2_method_gpw
integer, parameter, public ri_mp2_method_gpw
Definition input_constants.F:1046
input_constants::ri_mp2_laplace
integer, parameter, public ri_mp2_laplace
Definition input_constants.F:1046
input_constants::do_eri_gpw
integer, parameter, public do_eri_gpw
Definition input_constants.F:1226
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
libint_2c_3c
2- and 3-center electron repulsion integral routines based on libint2 Currently available operators: ...
Definition libint_2c_3c.F:14
libint_2c_3c::compare_potential_types
pure logical function, public compare_potential_types(potential1, potential2)
Helper function to compare libint_potential_types.
Definition libint_2c_3c.F:963
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
message_passing::mp_para_env_release
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)
Definition message_passing.F:2027
message_passing::mp_request_null
type(mp_request_type), parameter, public mp_request_null
Definition message_passing.F:814
mp2_eri_gpw
Routines to calculate 2c- and 3c-integrals for RI with GPW.
Definition mp2_eri_gpw.F:13
mp2_eri_gpw::integrate_potential_forces_3c_2c
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
mp2_eri_gpw::prepare_gpw
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
mp2_eri_gpw::integrate_potential_forces_2c
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
mp2_eri_gpw::integrate_potential_forces_3c_1c
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
mp2_eri_gpw::cleanup_gpw
subroutine, public cleanup_gpw(qs_env, e_cutoff_old, cutoff_old, relative_cutoff_old, para_env_sub, pw_env_sub, task_list_sub, auxbas_pw_pool, rho_r, rho_g, pot_g, psi_l)
Cleanup GPW integration for RI-MP2/RI-RPA.
Definition mp2_eri_gpw.F:1060
mp2_eri
Interface to direct methods for electron repulsion integrals for MP2.
Definition mp2_eri.F:12
mp2_eri::mp2_eri_2c_integrate
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
mp2_eri::mp2_eri_3c_integrate
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
mp2_eri::mp2_eri_deallocate_forces
pure subroutine, public mp2_eri_deallocate_forces(force)
...
Definition mp2_eri.F:5842
mp2_ri_grad
Routines to calculate gradients of RI-GPW-MP2 energy using pw.
Definition mp2_ri_grad.F:13
mp2_ri_grad::calc_ri_mp2_nonsep
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
mp2_types
Types needed for MP2 calculations.
Definition mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_pool_types
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
pw_types
Definition pw_types.F:24
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
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.
Definition qs_environment_types.F:502
qs_force_types
Definition qs_force_types.F:13
qs_force_types::zero_qs_force
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
Definition qs_force_types.F:250
qs_force_types::allocate_qs_force
subroutine, public allocate_qs_force(qs_force, natom_of_kind)
Allocate a Quickstep force data structure.
Definition qs_force_types.F:79
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_ks_types
Definition qs_ks_types.F:15
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
task_list_types
types for task lists
Definition task_list_types.F:14
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
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
virial_types
Definition virial_types.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_eri_mme_interface::cp_eri_mme_param
Definition cp_eri_mme_interface.F:79
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
message_passing::mp_request_type
Definition message_passing.F:623
mp2_eri::mp2_eri_force
Definition mp2_eri.F:83
mp2_types::integ_mat_buffer_type_2d
Definition mp2_types.F:385
mp2_types::integ_mat_buffer_type
Definition mp2_types.F:377
mp2_types::mp2_type
Definition mp2_types.F:317
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:122
pw_pool_types::pw_pool_type
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:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
task_list_types::task_list_type
Definition task_list_types.F:58
virial_types::virial_type
Definition virial_types.F:26