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