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