Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Routines to calculate MP2 energy
10 : !> \par History
11 : !> 05.2011 created [Mauro Del Ben]
12 : !> \author Mauro Del Ben
13 : ! **************************************************************************************************
14 : MODULE mp2
15 : USE admm_types, ONLY: admm_type
16 : USE admm_utils, ONLY: admm_correct_for_eigenvalues,&
17 : admm_uncorrect_for_eigenvalues
18 : USE atomic_kind_types, ONLY: atomic_kind_type,&
19 : get_atomic_kind_set
20 : USE bibliography, ONLY: Bussy2023,&
21 : DelBen2012,&
22 : DelBen2015b,&
23 : Rybkin2016,&
24 : Stein2022,&
25 : Stein2024,&
26 : cite_reference
27 : USE cp_blacs_env, ONLY: cp_blacs_env_type
28 : USE cp_control_types, ONLY: dft_control_type
29 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
30 : dbcsr_create,&
31 : dbcsr_get_info,&
32 : dbcsr_p_type
33 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34 : dbcsr_allocate_matrix_set
35 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
36 : cp_fm_syrk,&
37 : cp_fm_triangular_invert,&
38 : cp_fm_uplo_to_full
39 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose
40 : USE cp_fm_diag, ONLY: choose_eigv_solver
41 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
42 : cp_fm_struct_release,&
43 : cp_fm_struct_type
44 : USE cp_fm_types, ONLY: cp_fm_create,&
45 : cp_fm_get_submatrix,&
46 : cp_fm_release,&
47 : cp_fm_set_all,&
48 : cp_fm_to_fm,&
49 : cp_fm_type
50 : USE cp_log_handling, ONLY: cp_get_default_logger,&
51 : cp_logger_type
52 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
53 : cp_print_key_unit_nr
54 : USE exstates_types, ONLY: excited_energy_type
55 : USE hfx_exx, ONLY: calculate_exx
56 : USE hfx_types, ONLY: &
57 : alloc_containers, dealloc_containers, hfx_basis_info_type, hfx_basis_type, &
58 : hfx_container_type, hfx_create_basis_types, hfx_init_container, hfx_release_basis_types, &
59 : hfx_type
60 : USE input_constants, ONLY: cholesky_inverse,&
61 : cholesky_off,&
62 : do_eri_gpw,&
63 : do_eri_mme,&
64 : rpa_exchange_axk,&
65 : rpa_exchange_none,&
66 : rpa_exchange_sosex,&
67 : sigma_none
68 : USE input_section_types, ONLY: section_vals_get,&
69 : section_vals_get_subs_vals,&
70 : section_vals_type
71 : USE kinds, ONLY: dp,&
72 : int_8
73 : USE kpoint_types, ONLY: kpoint_type
74 : USE machine, ONLY: m_flush,&
75 : m_memory,&
76 : m_walltime
77 : USE message_passing, ONLY: mp_para_env_type
78 : USE mp2_direct_method, ONLY: mp2_direct_energy
79 : USE mp2_gpw, ONLY: mp2_gpw_main
80 : USE mp2_optimize_ri_basis, ONLY: optimize_ri_basis_main
81 : USE mp2_types, ONLY: mp2_biel_type,&
82 : mp2_method_direct,&
83 : mp2_method_gpw,&
84 : mp2_ri_optimize_basis,&
85 : mp2_type,&
86 : ri_mp2_laplace,&
87 : ri_mp2_method_gpw,&
88 : ri_rpa_method_gpw
89 : USE parallel_gemm_api, ONLY: parallel_gemm
90 : USE particle_types, ONLY: particle_type
91 : USE qs_energy_types, ONLY: qs_energy_type
92 : USE qs_environment_types, ONLY: get_qs_env,&
93 : qs_environment_type
94 : USE qs_kind_types, ONLY: qs_kind_type
95 : USE qs_mo_types, ONLY: allocate_mo_set,&
96 : deallocate_mo_set,&
97 : get_mo_set,&
98 : init_mo_set,&
99 : mo_set_type
100 : USE qs_scf_methods, ONLY: eigensolver,&
101 : eigensolver_symm
102 : USE qs_scf_types, ONLY: qs_scf_env_type
103 : USE rpa_gw_sigma_x, ONLY: compute_vec_Sigma_x_minus_vxc_gw
104 : USE scf_control_types, ONLY: scf_control_type
105 : USE virial_types, ONLY: virial_type
106 :
107 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
108 :
109 : #include "./base/base_uses.f90"
110 :
111 : IMPLICIT NONE
112 :
113 : PRIVATE
114 :
115 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2'
116 :
117 : PUBLIC :: mp2_main
118 :
119 : CONTAINS
120 :
121 : ! **************************************************************************************************
122 : !> \brief the main entry point for MP2 calculations
123 : !> \param qs_env ...
124 : !> \param calc_forces ...
125 : !> \author Mauro Del Ben
126 : ! **************************************************************************************************
127 696 : SUBROUTINE mp2_main(qs_env, calc_forces)
128 : TYPE(qs_environment_type), POINTER :: qs_env
129 : LOGICAL, INTENT(IN) :: calc_forces
130 :
131 : CHARACTER(len=*), PARAMETER :: routineN = 'mp2_main'
132 :
133 : INTEGER :: bin, cholesky_method, dimen, handle, handle2, i, i_thread, iatom, ii, ikind, &
134 : irep, ispin, max_nset, my_bin_size, n_rep_hf, n_threads, nao, natom, ndep, &
135 : nfullcols_total, nfullrows_total, nkind, nmo, nspins, unit_nr
136 : INTEGER(KIND=int_8) :: mem
137 696 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of, nelec
138 : LOGICAL :: calc_ex, do_admm, do_admm_rpa_exx, do_dynamic_load_balancing, do_exx, do_gw, &
139 : do_im_time, do_kpoints_cubic_RPA, free_hfx_buffer, reuse_hfx, update_xc_energy
140 : REAL(KIND=dp) :: E_admm_from_GW(2), E_ex_from_GW, Emp2, Emp2_AA, Emp2_AA_Cou, Emp2_AA_ex, &
141 : Emp2_AB, Emp2_AB_Cou, Emp2_AB_ex, Emp2_BB, Emp2_BB_Cou, Emp2_BB_ex, Emp2_Cou, Emp2_ex, &
142 : Emp2_S, Emp2_T, maxocc, mem_real, t1, t2, t3
143 696 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
144 696 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: Auto
145 696 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: C
146 696 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
147 : TYPE(admm_type), POINTER :: admm_env
148 696 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
149 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
150 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
151 : TYPE(cp_fm_type) :: evecs, fm_matrix_ks, fm_matrix_s, &
152 : fm_matrix_work
153 : TYPE(cp_fm_type), POINTER :: fm_matrix_ks_red, fm_matrix_s_red, &
154 : fm_work_red, mo_coeff
155 : TYPE(cp_logger_type), POINTER :: logger
156 696 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
157 696 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_transl, matrix_s_kp
158 : TYPE(dft_control_type), POINTER :: dft_control
159 : TYPE(excited_energy_type), POINTER :: ex_env
160 : TYPE(hfx_basis_info_type) :: basis_info
161 696 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
162 696 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
163 : TYPE(hfx_container_type), POINTER :: maxval_container
164 : TYPE(hfx_type), POINTER :: actual_x_data
165 : TYPE(kpoint_type), POINTER :: kpoints
166 696 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mos_mp2
167 696 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
168 696 : TYPE(mp2_biel_type) :: mp2_biel
169 : TYPE(mp2_type), POINTER :: mp2_env
170 : TYPE(mp_para_env_type), POINTER :: para_env
171 696 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
172 : TYPE(qs_energy_type), POINTER :: energy
173 696 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
174 : TYPE(qs_scf_env_type), POINTER :: scf_env
175 : TYPE(scf_control_type), POINTER :: scf_control
176 : TYPE(section_vals_type), POINTER :: hfx_sections, input
177 : TYPE(virial_type), POINTER :: virial
178 :
179 : ! If SCF has not converged we should abort MP2 calculation
180 696 : IF (qs_env%mp2_env%hf_fail) THEN
181 : CALL cp_abort(__LOCATION__, "SCF not converged: "// &
182 0 : "not possible to run MP2")
183 : END IF
184 :
185 696 : NULLIFY (virial, dft_control, blacs_env, kpoints, fm_matrix_s_red, fm_matrix_ks_red, fm_work_red)
186 696 : CALL timeset(routineN, handle)
187 696 : logger => cp_get_default_logger()
188 :
189 696 : CALL cite_reference(DelBen2012)
190 :
191 696 : do_kpoints_cubic_RPA = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
192 :
193 : ! for cubic RPA and GW, we have kpoints and therefore, we get other matrices later
194 696 : IF (do_kpoints_cubic_RPA) THEN
195 :
196 : CALL get_qs_env(qs_env, &
197 : input=input, &
198 : atomic_kind_set=atomic_kind_set, &
199 : qs_kind_set=qs_kind_set, &
200 : dft_control=dft_control, &
201 : particle_set=particle_set, &
202 : para_env=para_env, &
203 : blacs_env=blacs_env, &
204 : energy=energy, &
205 : kpoints=kpoints, &
206 : scf_env=scf_env, &
207 : scf_control=scf_control, &
208 : matrix_ks_kp=matrix_ks_transl, &
209 : matrix_s_kp=matrix_s_kp, &
210 4 : mp2_env=mp2_env)
211 :
212 : CALL get_gamma(matrix_s, matrix_ks, mos, &
213 4 : matrix_s_kp, matrix_ks_transl, kpoints)
214 :
215 : ELSE
216 :
217 : CALL get_qs_env(qs_env, &
218 : input=input, &
219 : atomic_kind_set=atomic_kind_set, &
220 : qs_kind_set=qs_kind_set, &
221 : dft_control=dft_control, &
222 : particle_set=particle_set, &
223 : para_env=para_env, &
224 : blacs_env=blacs_env, &
225 : energy=energy, &
226 : mos=mos, &
227 : scf_env=scf_env, &
228 : scf_control=scf_control, &
229 : virial=virial, &
230 : matrix_ks=matrix_ks, &
231 : matrix_s=matrix_s, &
232 : mp2_env=mp2_env, &
233 692 : admm_env=admm_env)
234 :
235 : END IF
236 :
237 : ! IF DO_BSE In TDDFT, SAVE ks_matrix to ex_env
238 696 : IF (dft_control%tddfpt2_control%do_bse .OR. dft_control%tddfpt2_control%do_bse_w_only .OR. &
239 : dft_control%tddfpt2_control%do_bse_gw_only) THEN
240 4 : NULLIFY (ex_env)
241 4 : CALL get_qs_env(qs_env, exstate_env=ex_env)
242 4 : nspins = 1 ! for now only open-shell
243 4 : CALL dbcsr_allocate_matrix_set(ex_env%matrix_ks, nspins)
244 8 : DO ispin = 1, nspins
245 4 : ALLOCATE (ex_env%matrix_ks(ispin)%matrix)
246 4 : CALL dbcsr_create(ex_env%matrix_ks(ispin)%matrix, template=matrix_s(1)%matrix)
247 700 : CALL dbcsr_copy(ex_env%matrix_ks(ispin)%matrix, matrix_ks(ispin)%matrix)
248 : END DO
249 : END IF
250 :
251 : unit_nr = cp_print_key_unit_nr(logger, input, "DFT%XC%WF_CORRELATION%PRINT", &
252 696 : extension=".mp2Log")
253 :
254 696 : IF (unit_nr > 0) THEN
255 348 : IF (mp2_env%method /= ri_rpa_method_gpw) THEN
256 225 : WRITE (unit_nr, *)
257 225 : WRITE (unit_nr, *)
258 225 : WRITE (unit_nr, '(T2,A)') 'MP2 section'
259 225 : WRITE (unit_nr, '(T2,A)') '-----------'
260 225 : WRITE (unit_nr, *)
261 : ELSE
262 123 : WRITE (unit_nr, *)
263 123 : WRITE (unit_nr, *)
264 123 : WRITE (unit_nr, '(T2,A)') 'RI-RPA section'
265 123 : WRITE (unit_nr, '(T2,A)') '--------------'
266 123 : WRITE (unit_nr, *)
267 : END IF
268 : END IF
269 :
270 696 : IF (calc_forces) THEN
271 322 : CALL cite_reference(DelBen2015b)
272 322 : CALL cite_reference(Rybkin2016)
273 322 : CALL cite_reference(Stein2022)
274 322 : CALL cite_reference(Bussy2023)
275 322 : CALL cite_reference(Stein2024)
276 : !Gradients available for RI-MP2, and low-scaling Laplace MP2/RPA
277 322 : IF (.NOT. (mp2_env%method == ri_mp2_method_gpw .OR. &
278 : mp2_env%method == ri_rpa_method_gpw .OR. mp2_env%method == ri_mp2_laplace)) THEN
279 0 : CPABORT("No forces/gradients for the selected method.")
280 : END IF
281 : END IF
282 :
283 696 : IF (.NOT. do_kpoints_cubic_RPA) THEN
284 692 : IF (virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. mp2_env%eri_method == do_eri_mme) THEN
285 0 : CPABORT("analytical stress not implemented with ERI_METHOD = MME")
286 : END IF
287 : END IF
288 :
289 696 : IF (mp2_env%do_im_time .AND. mp2_env%eri_method /= do_eri_gpw) THEN
290 122 : mp2_env%mp2_num_proc = 1
291 : END IF
292 :
293 696 : IF (mp2_env%mp2_num_proc < 1 .OR. mp2_env%mp2_num_proc > para_env%num_pe) THEN
294 0 : CPWARN("GROUP_SIZE is reset because of a too small or too large value.")
295 0 : mp2_env%mp2_num_proc = MAX(MIN(para_env%num_pe, mp2_env%mp2_num_proc), 1)
296 : END IF
297 :
298 696 : IF (MOD(para_env%num_pe, mp2_env%mp2_num_proc) /= 0) THEN
299 0 : CPABORT("GROUP_SIZE must be a divisor of the total number of MPI ranks!")
300 : END IF
301 :
302 696 : IF (.NOT. mp2_env%do_im_time) THEN
303 562 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T76,I5)') 'Used number of processes per group:', mp2_env%mp2_num_proc
304 843 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'Maximum allowed memory usage per MPI process:', &
305 562 : mp2_env%mp2_memory, ' MiB'
306 : END IF
307 :
308 : IF ((mp2_env%method /= mp2_method_gpw) .AND. &
309 : (mp2_env%method /= ri_mp2_method_gpw) .AND. &
310 696 : (mp2_env%method /= ri_rpa_method_gpw) .AND. &
311 : (mp2_env%method /= ri_mp2_laplace)) THEN
312 24 : CALL m_memory(mem)
313 24 : mem_real = (mem + 1024*1024 - 1)/(1024*1024)
314 24 : CALL para_env%max(mem_real)
315 24 : mp2_env%mp2_memory = MIN(mp2_env%mp2_memory, mem_real)
316 24 : IF (mp2_env%mp2_memory < 0.0_dp) mp2_env%mp2_memory = 1.0_dp
317 :
318 36 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T68,F9.2,A4)') 'Available memory per MPI process for MP2:', &
319 24 : mp2_env%mp2_memory, ' MiB'
320 : END IF
321 :
322 696 : IF (unit_nr > 0) CALL m_flush(unit_nr)
323 :
324 696 : nspins = dft_control%nspins
325 696 : natom = SIZE(particle_set, 1)
326 :
327 696 : CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
328 696 : nkind = SIZE(atomic_kind_set, 1)
329 :
330 696 : do_admm_rpa_exx = mp2_env%ri_rpa%do_admm
331 696 : IF (do_admm_rpa_exx .AND. .NOT. dft_control%do_admm) THEN
332 0 : CPABORT("Need an ADMM input section for ADMM RI_RPA EXX to work")
333 : END IF
334 : IF (do_admm_rpa_exx) THEN
335 18 : CALL hfx_create_basis_types(basis_parameter, basis_info, qs_kind_set, "AUX_FIT")
336 : ELSE
337 678 : CALL hfx_create_basis_types(basis_parameter, basis_info, qs_kind_set, "ORB")
338 : END IF
339 :
340 696 : dimen = 0
341 696 : max_nset = 0
342 2722 : DO iatom = 1, natom
343 2026 : ikind = kind_of(iatom)
344 7844 : dimen = dimen + SUM(basis_parameter(ikind)%nsgf)
345 2722 : max_nset = MAX(max_nset, basis_parameter(ikind)%nset)
346 : END DO
347 :
348 696 : CALL get_mo_set(mo_set=mos(1), nao=nao)
349 :
350 : ! diagonalize the KS matrix in order to have the full set of MO's
351 : ! get S and KS matrices in fm_type (create also a working array)
352 696 : NULLIFY (fm_struct)
353 696 : CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
354 : CALL cp_fm_struct_create(fm_struct, context=blacs_env, nrow_global=nfullrows_total, &
355 696 : ncol_global=nfullcols_total, para_env=para_env)
356 696 : CALL cp_fm_create(fm_matrix_s, fm_struct, name="fm_matrix_s")
357 696 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, fm_matrix_s)
358 :
359 696 : CALL cp_fm_create(fm_matrix_ks, fm_struct, name="fm_matrix_ks")
360 :
361 696 : CALL cp_fm_create(fm_matrix_work, fm_struct, name="fm_matrix_work")
362 696 : CALL cp_fm_set_all(matrix=fm_matrix_work, alpha=0.0_dp)
363 :
364 696 : CALL cp_fm_struct_release(fm_struct)
365 :
366 696 : nmo = nao
367 2088 : ALLOCATE (nelec(nspins))
368 696 : IF (scf_env%cholesky_method == cholesky_off) THEN
369 42 : ALLOCATE (evals(nao))
370 296 : evals = 0
371 :
372 14 : CALL cp_fm_create(evecs, fm_matrix_s%matrix_struct)
373 :
374 : ! Perform an EVD
375 14 : CALL choose_eigv_solver(fm_matrix_s, evecs, evals)
376 :
377 : ! Determine the number of neglectable eigenvalues assuming that the eigenvalues are in ascending order
378 : ! (Required by Lapack)
379 14 : ndep = 0
380 40 : DO ii = 1, nao
381 40 : IF (evals(ii) > scf_control%eps_eigval) THEN
382 14 : ndep = ii - 1
383 14 : EXIT
384 : END IF
385 : END DO
386 14 : nmo = nao - ndep
387 :
388 30 : DO ispin = 1, nspins
389 30 : CALL get_mo_set(mo_set=mos(ispin), nelectron=nelec(ispin))
390 : END DO
391 30 : IF (MAXVAL(nelec)/(3 - nspins) > nmo) THEN
392 : ! Should not happen as the following MO calculation is the same as during the SCF steps
393 0 : CPABORT("Not enough MOs found!")
394 : END IF
395 :
396 : ! Set the eigenvalue of the eigenvectors belonging to the linear subspace to zero
397 40 : evals(1:ndep) = 0.0_dp
398 : ! Determine the eigenvalues of the inverse square root
399 270 : evals(ndep + 1:nao) = 1.0_dp/SQRT(evals(ndep + 1:nao))
400 :
401 14 : IF (ndep > 0) THEN
402 14 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T76,I5)') 'Number of removed MOs:', ndep
403 14 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T76,I5)') 'Number of available MOs:', nmo
404 :
405 : ! Create reduced matrices
406 14 : NULLIFY (fm_struct)
407 14 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_matrix_s%matrix_struct, ncol_global=nmo)
408 :
409 14 : ALLOCATE (fm_matrix_s_red, fm_work_red)
410 14 : CALL cp_fm_create(fm_matrix_s_red, fm_struct)
411 14 : CALL cp_fm_create(fm_work_red, fm_struct)
412 14 : CALL cp_fm_struct_release(fm_struct)
413 :
414 14 : ALLOCATE (fm_matrix_ks_red)
415 : CALL cp_fm_struct_create(fm_struct, template_fmstruct=fm_matrix_s%matrix_struct, &
416 14 : nrow_global=nmo, ncol_global=nmo)
417 14 : CALL cp_fm_create(fm_matrix_ks_red, fm_struct)
418 14 : CALL cp_fm_struct_release(fm_struct)
419 :
420 : ! Scale the eigenvalues and copy them to
421 14 : CALL cp_fm_to_fm(evecs, fm_matrix_s_red, nmo, ndep + 1)
422 14 : CALL cp_fm_column_scale(fm_matrix_s_red, evals(ndep + 1:))
423 :
424 : ! Obtain ortho from (P)DGEMM, skip the linear dependent columns
425 : CALL parallel_gemm("N", "T", nao, nao, nmo, 1.0_dp, fm_matrix_s_red, evecs, &
426 14 : 0.0_dp, fm_matrix_s, b_first_col=ndep + 1)
427 : ELSE
428 : ! Take the square roots of the target values to allow application of SYRK
429 0 : evals = SQRT(evals)
430 0 : CALL cp_fm_column_scale(evecs, evals)
431 0 : CALL cp_fm_syrk("U", "N", nao, 1.0_dp, evecs, 1, 1, 0.0_dp, fm_matrix_s)
432 0 : CALL cp_fm_uplo_to_full(fm_matrix_s, fm_matrix_work)
433 : END IF
434 :
435 14 : CALL cp_fm_release(evecs)
436 28 : cholesky_method = cholesky_off
437 : ELSE
438 : ! calculate S^(-1/2) (cholesky decomposition)
439 682 : CALL cp_fm_cholesky_decompose(fm_matrix_s)
440 682 : CALL cp_fm_triangular_invert(fm_matrix_s)
441 682 : cholesky_method = cholesky_inverse
442 : END IF
443 :
444 2942 : ALLOCATE (mos_mp2(nspins))
445 1550 : DO ispin = 1, nspins
446 :
447 854 : CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc, nelectron=nelec(ispin))
448 :
449 : CALL allocate_mo_set(mo_set=mos_mp2(ispin), &
450 : nao=nao, &
451 : nmo=nmo, &
452 : nelectron=nelec(ispin), &
453 : n_el_f=REAL(nelec(ispin), dp), &
454 : maxocc=maxocc, &
455 854 : flexible_electron_count=dft_control%relax_multiplicity)
456 :
457 854 : CALL get_mo_set(mos_mp2(ispin), nao=nao)
458 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
459 : ncol_global=nmo, para_env=para_env, &
460 854 : context=blacs_env)
461 :
462 : CALL init_mo_set(mos_mp2(ispin), &
463 : fm_struct=fm_struct, &
464 854 : name="mp2_mos")
465 3258 : CALL cp_fm_struct_release(fm_struct)
466 : END DO
467 :
468 1550 : DO ispin = 1, nspins
469 :
470 : ! If ADMM we should make the ks matrix up-to-date
471 854 : IF (dft_control%do_admm) THEN
472 94 : CALL admm_correct_for_eigenvalues(ispin, admm_env, matrix_ks(ispin)%matrix)
473 : END IF
474 :
475 854 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, fm_matrix_ks)
476 :
477 854 : IF (dft_control%do_admm) THEN
478 94 : CALL admm_uncorrect_for_eigenvalues(ispin, admm_env, matrix_ks(ispin)%matrix)
479 : END IF
480 :
481 854 : IF (cholesky_method == cholesky_inverse) THEN
482 :
483 : ! diagonalize KS matrix
484 : CALL eigensolver(matrix_ks_fm=fm_matrix_ks, &
485 : mo_set=mos_mp2(ispin), &
486 : ortho=fm_matrix_s, &
487 : work=fm_matrix_work, &
488 : cholesky_method=cholesky_method, &
489 : do_level_shift=.FALSE., &
490 : level_shift=0.0_dp, &
491 838 : use_jacobi=.FALSE.)
492 :
493 16 : ELSE IF (cholesky_method == cholesky_off) THEN
494 :
495 16 : IF (ASSOCIATED(fm_matrix_s_red)) THEN
496 : CALL eigensolver_symm(matrix_ks_fm=fm_matrix_ks, &
497 : mo_set=mos_mp2(ispin), &
498 : ortho=fm_matrix_s, &
499 : work=fm_matrix_work, &
500 : do_level_shift=.FALSE., &
501 : level_shift=0.0_dp, &
502 : use_jacobi=.FALSE., &
503 : jacobi_threshold=0.0_dp, &
504 : ortho_red=fm_matrix_s_red, &
505 : matrix_ks_fm_red=fm_matrix_ks_red, &
506 16 : work_red=fm_work_red)
507 : ELSE
508 : CALL eigensolver_symm(matrix_ks_fm=fm_matrix_ks, &
509 : mo_set=mos_mp2(ispin), &
510 : ortho=fm_matrix_s, &
511 : work=fm_matrix_work, &
512 : do_level_shift=.FALSE., &
513 : level_shift=0.0_dp, &
514 : use_jacobi=.FALSE., &
515 0 : jacobi_threshold=0.0_dp)
516 : END IF
517 : END IF
518 :
519 1550 : CALL get_mo_set(mos_mp2(ispin), mo_coeff=mo_coeff)
520 : END DO
521 :
522 696 : CALL cp_fm_release(fm_matrix_s)
523 696 : CALL cp_fm_release(fm_matrix_ks)
524 696 : CALL cp_fm_release(fm_matrix_work)
525 696 : IF (ASSOCIATED(fm_matrix_s_red)) THEN
526 14 : CALL cp_fm_release(fm_matrix_s_red)
527 14 : DEALLOCATE (fm_matrix_s_red)
528 : END IF
529 696 : IF (ASSOCIATED(fm_matrix_ks_red)) THEN
530 14 : CALL cp_fm_release(fm_matrix_ks_red)
531 14 : DEALLOCATE (fm_matrix_ks_red)
532 : END IF
533 696 : IF (ASSOCIATED(fm_work_red)) THEN
534 14 : CALL cp_fm_release(fm_work_red)
535 14 : DEALLOCATE (fm_work_red)
536 : END IF
537 :
538 696 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
539 :
540 : ! build the table of index
541 696 : t1 = m_walltime()
542 2784 : ALLOCATE (mp2_biel%index_table(natom, max_nset))
543 :
544 696 : CALL build_index_table(natom, max_nset, mp2_biel%index_table, basis_parameter, kind_of)
545 :
546 : ! free the hfx_container (for now if forces are required we don't release the HFX stuff)
547 696 : free_hfx_buffer = .FALSE.
548 696 : IF (ASSOCIATED(qs_env%x_data)) THEN
549 428 : free_hfx_buffer = .TRUE.
550 428 : IF (calc_forces .AND. (.NOT. mp2_env%ri_grad%free_hfx_buffer)) free_hfx_buffer = .FALSE.
551 428 : IF (qs_env%x_data(1, 1)%do_hfx_ri) free_hfx_buffer = .FALSE.
552 428 : IF (calc_forces .AND. mp2_env%do_im_time) free_hfx_buffer = .FALSE.
553 428 : IF (mp2_env%ri_rpa%reuse_hfx) free_hfx_buffer = .FALSE.
554 : END IF
555 :
556 696 : IF (.NOT. do_kpoints_cubic_RPA) THEN
557 692 : IF (virial%pv_numer) THEN
558 : ! in the case of numerical stress we don't have to free the HFX integrals
559 72 : free_hfx_buffer = .FALSE.
560 72 : mp2_env%ri_grad%free_hfx_buffer = free_hfx_buffer
561 : END IF
562 : END IF
563 :
564 : ! calculate the matrix sigma_x - vxc for G0W0
565 696 : t3 = 0
566 696 : IF (mp2_env%ri_rpa%do_ri_g0w0) THEN
567 108 : CALL compute_vec_Sigma_x_minus_vxc_gw(qs_env, mp2_env, mos_mp2, E_ex_from_GW, E_admm_from_GW, t3, unit_nr)
568 : END IF
569 :
570 696 : IF (free_hfx_buffer) THEN
571 224 : CALL timeset(routineN//"_free_hfx", handle2)
572 224 : CALL section_vals_get(hfx_sections, n_repetition=n_rep_hf)
573 224 : n_threads = 1
574 224 : !$ n_threads = omp_get_max_threads()
575 :
576 448 : DO irep = 1, n_rep_hf
577 672 : DO i_thread = 0, n_threads - 1
578 224 : actual_x_data => qs_env%x_data(irep, i_thread + 1)
579 :
580 224 : do_dynamic_load_balancing = .TRUE.
581 224 : IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.
582 :
583 : IF (do_dynamic_load_balancing) THEN
584 0 : my_bin_size = SIZE(actual_x_data%distribution_energy)
585 : ELSE
586 224 : my_bin_size = 1
587 : END IF
588 :
589 448 : IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly) THEN
590 222 : CALL dealloc_containers(actual_x_data%store_ints, actual_x_data%memory_parameter%actual_memory_usage)
591 : END IF
592 : END DO
593 : END DO
594 224 : CALL timestop(handle2)
595 : END IF
596 :
597 696 : Emp2 = 0.D+00
598 696 : Emp2_Cou = 0.D+00
599 696 : Emp2_ex = 0.D+00
600 696 : calc_ex = .TRUE.
601 :
602 696 : t1 = m_walltime()
603 714 : SELECT CASE (mp2_env%method)
604 : CASE (mp2_method_direct)
605 18 : IF (unit_nr > 0) WRITE (unit_nr, *)
606 :
607 72 : ALLOCATE (Auto(dimen, nspins))
608 90 : ALLOCATE (C(dimen, dimen, nspins))
609 :
610 40 : DO ispin = 1, nspins
611 : ! get the alpha coeff and eigenvalues
612 : CALL get_mo_set(mo_set=mos_mp2(ispin), &
613 : eigenvalues=mo_eigenvalues, &
614 22 : mo_coeff=mo_coeff)
615 :
616 22 : CALL cp_fm_get_submatrix(mo_coeff, C(:, :, ispin), 1, 1, dimen, dimen, .FALSE.)
617 1072 : Auto(:, ispin) = mo_eigenvalues(:)
618 : END DO
619 :
620 18 : IF (nspins == 2) THEN
621 4 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Unrestricted Canonical Direct Methods:'
622 : ! for now, require the mos to be always present
623 :
624 : ! calculate the alpha-alpha MP2
625 4 : Emp2_AA = 0.0_dp
626 4 : Emp2_AA_Cou = 0.0_dp
627 4 : Emp2_AA_ex = 0.0_dp
628 : CALL mp2_direct_energy(dimen, nelec(1), nelec(1), mp2_biel, &
629 : mp2_env, C(:, :, 1), Auto(:, 1), Emp2_AA, Emp2_AA_Cou, Emp2_AA_ex, &
630 4 : qs_env, para_env, unit_nr)
631 4 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Energy Alpha-Alpha = ', Emp2_AA
632 4 : IF (unit_nr > 0) WRITE (unit_nr, *)
633 :
634 4 : Emp2_BB = 0.0_dp
635 4 : Emp2_BB_Cou = 0.0_dp
636 4 : Emp2_BB_ex = 0.0_dp
637 : CALL mp2_direct_energy(dimen, nelec(2), nelec(2), mp2_biel, mp2_env, &
638 : C(:, :, 2), Auto(:, 2), Emp2_BB, Emp2_BB_Cou, Emp2_BB_ex, &
639 4 : qs_env, para_env, unit_nr)
640 4 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Energy Beta-Beta= ', Emp2_BB
641 4 : IF (unit_nr > 0) WRITE (unit_nr, *)
642 :
643 4 : Emp2_AB = 0.0_dp
644 4 : Emp2_AB_Cou = 0.0_dp
645 4 : Emp2_AB_ex = 0.0_dp
646 : CALL mp2_direct_energy(dimen, nelec(1), nelec(2), mp2_biel, mp2_env, C(:, :, 1), &
647 : Auto(:, 1), Emp2_AB, Emp2_AB_Cou, Emp2_AB_ex, &
648 4 : qs_env, para_env, unit_nr, C(:, :, 2), Auto(:, 2))
649 4 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Energy Alpha-Beta= ', Emp2_AB
650 4 : IF (unit_nr > 0) WRITE (unit_nr, *)
651 :
652 4 : Emp2 = Emp2_AA + Emp2_BB + Emp2_AB*2.0_dp !+Emp2_BA
653 4 : Emp2_Cou = Emp2_AA_Cou + Emp2_BB_Cou + Emp2_AB_Cou*2.0_dp !+Emp2_BA
654 4 : Emp2_ex = Emp2_AA_ex + Emp2_BB_ex + Emp2_AB_ex*2.0_dp !+Emp2_BA
655 :
656 4 : Emp2_S = Emp2_AB*2.0_dp
657 4 : Emp2_T = Emp2_AA + Emp2_BB
658 :
659 : ELSE
660 :
661 14 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A)') 'Canonical Direct Methods:'
662 :
663 : CALL mp2_direct_energy(dimen, nelec(1)/2, nelec(1)/2, mp2_biel, mp2_env, &
664 : C(:, :, 1), Auto(:, 1), Emp2, Emp2_Cou, Emp2_ex, &
665 14 : qs_env, para_env, unit_nr)
666 :
667 : END IF
668 :
669 18 : DEALLOCATE (C, Auto)
670 :
671 : CASE (mp2_ri_optimize_basis)
672 : ! optimize ri basis set or tests for RI-MP2 gradients
673 6 : IF (unit_nr > 0) THEN
674 3 : WRITE (unit_nr, *)
675 3 : WRITE (unit_nr, '(T3,A)') 'Optimization of the auxiliary RI-MP2 basis'
676 3 : WRITE (unit_nr, *)
677 : END IF
678 :
679 24 : ALLOCATE (Auto(dimen, nspins))
680 30 : ALLOCATE (C(dimen, dimen, nspins))
681 :
682 18 : DO ispin = 1, nspins
683 : ! get the alpha coeff and eigenvalues
684 : CALL get_mo_set(mo_set=mos_mp2(ispin), &
685 : eigenvalues=mo_eigenvalues, &
686 12 : mo_coeff=mo_coeff)
687 :
688 12 : CALL cp_fm_get_submatrix(mo_coeff, C(:, :, ispin), 1, 1, dimen, dimen, .FALSE.)
689 174 : Auto(:, ispin) = mo_eigenvalues(:)
690 : END DO
691 :
692 : ! optimize basis
693 6 : IF (nspins == 2) THEN
694 : CALL optimize_ri_basis_main(Emp2, Emp2_Cou, Emp2_ex, Emp2_S, Emp2_T, dimen, natom, nelec(1), &
695 : mp2_biel, mp2_env, C(:, :, 1), Auto(:, 1), &
696 : kind_of, qs_env, para_env, unit_nr, &
697 6 : nelec(2), C(:, :, 2), Auto(:, 2))
698 :
699 : ELSE
700 : CALL optimize_ri_basis_main(Emp2, Emp2_Cou, Emp2_ex, Emp2_S, Emp2_T, dimen, natom, nelec(1)/2, &
701 : mp2_biel, mp2_env, C(:, :, 1), Auto(:, 1), &
702 0 : kind_of, qs_env, para_env, unit_nr)
703 : END IF
704 :
705 6 : DEALLOCATE (Auto, C)
706 :
707 : CASE (mp2_method_gpw)
708 : ! check if calculate the exchange contribution
709 14 : IF (mp2_env%scale_T == 0.0_dp .AND. (nspins == 2)) calc_ex = .FALSE.
710 :
711 : ! go with mp2_gpw
712 : CALL mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T, &
713 368 : mos_mp2, para_env, unit_nr, calc_forces, calc_ex)
714 :
715 : CASE (ri_mp2_method_gpw)
716 : ! check if calculate the exchange contribution
717 354 : IF (mp2_env%scale_T == 0.0_dp .AND. (nspins == 2)) calc_ex = .FALSE.
718 :
719 : ! go with mp2_gpw
720 : CALL mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T, &
721 354 : mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_mp2=.TRUE.)
722 :
723 : CASE (ri_rpa_method_gpw)
724 : ! perform RI-RPA energy calculation (since most part of the calculation
725 : ! is actually equal to the RI-MP2-GPW we decided to put RPA in the MP2
726 : ! section to avoid code replication)
727 :
728 246 : calc_ex = .FALSE.
729 :
730 : ! go with ri_rpa_gpw
731 : CALL mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T, &
732 246 : mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_rpa=.TRUE.)
733 : ! Scale energy contributions
734 246 : Emp2 = Emp2*mp2_env%ri_rpa%scale_rpa
735 246 : mp2_env%ri_rpa%ener_exchange = mp2_env%ri_rpa%ener_exchange*mp2_env%ri_rpa%scale_rpa
736 :
737 : CASE (ri_mp2_laplace)
738 : ! perform RI-SOS-Laplace-MP2 energy calculation, most part of the code in common
739 : ! with the RI-RPA part
740 :
741 : ! In SOS-MP2 only the coulomb-like contribution of the MP2 energy is computed
742 58 : calc_ex = .FALSE.
743 :
744 : ! go with sos_laplace_mp2_gpw
745 : CALL mp2_gpw_main(qs_env, mp2_env, Emp2, Emp2_Cou, Emp2_EX, Emp2_S, Emp2_T, &
746 58 : mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_sos_laplace_mp2=.TRUE.)
747 :
748 : CASE DEFAULT
749 710 : CPABORT("")
750 : END SELECT
751 :
752 696 : t2 = m_walltime()
753 696 : IF (unit_nr > 0) WRITE (unit_nr, *)
754 696 : IF (mp2_env%method /= ri_rpa_method_gpw) THEN
755 450 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total MP2 Time=', t2 - t1
756 450 : IF (mp2_env%method == ri_mp2_laplace) THEN
757 58 : Emp2_S = Emp2
758 58 : Emp2_T = 0.0_dp
759 58 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Energy SO component (singlet) = ', Emp2_S
760 58 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Scaling factor SO = ', mp2_env%scale_S
761 : ELSE
762 392 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Coulomb Energy = ', Emp2_Cou/2.0_dp
763 392 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Exchange Energy = ', Emp2_ex
764 392 : IF (nspins == 1) THEN
765 : ! valid only in the closed shell case
766 292 : Emp2_S = Emp2_Cou/2.0_dp
767 292 : IF (calc_ex) THEN
768 292 : Emp2_T = Emp2_ex + Emp2_Cou/2.0_dp
769 : ELSE
770 : ! unknown if Emp2_ex is not computed
771 0 : Emp2_T = 0.0_dp
772 : END IF
773 : END IF
774 392 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Energy SO component (singlet) = ', Emp2_S
775 392 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'MP2 Energy SS component (triplet) = ', Emp2_T
776 392 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Scaling factor SO = ', mp2_env%scale_S
777 392 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Scaling factor SS = ', mp2_env%scale_T
778 : END IF
779 450 : Emp2_S = Emp2_S*mp2_env%scale_S
780 450 : Emp2_T = Emp2_T*mp2_env%scale_T
781 450 : Emp2 = Emp2_S + Emp2_T
782 450 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Second order perturbation energy = ', Emp2
783 : ELSE
784 246 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.6)') 'Total RI-RPA Time=', t2 - t1
785 :
786 246 : update_xc_energy = .TRUE.
787 246 : IF (mp2_env%ri_rpa%do_ri_g0w0 .AND. .NOT. mp2_env%ri_g0w0%update_xc_energy) update_xc_energy = .FALSE.
788 82 : IF (.NOT. update_xc_energy) Emp2 = 0.0_dp
789 :
790 246 : IF (unit_nr > 0 .AND. update_xc_energy) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA energy = ', Emp2
791 246 : IF (unit_nr > 0 .AND. mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
792 5 : WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Sigma corr. to RI-RPA energy = ', &
793 10 : mp2_env%ri_rpa%e_sigma_corr
794 : END IF
795 246 : IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_axk) THEN
796 10 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA-AXK energy=', mp2_env%ri_rpa%ener_exchange
797 236 : ELSE IF (mp2_env%ri_rpa%exchange_correction == rpa_exchange_sosex) THEN
798 2 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'RI-RPA-SOSEX energy=', mp2_env%ri_rpa%ener_exchange
799 : END IF
800 246 : IF (mp2_env%ri_rpa%do_rse) THEN
801 9 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Diagonal singles correction (dRSE) = ', &
802 6 : mp2_env%ri_rpa%rse_corr_diag
803 9 : IF (unit_nr > 0) WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Full singles correction (RSE) =', &
804 6 : mp2_env%ri_rpa%rse_corr
805 6 : IF (dft_control%do_admm) CPABORT("RPA RSE not implemented with RI_RPA%ADMM on")
806 : END IF
807 : END IF
808 696 : IF (unit_nr > 0) WRITE (unit_nr, *)
809 :
810 : ! we have it !!!!
811 696 : IF (mp2_env%ri_rpa%exchange_correction /= rpa_exchange_none) THEN
812 12 : Emp2 = Emp2 + mp2_env%ri_rpa%ener_exchange
813 : END IF
814 696 : IF (mp2_env%ri_rpa%do_rse) THEN
815 6 : Emp2 = Emp2 + mp2_env%ri_rpa%rse_corr
816 : END IF
817 696 : IF (mp2_env%ri_rpa%sigma_param /= sigma_none) THEN
818 : !WRITE (unit_nr, '(T3,A,T56,F25.14)') 'Sigma corr. to RI-RPA energy = ',&
819 10 : Emp2 = Emp2 + mp2_env%ri_rpa%e_sigma_corr
820 : END IF
821 696 : energy%mp2 = Emp2
822 696 : energy%total = energy%total + Emp2
823 :
824 1550 : DO ispin = 1, nspins
825 1550 : CALL deallocate_mo_set(mo_set=mos_mp2(ispin))
826 : END DO
827 696 : DEALLOCATE (mos_mp2)
828 :
829 : ! if necessary reallocate hfx buffer
830 696 : IF (free_hfx_buffer .AND. (.NOT. calc_forces) .AND. &
831 : (mp2_env%ri_g0w0%do_ri_Sigma_x .OR. .NOT. mp2_env%ri_rpa_im_time%do_kpoints_from_Gamma)) THEN
832 120 : CALL timeset(routineN//"_alloc_hfx", handle2)
833 240 : DO irep = 1, n_rep_hf
834 360 : DO i_thread = 0, n_threads - 1
835 120 : actual_x_data => qs_env%x_data(irep, i_thread + 1)
836 :
837 120 : do_dynamic_load_balancing = .TRUE.
838 120 : IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.
839 :
840 : IF (do_dynamic_load_balancing) THEN
841 0 : my_bin_size = SIZE(actual_x_data%distribution_energy)
842 : ELSE
843 120 : my_bin_size = 1
844 : END IF
845 :
846 240 : IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly) THEN
847 120 : CALL alloc_containers(actual_x_data%store_ints, my_bin_size)
848 :
849 240 : DO bin = 1, my_bin_size
850 120 : maxval_container => actual_x_data%store_ints%maxval_container(bin)
851 120 : integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
852 120 : CALL hfx_init_container(maxval_container, actual_x_data%memory_parameter%actual_memory_usage, .FALSE.)
853 7920 : DO i = 1, 64
854 7800 : CALL hfx_init_container(integral_containers(i), actual_x_data%memory_parameter%actual_memory_usage, .FALSE.)
855 : END DO
856 : END DO
857 : END IF
858 : END DO
859 : END DO
860 120 : CALL timestop(handle2)
861 : END IF
862 :
863 696 : CALL hfx_release_basis_types(basis_parameter)
864 :
865 : ! if required calculate the EXX contribution from the DFT density
866 696 : IF (mp2_env%method == ri_rpa_method_gpw .AND. .NOT. calc_forces) THEN
867 : do_exx = .FALSE.
868 194 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
869 194 : CALL section_vals_get(hfx_sections, explicit=do_exx)
870 194 : IF (do_exx) THEN
871 132 : do_gw = mp2_env%ri_rpa%do_ri_g0w0
872 132 : do_admm = mp2_env%ri_rpa%do_admm
873 132 : reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
874 132 : do_im_time = qs_env%mp2_env%do_im_time
875 :
876 : CALL calculate_exx(qs_env=qs_env, &
877 : unit_nr=unit_nr, &
878 : hfx_sections=hfx_sections, &
879 : x_data=qs_env%mp2_env%ri_rpa%x_data, &
880 : do_gw=do_gw, &
881 : do_admm=do_admm, &
882 : calc_forces=.FALSE., &
883 : reuse_hfx=reuse_hfx, &
884 : do_im_time=do_im_time, &
885 : E_ex_from_GW=E_ex_from_GW, &
886 : E_admm_from_GW=E_admm_from_GW, &
887 132 : t3=t3)
888 :
889 : END IF
890 : END IF
891 :
892 : CALL cp_print_key_finished_output(unit_nr, logger, input, &
893 696 : "DFT%XC%WF_CORRELATION%PRINT")
894 :
895 696 : CALL timestop(handle)
896 :
897 3480 : END SUBROUTINE mp2_main
898 :
899 : ! **************************************************************************************************
900 : !> \brief ...
901 : !> \param natom ...
902 : !> \param max_nset ...
903 : !> \param index_table ...
904 : !> \param basis_parameter ...
905 : !> \param kind_of ...
906 : ! **************************************************************************************************
907 696 : PURE SUBROUTINE build_index_table(natom, max_nset, index_table, basis_parameter, kind_of)
908 : INTEGER, INTENT(IN) :: natom, max_nset
909 : INTEGER, DIMENSION(natom, max_nset), INTENT(OUT) :: index_table
910 : TYPE(hfx_basis_type), DIMENSION(:), POINTER :: basis_parameter
911 : INTEGER, DIMENSION(natom), INTENT(IN) :: kind_of
912 :
913 : INTEGER :: counter, iatom, ikind, iset, nset
914 :
915 8784 : index_table = -HUGE(0)
916 : counter = 0
917 2722 : DO iatom = 1, natom
918 2026 : ikind = kind_of(iatom)
919 2026 : nset = basis_parameter(ikind)%nset
920 8540 : DO iset = 1, nset
921 5818 : index_table(iatom, iset) = counter + 1
922 7844 : counter = counter + basis_parameter(ikind)%nsgf(iset)
923 : END DO
924 : END DO
925 :
926 696 : END SUBROUTINE build_index_table
927 :
928 : ! **************************************************************************************************
929 : !> \brief ...
930 : !> \param matrix_s ...
931 : !> \param matrix_ks ...
932 : !> \param mos ...
933 : !> \param matrix_s_kp ...
934 : !> \param matrix_ks_transl ...
935 : !> \param kpoints ...
936 : ! **************************************************************************************************
937 4 : PURE SUBROUTINE get_gamma(matrix_s, matrix_ks, mos, matrix_s_kp, matrix_ks_transl, kpoints)
938 :
939 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_ks
940 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
941 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrix_ks_transl
942 : TYPE(kpoint_type), POINTER :: kpoints
943 :
944 : INTEGER :: nspins
945 :
946 4 : nspins = SIZE(matrix_ks_transl, 1)
947 :
948 4 : matrix_ks(1:nspins) => matrix_ks_transl(1:nspins, 1)
949 4 : matrix_s(1:1) => matrix_s_kp(1:1, 1)
950 4 : mos(1:nspins) => kpoints%kp_env(1)%kpoint_env%mos(1:nspins, 1)
951 :
952 4 : END SUBROUTINE get_gamma
953 :
954 : END MODULE mp2
955 :
|
