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