LCOV - code coverage report
Current view: top level - src - mp2_cphf.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 97.2 % 883 858
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 8 8

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Routines to calculate CPHF like update and solve Z-vector equation
      10              : !>        for MP2 gradients (only GPW)
      11              : !> \par History
      12              : !>      11.2013 created [Mauro Del Ben]
      13              : ! **************************************************************************************************
      14              : MODULE mp2_cphf
      15              :    USE admm_methods,                    ONLY: admm_projection_derivative
      16              :    USE admm_types,                      ONLY: admm_type,&
      17              :                                               get_admm_env
      18              :    USE atomic_kind_types,               ONLY: atomic_kind_type
      19              :    USE cell_types,                      ONLY: cell_type
      20              :    USE core_ae,                         ONLY: build_core_ae
      21              :    USE core_ppl,                        ONLY: build_core_ppl
      22              :    USE core_ppnl,                       ONLY: build_core_ppnl
      23              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      24              :    USE cp_control_types,                ONLY: dft_control_type
      25              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      26              :                                               dbcsr_copy,&
      27              :                                               dbcsr_p_type,&
      28              :                                               dbcsr_release,&
      29              :                                               dbcsr_scale,&
      30              :                                               dbcsr_set
      31              :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      32              :                                               cp_dbcsr_plus_fm_fm_t,&
      33              :                                               dbcsr_allocate_matrix_set,&
      34              :                                               dbcsr_deallocate_matrix_set
      35              :    USE cp_fm_basic_linalg,              ONLY: cp_fm_uplo_to_full
      36              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      37              :                                               cp_fm_struct_p_type,&
      38              :                                               cp_fm_struct_release,&
      39              :                                               cp_fm_struct_type
      40              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      41              :                                               cp_fm_get_info,&
      42              :                                               cp_fm_release,&
      43              :                                               cp_fm_set_all,&
      44              :                                               cp_fm_to_fm_submat,&
      45              :                                               cp_fm_type
      46              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      47              :                                               cp_logger_type
      48              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      49              :                                               cp_print_key_unit_nr
      50              :    USE hfx_admm_utils,                  ONLY: tddft_hfx_matrix
      51              :    USE hfx_derivatives,                 ONLY: derivatives_four_center
      52              :    USE hfx_exx,                         ONLY: add_exx_to_rhs
      53              :    USE hfx_ri,                          ONLY: hfx_ri_update_forces
      54              :    USE hfx_types,                       ONLY: alloc_containers,&
      55              :                                               hfx_container_type,&
      56              :                                               hfx_init_container,&
      57              :                                               hfx_type
      58              :    USE input_constants,                 ONLY: do_admm_aux_exch_func_none,&
      59              :                                               ot_precond_full_all,&
      60              :                                               z_solver_cg,&
      61              :                                               z_solver_pople,&
      62              :                                               z_solver_richardson,&
      63              :                                               z_solver_sd
      64              :    USE input_section_types,             ONLY: section_vals_get,&
      65              :                                               section_vals_get_subs_vals,&
      66              :                                               section_vals_type
      67              :    USE kahan_sum,                       ONLY: accurate_dot_product
      68              :    USE kinds,                           ONLY: dp
      69              :    USE linear_systems,                  ONLY: solve_system
      70              :    USE machine,                         ONLY: m_flush,&
      71              :                                               m_walltime
      72              :    USE mathconstants,                   ONLY: fourpi
      73              :    USE message_passing,                 ONLY: mp_para_env_type
      74              :    USE mp2_types,                       ONLY: mp2_type,&
      75              :                                               ri_rpa_method_gpw
      76              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      77              :    USE particle_types,                  ONLY: particle_type
      78              :    USE pw_env_types,                    ONLY: pw_env_get,&
      79              :                                               pw_env_type
      80              :    USE pw_methods,                      ONLY: pw_axpy,&
      81              :                                               pw_copy,&
      82              :                                               pw_derive,&
      83              :                                               pw_integral_ab,&
      84              :                                               pw_scale,&
      85              :                                               pw_transfer,&
      86              :                                               pw_zero
      87              :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      88              :    USE pw_poisson_types,                ONLY: pw_poisson_type
      89              :    USE pw_pool_types,                   ONLY: pw_pool_type
      90              :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      91              :                                               pw_r3d_rs_type
      92              :    USE qs_2nd_kernel_ao,                ONLY: apply_2nd_order_kernel
      93              :    USE qs_density_matrices,             ONLY: calculate_whz_matrix
      94              :    USE qs_dispersion_pairpot,           ONLY: calculate_dispersion_pairpot
      95              :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      96              :    USE qs_energy_types,                 ONLY: qs_energy_type
      97              :    USE qs_environment_types,            ONLY: get_qs_env,&
      98              :                                               qs_environment_type,&
      99              :                                               set_qs_env
     100              :    USE qs_force_types,                  ONLY: deallocate_qs_force,&
     101              :                                               qs_force_type,&
     102              :                                               sum_qs_force,&
     103              :                                               zero_qs_force
     104              :    USE qs_integrate_potential,          ONLY: integrate_v_core_rspace,&
     105              :                                               integrate_v_rspace
     106              :    USE qs_kind_types,                   ONLY: qs_kind_type
     107              :    USE qs_kinetic,                      ONLY: build_kinetic_matrix
     108              :    USE qs_ks_reference,                 ONLY: ks_ref_potential
     109              :    USE qs_ks_types,                     ONLY: qs_ks_env_type,&
     110              :                                               set_ks_env
     111              :    USE qs_linres_types,                 ONLY: linres_control_type
     112              :    USE qs_mo_types,                     ONLY: get_mo_set,&
     113              :                                               mo_set_type
     114              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
     115              :    USE qs_overlap,                      ONLY: build_overlap_matrix
     116              :    USE qs_p_env_methods,                ONLY: p_env_check_i_alloc,&
     117              :                                               p_env_create,&
     118              :                                               p_env_psi0_changed,&
     119              :                                               p_env_update_rho
     120              :    USE qs_p_env_types,                  ONLY: p_env_release,&
     121              :                                               qs_p_env_type
     122              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     123              :                                               qs_rho_type
     124              :    USE task_list_types,                 ONLY: task_list_type
     125              :    USE virial_types,                    ONLY: virial_type,&
     126              :                                               zero_virial
     127              : 
     128              : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
     129              : 
     130              : #include "./base/base_uses.f90"
     131              : 
     132              :    IMPLICIT NONE
     133              : 
     134              :    PRIVATE
     135              : 
     136              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_cphf'
     137              :    LOGICAL, PARAMETER, PRIVATE :: debug_forces = .TRUE.
     138              : 
     139              :    PUBLIC :: solve_z_vector_eq, update_mp2_forces
     140              : 
     141              : CONTAINS
     142              : 
     143              : ! **************************************************************************************************
     144              : !> \brief Solve Z-vector equations necessary for the calculation of the MP2
     145              : !>        gradients, in order to be consistent here the parameters for the
     146              : !>        calculation of the CPHF like updats have to be exactly equal to the
     147              : !>        SCF case
     148              : !> \param qs_env ...
     149              : !> \param mp2_env ...
     150              : !> \param para_env ...
     151              : !> \param dft_control ...
     152              : !> \param mo_coeff ...
     153              : !> \param homo ...
     154              : !> \param Eigenval ...
     155              : !> \param unit_nr ...
     156              : !> \author Mauro Del Ben, Vladimir Rybkin
     157              : ! **************************************************************************************************
     158          272 :    SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
     159          272 :                                 mo_coeff, homo, Eigenval, unit_nr)
     160              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     161              :       TYPE(mp2_type), INTENT(INOUT)                      :: mp2_env
     162              :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
     163              :       TYPE(dft_control_type), INTENT(IN), POINTER        :: dft_control
     164              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff
     165              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: homo
     166              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
     167              :       INTEGER, INTENT(IN)                                :: unit_nr
     168              : 
     169              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'solve_z_vector_eq'
     170              : 
     171              :       INTEGER :: bin, handle, handle2, i, i_global, i_thread, iiB, irep, ispin, j_global, jjB, &
     172              :          my_bin_size, n_rep_hf, n_threads, nao, ncol_local, nmo, nrow_local, nspins, transf_type_in
     173          272 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: virtual
     174          272 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     175              :       LOGICAL                                            :: alpha_beta, do_dynamic_load_balancing, &
     176              :                                                             do_exx, do_hfx, restore_p_screen
     177              :       REAL(KIND=dp)                                      :: focc
     178              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     179              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     180              :       TYPE(cp_fm_type)                                   :: fm_back, fm_G_mu_nu, fm_mo_mo
     181          272 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: L_jb, mo_coeff_o, mo_coeff_v, P_ia, &
     182          272 :                                                             P_mo, W_Mo
     183          272 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_p_mp2, &
     184          272 :                                                             matrix_p_mp2_admm, matrix_s, P_mu_nu, &
     185          272 :                                                             rho1_ao, rho_ao, rho_ao_aux_fit
     186          272 :       TYPE(hfx_container_type), DIMENSION(:), POINTER    :: integral_containers
     187              :       TYPE(hfx_container_type), POINTER                  :: maxval_container
     188              :       TYPE(hfx_type), POINTER                            :: actual_x_data
     189              :       TYPE(linres_control_type), POINTER                 :: linres_control
     190              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     191              :       TYPE(qs_p_env_type), POINTER                       :: p_env
     192              :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit
     193              :       TYPE(section_vals_type), POINTER                   :: hfx_section, hfx_sections, input
     194              : 
     195          272 :       CALL timeset(routineN, handle)
     196              : 
     197              :       ! start collecting stuff
     198          272 :       CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
     199          272 :       CPASSERT(SIZE(eigenval, 1) == nmo)
     200          272 :       CPASSERT(nao >= nmo)
     201          272 :       NULLIFY (input, matrix_s, blacs_env, rho, &
     202          272 :                matrix_p_mp2, matrix_p_mp2_admm, matrix_ks)
     203              :       CALL get_qs_env(qs_env, &
     204              :                       ks_env=ks_env, &
     205              :                       input=input, &
     206              :                       matrix_s=matrix_s, &
     207              :                       matrix_ks=matrix_ks, &
     208              :                       matrix_p_mp2=matrix_p_mp2, &
     209              :                       matrix_p_mp2_admm=matrix_p_mp2_admm, &
     210              :                       blacs_env=blacs_env, &
     211          272 :                       rho=rho)
     212              : 
     213          272 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     214              : 
     215              :       ! Get number of relevant spin states
     216          272 :       nspins = dft_control%nspins
     217          272 :       alpha_beta = (nspins == 2)
     218              : 
     219          272 :       CALL MOVE_ALLOC(mp2_env%ri_grad%P_mo, P_mo)
     220          272 :       CALL MOVE_ALLOC(mp2_env%ri_grad%W_mo, W_mo)
     221          272 :       CALL MOVE_ALLOC(mp2_env%ri_grad%L_jb, L_jb)
     222              : 
     223          816 :       ALLOCATE (virtual(nspins))
     224          630 :       virtual(:) = nmo - homo(:)
     225              : 
     226          272 :       NULLIFY (P_mu_nu)
     227          272 :       CALL dbcsr_allocate_matrix_set(P_mu_nu, nspins)
     228          630 :       DO ispin = 1, nspins
     229          358 :          ALLOCATE (P_mu_nu(ispin)%matrix)
     230          358 :          CALL dbcsr_copy(P_mu_nu(ispin)%matrix, rho_ao(1)%matrix, name="P_mu_nu")
     231          630 :          CALL dbcsr_set(P_mu_nu(ispin)%matrix, 0.0_dp)
     232              :       END DO
     233              : 
     234          272 :       NULLIFY (fm_struct_tmp)
     235              :       CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     236          272 :                                nrow_global=nao, ncol_global=nao)
     237          272 :       CALL cp_fm_create(fm_G_mu_nu, fm_struct_tmp, name="G_mu_nu")
     238          272 :       CALL cp_fm_create(fm_back, fm_struct_tmp, name="fm_back")
     239          272 :       CALL cp_fm_struct_release(fm_struct_tmp)
     240          272 :       CALL cp_fm_set_all(fm_G_mu_nu, 0.0_dp)
     241          272 :       CALL cp_fm_set_all(fm_back, 0.0_dp)
     242              : 
     243         1804 :       ALLOCATE (mo_coeff_o(nspins), mo_coeff_v(nspins))
     244          630 :       DO ispin = 1, nspins
     245          358 :          NULLIFY (fm_struct_tmp)
     246              :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     247          358 :                                   nrow_global=nao, ncol_global=homo(ispin))
     248          358 :          CALL cp_fm_create(mo_coeff_o(ispin), fm_struct_tmp, name="mo_coeff_o")
     249          358 :          CALL cp_fm_struct_release(fm_struct_tmp)
     250          358 :          CALL cp_fm_set_all(mo_coeff_o(ispin), 0.0_dp)
     251              :          CALL cp_fm_to_fm_submat(msource=mo_coeff(ispin), mtarget=mo_coeff_o(ispin), &
     252              :                                  nrow=nao, ncol=homo(ispin), &
     253              :                                  s_firstrow=1, s_firstcol=1, &
     254          358 :                                  t_firstrow=1, t_firstcol=1)
     255              : 
     256          358 :          NULLIFY (fm_struct_tmp)
     257              :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     258          358 :                                   nrow_global=nao, ncol_global=virtual(ispin))
     259          358 :          CALL cp_fm_create(mo_coeff_v(ispin), fm_struct_tmp, name="mo_coeff_v")
     260          358 :          CALL cp_fm_struct_release(fm_struct_tmp)
     261          358 :          CALL cp_fm_set_all(mo_coeff_v(ispin), 0.0_dp)
     262              :          CALL cp_fm_to_fm_submat(msource=mo_coeff(ispin), mtarget=mo_coeff_v(ispin), &
     263              :                                  nrow=nao, ncol=virtual(ispin), &
     264              :                                  s_firstrow=1, s_firstcol=homo(ispin) + 1, &
     265          630 :                                  t_firstrow=1, t_firstcol=1)
     266              :       END DO
     267              : 
     268              :       ! hfx section
     269          272 :       NULLIFY (hfx_sections)
     270          272 :       hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
     271          272 :       CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
     272          272 :       IF (do_hfx) THEN
     273              :          ! here we check if we have to reallocate the HFX container
     274          192 :          IF (mp2_env%ri_grad%free_hfx_buffer .AND. (.NOT. qs_env%x_data(1, 1)%do_hfx_ri)) THEN
     275          102 :             CALL timeset(routineN//"_alloc_hfx", handle2)
     276          102 :             n_threads = 1
     277          102 : !$          n_threads = omp_get_max_threads()
     278              : 
     279          204 :             DO irep = 1, n_rep_hf
     280          306 :                DO i_thread = 0, n_threads - 1
     281          102 :                   actual_x_data => qs_env%x_data(irep, i_thread + 1)
     282              : 
     283          102 :                   do_dynamic_load_balancing = .TRUE.
     284          102 :                   IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.
     285              : 
     286              :                   IF (do_dynamic_load_balancing) THEN
     287            0 :                      my_bin_size = SIZE(actual_x_data%distribution_energy)
     288              :                   ELSE
     289          102 :                      my_bin_size = 1
     290              :                   END IF
     291              : 
     292          204 :                   IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly) THEN
     293          102 :                      CALL alloc_containers(actual_x_data%store_ints, my_bin_size)
     294              : 
     295          204 :                      DO bin = 1, my_bin_size
     296          102 :                         maxval_container => actual_x_data%store_ints%maxval_container(bin)
     297          102 :                         integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
     298          102 :                         CALL hfx_init_container(maxval_container, actual_x_data%memory_parameter%actual_memory_usage, .FALSE.)
     299         6732 :                         DO i = 1, 64
     300              :                            CALL hfx_init_container(integral_containers(i), &
     301         6630 :                                                    actual_x_data%memory_parameter%actual_memory_usage, .FALSE.)
     302              :                         END DO
     303              :                      END DO
     304              :                   END IF
     305              :                END DO
     306              :             END DO
     307          102 :             CALL timestop(handle2)
     308              :          END IF
     309              : 
     310              :          ! set up parameters for P_screening
     311          192 :          restore_p_screen = qs_env%x_data(1, 1)%screening_parameter%do_initial_p_screening
     312          192 :          IF (qs_env%x_data(1, 1)%screening_parameter%do_initial_p_screening) THEN
     313            8 :             IF (mp2_env%ri_grad%free_hfx_buffer) THEN
     314            8 :                mp2_env%p_screen = .FALSE.
     315              :             ELSE
     316            0 :                mp2_env%p_screen = .TRUE.
     317              :             END IF
     318              :          END IF
     319              :       END IF
     320              : 
     321              :       ! Add exx part for RPA
     322          272 :       do_exx = .FALSE.
     323          272 :       IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
     324           24 :          hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
     325           24 :          CALL section_vals_get(hfx_section, explicit=do_exx)
     326              :       END IF
     327          272 :       IF (do_exx) THEN
     328              :          CALL add_exx_to_rhs(rhs=P_mu_nu, &
     329              :                              qs_env=qs_env, &
     330              :                              ext_hfx_section=hfx_section, &
     331              :                              x_data=mp2_env%ri_rpa%x_data, &
     332              :                              recalc_integrals=.FALSE., &
     333              :                              do_admm=mp2_env%ri_rpa%do_admm, &
     334              :                              do_exx=do_exx, &
     335            8 :                              reuse_hfx=mp2_env%ri_rpa%reuse_hfx)
     336              : 
     337            8 :          focc = 1.0_dp
     338            8 :          IF (nspins == 1) focc = 2.0_dp
     339              :          !focc = 0.0_dp
     340           16 :          DO ispin = 1, nspins
     341            8 :             CALL dbcsr_add(P_mu_nu(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, -1.0_dp)
     342            8 :             CALL copy_dbcsr_to_fm(matrix=P_mu_nu(ispin)%matrix, fm=fm_G_mu_nu)
     343              :             CALL parallel_gemm("N", "N", nao, homo(ispin), nao, 1.0_dp, &
     344            8 :                                fm_G_mu_nu, mo_coeff_o(ispin), 0.0_dp, fm_back)
     345              :             CALL parallel_gemm("T", "N", homo(ispin), virtual(ispin), nao, focc, &
     346            8 :                                fm_back, mo_coeff_v(ispin), 1.0_dp, L_jb(ispin))
     347              :             CALL parallel_gemm("T", "N", homo(ispin), homo(ispin), nao, -focc, &
     348           16 :                                fm_back, mo_coeff_o(ispin), 1.0_dp, W_mo(ispin))
     349              :          END DO
     350              :       END IF
     351              : 
     352              :       ! Prepare arrays for linres code
     353              :       NULLIFY (linres_control)
     354          272 :       ALLOCATE (linres_control)
     355          272 :       linres_control%do_kernel = .TRUE.
     356              :       linres_control%lr_triplet = .FALSE.
     357              :       linres_control%linres_restart = .FALSE.
     358          272 :       linres_control%max_iter = mp2_env%ri_grad%cphf_max_num_iter
     359          272 :       linres_control%eps = mp2_env%ri_grad%cphf_eps_conv
     360          272 :       linres_control%eps_filter = mp2_env%mp2_gpw%eps_filter
     361          272 :       linres_control%restart_every = 50
     362          272 :       linres_control%preconditioner_type = ot_precond_full_all
     363          272 :       linres_control%energy_gap = 0.02_dp
     364              : 
     365              :       NULLIFY (p_env)
     366         1904 :       ALLOCATE (p_env)
     367          272 :       CALL p_env_create(p_env, qs_env, p1_option=P_mu_nu, orthogonal_orbitals=.TRUE., linres_control=linres_control)
     368          272 :       CALL set_qs_env(qs_env, linres_control=linres_control)
     369          272 :       CALL p_env_psi0_changed(p_env, qs_env)
     370          272 :       p_env%new_preconditioner = .TRUE.
     371          272 :       CALL p_env_check_i_alloc(p_env, qs_env)
     372          272 :       mp2_env%ri_grad%p_env => p_env
     373              : 
     374              :       ! update Lagrangian with the CPHF like update, occ-occ block, first call (recompute hfx integrals if needed)
     375          272 :       transf_type_in = 1
     376              :       ! In alpha-beta case, L_bj_alpha has Coulomb and XC alpha-alpha part
     377              :       ! and (only) Coulomb alpha-beta part and vice versa.
     378              : 
     379              :       ! Complete in closed shell case, alpha-alpha (Coulomb and XC)
     380              :       ! part of L_bj(alpha) for open shell
     381              : 
     382              :       CALL cphf_like_update(qs_env, mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     383              :                             P_mo, fm_G_mu_nu, fm_back, transf_type_in, L_jb, &
     384              :                             recalc_hfx_integrals=(.NOT. do_exx .AND. mp2_env%ri_grad%free_hfx_buffer) &
     385          364 :                             .OR. (do_exx .AND. .NOT. mp2_env%ri_rpa%reuse_hfx))
     386              : 
     387              :       ! at this point Lagrangian is completed ready to solve the Z-vector equations
     388              :       ! P_ia will contain the solution of these equations
     389          902 :       ALLOCATE (P_ia(nspins))
     390          630 :       DO ispin = 1, nspins
     391          358 :          NULLIFY (fm_struct_tmp)
     392              :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     393          358 :                                   nrow_global=homo(ispin), ncol_global=virtual(ispin))
     394          358 :          CALL cp_fm_create(P_ia(ispin), fm_struct_tmp, name="P_ia")
     395          358 :          CALL cp_fm_struct_release(fm_struct_tmp)
     396          630 :          CALL cp_fm_set_all(P_ia(ispin), 0.0_dp)
     397              :       END DO
     398              : 
     399              :       CALL solve_z_vector_eq_low(qs_env, mp2_env, unit_nr, &
     400              :                                  mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     401          272 :                                  L_jb, fm_G_mu_nu, fm_back, P_ia)
     402              : 
     403              :       ! release fm stuff
     404          272 :       CALL cp_fm_release(fm_G_mu_nu)
     405          272 :       CALL cp_fm_release(L_jb)
     406          272 :       CALL cp_fm_release(mo_coeff_o)
     407          272 :       CALL cp_fm_release(mo_coeff_v)
     408              : 
     409          272 :       CALL cp_fm_create(fm_mo_mo, P_mo(1)%matrix_struct)
     410              : 
     411          630 :       DO ispin = 1, nspins
     412              :          ! update the MP2-MO density matrix with the occ-virt block
     413              :          CALL cp_fm_to_fm_submat(msource=P_ia(ispin), mtarget=P_mo(ispin), &
     414              :                                  nrow=homo(ispin), ncol=virtual(ispin), &
     415              :                                  s_firstrow=1, s_firstcol=1, &
     416          358 :                                  t_firstrow=1, t_firstcol=homo(ispin) + 1)
     417              :          ! transpose P_MO matrix (easy way to symmetrize)
     418          358 :          CALL cp_fm_set_all(fm_mo_mo, 0.0_dp)
     419              :          ! P_mo now is ready
     420          630 :          CALL cp_fm_uplo_to_full(matrix=P_mo(ispin), work=fm_mo_mo)
     421              :       END DO
     422          272 :       CALL cp_fm_release(P_ia)
     423          272 :       CALL cp_fm_release(fm_mo_mo)
     424              : 
     425              :       ! do the final update to the energy weighted matrix W_MO
     426          630 :       DO ispin = 1, nspins
     427              :          CALL cp_fm_get_info(matrix=W_mo(ispin), &
     428              :                              nrow_local=nrow_local, &
     429              :                              ncol_local=ncol_local, &
     430              :                              row_indices=row_indices, &
     431          358 :                              col_indices=col_indices)
     432         7402 :          DO jjB = 1, ncol_local
     433         6772 :             j_global = col_indices(jjB)
     434         7130 :             IF (j_global <= homo(ispin)) THEN
     435        13994 :                DO iiB = 1, nrow_local
     436        12692 :                   i_global = row_indices(iiB)
     437              :                   W_mo(ispin)%local_data(iiB, jjB) = W_mo(ispin)%local_data(iiB, jjB) &
     438        12692 :                                                      - P_mo(ispin)%local_data(iiB, jjB)*Eigenval(j_global, ispin)
     439        12692 :                   IF (i_global == j_global .AND. nspins == 1) W_mo(ispin)%local_data(iiB, jjB) = &
     440          359 :                      W_mo(ispin)%local_data(iiB, jjB) - 2.0_dp*Eigenval(j_global, ispin)
     441        12692 :                   IF (i_global == j_global .AND. nspins == 2) W_mo(ispin)%local_data(iiB, jjB) = &
     442         1594 :                      W_mo(ispin)%local_data(iiB, jjB) - Eigenval(j_global, ispin)
     443              :                END DO
     444              :             ELSE
     445        67764 :                DO iiB = 1, nrow_local
     446        62294 :                   i_global = row_indices(iiB)
     447        67764 :                   IF (i_global <= homo(ispin)) THEN
     448              :                      ! virt-occ
     449              :                      W_mo(ispin)%local_data(iiB, jjB) = W_mo(ispin)%local_data(iiB, jjB) &
     450        10193 :                                                         - P_mo(ispin)%local_data(iiB, jjB)*Eigenval(i_global, ispin)
     451              :                   ELSE
     452              :                      ! virt-virt
     453              :                      W_mo(ispin)%local_data(iiB, jjB) = W_mo(ispin)%local_data(iiB, jjB) &
     454        52101 :                                                         - P_mo(ispin)%local_data(iiB, jjB)*Eigenval(j_global, ispin)
     455              :                   END IF
     456              :                END DO
     457              :             END IF
     458              :          END DO
     459              :       END DO
     460              : 
     461              :       ! create the MP2 energy weighted density matrix
     462          272 :       NULLIFY (p_env%w1)
     463          272 :       CALL dbcsr_allocate_matrix_set(p_env%w1, 1)
     464          272 :       ALLOCATE (p_env%w1(1)%matrix)
     465              :       CALL dbcsr_copy(p_env%w1(1)%matrix, matrix_s(1)%matrix, &
     466          272 :                       name="W MATRIX MP2")
     467          272 :       CALL dbcsr_set(p_env%w1(1)%matrix, 0.0_dp)
     468              : 
     469              :       ! backtnsform the collected parts of the energy-weighted density matrix into AO basis
     470          630 :       DO ispin = 1, nspins
     471              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, 1.0_dp, &
     472          358 :                             mo_coeff(ispin), W_mo(ispin), 0.0_dp, fm_back)
     473          630 :          CALL cp_dbcsr_plus_fm_fm_t(p_env%w1(1)%matrix, fm_back, mo_coeff(ispin), nmo, 1.0_dp, .TRUE., 1)
     474              :       END DO
     475          272 :       CALL cp_fm_release(W_mo)
     476              : 
     477          272 :       CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
     478              : 
     479          630 :       DO ispin = 1, nspins
     480          358 :          CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
     481              : 
     482              :          CALL parallel_gemm('N', 'N', nao, nmo, nmo, 1.0_dp, &
     483          358 :                             mo_coeff(ispin), P_mo(ispin), 0.0_dp, fm_back)
     484          358 :          CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff(ispin), nmo, 1.0_dp, .TRUE.)
     485              : 
     486          630 :          CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
     487              :       END DO
     488          272 :       CALL cp_fm_release(P_mo)
     489          272 :       CALL cp_fm_release(fm_back)
     490              : 
     491          272 :       CALL p_env_update_rho(p_env, qs_env)
     492              : 
     493              :       ! create mp2 DBCSR density
     494          272 :       CALL dbcsr_allocate_matrix_set(matrix_p_mp2, nspins)
     495          630 :       DO ispin = 1, nspins
     496          358 :          ALLOCATE (matrix_p_mp2(ispin)%matrix)
     497              :          CALL dbcsr_copy(matrix_p_mp2(ispin)%matrix, p_env%p1(ispin)%matrix, &
     498          630 :                          name="P MATRIX MP2")
     499              :       END DO
     500              : 
     501          272 :       IF (dft_control%do_admm) THEN
     502           40 :          CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux_fit)
     503           40 :          CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
     504              : 
     505              :          ! create mp2 DBCSR density in auxiliary basis
     506           40 :          CALL dbcsr_allocate_matrix_set(matrix_p_mp2_admm, nspins)
     507           98 :          DO ispin = 1, nspins
     508           58 :             ALLOCATE (matrix_p_mp2_admm(ispin)%matrix)
     509              :             CALL dbcsr_copy(matrix_p_mp2_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, &
     510           98 :                             name="P MATRIX MP2 ADMM")
     511              :          END DO
     512              :       END IF
     513              : 
     514          272 :       CALL set_ks_env(ks_env, matrix_p_mp2=matrix_p_mp2, matrix_p_mp2_admm=matrix_p_mp2_admm)
     515              : 
     516              :       ! We will need one more hfx calculation for HF gradient part
     517          272 :       mp2_env%not_last_hfx = .FALSE.
     518          272 :       mp2_env%p_screen = restore_p_screen
     519              : 
     520          272 :       CALL timestop(handle)
     521              : 
     522         1360 :    END SUBROUTINE solve_z_vector_eq
     523              : 
     524              : ! **************************************************************************************************
     525              : !> \brief Here we performe the CPHF like update using GPW,
     526              : !>        transf_type_in  defines the type of transformation for the matrix in input
     527              : !>        transf_type_in = 1 -> occ-occ back transformation
     528              : !>        transf_type_in = 2 -> virt-virt back transformation
     529              : !>        transf_type_in = 3 -> occ-virt back transformation including the
     530              : !>                              eigenvalues energy differences for the diagonal elements
     531              : !> \param qs_env ...
     532              : !> \param mo_coeff_o ...
     533              : !> \param mo_coeff_v ...
     534              : !> \param Eigenval ...
     535              : !> \param p_env ...
     536              : !> \param fm_mo ...
     537              : !> \param fm_ao ...
     538              : !> \param fm_back ...
     539              : !> \param transf_type_in ...
     540              : !> \param fm_mo_out ...
     541              : !> \param recalc_hfx_integrals ...
     542              : !> \author Mauro Del Ben, Vladimir Rybkin
     543              : ! **************************************************************************************************
     544         1512 :    SUBROUTINE cphf_like_update(qs_env, mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     545         1512 :                                fm_mo, fm_ao, fm_back, transf_type_in, &
     546         1512 :                                fm_mo_out, recalc_hfx_integrals)
     547              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     548              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff_o, mo_coeff_v
     549              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
     550              :       TYPE(qs_p_env_type)                                :: p_env
     551              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_mo
     552              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_ao
     553              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_back
     554              :       INTEGER, INTENT(IN)                                :: transf_type_in
     555              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: fm_mo_out
     556              :       LOGICAL, INTENT(IN), OPTIONAL                      :: recalc_hfx_integrals
     557              : 
     558              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'cphf_like_update'
     559              : 
     560              :       INTEGER                                            :: handle, homo, i_global, iiB, ispin, &
     561              :                                                             j_global, jjB, nao, ncol_local, &
     562              :                                                             nrow_local, nspins, virtual
     563         1512 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     564         1512 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho1_ao
     565              : 
     566         1512 :       CALL timeset(routineN, handle)
     567              : 
     568         1512 :       nspins = SIZE(Eigenval, 2)
     569              : 
     570         1512 :       CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
     571              : 
     572              :       ! Determine the first-order density matrices in AO basis
     573         3488 :       DO ispin = 1, nspins
     574         1976 :          CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
     575              : 
     576         1976 :          CALL cp_fm_get_info(mo_coeff_o(ispin), nrow_global=nao, ncol_global=homo)
     577         1976 :          CALL cp_fm_get_info(mo_coeff_v(ispin), nrow_global=nao, ncol_global=virtual)
     578              : 
     579              :          ASSOCIATE (mat_in => fm_mo(ispin))
     580              : 
     581              :             ! perform back transformation
     582         2334 :             SELECT CASE (transf_type_in)
     583              :             CASE (1)
     584              :                ! occ-occ block
     585              :                CALL parallel_gemm('N', 'N', nao, homo, homo, 1.0_dp, &
     586          358 :                                   mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
     587          358 :                CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_o(ispin), homo, 1.0_dp, .TRUE.)
     588              :                ! virt-virt block
     589              :                CALL parallel_gemm('N', 'N', nao, virtual, virtual, 1.0_dp, &
     590              :                                   mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back, &
     591              :                                   b_first_col=homo + 1, &
     592          358 :                                   b_first_row=homo + 1)
     593          358 :                CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_v(ispin), virtual, 1.0_dp, .TRUE.)
     594              : 
     595              :             CASE (3)
     596              :                ! virt-occ blocks
     597              :                CALL parallel_gemm('N', 'N', nao, virtual, homo, 1.0_dp, &
     598         1618 :                                   mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
     599         1618 :                CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_v(ispin), virtual, 1.0_dp, .TRUE.)
     600              :                ! and symmetrize (here again multiply instead of transposing)
     601              :                CALL parallel_gemm('N', 'T', nao, homo, virtual, 1.0_dp, &
     602         1618 :                                   mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back)
     603         3594 :                CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_o(ispin), homo, 1.0_dp, .TRUE.)
     604              : 
     605              :             CASE DEFAULT
     606              :                ! nothing
     607              :             END SELECT
     608              :          END ASSOCIATE
     609              : 
     610         5464 :          CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
     611              :       END DO
     612              : 
     613         1512 :       CALL p_env_update_rho(p_env, qs_env)
     614              : 
     615         1512 :       CALL apply_2nd_order_kernel(qs_env, p_env, recalc_hfx_integrals)
     616              : 
     617         3488 :       DO ispin = 1, nspins
     618         1976 :          CALL cp_fm_get_info(mo_coeff_o(ispin), nrow_global=nao, ncol_global=homo)
     619         1976 :          CALL cp_fm_get_info(mo_coeff_v(ispin), nrow_global=nao, ncol_global=virtual)
     620              : 
     621         1976 :          IF (transf_type_in == 3) THEN
     622              : 
     623              :             ! scale for the orbital energy differences for the diagonal elements
     624              :             CALL cp_fm_get_info(matrix=fm_mo_out(ispin), &
     625              :                                 nrow_local=nrow_local, &
     626              :                                 ncol_local=ncol_local, &
     627              :                                 row_indices=row_indices, &
     628         1618 :                                 col_indices=col_indices)
     629        28110 :             DO jjB = 1, ncol_local
     630        26492 :                j_global = col_indices(jjB)
     631        79179 :                DO iiB = 1, nrow_local
     632        51069 :                   i_global = row_indices(iiB)
     633              :                   fm_mo_out(ispin)%local_data(iiB, jjB) = fm_mo(ispin)%local_data(iiB, jjB)* &
     634        77561 :                                                           (Eigenval(j_global + homo, ispin) - Eigenval(i_global, ispin))
     635              :                END DO
     636              :             END DO
     637              :          END IF
     638              : 
     639              :          ! copy back to fm
     640         1976 :          CALL cp_fm_set_all(fm_ao, 0.0_dp)
     641         1976 :          CALL copy_dbcsr_to_fm(matrix=p_env%kpp1(ispin)%matrix, fm=fm_ao)
     642         1976 :          CALL cp_fm_set_all(fm_back, 0.0_dp)
     643         1976 :          CALL cp_fm_uplo_to_full(fm_ao, fm_back)
     644              : 
     645         3488 :          ASSOCIATE (mat_out => fm_mo_out(ispin))
     646              : 
     647              :             ! transform to MO basis, here we always sum the result into the input matrix
     648              : 
     649              :             ! occ-virt block
     650              :             CALL parallel_gemm('T', 'N', homo, nao, nao, 1.0_dp, &
     651         1976 :                                mo_coeff_o(ispin), fm_ao, 0.0_dp, fm_back)
     652              :             CALL parallel_gemm('N', 'N', homo, virtual, nao, 1.0_dp, &
     653         1976 :                                fm_back, mo_coeff_v(ispin), 1.0_dp, mat_out)
     654              :          END ASSOCIATE
     655              :       END DO
     656              : 
     657         1512 :       CALL timestop(handle)
     658              : 
     659         1512 :    END SUBROUTINE cphf_like_update
     660              : 
     661              : ! **************************************************************************************************
     662              : !> \brief Low level subroutine for the iterative solution of a large
     663              : !>        system of linear equation
     664              : !> \param qs_env ...
     665              : !> \param mp2_env ...
     666              : !> \param unit_nr ...
     667              : !> \param mo_coeff_o ...
     668              : !> \param mo_coeff_v ...
     669              : !> \param Eigenval ...
     670              : !> \param p_env ...
     671              : !> \param L_jb ...
     672              : !> \param fm_G_mu_nu ...
     673              : !> \param fm_back ...
     674              : !> \param P_ia ...
     675              : !> \author Mauro Del Ben, Vladimir Rybkin
     676              : ! **************************************************************************************************
     677          272 :    SUBROUTINE solve_z_vector_eq_low(qs_env, mp2_env, unit_nr, &
     678          544 :                                     mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     679          272 :                                     L_jb, fm_G_mu_nu, fm_back, P_ia)
     680              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     681              :       TYPE(mp2_type), INTENT(IN)                         :: mp2_env
     682              :       INTEGER, INTENT(IN)                                :: unit_nr
     683              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff_o, mo_coeff_v
     684              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
     685              :       TYPE(qs_p_env_type), INTENT(IN), POINTER           :: p_env
     686              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: L_jb
     687              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_G_mu_nu
     688              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_back
     689              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: P_ia
     690              : 
     691              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_eq_low'
     692              : 
     693              :       INTEGER                                            :: handle, i_global, iiB, ispin, j_global, &
     694              :                                                             jjB, ncol_local, nmo, nrow_local, &
     695              :                                                             nspins, virtual
     696          272 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     697          272 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: precond
     698              : 
     699          272 :       CALL timeset(routineN, handle)
     700              : 
     701          272 :       nmo = SIZE(eigenval, 1)
     702          272 :       nspins = SIZE(eigenval, 2)
     703              : 
     704              :       ! Pople method
     705              :       ! change sign to L_jb
     706          630 :       DO ispin = 1, nspins
     707        16320 :          L_jb(ispin)%local_data(:, :) = -L_jb(ispin)%local_data(:, :)
     708              :       END DO
     709              : 
     710              :       ! create fm structure
     711         1174 :       ALLOCATE (precond(nspins))
     712          630 :       DO ispin = 1, nspins
     713              :          ! create preconditioner (for now only orbital energy differences)
     714          358 :          CALL cp_fm_create(precond(ispin), P_ia(ispin)%matrix_struct, name="precond")
     715          358 :          CALL cp_fm_set_all(precond(ispin), 1.0_dp)
     716              :          CALL cp_fm_get_info(matrix=precond(ispin), &
     717              :                              nrow_local=nrow_local, &
     718              :                              ncol_local=ncol_local, &
     719              :                              row_indices=row_indices, &
     720              :                              col_indices=col_indices, &
     721          358 :                              ncol_global=virtual)
     722         6100 :          DO jjB = 1, ncol_local
     723         5470 :             j_global = col_indices(jjB)
     724        16021 :             DO iiB = 1, nrow_local
     725        10193 :                i_global = row_indices(iiB)
     726        15663 :                precond(ispin)%local_data(iiB, jjB) = 1.0_dp/(Eigenval(j_global + nmo - virtual, ispin) - Eigenval(i_global, ispin))
     727              :             END DO
     728              :          END DO
     729              :       END DO
     730              : 
     731          520 :       SELECT CASE (mp2_env%ri_grad%z_solver_method)
     732              :       CASE (z_solver_pople)
     733              :          CALL solve_z_vector_pople(qs_env, mp2_env, unit_nr, &
     734              :                                    mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     735          248 :                                    L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
     736              :       CASE (z_solver_cg, z_solver_richardson, z_solver_sd)
     737              :          CALL solve_z_vector_cg(qs_env, mp2_env, unit_nr, &
     738              :                                 mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     739           24 :                                 L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
     740              :       CASE DEFAULT
     741          272 :          CPABORT("Unknown solver")
     742              :       END SELECT
     743              : 
     744          272 :       CALL cp_fm_release(precond)
     745              : 
     746          272 :       CALL timestop(handle)
     747              : 
     748          544 :    END SUBROUTINE solve_z_vector_eq_low
     749              : 
     750              : ! **************************************************************************************************
     751              : !> \brief ...
     752              : !> \param qs_env ...
     753              : !> \param mp2_env ...
     754              : !> \param unit_nr ...
     755              : !> \param mo_coeff_o ...
     756              : !> \param mo_coeff_v ...
     757              : !> \param Eigenval ...
     758              : !> \param p_env ...
     759              : !> \param L_jb ...
     760              : !> \param fm_G_mu_nu ...
     761              : !> \param fm_back ...
     762              : !> \param P_ia ...
     763              : !> \param precond ...
     764              : ! **************************************************************************************************
     765          248 :    SUBROUTINE solve_z_vector_pople(qs_env, mp2_env, unit_nr, &
     766          496 :                                    mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
     767          248 :                                    L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
     768              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
     769              :       TYPE(mp2_type), INTENT(IN)                         :: mp2_env
     770              :       INTEGER, INTENT(IN)                                :: unit_nr
     771              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff_o, mo_coeff_v
     772              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
     773              :       TYPE(qs_p_env_type), INTENT(IN), POINTER           :: p_env
     774              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: L_jb
     775              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_G_mu_nu
     776              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_back
     777              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: P_ia, precond
     778              : 
     779              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_pople'
     780              : 
     781              :       INTEGER                                            :: cycle_counter, handle, iiB, iiter, &
     782              :                                                             ispin, max_num_iter, nspins, &
     783              :                                                             transf_type_in
     784              :       LOGICAL                                            :: converged
     785              :       REAL(KIND=dp)                                      :: conv, eps_conv, scale_cphf, t1, t2
     786          248 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: proj_bi_xj, temp_vals, x_norm, xi_b
     787          248 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: A_small, b_small, xi_Axi
     788              :       TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
     789          248 :          DIMENSION(:)                                    :: fm_struct_tmp
     790          248 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: b_i, residual
     791          248 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: Ax, xn
     792              : 
     793          248 :       CALL timeset(routineN, handle)
     794              : 
     795          248 :       nspins = SIZE(eigenval, 2)
     796              : 
     797          248 :       eps_conv = mp2_env%ri_grad%cphf_eps_conv
     798              : 
     799          248 :       IF (accurate_dot_product_spin(L_jb, L_jb) >= (eps_conv*eps_conv)) THEN
     800              : 
     801          248 :          max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
     802          248 :          scale_cphf = mp2_env%ri_grad%scale_step_size
     803              : 
     804              :          ! set the transformation type (equal for all methods all updates)
     805          248 :          transf_type_in = 3
     806              : 
     807              :          ! set convergence flag
     808          248 :          converged = .FALSE.
     809              : 
     810         2490 :          ALLOCATE (fm_struct_tmp(nspins), b_i(nspins), residual(nspins))
     811          582 :          DO ispin = 1, nspins
     812          334 :             fm_struct_tmp(ispin)%struct => P_ia(ispin)%matrix_struct
     813              : 
     814          334 :             CALL cp_fm_create(b_i(ispin), fm_struct_tmp(ispin)%struct, name="b_i")
     815          334 :             CALL cp_fm_set_all(b_i(ispin), 0.0_dp)
     816        14656 :             b_i(ispin)%local_data(:, :) = precond(ispin)%local_data(:, :)*L_jb(ispin)%local_data(:, :)
     817              : 
     818              :             ! create the residual vector (r), we check convergence on the norm of
     819              :             ! this vector r=(Ax-b)
     820          334 :             CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
     821          582 :             CALL cp_fm_set_all(residual(ispin), 0.0_dp)
     822              :          END DO
     823              : 
     824          248 :          IF (unit_nr > 0) THEN
     825          124 :             WRITE (unit_nr, *)
     826          124 :             WRITE (unit_nr, '(T3,A)') "MP2_CPHF| Iterative solution of Z-Vector equations (Pople's method)"
     827          124 :             WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
     828          124 :             WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
     829          124 :             WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling of initial guess: ', scale_cphf
     830          124 :             WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
     831          124 :             WRITE (unit_nr, '(T4,A,T15,A,T33,A)') 'Step', 'Time', 'Convergence'
     832          124 :             WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
     833              :          END IF
     834              : 
     835        11882 :          ALLOCATE (xn(nspins, max_num_iter))
     836        11634 :          ALLOCATE (Ax(nspins, max_num_iter))
     837          744 :          ALLOCATE (x_norm(max_num_iter))
     838          496 :          ALLOCATE (xi_b(max_num_iter))
     839          992 :          ALLOCATE (xi_Axi(max_num_iter, 0:max_num_iter))
     840         4946 :          x_norm = 0.0_dp
     841         4946 :          xi_b = 0.0_dp
     842       139110 :          xi_Axi = 0.0_dp
     843              : 
     844          248 :          cycle_counter = 0
     845         1056 :          DO iiter = 1, max_num_iter
     846         1034 :             cycle_counter = cycle_counter + 1
     847              : 
     848         1034 :             t1 = m_walltime()
     849              : 
     850              :             ! create and update x_i (orthogonalization with previous vectors)
     851         2446 :             DO ispin = 1, nspins
     852         1412 :                CALL cp_fm_create(xn(ispin, iiter), fm_struct_tmp(ispin)%struct, name="xi")
     853         2446 :                CALL cp_fm_set_all(xn(ispin, iiter), 0.0_dp)
     854              :             END DO
     855              : 
     856         2854 :             ALLOCATE (proj_bi_xj(iiter - 1))
     857         3028 :             proj_bi_xj = 0.0_dp
     858              :             ! first compute the projection of the actual b_i into all previous x_i
     859              :             ! already scaled with the norm of each x_i
     860         3028 :             DO iiB = 1, iiter - 1
     861         3028 :                proj_bi_xj(iiB) = proj_bi_xj(iiB) + accurate_dot_product_spin(b_i, xn(:, iiB))/x_norm(iiB)
     862              :             END DO
     863              : 
     864              :             ! update actual x_i
     865         2446 :             DO ispin = 1, nspins
     866        67285 :                xn(ispin, iiter)%local_data(:, :) = scale_cphf*b_i(ispin)%local_data(:, :)
     867         5338 :                DO iiB = 1, iiter - 1
     868              :                   xn(ispin, iiter)%local_data(:, :) = xn(ispin, iiter)%local_data(:, :) - &
     869       151476 :                                                       xn(ispin, iiB)%local_data(:, :)*proj_bi_xj(iiB)
     870              :                END DO
     871              :             END DO
     872         1034 :             DEALLOCATE (proj_bi_xj)
     873              : 
     874              :             ! create Ax(iiter) that will store the matrix vector product for this cycle
     875         2446 :             DO ispin = 1, nspins
     876         1412 :                CALL cp_fm_create(Ax(ispin, iiter), fm_struct_tmp(ispin)%struct, name="Ai")
     877         2446 :                CALL cp_fm_set_all(Ax(ispin, iiter), 0.0_dp)
     878              :             END DO
     879              : 
     880              :             CALL cphf_like_update(qs_env, mo_coeff_o, &
     881              :                                   mo_coeff_v, Eigenval, p_env, &
     882              :                                   xn(:, iiter), fm_G_mu_nu, fm_back, transf_type_in, &
     883         1034 :                                   Ax(:, iiter))
     884              : 
     885              :             ! in order to reduce the number of  parallel sums here we
     886              :             ! cluster all necessary scalar products into a single vector
     887              :             ! temp_vals contains:
     888              :             ! 1:iiter -> <Ax_i|x_j>
     889              :             ! iiter+1 -> <x_i|b>
     890              :             ! iiter+2 -> <x_i|x_i>
     891              : 
     892         3102 :             ALLOCATE (temp_vals(iiter + 2))
     893         6130 :             temp_vals = 0.0_dp
     894              :             ! <Ax_i|x_j>
     895         4062 :             DO iiB = 1, iiter
     896         4062 :                temp_vals(iiB) = temp_vals(iiB) + accurate_dot_product_spin(Ax(:, iiter), xn(:, iiB))
     897              :             END DO
     898              :             ! <x_i|b>
     899         1034 :             temp_vals(iiter + 1) = temp_vals(iiter + 1) + accurate_dot_product_spin(xn(:, iiter), L_jb)
     900              :             ! norm
     901         1034 :             temp_vals(iiter + 2) = temp_vals(iiter + 2) + accurate_dot_product_spin(xn(:, iiter), xn(:, iiter))
     902              :             ! update <Ax_i|x_j>,  <x_i|b> and norm <x_i|x_i>
     903         4062 :             xi_Axi(iiter, 1:iiter) = temp_vals(1:iiter)
     904         4062 :             xi_Axi(1:iiter, iiter) = temp_vals(1:iiter)
     905         1034 :             xi_b(iiter) = temp_vals(iiter + 1)
     906         1034 :             x_norm(iiter) = temp_vals(iiter + 2)
     907         1034 :             DEALLOCATE (temp_vals)
     908              : 
     909              :             ! solve reduced system
     910         1034 :             IF (ALLOCATED(A_small)) DEALLOCATE (A_small)
     911         1034 :             IF (ALLOCATED(b_small)) DEALLOCATE (b_small)
     912         4136 :             ALLOCATE (A_small(iiter, iiter))
     913         3102 :             ALLOCATE (b_small(iiter, 1))
     914        16110 :             A_small(1:iiter, 1:iiter) = xi_Axi(1:iiter, 1:iiter)
     915         4062 :             b_small(1:iiter, 1) = xi_b(1:iiter)
     916              : 
     917         1034 :             CALL solve_system(matrix=A_small, mysize=iiter, eigenvectors=b_small)
     918              : 
     919              :             ! check for convergence
     920         2446 :             DO ispin = 1, nspins
     921         1412 :                CALL cp_fm_set_all(residual(ispin), 0.0_dp)
     922         5716 :                DO iiB = 1, iiter
     923              :                   residual(ispin)%local_data(:, :) = &
     924              :                      residual(ispin)%local_data(:, :) + &
     925       218761 :                      b_small(iiB, 1)*Ax(ispin, iiB)%local_data(:, :)
     926              :                END DO
     927              : 
     928              :                residual(ispin)%local_data(:, :) = &
     929              :                   residual(ispin)%local_data(:, :) - &
     930        68319 :                   L_jb(ispin)%local_data(:, :)
     931              :             END DO
     932              : 
     933         1034 :             conv = SQRT(accurate_dot_product_spin(residual, residual))
     934              : 
     935         1034 :             t2 = m_walltime()
     936              : 
     937         1034 :             IF (unit_nr > 0) THEN
     938          517 :                WRITE (unit_nr, '(T3,I5,T13,F6.1,11X,F14.8)') iiter, t2 - t1, conv
     939          517 :                CALL m_flush(unit_nr)
     940              :             END IF
     941              : 
     942         1034 :             IF (conv <= eps_conv) THEN
     943              :                converged = .TRUE.
     944              :                EXIT
     945              :             END IF
     946              : 
     947              :             ! update b_i for the next round
     948         1926 :             DO ispin = 1, nspins
     949              :                b_i(ispin)%local_data(:, :) = b_i(ispin)%local_data(:, :) &
     950              :                                              + precond(ispin)%local_data(:, :) &
     951        54956 :                                              *Ax(ispin, iiter)%local_data(:, :)
     952              :             END DO
     953              : 
     954          830 :             scale_cphf = 1.0_dp
     955              : 
     956              :          END DO
     957              : 
     958          248 :          IF (unit_nr > 0) THEN
     959          124 :             WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
     960          124 :             IF (converged) THEN
     961          113 :                WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycle_counter, ' steps'
     962              :             ELSE
     963           11 :                WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations NOT converged in', cycle_counter, ' steps'
     964              :             END IF
     965              :          END IF
     966              : 
     967              :          ! store solution into P_ia
     968         1282 :          DO iiter = 1, cycle_counter
     969         2694 :             DO ispin = 1, nspins
     970              :                P_ia(ispin)%local_data(:, :) = P_ia(ispin)%local_data(:, :) + &
     971        68319 :                                               b_small(iiter, 1)*xn(ispin, iiter)%local_data(:, :)
     972              :             END DO
     973              :          END DO
     974              : 
     975              :          ! Release arrays
     976          248 :          DEALLOCATE (x_norm)
     977          248 :          DEALLOCATE (xi_b)
     978          248 :          DEALLOCATE (xi_Axi)
     979              : 
     980          248 :          CALL cp_fm_release(b_i)
     981          248 :          CALL cp_fm_release(residual)
     982          248 :          CALL cp_fm_release(Ax)
     983          248 :          CALL cp_fm_release(xn)
     984              : 
     985              :       ELSE
     986            0 :          IF (unit_nr > 0) THEN
     987            0 :             WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
     988            0 :             WRITE (unit_nr, '(T3,A)') 'Residual smaller than EPS_CONV. Skip solution of Z-vector equation.'
     989              :          END IF
     990              :       END IF
     991              : 
     992          248 :       CALL timestop(handle)
     993              : 
     994          248 :    END SUBROUTINE solve_z_vector_pople
     995              : 
     996              : ! **************************************************************************************************
     997              : !> \brief ...
     998              : !> \param qs_env ...
     999              : !> \param mp2_env ...
    1000              : !> \param unit_nr ...
    1001              : !> \param mo_coeff_o ...
    1002              : !> \param mo_coeff_v ...
    1003              : !> \param Eigenval ...
    1004              : !> \param p_env ...
    1005              : !> \param L_jb ...
    1006              : !> \param fm_G_mu_nu ...
    1007              : !> \param fm_back ...
    1008              : !> \param P_ia ...
    1009              : !> \param precond ...
    1010              : ! **************************************************************************************************
    1011           24 :    SUBROUTINE solve_z_vector_cg(qs_env, mp2_env, unit_nr, &
    1012           48 :                                 mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
    1013           48 :                                 L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
    1014              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
    1015              :       TYPE(mp2_type), INTENT(IN)                         :: mp2_env
    1016              :       INTEGER, INTENT(IN)                                :: unit_nr
    1017              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: mo_coeff_o, mo_coeff_v
    1018              :       REAL(KIND=dp), DIMENSION(:, :), INTENT(IN)         :: Eigenval
    1019              :       TYPE(qs_p_env_type), INTENT(IN), POINTER           :: p_env
    1020              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: L_jb
    1021              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: fm_G_mu_nu
    1022              :       TYPE(cp_fm_type), INTENT(IN)                       :: fm_back
    1023              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: P_ia, precond
    1024              : 
    1025              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'solve_z_vector_cg'
    1026              : 
    1027              :       INTEGER :: cycles_passed, handle, iiter, ispin, max_num_iter, nspins, restart_counter, &
    1028              :          restart_every, transf_type_in, z_solver_method
    1029              :       LOGICAL                                            :: converged, do_restart
    1030              :       REAL(KIND=dp) :: eps_conv, norm_residual, norm_residual_old, &
    1031              :          residual_dot_diff_search_vec_new, residual_dot_diff_search_vec_old, &
    1032              :          residual_dot_search_vec, residual_new_dot_diff_search_vec_old, scale_result, &
    1033              :          scale_search, scale_step_size, search_vec_dot_A_search_vec, t1, t2
    1034              :       TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
    1035           24 :          DIMENSION(:)                                    :: fm_struct_tmp
    1036           24 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: A_dot_search_vector, diff_search_vector, &
    1037           24 :                                                             residual, search_vector
    1038              : 
    1039           24 :       CALL timeset(routineN, handle)
    1040              : 
    1041           24 :       max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
    1042           24 :       eps_conv = mp2_env%ri_grad%cphf_eps_conv
    1043           24 :       z_solver_method = mp2_env%ri_grad%z_solver_method
    1044           24 :       restart_every = mp2_env%ri_grad%cphf_restart
    1045           24 :       scale_step_size = mp2_env%ri_grad%scale_step_size
    1046           24 :       transf_type_in = 3
    1047           24 :       nspins = SIZE(eigenval, 2)
    1048              : 
    1049           24 :       IF (unit_nr > 0) THEN
    1050           12 :          WRITE (unit_nr, *)
    1051            8 :          SELECT CASE (z_solver_method)
    1052              :          CASE (z_solver_cg)
    1053            8 :             IF (mp2_env%ri_grad%polak_ribiere) THEN
    1054            2 :                WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Polak-Ribiere step)'
    1055              :             ELSE
    1056            6 :                WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Fletcher-Reeves step)'
    1057              :             END IF
    1058              :          CASE (z_solver_richardson)
    1059            2 :             WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Richardson method)'
    1060              :          CASE (z_solver_sd)
    1061            2 :             WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Steepest Descent method)'
    1062              :          CASE DEFAULT
    1063           12 :             CPABORT("Unknown solver")
    1064              :          END SELECT
    1065           12 :          WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
    1066           12 :          WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
    1067           12 :          WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Number of steps for restart: ', restart_every
    1068           12 :          WRITE (unit_nr, '(T3, A)') 'MP2_CPHF| Restart after no decrease'
    1069           12 :          WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling factor of each step: ', scale_step_size
    1070           12 :          WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
    1071           12 :          WRITE (unit_nr, '(T4,A,T13,A,T28,A,T43,A)') 'Step', 'Restart', 'Time', 'Convergence'
    1072           12 :          WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
    1073              :       END IF
    1074              : 
    1075            0 :       ALLOCATE (fm_struct_tmp(nspins), residual(nspins), diff_search_vector(nspins), &
    1076          312 :                 search_vector(nspins), A_dot_search_vector(nspins))
    1077           48 :       DO ispin = 1, nspins
    1078           24 :          fm_struct_tmp(ispin)%struct => P_ia(ispin)%matrix_struct
    1079              : 
    1080           24 :          CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
    1081           24 :          CALL cp_fm_set_all(residual(ispin), 0.0_dp)
    1082              : 
    1083           24 :          CALL cp_fm_create(diff_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="difference search vector")
    1084           24 :          CALL cp_fm_set_all(diff_search_vector(ispin), 0.0_dp)
    1085              : 
    1086           24 :          CALL cp_fm_create(search_vector(ispin), fm_struct_tmp(ispin)%struct, name="search vector")
    1087           24 :          CALL cp_fm_set_all(search_vector(ispin), 0.0_dp)
    1088              : 
    1089           24 :          CALL cp_fm_create(A_dot_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="A times search vector")
    1090           48 :          CALL cp_fm_set_all(A_dot_search_vector(ispin), 0.0_dp)
    1091              :       END DO
    1092              : 
    1093           24 :       converged = .FALSE.
    1094           24 :       cycles_passed = max_num_iter
    1095              :       ! By that, we enforce the setup of the matrices
    1096           24 :       do_restart = .TRUE.
    1097              : 
    1098           24 :       t1 = m_walltime()
    1099              : 
    1100          164 :       DO iiter = 1, max_num_iter
    1101              : 
    1102              :          ! During the first iteration, P_ia=0 such that the application of the 2nd order matrix is zero
    1103          160 :          IF (do_restart) THEN
    1104              :             ! We do not consider the first step to be a restart
    1105              :             ! Do not recalculate residual if it is already enforced to save FLOPs
    1106           52 :             IF (.NOT. mp2_env%ri_grad%recalc_residual .OR. (iiter == 1)) THEN
    1107           52 :                IF (iiter > 1) THEN
    1108              :                   CALL cphf_like_update(qs_env, mo_coeff_o, &
    1109              :                                         mo_coeff_v, Eigenval, p_env, &
    1110              :                                         P_ia, fm_G_mu_nu, fm_back, transf_type_in, &
    1111           28 :                                         residual)
    1112              :                ELSE
    1113           24 :                   do_restart = .FALSE.
    1114              : 
    1115           48 :                   DO ispin = 1, nspins
    1116           48 :                      CALL cp_fm_set_all(residual(ispin), 0.0_dp)
    1117              :                   END DO
    1118              :                END IF
    1119              : 
    1120          104 :                DO ispin = 1, nspins
    1121              :                   residual(ispin)%local_data(:, :) = L_jb(ispin)%local_data(:, :) &
    1122         3068 :                                                      - residual(ispin)%local_data(:, :)
    1123              :                END DO
    1124              :             END IF
    1125              : 
    1126          104 :             DO ispin = 1, nspins
    1127              :                diff_search_vector(ispin)%local_data(:, :) = &
    1128         3016 :                   precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
    1129         3068 :                search_vector(ispin)%local_data(:, :) = diff_search_vector(ispin)%local_data(:, :)
    1130              :             END DO
    1131              : 
    1132              :             restart_counter = 1
    1133              :          END IF
    1134              : 
    1135          160 :          norm_residual_old = SQRT(accurate_dot_product_spin(residual, residual))
    1136              : 
    1137          160 :          residual_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
    1138              : 
    1139              :          CALL cphf_like_update(qs_env, mo_coeff_o, &
    1140              :                                mo_coeff_v, Eigenval, p_env, &
    1141              :                                search_vector, fm_G_mu_nu, fm_back, transf_type_in, &
    1142          160 :                                A_dot_search_vector)
    1143              : 
    1144          160 :          IF (z_solver_method /= z_solver_richardson) THEN
    1145          120 :             search_vec_dot_A_search_vec = accurate_dot_product_spin(search_vector, A_dot_search_vector)
    1146              : 
    1147          120 :             IF (z_solver_method == z_solver_cg) THEN
    1148           94 :                scale_result = residual_dot_diff_search_vec_old/search_vec_dot_A_search_vec
    1149              :             ELSE
    1150           26 :                residual_dot_search_vec = accurate_dot_product_spin(residual, search_vector)
    1151           26 :                scale_result = residual_dot_search_vec/search_vec_dot_A_search_vec
    1152              :             END IF
    1153              : 
    1154          120 :             scale_result = scale_result*scale_step_size
    1155              : 
    1156              :          ELSE
    1157              : 
    1158              :             scale_result = scale_step_size
    1159              : 
    1160              :          END IF
    1161              : 
    1162          320 :          DO ispin = 1, nspins
    1163              :             P_ia(ispin)%local_data(:, :) = P_ia(ispin)%local_data(:, :) &
    1164         9440 :                                            + scale_result*search_vector(ispin)%local_data(:, :)
    1165              :          END DO
    1166              : 
    1167          160 :          IF (.NOT. mp2_env%ri_grad%recalc_residual) THEN
    1168              : 
    1169          284 :             DO ispin = 1, nspins
    1170              :                residual(ispin)%local_data(:, :) = residual(ispin)%local_data(:, :) &
    1171         8378 :                                                   - scale_result*A_dot_search_vector(ispin)%local_data(:, :)
    1172              :             END DO
    1173              :          ELSE
    1174              :             CALL cphf_like_update(qs_env, mo_coeff_o, &
    1175              :                                   mo_coeff_v, Eigenval, p_env, &
    1176              :                                   P_ia, fm_G_mu_nu, fm_back, transf_type_in, &
    1177           18 :                                   residual)
    1178              : 
    1179           36 :             DO ispin = 1, nspins
    1180         1062 :                residual(ispin)%local_data(:, :) = L_jb(ispin)%local_data(:, :) - residual(ispin)%local_data(:, :)
    1181              :             END DO
    1182              :          END IF
    1183              : 
    1184          160 :          norm_residual = SQRT(accurate_dot_product_spin(residual, residual))
    1185              : 
    1186          160 :          t2 = m_walltime()
    1187              : 
    1188          160 :          IF (unit_nr > 0) THEN
    1189           80 :             WRITE (unit_nr, '(T3,I4,T16,L1,T26,F6.1,8X,F14.8)') iiter, do_restart, t2 - t1, norm_residual
    1190           80 :             CALL m_flush(unit_nr)
    1191              :          END IF
    1192              : 
    1193          160 :          IF (norm_residual <= eps_conv) THEN
    1194           20 :             converged = .TRUE.
    1195           20 :             cycles_passed = iiter
    1196           20 :             EXIT
    1197              :          END IF
    1198              : 
    1199          140 :          t1 = m_walltime()
    1200              : 
    1201          140 :          IF (z_solver_method == z_solver_richardson) THEN
    1202           80 :             DO ispin = 1, nspins
    1203              :                search_vector(ispin)%local_data(:, :) = &
    1204         2360 :                   scale_step_size*precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
    1205              :             END DO
    1206          100 :          ELSE IF (z_solver_method == z_solver_sd) THEN
    1207           44 :             DO ispin = 1, nspins
    1208              :                search_vector(ispin)%local_data(:, :) = &
    1209         1298 :                   precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
    1210              :             END DO
    1211              :          ELSE
    1212           78 :             IF (mp2_env%ri_grad%polak_ribiere) &
    1213           28 :                residual_new_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
    1214              : 
    1215          156 :             DO ispin = 1, nspins
    1216              :                diff_search_vector(ispin)%local_data(:, :) = &
    1217         4602 :                   precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
    1218              :             END DO
    1219              : 
    1220           78 :             residual_dot_diff_search_vec_new = accurate_dot_product_spin(residual, diff_search_vector)
    1221              : 
    1222           78 :             scale_search = residual_dot_diff_search_vec_new/residual_dot_diff_search_vec_old
    1223           78 :             IF (mp2_env%ri_grad%polak_ribiere) scale_search = &
    1224           28 :                scale_search - residual_new_dot_diff_search_vec_old/residual_dot_diff_search_vec_old
    1225              : 
    1226          156 :             DO ispin = 1, nspins
    1227              :                search_vector(ispin)%local_data(:, :) = scale_search*search_vector(ispin)%local_data(:, :) &
    1228         4602 :                                                        + diff_search_vector(ispin)%local_data(:, :)
    1229              :             END DO
    1230              : 
    1231              :             ! Make new to old
    1232          140 :             residual_dot_diff_search_vec_old = residual_dot_diff_search_vec_new
    1233              :          END IF
    1234              : 
    1235              :          ! Check whether the residual decrease or restart is enforced and ask for restart
    1236          140 :          do_restart = (norm_residual >= norm_residual_old .OR. (MOD(restart_counter, restart_every) == 0))
    1237              : 
    1238          140 :          restart_counter = restart_counter + 1
    1239          144 :          norm_residual_old = norm_residual
    1240              : 
    1241              :       END DO
    1242              : 
    1243           24 :       IF (unit_nr > 0) THEN
    1244           12 :          WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
    1245           12 :          IF (converged) THEN
    1246           10 :             WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycles_passed, ' steps'
    1247              :          ELSE
    1248            2 :             WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations NOT converged in', max_num_iter, ' steps'
    1249              :          END IF
    1250              :       END IF
    1251              : 
    1252           24 :       DEALLOCATE (fm_struct_tmp)
    1253           24 :       CALL cp_fm_release(residual)
    1254           24 :       CALL cp_fm_release(diff_search_vector)
    1255           24 :       CALL cp_fm_release(search_vector)
    1256           24 :       CALL cp_fm_release(A_dot_search_vector)
    1257              : 
    1258           24 :       CALL timestop(handle)
    1259              : 
    1260           48 :    END SUBROUTINE solve_z_vector_cg
    1261              : 
    1262              : ! **************************************************************************************************
    1263              : !> \brief ...
    1264              : !> \param matrix1 ...
    1265              : !> \param matrix2 ...
    1266              : !> \return ...
    1267              : ! **************************************************************************************************
    1268         9104 :    FUNCTION accurate_dot_product_spin(matrix1, matrix2) RESULT(dotproduct)
    1269              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: matrix1, matrix2
    1270              :       REAL(KIND=dp)                                      :: dotproduct
    1271              : 
    1272              :       INTEGER                                            :: ispin
    1273              : 
    1274         9104 :       dotproduct = 0.0_dp
    1275        21602 :       DO ispin = 1, SIZE(matrix1)
    1276        21602 :          dotproduct = dotproduct + accurate_dot_product(matrix1(ispin)%local_data, matrix2(ispin)%local_data)
    1277              :       END DO
    1278         9104 :       CALL matrix1(1)%matrix_struct%para_env%sum(dotproduct)
    1279              : 
    1280         9104 :    END FUNCTION accurate_dot_product_spin
    1281              : 
    1282              : ! **************************************************************************************************
    1283              : !> \brief ...
    1284              : !> \param qs_env ...
    1285              : ! **************************************************************************************************
    1286          272 :    SUBROUTINE update_mp2_forces(qs_env)
    1287              :       TYPE(qs_environment_type), INTENT(IN), POINTER     :: qs_env
    1288              : 
    1289              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'update_mp2_forces'
    1290              : 
    1291              :       INTEGER                                            :: alpha, beta, handle, idir, iounit, &
    1292              :                                                             ispin, nimages, nocc, nspins
    1293              :       INTEGER, DIMENSION(3)                              :: comp
    1294          272 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
    1295              :       LOGICAL                                            :: do_exx, do_hfx, use_virial
    1296              :       REAL(KIND=dp)                                      :: e_dummy, e_hartree, e_xc, ehartree, exc
    1297              :       REAL(KIND=dp), DIMENSION(3)                        :: deb
    1298              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: h_stress, pv_virial
    1299              :       TYPE(admm_type), POINTER                           :: admm_env
    1300          272 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1301              :       TYPE(cell_type), POINTER                           :: cell
    1302              :       TYPE(cp_logger_type), POINTER                      :: logger
    1303              :       TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
    1304          272 :          TARGET                                          :: matrix_ks_aux
    1305          272 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_p_mp2, &
    1306          272 :                                                             matrix_p_mp2_admm, matrix_s, rho1, &
    1307          272 :                                                             rho_ao, rho_ao_aux, scrm
    1308          272 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: rho_ao_kp, scrm_kp
    1309              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1310          272 :       TYPE(hfx_type), DIMENSION(:, :), POINTER           :: x_data
    1311              :       TYPE(linres_control_type), POINTER                 :: linres_control
    1312          272 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
    1313              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1314              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1315          272 :          POINTER                                         :: sab_orb, sac_ae, sac_ppl, sap_ppnl
    1316          272 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1317              :       TYPE(pw_c1d_gs_type)                               :: pot_g, rho_tot_g, temp_pw_g
    1318          272 :       TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:)    :: dvg
    1319          272 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g, rho_mp2_g
    1320              :       TYPE(pw_c1d_gs_type), POINTER                      :: rho_core
    1321              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1322              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
    1323              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1324              :       TYPE(pw_r3d_rs_type)                               :: pot_r, vh_rspace, vhxc_rspace
    1325          272 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_mp2_r, rho_mp2_r_aux, rho_r, &
    1326          272 :                                                             tau_mp2_r, vadmm_rspace, vtau_rspace, &
    1327          272 :                                                             vxc_rspace
    1328              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
    1329              :       TYPE(qs_energy_type), POINTER                      :: energy
    1330          272 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
    1331          272 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1332              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1333              :       TYPE(qs_p_env_type), POINTER                       :: p_env
    1334              :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux
    1335              :       TYPE(section_vals_type), POINTER                   :: hfx_section, hfx_sections, input, &
    1336              :                                                             xc_section
    1337              :       TYPE(task_list_type), POINTER                      :: task_list_aux_fit
    1338              :       TYPE(virial_type), POINTER                         :: virial
    1339              : 
    1340          272 :       CALL timeset(routineN, handle)
    1341              : 
    1342          272 :       NULLIFY (input, pw_env, matrix_s, rho, energy, force, virial, &
    1343          272 :                matrix_p_mp2, matrix_p_mp2_admm, matrix_ks, rho_core)
    1344              :       CALL get_qs_env(qs_env, &
    1345              :                       ks_env=ks_env, &
    1346              :                       dft_control=dft_control, &
    1347              :                       pw_env=pw_env, &
    1348              :                       input=input, &
    1349              :                       mos=mos, &
    1350              :                       para_env=para_env, &
    1351              :                       matrix_s=matrix_s, &
    1352              :                       matrix_ks=matrix_ks, &
    1353              :                       matrix_p_mp2=matrix_p_mp2, &
    1354              :                       matrix_p_mp2_admm=matrix_p_mp2_admm, &
    1355              :                       rho=rho, &
    1356              :                       cell=cell, &
    1357              :                       force=force, &
    1358              :                       virial=virial, &
    1359              :                       sab_orb=sab_orb, &
    1360              :                       energy=energy, &
    1361              :                       rho_core=rho_core, &
    1362          272 :                       x_data=x_data)
    1363              : 
    1364          272 :       logger => cp_get_default_logger()
    1365              :       iounit = cp_print_key_unit_nr(logger, input, "DFT%XC%WF_CORRELATION%PRINT", &
    1366          272 :                                     extension=".mp2Log")
    1367              : 
    1368          272 :       do_exx = .FALSE.
    1369          272 :       IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
    1370           24 :          hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
    1371           24 :          CALL section_vals_get(hfx_section, explicit=do_exx)
    1372              :       END IF
    1373              : 
    1374          272 :       nimages = dft_control%nimages
    1375          272 :       CPASSERT(nimages == 1)
    1376          272 :       NULLIFY (cell_to_index)
    1377              : 
    1378          272 :       p_env => qs_env%mp2_env%ri_grad%p_env
    1379              : 
    1380          272 :       CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp, rho_r=rho_r, rho_g=rho_g)
    1381          272 :       nspins = SIZE(rho_ao)
    1382              : 
    1383              :       ! check if we have to calculate the virial
    1384          272 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
    1385          272 :       IF (use_virial) virial%pv_calculate = .TRUE.
    1386              : 
    1387          272 :       CALL zero_qs_force(force)
    1388          272 :       IF (use_virial) CALL zero_virial(virial, .FALSE.)
    1389              : 
    1390          630 :       DO ispin = 1, nspins
    1391          630 :          CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
    1392              :       END DO
    1393              : 
    1394          272 :       IF (nspins == 2) THEN
    1395           86 :          CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, 1.0_dp)
    1396              :       END IF
    1397              : 
    1398              :       ! Kinetic energy matrix
    1399          272 :       NULLIFY (scrm)
    1400              :       IF (debug_forces) THEN
    1401         1088 :          deb(1:3) = force(1)%kinetic(1:3, 1)
    1402          272 :          IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1403              :       END IF
    1404              :       CALL build_kinetic_matrix(ks_env, matrix_t=scrm, &
    1405              :                                 matrix_name="KINETIC ENERGY MATRIX", &
    1406              :                                 basis_type="ORB", &
    1407              :                                 sab_nl=sab_orb, calculate_forces=.TRUE., &
    1408          272 :                                 matrix_p=rho_ao(1)%matrix)
    1409              :       IF (debug_forces) THEN
    1410         1088 :          deb(1:3) = force(1)%kinetic(1:3, 1) - deb(1:3)
    1411          272 :          CALL para_env%sum(deb)
    1412          272 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT       ", deb
    1413          272 :          IF (use_virial) THEN
    1414           50 :             e_dummy = third_tr(virial%pv_virial) - e_dummy
    1415           50 :             CALL para_env%sum(e_dummy)
    1416           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: T          ", e_dummy
    1417              :          END IF
    1418              :       END IF
    1419          272 :       CALL dbcsr_deallocate_matrix_set(scrm)
    1420              : 
    1421          272 :       IF (nspins == 2) THEN
    1422           86 :          CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, -1.0_dp)
    1423              :       END IF
    1424              : 
    1425              :       ! Add pseudo potential terms
    1426          272 :       scrm_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
    1427              :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
    1428          272 :                       atomic_kind_set=atomic_kind_set, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
    1429          272 :       IF (ASSOCIATED(sac_ae)) THEN
    1430              :          IF (debug_forces) THEN
    1431            0 :             deb(1:3) = force(1)%all_potential(1:3, 1)
    1432            0 :             IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1433              :          END IF
    1434              :          CALL build_core_ae(scrm_kp, rho_ao_kp, force, &
    1435              :                             virial, .TRUE., use_virial, 1, &
    1436              :                             qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
    1437            0 :                             nimages, cell_to_index, basis_type="ORB")
    1438              :          IF (debug_forces) THEN
    1439            0 :             deb(1:3) = force(1)%all_potential(1:3, 1) - deb(1:3)
    1440            0 :             CALL para_env%sum(deb)
    1441            0 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae     ", deb
    1442            0 :             IF (use_virial) THEN
    1443            0 :                e_dummy = third_tr(virial%pv_virial) - e_dummy
    1444            0 :                CALL para_env%sum(e_dummy)
    1445            0 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hae        ", e_dummy
    1446              :             END IF
    1447              :          END IF
    1448              :       END IF
    1449          272 :       IF (ASSOCIATED(sac_ppl)) THEN
    1450              :          IF (debug_forces) THEN
    1451         1088 :             deb(1:3) = force(1)%gth_ppl(1:3, 1)
    1452          272 :             IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1453              :          END IF
    1454              :          CALL build_core_ppl(scrm_kp, rho_ao_kp, force, &
    1455              :                              virial, .TRUE., use_virial, 1, &
    1456              :                              qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
    1457          272 :                              nimages=1, basis_type="ORB")
    1458              :          IF (debug_forces) THEN
    1459         1088 :             deb(1:3) = force(1)%gth_ppl(1:3, 1) - deb(1:3)
    1460          272 :             CALL para_env%sum(deb)
    1461          272 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl    ", deb
    1462          272 :             IF (use_virial) THEN
    1463           50 :                e_dummy = third_tr(virial%pv_virial) - e_dummy
    1464           50 :                CALL para_env%sum(e_dummy)
    1465           50 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppl       ", e_dummy
    1466              :             END IF
    1467              :          END IF
    1468              :       END IF
    1469          272 :       IF (ASSOCIATED(sap_ppnl)) THEN
    1470              :          IF (debug_forces) THEN
    1471         1048 :             deb(1:3) = force(1)%gth_ppnl(1:3, 1)
    1472          262 :             IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1473              :          END IF
    1474              :          CALL build_core_ppnl(scrm_kp, rho_ao_kp, force, &
    1475              :                               virial, .TRUE., use_virial, 1, &
    1476              :                               qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
    1477          262 :                               nimages=1, basis_type="ORB")
    1478              :          IF (debug_forces) THEN
    1479         1048 :             deb(1:3) = force(1)%gth_ppnl(1:3, 1) - deb(1:3)
    1480          262 :             CALL para_env%sum(deb)
    1481          262 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl   ", deb
    1482          262 :             IF (use_virial) THEN
    1483           50 :                e_dummy = third_tr(virial%pv_virial) - e_dummy
    1484           50 :                CALL para_env%sum(e_dummy)
    1485           50 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppnl      ", e_dummy
    1486              :             END IF
    1487              :          END IF
    1488              :       END IF
    1489              : 
    1490              :       ! Get the different components of the KS potential
    1491          272 :       NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
    1492          272 :       IF (use_virial) THEN
    1493           50 :          h_stress = 0.0_dp
    1494           50 :          CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
    1495              :          ! Update virial
    1496          650 :          virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
    1497          650 :          virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
    1498           50 :          IF (.NOT. do_exx) THEN
    1499          650 :             virial%pv_exc = virial%pv_exc - virial%pv_xc
    1500          650 :             virial%pv_virial = virial%pv_virial - virial%pv_xc
    1501              :          ELSE
    1502            0 :             virial%pv_xc = 0.0_dp
    1503              :          END IF
    1504              :          IF (debug_forces) THEN
    1505           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: xc       ", third_tr(h_stress)
    1506           50 :             CALL para_env%sum(virial%pv_xc(1, 1))
    1507           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Corr xc   ", third_tr(virial%pv_xc)
    1508              :          END IF
    1509              :       ELSE
    1510          222 :          CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc)
    1511              :       END IF
    1512              : 
    1513              :       ! Vhxc
    1514          272 :       CALL get_qs_env(qs_env, pw_env=pw_env)
    1515              :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
    1516          272 :                       poisson_env=poisson_env)
    1517          272 :       CALL auxbas_pw_pool%create_pw(vhxc_rspace)
    1518          872 :       IF (use_virial) h_stress = virial%pv_virial
    1519              :       IF (debug_forces) THEN
    1520         1088 :          deb(1:3) = force(1)%rho_elec(1:3, 1)
    1521          272 :          IF (use_virial) e_dummy = third_tr(h_stress)
    1522              :       END IF
    1523          272 :       IF (do_exx) THEN
    1524           16 :          DO ispin = 1, nspins
    1525            8 :             CALL pw_transfer(vh_rspace, vhxc_rspace)
    1526            8 :             CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
    1527              :             CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
    1528              :                                     hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
    1529            8 :                                     qs_env=qs_env, calculate_forces=.TRUE.)
    1530            8 :             CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
    1531            8 :             CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
    1532              :             CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
    1533              :                                     hmat=matrix_ks(ispin), pmat=matrix_p_mp2(ispin), &
    1534            8 :                                     qs_env=qs_env, calculate_forces=.TRUE.)
    1535           16 :             IF (ASSOCIATED(vtau_rspace)) THEN
    1536              :                CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
    1537              :                                        hmat=matrix_ks(ispin), pmat=matrix_p_mp2(ispin), &
    1538            0 :                                        qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE.)
    1539              :             END IF
    1540              :          END DO
    1541              :       ELSE
    1542          614 :          DO ispin = 1, nspins
    1543          350 :             CALL pw_transfer(vh_rspace, vhxc_rspace)
    1544          350 :             CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
    1545              :             CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
    1546              :                                     hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
    1547          350 :                                     qs_env=qs_env, calculate_forces=.TRUE.)
    1548          614 :             IF (ASSOCIATED(vtau_rspace)) THEN
    1549              :                CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
    1550              :                                        hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
    1551           60 :                                        qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE.)
    1552              :             END IF
    1553              :          END DO
    1554              :       END IF
    1555              :       IF (debug_forces) THEN
    1556         1088 :          deb(1:3) = force(1)%rho_elec(1:3, 1) - deb(1:3)
    1557          272 :          CALL para_env%sum(deb)
    1558          272 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dVhxc    ", deb
    1559          272 :          IF (use_virial) THEN
    1560           50 :             e_dummy = third_tr(virial%pv_virial) - e_dummy
    1561           50 :             CALL para_env%sum(e_dummy)
    1562           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Vhxc      ", e_dummy
    1563              :          END IF
    1564              :       END IF
    1565          272 :       IF (use_virial) THEN
    1566          650 :          h_stress = virial%pv_virial - h_stress
    1567          650 :          virial%pv_ehartree = virial%pv_ehartree + h_stress
    1568              : 
    1569           50 :          CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, tau_r=tau_mp2_r)
    1570           50 :          e_xc = 0.0_dp
    1571          124 :          DO ispin = 1, nspins
    1572              :             ! The potentials have been scaled in ks_ref_potential, but for pw_integral_ab, we need the unscaled potentials
    1573           74 :             CALL pw_scale(vxc_rspace(ispin), 1.0_dp/vxc_rspace(ispin)%pw_grid%dvol)
    1574           74 :             e_xc = e_xc + pw_integral_ab(rho_mp2_r(ispin), vxc_rspace(ispin))
    1575           74 :             IF (ASSOCIATED(vtau_rspace)) CALL pw_scale(vtau_rspace(ispin), 1.0_dp/vtau_rspace(ispin)%pw_grid%dvol)
    1576          124 :             IF (ASSOCIATED(vtau_rspace)) e_xc = e_xc + pw_integral_ab(tau_mp2_r(ispin), vtau_rspace(ispin))
    1577              :          END DO
    1578           50 :          IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG VIRIAL:: vxc*d1   ", e_xc
    1579          200 :          DO alpha = 1, 3
    1580          150 :             virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) - e_xc/REAL(para_env%num_pe, dp)
    1581          200 :             virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) - e_xc/REAL(para_env%num_pe, dp)
    1582              :          END DO
    1583              :       END IF
    1584          630 :       DO ispin = 1, nspins
    1585          358 :          CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
    1586          630 :          IF (ASSOCIATED(vtau_rspace)) THEN
    1587           60 :             CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
    1588              :          END IF
    1589              :       END DO
    1590          272 :       DEALLOCATE (vxc_rspace)
    1591          272 :       CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
    1592          272 :       IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
    1593              : 
    1594          630 :       DO ispin = 1, nspins
    1595          630 :          CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
    1596              :       END DO
    1597              : 
    1598              :       ! pw stuff
    1599          272 :       NULLIFY (poisson_env, auxbas_pw_pool)
    1600              :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
    1601          272 :                       poisson_env=poisson_env)
    1602              : 
    1603              :       ! get some of the grids ready
    1604          272 :       CALL auxbas_pw_pool%create_pw(pot_r)
    1605          272 :       CALL auxbas_pw_pool%create_pw(pot_g)
    1606          272 :       CALL auxbas_pw_pool%create_pw(rho_tot_g)
    1607              : 
    1608          272 :       CALL pw_zero(rho_tot_g)
    1609              : 
    1610          272 :       CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, rho_g=rho_mp2_g)
    1611          630 :       DO ispin = 1, nspins
    1612          630 :          CALL pw_axpy(rho_mp2_g(ispin), rho_tot_g)
    1613              :       END DO
    1614              : 
    1615          272 :       IF (use_virial) THEN
    1616          200 :          ALLOCATE (dvg(3))
    1617          200 :          DO idir = 1, 3
    1618          200 :             CALL auxbas_pw_pool%create_pw(dvg(idir))
    1619              :          END DO
    1620           50 :          CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g, dvhartree=dvg)
    1621              :       ELSE
    1622          222 :          CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
    1623              :       END IF
    1624              : 
    1625          272 :       CALL pw_transfer(pot_g, pot_r)
    1626          272 :       CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
    1627          272 :       CALL pw_axpy(pot_r, vh_rspace)
    1628              : 
    1629              :       ! calculate core forces
    1630          272 :       CALL integrate_v_core_rspace(vh_rspace, qs_env)
    1631              :       IF (debug_forces) THEN
    1632         1088 :          deb(:) = force(1)%rho_core(:, 1)
    1633          272 :          CALL para_env%sum(deb)
    1634          272 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: core      ", deb
    1635          272 :          IF (use_virial) THEN
    1636           50 :             e_dummy = third_tr(virial%pv_virial) - e_dummy
    1637           50 :             CALL para_env%sum(e_dummy)
    1638           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: core      ", e_dummy
    1639              :          END IF
    1640              :       END IF
    1641          272 :       CALL auxbas_pw_pool%give_back_pw(vh_rspace)
    1642              : 
    1643          272 :       IF (use_virial) THEN
    1644              :          ! update virial if necessary with the volume term
    1645              :          ! first create pw auxiliary stuff
    1646           50 :          CALL auxbas_pw_pool%create_pw(temp_pw_g)
    1647              : 
    1648              :          ! make a copy of the MP2 density in G space
    1649           50 :          CALL pw_copy(rho_tot_g, temp_pw_g)
    1650              : 
    1651              :          ! calculate total SCF density and potential
    1652           50 :          CALL pw_copy(rho_g(1), rho_tot_g)
    1653           50 :          IF (nspins == 2) CALL pw_axpy(rho_g(2), rho_tot_g)
    1654           50 :          CALL pw_axpy(rho_core, rho_tot_g)
    1655           50 :          CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
    1656              : 
    1657              :          ! finally update virial with the volume contribution
    1658           50 :          e_hartree = pw_integral_ab(temp_pw_g, pot_g)
    1659           50 :          IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: vh1*d0   ", e_hartree
    1660              : 
    1661           50 :          h_stress = 0.0_dp
    1662          200 :          DO alpha = 1, 3
    1663          150 :             comp = 0
    1664          150 :             comp(alpha) = 1
    1665          150 :             CALL pw_copy(pot_g, rho_tot_g)
    1666          150 :             CALL pw_derive(rho_tot_g, comp)
    1667          150 :             h_stress(alpha, alpha) = -e_hartree
    1668          500 :             DO beta = alpha, 3
    1669              :                h_stress(alpha, beta) = h_stress(alpha, beta) &
    1670          300 :                                        - 2.0_dp*pw_integral_ab(rho_tot_g, dvg(beta))/fourpi
    1671          450 :                h_stress(beta, alpha) = h_stress(alpha, beta)
    1672              :             END DO
    1673              :          END DO
    1674           50 :          IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hartree  ", third_tr(h_stress)
    1675              : 
    1676              :          ! free stuff
    1677           50 :          CALL auxbas_pw_pool%give_back_pw(temp_pw_g)
    1678          200 :          DO idir = 1, 3
    1679          200 :             CALL auxbas_pw_pool%give_back_pw(dvg(idir))
    1680              :          END DO
    1681           50 :          DEALLOCATE (dvg)
    1682              : 
    1683          650 :          virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
    1684          650 :          virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
    1685              : 
    1686              :       END IF
    1687              : 
    1688          630 :       DO ispin = 1, nspins
    1689          358 :          CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
    1690          630 :          IF (dft_control%do_admm) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
    1691              :       END DO
    1692              : 
    1693          272 :       CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
    1694              : 
    1695          272 :       IF (dft_control%do_admm) THEN
    1696           40 :          CALL get_qs_env(qs_env, admm_env=admm_env)
    1697           40 :          xc_section => admm_env%xc_section_primary
    1698              :       ELSE
    1699          232 :          xc_section => section_vals_get_subs_vals(input, "DFT%XC")
    1700              :       END IF
    1701              : 
    1702          272 :       IF (use_virial) THEN
    1703           50 :          h_stress = 0.0_dp
    1704          650 :          pv_virial = virial%pv_virial
    1705              :       END IF
    1706              :       IF (debug_forces) THEN
    1707         1088 :          deb = force(1)%rho_elec(1:3, 1)
    1708          272 :          IF (use_virial) e_dummy = third_tr(pv_virial)
    1709              :       END IF
    1710          272 :       CALL apply_2nd_order_kernel(qs_env, p_env, .FALSE., .TRUE., use_virial, h_stress)
    1711          272 :       IF (use_virial) THEN
    1712          650 :          virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_virial)
    1713              :          IF (debug_forces) THEN
    1714          650 :             e_dummy = third_tr(virial%pv_virial - pv_virial)
    1715           50 :             CALL para_env%sum(e_dummy)
    1716           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kh       ", e_dummy
    1717              :          END IF
    1718          650 :          virial%pv_exc = virial%pv_exc + h_stress
    1719          650 :          virial%pv_virial = virial%pv_virial + h_stress
    1720              :          IF (debug_forces) THEN
    1721           50 :             e_dummy = third_tr(h_stress)
    1722           50 :             CALL para_env%sum(e_dummy)
    1723           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kxc       ", e_dummy
    1724              :          END IF
    1725              :       END IF
    1726              :       IF (debug_forces) THEN
    1727         1088 :          deb(:) = force(1)%rho_elec(:, 1) - deb
    1728          272 :          CALL para_env%sum(deb)
    1729          272 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P0*Khxc    ", deb
    1730          272 :          IF (use_virial) THEN
    1731           50 :             e_dummy = third_tr(virial%pv_virial) - e_dummy
    1732           50 :             CALL para_env%sum(e_dummy)
    1733           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Khxc      ", e_dummy
    1734              :          END IF
    1735              :       END IF
    1736              : 
    1737              :       ! hfx section
    1738          272 :       NULLIFY (hfx_sections)
    1739          272 :       hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
    1740          272 :       CALL section_vals_get(hfx_sections, explicit=do_hfx)
    1741          272 :       IF (do_hfx) THEN
    1742          192 :          IF (do_exx) THEN
    1743            0 :             IF (dft_control%do_admm) THEN
    1744            0 :                CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
    1745            0 :                CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
    1746            0 :                rho1 => p_env%p1_admm
    1747              :             ELSE
    1748            0 :                rho1 => p_env%p1
    1749              :             END IF
    1750              :          ELSE
    1751          192 :             IF (dft_control%do_admm) THEN
    1752           36 :                CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
    1753           36 :                CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
    1754           90 :                DO ispin = 1, nspins
    1755           90 :                   CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
    1756              :                END DO
    1757           36 :                rho1 => p_env%p1_admm
    1758              :             ELSE
    1759          344 :                DO ispin = 1, nspins
    1760          344 :                   CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
    1761              :                END DO
    1762          156 :                rho1 => p_env%p1
    1763              :             END IF
    1764              :          END IF
    1765              : 
    1766          192 :          IF (x_data(1, 1)%do_hfx_ri) THEN
    1767              : 
    1768              :             CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
    1769              :                                       x_data(1, 1)%general_parameter%fraction, &
    1770              :                                       rho_ao=rho_ao_kp, rho_ao_resp=rho1, &
    1771            6 :                                       use_virial=use_virial, resp_only=do_exx)
    1772              : 
    1773              :          ELSE
    1774              :             CALL derivatives_four_center(qs_env, rho_ao_kp, rho1, hfx_sections, para_env, &
    1775          186 :                                          1, use_virial, resp_only=do_exx)
    1776              :          END IF
    1777              : 
    1778          192 :          IF (use_virial) THEN
    1779          338 :             virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
    1780          338 :             virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
    1781              :          END IF
    1782              :          IF (debug_forces) THEN
    1783          768 :             deb(1:3) = force(1)%fock_4c(1:3, 1) - deb(1:3)
    1784          192 :             CALL para_env%sum(deb)
    1785          192 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx  ", deb
    1786          192 :             IF (use_virial) THEN
    1787           26 :                e_dummy = third_tr(virial%pv_fock_4c)
    1788           26 :                CALL para_env%sum(e_dummy)
    1789           26 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: hfx    ", e_dummy
    1790              :             END IF
    1791              :          END IF
    1792              : 
    1793          192 :          IF (.NOT. do_exx) THEN
    1794          192 :          IF (dft_control%do_admm) THEN
    1795           36 :             CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
    1796           90 :             DO ispin = 1, nspins
    1797           90 :                CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
    1798              :             END DO
    1799              :          ELSE
    1800          344 :             DO ispin = 1, nspins
    1801          344 :                CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
    1802              :             END DO
    1803              :          END IF
    1804              :          END IF
    1805              : 
    1806          192 :          IF (dft_control%do_admm) THEN
    1807              :             IF (debug_forces) THEN
    1808          144 :                deb = force(1)%overlap_admm(1:3, 1)
    1809           36 :                IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1810              :             END IF
    1811              :             ! The 2nd order kernel contains a factor of two in apply_xc_admm_ao which we don't need for the projection derivatives
    1812           36 :             IF (nspins == 1) CALL dbcsr_scale(p_env%kpp1_admm(1)%matrix, 0.5_dp)
    1813           36 :             CALL admm_projection_derivative(qs_env, p_env%kpp1_admm, rho_ao)
    1814              :             IF (debug_forces) THEN
    1815          144 :                deb(:) = force(1)%overlap_admm(:, 1) - deb
    1816           36 :                CALL para_env%sum(deb)
    1817           36 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*KADMM*dS'", deb
    1818           36 :                IF (use_virial) THEN
    1819           12 :                   e_dummy = third_tr(virial%pv_virial) - e_dummy
    1820           12 :                   CALL para_env%sum(e_dummy)
    1821           12 :                   IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: KADMM*S'  ", e_dummy
    1822              :                END IF
    1823              :             END IF
    1824              : 
    1825          162 :             ALLOCATE (matrix_ks_aux(nspins))
    1826           90 :             DO ispin = 1, nspins
    1827           54 :                NULLIFY (matrix_ks_aux(ispin)%matrix)
    1828           54 :                ALLOCATE (matrix_ks_aux(ispin)%matrix)
    1829           54 :                CALL dbcsr_copy(matrix_ks_aux(ispin)%matrix, p_env%kpp1_admm(ispin)%matrix)
    1830           90 :                CALL dbcsr_set(matrix_ks_aux(ispin)%matrix, 0.0_dp)
    1831              :             END DO
    1832              : 
    1833              :             ! Calculate kernel
    1834           36 :             CALL tddft_hfx_matrix(matrix_ks_aux, rho_ao_aux, qs_env, .FALSE.)
    1835              : 
    1836           36 :             IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
    1837           24 :                CALL get_qs_env(qs_env, ks_env=ks_env)
    1838           24 :                CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list_aux_fit)
    1839              : 
    1840           60 :                DO ispin = 1, nspins
    1841           60 :                   CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
    1842              :                END DO
    1843              : 
    1844           24 :                IF (use_virial) THEN
    1845            8 :                   CALL qs_rho_get(p_env%rho1_admm, rho_r=rho_mp2_r_aux)
    1846            8 :                   e_xc = 0.0_dp
    1847           20 :                   DO ispin = 1, nspins
    1848           20 :                      e_xc = e_xc + pw_integral_ab(rho_mp2_r_aux(ispin), vadmm_rspace(ispin))
    1849              :                   END DO
    1850              : 
    1851            8 :                   e_xc = -e_xc/vadmm_rspace(1)%pw_grid%dvol/REAL(para_env%num_pe, dp)
    1852              : 
    1853              :                   ! Update the virial
    1854           32 :                   DO alpha = 1, 3
    1855           24 :                      virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) + e_xc
    1856           32 :                      virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) + e_xc
    1857              :                   END DO
    1858              :                   IF (debug_forces) THEN
    1859            8 :                      IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: P1*VADMM  ", e_xc
    1860              :                   END IF
    1861              :                END IF
    1862              : 
    1863          120 :                IF (use_virial) h_stress = virial%pv_virial
    1864              :                IF (debug_forces) THEN
    1865           96 :                   deb = force(1)%rho_elec(1:3, 1)
    1866           24 :                   IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1867              :                END IF
    1868           60 :                DO ispin = 1, nspins
    1869           36 :                   IF (do_exx) THEN
    1870              :                      CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=matrix_ks_aux(ispin), qs_env=qs_env, &
    1871              :                                              calculate_forces=.TRUE., basis_type="AUX_FIT", &
    1872            0 :                                              task_list_external=task_list_aux_fit, pmat=matrix_p_mp2_admm(ispin))
    1873              :                   ELSE
    1874              :                      CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=matrix_ks_aux(ispin), qs_env=qs_env, &
    1875              :                                              calculate_forces=.TRUE., basis_type="AUX_FIT", &
    1876           36 :                                              task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
    1877              :                   END IF
    1878           60 :                   CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
    1879              :                END DO
    1880          120 :                IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - h_stress)
    1881           24 :                DEALLOCATE (vadmm_rspace)
    1882              :                IF (debug_forces) THEN
    1883           96 :                   deb(:) = force(1)%rho_elec(:, 1) - deb
    1884           24 :                   CALL para_env%sum(deb)
    1885           24 :                   IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM' ", deb
    1886           24 :                   IF (use_virial) THEN
    1887            8 :                      e_dummy = third_tr(virial%pv_virial) - e_dummy
    1888            8 :                      CALL para_env%sum(e_dummy)
    1889            8 :                      IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM'   ", e_dummy
    1890              :                   END IF
    1891              :                END IF
    1892              : 
    1893           60 :                DO ispin = 1, nspins
    1894           60 :                   CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
    1895              :                END DO
    1896              : 
    1897              :             END IF
    1898              : 
    1899           90 :             DO ispin = 1, nspins
    1900           90 :                CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
    1901              :             END DO
    1902              : 
    1903              :             IF (debug_forces) THEN
    1904          144 :                deb = force(1)%overlap_admm(1:3, 1)
    1905           36 :                IF (use_virial) e_dummy = third_tr(virial%pv_virial)
    1906              :             END IF
    1907              :             ! Add the second half of the projector deriatives contracting the first order density matrix with the fockian in the auxiliary basis
    1908           36 :             IF (do_exx) THEN
    1909            0 :                CALL admm_projection_derivative(qs_env, matrix_ks_aux, matrix_p_mp2)
    1910              :             ELSE
    1911           36 :                CALL admm_projection_derivative(qs_env, matrix_ks_aux, rho_ao)
    1912              :             END IF
    1913              :             IF (debug_forces) THEN
    1914          144 :                deb(:) = force(1)%overlap_admm(:, 1) - deb
    1915           36 :                CALL para_env%sum(deb)
    1916           36 :                IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM*dS'", deb
    1917           36 :                IF (use_virial) THEN
    1918           12 :                   e_dummy = third_tr(virial%pv_virial) - e_dummy
    1919           12 :                   CALL para_env%sum(e_dummy)
    1920           12 :                   IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM*S'  ", e_dummy
    1921              :                END IF
    1922              :             END IF
    1923              : 
    1924           90 :             DO ispin = 1, nspins
    1925           90 :                CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
    1926              :             END DO
    1927              : 
    1928           90 :             DO ispin = 1, nspins
    1929           54 :                CALL dbcsr_release(matrix_ks_aux(ispin)%matrix)
    1930           90 :                DEALLOCATE (matrix_ks_aux(ispin)%matrix)
    1931              :             END DO
    1932           36 :             DEALLOCATE (matrix_ks_aux)
    1933              :          END IF
    1934              :       END IF
    1935              : 
    1936          272 :       CALL dbcsr_scale(p_env%w1(1)%matrix, -1.0_dp)
    1937              : 
    1938              :       ! Finish matrix_w_mp2 with occ-occ block
    1939          630 :       DO ispin = 1, nspins
    1940          358 :          CALL get_mo_set(mo_set=mos(ispin), homo=nocc, nmo=alpha)
    1941              :          CALL calculate_whz_matrix(mos(ispin)%mo_coeff, p_env%kpp1(ispin)%matrix, &
    1942          630 :                                    p_env%w1(1)%matrix, 1.0_dp, nocc)
    1943              :       END DO
    1944              : 
    1945          272 :       IF (debug_forces .AND. use_virial) e_dummy = third_tr(virial%pv_virial)
    1946              : 
    1947          272 :       NULLIFY (scrm)
    1948              :       CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
    1949              :                                 matrix_name="OVERLAP MATRIX", &
    1950              :                                 basis_type_a="ORB", basis_type_b="ORB", &
    1951              :                                 sab_nl=sab_orb, calculate_forces=.TRUE., &
    1952          272 :                                 matrix_p=p_env%w1(1)%matrix)
    1953          272 :       CALL dbcsr_deallocate_matrix_set(scrm)
    1954              : 
    1955              :       IF (debug_forces) THEN
    1956         1088 :          deb = force(1)%overlap(1:3, 1)
    1957          272 :          CALL para_env%sum(deb)
    1958          272 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: W*dS     ", deb
    1959          272 :          IF (use_virial) THEN
    1960           50 :             e_dummy = third_tr(virial%pv_virial) - e_dummy
    1961           50 :             CALL para_env%sum(e_dummy)
    1962           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: S         ", e_dummy
    1963              :          END IF
    1964              :       END IF
    1965              : 
    1966          272 :       CALL auxbas_pw_pool%give_back_pw(pot_r)
    1967          272 :       CALL auxbas_pw_pool%give_back_pw(pot_g)
    1968          272 :       CALL auxbas_pw_pool%give_back_pw(rho_tot_g)
    1969              : 
    1970              :       ! Release linres stuff
    1971          272 :       CALL p_env_release(p_env)
    1972          272 :       DEALLOCATE (p_env)
    1973          272 :       NULLIFY (qs_env%mp2_env%ri_grad%p_env)
    1974              : 
    1975          272 :       CALL sum_qs_force(force, qs_env%mp2_env%ri_grad%mp2_force)
    1976          272 :       CALL deallocate_qs_force(qs_env%mp2_env%ri_grad%mp2_force)
    1977              : 
    1978          272 :       IF (use_virial) THEN
    1979         1250 :          virial%pv_mp2 = qs_env%mp2_env%ri_grad%mp2_virial
    1980         1250 :          virial%pv_virial = virial%pv_virial + qs_env%mp2_env%ri_grad%mp2_virial
    1981              :          IF (debug_forces) THEN
    1982           50 :             e_dummy = third_tr(virial%pv_mp2)
    1983           50 :             CALL para_env%sum(e_dummy)
    1984           50 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: MP2nonsep  ", e_dummy
    1985              :          END IF
    1986              :       END IF
    1987              :       ! Rewind the change from the beginning
    1988          272 :       IF (use_virial) virial%pv_calculate = .FALSE.
    1989              : 
    1990              :       ! Add the dispersion forces and virials
    1991          272 :       CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
    1992          272 :       CALL calculate_dispersion_pairpot(qs_env, dispersion_env, e_dummy, .TRUE.)
    1993              : 
    1994              :       CALL cp_print_key_finished_output(iounit, logger, input, &
    1995          272 :                                         "DFT%XC%WF_CORRELATION%PRINT")
    1996              : 
    1997          272 :       CALL timestop(handle)
    1998              : 
    1999          544 :    END SUBROUTINE update_mp2_forces
    2000              : 
    2001              : ! **************************************************************************************************
    2002              : !> \brief Calculates the third of the trace of a 3x3 matrix, for debugging purposes
    2003              : !> \param matrix ...
    2004              : !> \return ...
    2005              : ! **************************************************************************************************
    2006          965 :    PURE FUNCTION third_tr(matrix)
    2007              :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: matrix
    2008              :       REAL(KIND=dp)                                      :: third_tr
    2009              : 
    2010          965 :       third_tr = (matrix(1, 1) + matrix(2, 2) + matrix(3, 3))/3.0_dp
    2011              : 
    2012          965 :    END FUNCTION third_tr
    2013              : 
    2014              : END MODULE mp2_cphf
        

Generated by: LCOV version 2.0-1