114 LOGICAL,
INTENT(IN) :: calc_forces
116 CHARACTER(len=*),
PARAMETER :: routinen =
'mp2_main'
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 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 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: auto
129 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: c
130 REAL(kind=
dp),
DIMENSION(:),
POINTER :: mo_eigenvalues
135 TYPE(
cp_fm_type) :: fm_matrix_ks, fm_matrix_s, fm_matrix_work
138 TYPE(dbcsr_p_type),
DIMENSION(:),
POINTER :: matrix_ks, matrix_s
139 TYPE(dbcsr_p_type),
DIMENSION(:, :),
POINTER :: matrix_ks_transl, matrix_s_kp
145 TYPE(
hfx_type),
POINTER :: actual_x_data
147 TYPE(
mo_set_type),
ALLOCATABLE,
DIMENSION(:) :: mos_mp2
154 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
161 IF (qs_env%mp2_env%hf_fail)
THEN
162 CALL cp_abort(__location__,
"SCF not converged: "// &
163 "not possible to run MP2")
166 NULLIFY (virial, dft_control, blacs_env, kpoints)
167 CALL timeset(routinen, handle)
172 do_kpoints_cubic_rpa = qs_env%mp2_env%ri_rpa_im_time%do_im_time_kpoints
175 IF (do_kpoints_cubic_rpa)
THEN
179 atomic_kind_set=atomic_kind_set, &
180 qs_kind_set=qs_kind_set, &
181 dft_control=dft_control, &
182 particle_set=particle_set, &
184 blacs_env=blacs_env, &
188 scf_control=scf_control, &
189 matrix_ks_kp=matrix_ks_transl, &
190 matrix_s_kp=matrix_s_kp, &
193 CALL get_gamma(matrix_s, matrix_ks, mos, &
194 matrix_s_kp, matrix_ks_transl, kpoints)
200 atomic_kind_set=atomic_kind_set, &
201 qs_kind_set=qs_kind_set, &
202 dft_control=dft_control, &
203 particle_set=particle_set, &
205 blacs_env=blacs_env, &
209 scf_control=scf_control, &
211 matrix_ks=matrix_ks, &
221 IF (unit_nr > 0)
THEN
225 WRITE (unit_nr,
'(T2,A)')
'MP2 section'
226 WRITE (unit_nr,
'(T2,A)')
'-----------'
231 WRITE (unit_nr,
'(T2,A)')
'RI-RPA section'
232 WRITE (unit_nr,
'(T2,A)')
'--------------'
237 IF (calc_forces)
THEN
242 cpabort(
"No forces/gradients for the selected method.")
246 IF (.NOT. do_kpoints_cubic_rpa)
THEN
247 IF (virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. mp2_env%eri_method ==
do_eri_mme)
THEN
248 cpabort(
"analytical stress not implemented with ERI_METHOD = MME")
252 IF (mp2_env%do_im_time .AND. mp2_env%eri_method .NE.
do_eri_gpw)
THEN
253 mp2_env%mp2_num_proc = 1
256 IF (mp2_env%mp2_num_proc < 1 .OR. mp2_env%mp2_num_proc > para_env%num_pe)
THEN
257 cpwarn(
"GROUP_SIZE is reset because of a too small or too large value.")
258 mp2_env%mp2_num_proc = max(min(para_env%num_pe, mp2_env%mp2_num_proc), 1)
261 IF (mod(para_env%num_pe, mp2_env%mp2_num_proc) /= 0)
THEN
262 cpabort(
"GROUP_SIZE must be a divisor of the total number of MPI ranks!")
265 IF (.NOT. mp2_env%do_im_time)
THEN
266 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T76,I5)')
'Used number of processes per group:', mp2_env%mp2_num_proc
267 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T68,F9.2,A4)')
'Maximum allowed memory usage per MPI process:', &
268 mp2_env%mp2_memory,
' MiB'
276 mem_real = (mem + 1024*1024 - 1)/(1024*1024)
277 CALL para_env%max(mem_real)
278 mp2_env%mp2_memory = mp2_env%mp2_memory - mem_real
279 IF (mp2_env%mp2_memory < 0.0_dp) mp2_env%mp2_memory = 1.0_dp
281 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T68,F9.2,A4)')
'Available memory per MPI process for MP2:', &
282 mp2_env%mp2_memory,
' MiB'
285 IF (unit_nr > 0)
CALL m_flush(unit_nr)
287 nspins = dft_control%nspins
288 natom =
SIZE(particle_set, 1)
291 nkind =
SIZE(atomic_kind_set, 1)
293 do_admm_rpa_exx = mp2_env%ri_rpa%do_admm
294 IF (do_admm_rpa_exx .AND. .NOT. dft_control%do_admm)
THEN
295 cpabort(
"Need an ADMM input section for ADMM RI_RPA EXX to work")
297 IF (do_admm_rpa_exx)
THEN
306 ikind = kind_of(iatom)
307 dimen = dimen + sum(basis_parameter(ikind)%nsgf)
308 max_nset = max(max_nset, basis_parameter(ikind)%nset)
316 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nfullrows_total, nfullcols_total=nfullcols_total)
318 ncol_global=nfullcols_total, para_env=para_env)
319 CALL cp_fm_create(fm_matrix_s, fm_struct, name=
"fm_matrix_s")
322 CALL cp_fm_create(fm_matrix_ks, fm_struct, name=
"fm_matrix_ks")
324 CALL cp_fm_create(fm_matrix_work, fm_struct, name=
"fm_matrix_work")
329 CALL cp_fm_get_info(matrix=fm_matrix_ks, nrow_block=nrow_block, ncol_block=ncol_block)
332 CALL cp_fm_power(fm_matrix_s, fm_matrix_work, -0.5_dp, scf_control%eps_eigval, ndep)
341 ALLOCATE (mos_mp2(nspins))
342 ALLOCATE (nelec(nspins))
346 IF (dft_control%do_admm)
THEN
352 IF (dft_control%do_admm)
THEN
356 CALL get_mo_set(mo_set=mos(ispin), maxocc=maxocc, nelectron=nelec(ispin))
361 nelectron=nelec(ispin), &
362 n_el_f=real(nelec(ispin),
dp), &
364 flexible_electron_count=dft_control%relax_multiplicity)
368 ncol_global=nao, para_env=para_env, &
372 fm_struct=fm_struct, &
380 mo_set=mos_mp2(ispin), &
382 work=fm_matrix_work, &
383 cholesky_method=cholesky_method, &
384 do_level_shift=.false., &
385 level_shift=0.0_dp, &
391 mo_set=mos_mp2(ispin), &
393 work=fm_matrix_work, &
394 do_level_shift=.false., &
395 level_shift=0.0_dp, &
396 use_jacobi=.false., &
397 jacobi_threshold=0.0_dp)
411 ALLOCATE (mp2_biel%index_table(natom, max_nset))
413 CALL build_index_table(natom, max_nset, mp2_biel%index_table, basis_parameter, kind_of)
416 free_hfx_buffer = .false.
417 IF (
ASSOCIATED(qs_env%x_data))
THEN
418 free_hfx_buffer = .true.
419 IF (calc_forces .AND. (.NOT. mp2_env%ri_grad%free_hfx_buffer)) free_hfx_buffer = .false.
420 IF (qs_env%x_data(1, 1)%do_hfx_ri) free_hfx_buffer = .false.
421 IF (calc_forces .AND. mp2_env%do_im_time) free_hfx_buffer = .false.
422 IF (mp2_env%ri_rpa%reuse_hfx) free_hfx_buffer = .false.
425 IF (.NOT. do_kpoints_cubic_rpa)
THEN
426 IF (virial%pv_numer)
THEN
428 free_hfx_buffer = .false.
429 mp2_env%ri_grad%free_hfx_buffer = free_hfx_buffer
435 IF (mp2_env%ri_rpa%do_ri_g0w0)
THEN
439 IF (free_hfx_buffer)
THEN
440 CALL timeset(routinen//
"_free_hfx", handle2)
445 DO irep = 1, n_rep_hf
446 DO i_thread = 0, n_threads - 1
447 actual_x_data => qs_env%x_data(irep, i_thread + 1)
449 do_dynamic_load_balancing = .true.
450 IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .false.
452 IF (do_dynamic_load_balancing)
THEN
453 my_bin_size =
SIZE(actual_x_data%distribution_energy)
458 IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly)
THEN
459 CALL dealloc_containers(actual_x_data%store_ints, actual_x_data%memory_parameter%actual_memory_usage)
463 CALL timestop(handle2)
472 SELECT CASE (mp2_env%method)
474 IF (unit_nr > 0)
WRITE (unit_nr, *)
476 ALLOCATE (auto(dimen, nspins))
477 ALLOCATE (c(dimen, dimen, nspins))
482 eigenvalues=mo_eigenvalues, &
486 auto(:, ispin) = mo_eigenvalues(:)
489 IF (nspins == 2)
THEN
490 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A)')
'Unrestricted Canonical Direct Methods:'
498 mp2_env, c(:, :, 1), auto(:, 1), emp2_aa, emp2_aa_cou, emp2_aa_ex, &
499 qs_env, para_env, unit_nr)
500 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Energy Alpha-Alpha = ', emp2_aa
501 IF (unit_nr > 0)
WRITE (unit_nr, *)
507 c(:, :, 2), auto(:, 2), emp2_bb, emp2_bb_cou, emp2_bb_ex, &
508 qs_env, para_env, unit_nr)
509 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Energy Beta-Beta= ', emp2_bb
510 IF (unit_nr > 0)
WRITE (unit_nr, *)
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 qs_env, para_env, unit_nr, c(:, :, 2), auto(:, 2))
518 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Energy Alpha-Beta= ', emp2_ab
519 IF (unit_nr > 0)
WRITE (unit_nr, *)
521 emp2 = emp2_aa + emp2_bb + emp2_ab*2.0_dp
522 emp2_cou = emp2_aa_cou + emp2_bb_cou + emp2_ab_cou*2.0_dp
523 emp2_ex = emp2_aa_ex + emp2_bb_ex + emp2_ab_ex*2.0_dp
525 emp2_s = emp2_ab*2.0_dp
526 emp2_t = emp2_aa + emp2_bb
530 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A)')
'Canonical Direct Methods:'
533 c(:, :, 1), auto(:, 1), emp2, emp2_cou, emp2_ex, &
534 qs_env, para_env, unit_nr)
542 IF (unit_nr > 0)
THEN
544 WRITE (unit_nr,
'(T3,A)')
'Optimization of the auxiliary RI-MP2 basis'
548 ALLOCATE (auto(dimen, nspins))
549 ALLOCATE (c(dimen, dimen, nspins))
554 eigenvalues=mo_eigenvalues, &
558 auto(:, ispin) = mo_eigenvalues(:)
562 IF (nspins == 2)
THEN
564 mp2_biel, mp2_env, c(:, :, 1), auto(:, 1), &
565 kind_of, qs_env, para_env, unit_nr, &
566 nelec(2), c(:, :, 2), auto(:, 2))
570 mp2_biel, mp2_env, c(:, :, 1), auto(:, 1), &
571 kind_of, qs_env, para_env, unit_nr)
578 IF (mp2_env%scale_T == 0.0_dp .AND. (nspins == 2)) calc_ex = .false.
581 CALL mp2_gpw_main(qs_env, mp2_env, emp2, emp2_cou, emp2_ex, emp2_s, emp2_t, &
582 mos_mp2, para_env, unit_nr, calc_forces, calc_ex)
586 IF (mp2_env%scale_T == 0.0_dp .AND. (nspins == 2)) calc_ex = .false.
589 CALL mp2_gpw_main(qs_env, mp2_env, emp2, emp2_cou, emp2_ex, emp2_s, emp2_t, &
590 mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_mp2=.true.)
600 CALL mp2_gpw_main(qs_env, mp2_env, emp2, emp2_cou, emp2_ex, emp2_s, emp2_t, &
601 mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_rpa=.true.)
604 emp2 = emp2*mp2_env%ri_rpa%scale_rpa
605 mp2_env%ri_rpa%ener_axk = mp2_env%ri_rpa%ener_axk*mp2_env%ri_rpa%scale_rpa
615 CALL mp2_gpw_main(qs_env, mp2_env, emp2, emp2_cou, emp2_ex, emp2_s, emp2_t, &
616 mos_mp2, para_env, unit_nr, calc_forces, calc_ex, do_ri_sos_laplace_mp2=.true.)
623 IF (unit_nr > 0)
WRITE (unit_nr, *)
625 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.6)')
'Total MP2 Time=', t2 - t1
629 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Energy SO component (singlet) = ', emp2_s
630 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'Scaling factor SO = ', mp2_env%scale_S
632 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Coulomb Energy = ', emp2_cou/2.0_dp
633 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Exchange Energy = ', emp2_ex
634 IF (nspins == 1)
THEN
636 emp2_s = emp2_cou/2.0_dp
638 emp2_t = emp2_ex + emp2_cou/2.0_dp
644 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Energy SO component (singlet) = ', emp2_s
645 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'MP2 Energy SS component (triplet) = ', emp2_t
646 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'Scaling factor SO = ', mp2_env%scale_S
647 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'Scaling factor SS = ', mp2_env%scale_T
649 emp2_s = emp2_s*mp2_env%scale_S
650 emp2_t = emp2_t*mp2_env%scale_T
651 emp2 = emp2_s + emp2_t
652 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'Second order perturbation energy = ', emp2
654 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.6)')
'Total RI-RPA Time=', t2 - t1
656 update_xc_energy = .true.
657 IF (mp2_env%ri_rpa%do_ri_g0w0 .AND. .NOT. mp2_env%ri_g0w0%update_xc_energy) update_xc_energy = .false.
658 IF (.NOT. update_xc_energy) emp2 = 0.0_dp
660 IF (unit_nr > 0 .AND. update_xc_energy)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'RI-RPA energy = ', emp2
661 IF (mp2_env%ri_rpa%do_ri_axk)
THEN
662 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'RI-RPA-AXK energy=', mp2_env%ri_rpa%ener_axk
664 IF (mp2_env%ri_rpa%do_rse)
THEN
665 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'Diagonal singles correction (dRSE) = ', &
666 mp2_env%ri_rpa%rse_corr_diag
667 IF (unit_nr > 0)
WRITE (unit_nr,
'(T3,A,T56,F25.14)')
'Full singles correction (RSE) =', &
668 mp2_env%ri_rpa%rse_corr
669 IF (dft_control%do_admm) cpabort(
"RPA RSE not implemented with RI_RPA%ADMM on")
672 IF (unit_nr > 0)
WRITE (unit_nr, *)
675 IF (mp2_env%ri_rpa%do_ri_axk)
THEN
676 emp2 = emp2 + mp2_env%ri_rpa%ener_axk
678 IF (mp2_env%ri_rpa%do_rse)
THEN
679 emp2 = emp2 + mp2_env%ri_rpa%rse_corr
682 energy%total = energy%total + emp2
690 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 CALL timeset(routinen//
"_alloc_hfx", handle2)
693 DO irep = 1, n_rep_hf
694 DO i_thread = 0, n_threads - 1
695 actual_x_data => qs_env%x_data(irep, i_thread + 1)
697 do_dynamic_load_balancing = .true.
698 IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .false.
700 IF (do_dynamic_load_balancing)
THEN
701 my_bin_size =
SIZE(actual_x_data%distribution_energy)
706 IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly)
THEN
709 DO bin = 1, my_bin_size
710 maxval_container => actual_x_data%store_ints%maxval_container(bin)
711 integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
712 CALL hfx_init_container(maxval_container, actual_x_data%memory_parameter%actual_memory_usage, .false.)
714 CALL hfx_init_container(integral_containers(i), actual_x_data%memory_parameter%actual_memory_usage, .false.)
720 CALL timestop(handle2)
731 do_gw = mp2_env%ri_rpa%do_ri_g0w0
732 do_admm = mp2_env%ri_rpa%do_admm
733 reuse_hfx = qs_env%mp2_env%ri_rpa%reuse_hfx
734 do_im_time = qs_env%mp2_env%do_im_time
738 hfx_sections=hfx_sections, &
739 x_data=qs_env%mp2_env%ri_rpa%x_data, &
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, &
753 "DFT%XC%WF_CORRELATION%PRINT")
755 CALL timestop(handle)
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.