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