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