LCOV - code coverage report
Current view: top level - src - rtp_admm_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 134 146 91.8 %
Date: 2024-03-29 07:50:05 Functions: 7 7 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Utilities for rtp in combination with admm methods
      10             : !>        adapted routines from admm_method (author Manuel Guidon)
      11             : !>
      12             : !> \par History    Use new "force only" overlap routine [07.2014,JGH]
      13             : !> \author Florian Schiffmann
      14             : ! **************************************************************************************************
      15             : MODULE rtp_admm_methods
      16             :    USE admm_types,                      ONLY: admm_env_create,&
      17             :                                               admm_type,&
      18             :                                               get_admm_env
      19             :    USE cp_control_types,                ONLY: admm_control_type,&
      20             :                                               dft_control_type
      21             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      22             :                                               copy_fm_to_dbcsr,&
      23             :                                               cp_dbcsr_plus_fm_fm_t
      24             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_upper_to_full
      25             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose,&
      26             :                                               cp_fm_cholesky_invert
      27             :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      28             :                                               cp_fm_to_fm,&
      29             :                                               cp_fm_type
      30             :    USE dbcsr_api,                       ONLY: &
      31             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, &
      32             :         dbcsr_get_info, dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
      33             :    USE hfx_admm_utils,                  ONLY: create_admm_xc_section
      34             :    USE input_constants,                 ONLY: do_admm_basis_projection,&
      35             :                                               do_admm_purify_none
      36             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      37             :                                               section_vals_type
      38             :    USE kinds,                           ONLY: default_string_length,&
      39             :                                               dp
      40             :    USE mathconstants,                   ONLY: one,&
      41             :                                               zero
      42             :    USE message_passing,                 ONLY: mp_para_env_type
      43             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      44             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      45             :                                               pw_r3d_rs_type
      46             :    USE qs_collocate_density,            ONLY: calculate_rho_elec
      47             :    USE qs_environment_types,            ONLY: get_qs_env,&
      48             :                                               qs_environment_type,&
      49             :                                               set_qs_env
      50             :    USE qs_gapw_densities,               ONLY: prepare_gapw_den
      51             :    USE qs_kind_types,                   ONLY: get_qs_kind_set,&
      52             :                                               qs_kind_type
      53             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      54             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      55             :                                               mo_set_type
      56             :    USE qs_rho_atom_methods,             ONLY: calculate_rho_atom_coeff
      57             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      58             :                                               qs_rho_set,&
      59             :                                               qs_rho_type
      60             :    USE rt_propagation_types,            ONLY: get_rtp,&
      61             :                                               rt_prop_type
      62             :    USE task_list_types,                 ONLY: task_list_type
      63             : #include "./base/base_uses.f90"
      64             : 
      65             :    IMPLICIT NONE
      66             : 
      67             :    PRIVATE
      68             : 
      69             :    ! *** Public subroutines ***
      70             :    PUBLIC :: rtp_admm_calc_rho_aux, rtp_admm_merge_ks_matrix
      71             : 
      72             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rtp_admm_methods'
      73             : 
      74             : CONTAINS
      75             : 
      76             : ! **************************************************************************************************
      77             : !> \brief  Compute the ADMM density matrix in case of rtp (complex MO's)
      78             : !>
      79             : !> \param qs_env ...
      80             : !> \par History
      81             : ! **************************************************************************************************
      82          76 :    SUBROUTINE rtp_admm_calc_rho_aux(qs_env)
      83             : 
      84             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      85             : 
      86             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_calc_rho_aux'
      87             : 
      88             :       CHARACTER(LEN=default_string_length)               :: basis_type
      89             :       INTEGER                                            :: handle, ispin, nspins
      90             :       LOGICAL                                            :: gapw, s_mstruct_changed
      91          76 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tot_rho_r_aux
      92             :       TYPE(admm_type), POINTER                           :: admm_env
      93          76 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: rtp_coeff_aux_fit
      94          76 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p_aux, matrix_p_aux_im, &
      95          76 :                                                             matrix_s_aux_fit, &
      96          76 :                                                             matrix_s_aux_fit_vs_orb
      97             :       TYPE(dft_control_type), POINTER                    :: dft_control
      98          76 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos, mos_aux_fit
      99             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     100          76 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_g_aux
     101          76 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_r_aux
     102             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     103             :       TYPE(qs_rho_type), POINTER                         :: rho, rho_aux_fit
     104             :       TYPE(rt_prop_type), POINTER                        :: rtp
     105             :       TYPE(task_list_type), POINTER                      :: task_list_aux_fit
     106             : 
     107          76 :       CALL timeset(routineN, handle)
     108          76 :       NULLIFY (admm_env, matrix_p_aux, matrix_p_aux_im, mos, &
     109          76 :                mos_aux_fit, para_env, matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, rho, &
     110          76 :                ks_env, dft_control, tot_rho_r_aux, rho_r_aux, rho_g_aux, task_list_aux_fit)
     111             : 
     112             :       CALL get_qs_env(qs_env, &
     113             :                       admm_env=admm_env, &
     114             :                       ks_env=ks_env, &
     115             :                       dft_control=dft_control, &
     116             :                       para_env=para_env, &
     117             :                       mos=mos, &
     118             :                       rtp=rtp, &
     119             :                       rho=rho, &
     120          76 :                       s_mstruct_changed=s_mstruct_changed)
     121             :       CALL get_admm_env(admm_env, matrix_s_aux_fit=matrix_s_aux_fit, task_list_aux_fit=task_list_aux_fit, &
     122             :                         matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb, mos_aux_fit=mos_aux_fit, &
     123          76 :                         rho_aux_fit=rho_aux_fit)
     124          76 :       gapw = admm_env%do_gapw
     125             : 
     126          76 :       nspins = dft_control%nspins
     127             : 
     128          76 :       CALL get_rtp(rtp=rtp, admm_mos=rtp_coeff_aux_fit)
     129             :       CALL rtp_admm_fit_mo_coeffs(qs_env, admm_env, dft_control%admm_control, para_env, &
     130             :                                   matrix_s_aux_fit, matrix_s_aux_fit_vs_orb, &
     131             :                                   mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, &
     132          76 :                                   s_mstruct_changed)
     133             : 
     134         152 :       DO ispin = 1, nspins
     135             :          CALL qs_rho_get(rho_aux_fit, &
     136             :                          rho_ao=matrix_p_aux, &
     137             :                          rho_ao_im=matrix_p_aux_im, &
     138             :                          rho_r=rho_r_aux, &
     139             :                          rho_g=rho_g_aux, &
     140          76 :                          tot_rho_r=tot_rho_r_aux)
     141             : 
     142             :          CALL rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, &
     143             :                                     matrix_p_aux(ispin)%matrix, &
     144             :                                     matrix_p_aux_im(ispin)%matrix, &
     145          76 :                                     ispin)
     146             : 
     147             :          !IF GAPW, only do the soft basis with PW
     148          76 :          basis_type = "AUX_FIT"
     149          76 :          IF (gapw) THEN
     150          16 :             basis_type = "AUX_FIT_SOFT"
     151          16 :             task_list_aux_fit => admm_env%admm_gapw_env%task_list
     152             :          END IF
     153             : 
     154             :          CALL calculate_rho_elec(matrix_p=matrix_p_aux(ispin)%matrix, &
     155             :                                  rho=rho_r_aux(ispin), &
     156             :                                  rho_gspace=rho_g_aux(ispin), &
     157             :                                  total_rho=tot_rho_r_aux(ispin), &
     158             :                                  ks_env=ks_env, soft_valid=.FALSE., &
     159             :                                  basis_type="AUX_FIT", &
     160          76 :                                  task_list_external=task_list_aux_fit)
     161             : 
     162             :          !IF GAPW, also need to atomic densities
     163         152 :          IF (gapw) THEN
     164             :             CALL calculate_rho_atom_coeff(qs_env, matrix_p_aux, &
     165             :                                           rho_atom_set=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
     166             :                                           qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
     167             :                                           oce=admm_env%admm_gapw_env%oce, sab=admm_env%sab_aux_fit, &
     168          16 :                                           para_env=para_env)
     169             : 
     170             :             CALL prepare_gapw_den(qs_env, local_rho_set=admm_env%admm_gapw_env%local_rho_set, &
     171          16 :                                   do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set)
     172             :          END IF
     173             :       END DO
     174          76 :       CALL set_qs_env(qs_env, admm_env=admm_env)
     175          76 :       CALL qs_rho_set(rho_aux_fit, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
     176             : 
     177          76 :       CALL timestop(handle)
     178             : 
     179          76 :    END SUBROUTINE rtp_admm_calc_rho_aux
     180             : 
     181             : ! **************************************************************************************************
     182             : !> \brief ...
     183             : !> \param admm_env ...
     184             : !> \param rtp_coeff_aux_fit ...
     185             : !> \param density_matrix_aux ...
     186             : !> \param density_matrix_aux_im ...
     187             : !> \param ispin ...
     188             : ! **************************************************************************************************
     189          76 :    SUBROUTINE rtp_admm_calculate_dm(admm_env, rtp_coeff_aux_fit, density_matrix_aux, &
     190             :                                     density_matrix_aux_im, ispin)
     191             :       TYPE(admm_type), POINTER                           :: admm_env
     192             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: rtp_coeff_aux_fit
     193             :       TYPE(dbcsr_type), POINTER                          :: density_matrix_aux, density_matrix_aux_im
     194             :       INTEGER, INTENT(in)                                :: ispin
     195             : 
     196             :       CHARACTER(len=*), PARAMETER :: routineN = 'rtp_admm_calculate_dm'
     197             : 
     198             :       INTEGER                                            :: handle
     199             : 
     200          76 :       CALL timeset(routineN, handle)
     201             : 
     202         152 :       SELECT CASE (admm_env%purification_method)
     203             :       CASE (do_admm_purify_none)
     204             :          CALL calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
     205          76 :                                          rtp_coeff_aux_fit, ispin)
     206             :       CASE DEFAULT
     207          76 :          CPWARN("only purification NONE possible with RTP/EMD at the moment")
     208             :       END SELECT
     209             : 
     210          76 :       CALL timestop(handle)
     211             : 
     212          76 :    END SUBROUTINE rtp_admm_calculate_dm
     213             : 
     214             : ! **************************************************************************************************
     215             : !> \brief ...
     216             : !> \param qs_env ...
     217             : !> \param admm_env ...
     218             : !> \param admm_control ...
     219             : !> \param para_env ...
     220             : !> \param matrix_s_aux_fit ...
     221             : !> \param matrix_s_mixed ...
     222             : !> \param mos ...
     223             : !> \param mos_aux_fit ...
     224             : !> \param rtp ...
     225             : !> \param rtp_coeff_aux_fit ...
     226             : !> \param geometry_did_change ...
     227             : ! **************************************************************************************************
     228         152 :    SUBROUTINE rtp_admm_fit_mo_coeffs(qs_env, admm_env, admm_control, para_env, matrix_s_aux_fit, matrix_s_mixed, &
     229          76 :                                      mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
     230             : 
     231             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     232             :       TYPE(admm_type), POINTER                           :: admm_env
     233             :       TYPE(admm_control_type), POINTER                   :: admm_control
     234             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     235             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s_aux_fit, matrix_s_mixed
     236             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos, mos_aux_fit
     237             :       TYPE(rt_prop_type), POINTER                        :: rtp
     238             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: rtp_coeff_aux_fit
     239             :       LOGICAL, INTENT(IN)                                :: geometry_did_change
     240             : 
     241             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_fit_mo_coeffs'
     242             : 
     243             :       INTEGER                                            :: handle, nao_aux_fit, natoms
     244             :       LOGICAL                                            :: recalc_S
     245          76 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     246             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     247             : 
     248          76 :       CALL timeset(routineN, handle)
     249             : 
     250          76 :       NULLIFY (xc_section, qs_kind_set)
     251             : 
     252          76 :       IF (.NOT. (ASSOCIATED(admm_env))) THEN
     253             :          ! setup admm environment
     254           0 :          CALL get_qs_env(qs_env, input=input, natom=natoms, qs_kind_set=qs_kind_set)
     255           0 :          CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
     256           0 :          CALL admm_env_create(admm_env, admm_control, mos, para_env, natoms, nao_aux_fit)
     257           0 :          xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     258             :          CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
     259           0 :                                      admm_env=admm_env)
     260             : 
     261           0 :          IF (admm_control%method /= do_admm_basis_projection) THEN
     262           0 :             CPWARN("RTP requires BASIS_PROJECTION.")
     263             :          END IF
     264             :       END IF
     265             : 
     266          76 :       recalc_S = geometry_did_change .OR. (rtp%iter == 0 .AND. (rtp%istep == rtp%i_start))
     267             : 
     268         152 :       SELECT CASE (admm_env%purification_method)
     269             :       CASE (do_admm_purify_none)
     270             :          CALL rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
     271          76 :                                      mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, recalc_S)
     272             :       CASE DEFAULT
     273          76 :          CPWARN("Purification method not implemented in combination with RTP")
     274             :       END SELECT
     275             : 
     276          76 :       CALL timestop(handle)
     277             : 
     278          76 :    END SUBROUTINE rtp_admm_fit_mo_coeffs
     279             : ! **************************************************************************************************
     280             : !> \brief Calculates the MO coefficients for the auxiliary fitting basis set
     281             : !>        by minimizing int (psi_i - psi_aux_i)^2 using Lagrangian Multipliers
     282             : !>
     283             : !> \param qs_env ...
     284             : !> \param admm_env The ADMM env
     285             : !> \param para_env The parallel env
     286             : !> \param matrix_s_aux_fit the overlap matrix of the auxiliary fitting basis set
     287             : !> \param matrix_s_mixed the mixed overlap matrix of the auxiliary fitting basis
     288             : !>        set and the orbital basis set
     289             : !> \param mos the MO's of the orbital basis set
     290             : !> \param mos_aux_fit the MO's of the auxiliary fitting basis set
     291             : !> \param rtp ...
     292             : !> \param rtp_coeff_aux_fit ...
     293             : !> \param geometry_did_change flag to indicate if the geomtry changed
     294             : !> \par History
     295             : !>      05.2008 created [Manuel Guidon]
     296             : !> \author Manuel Guidon
     297             : ! **************************************************************************************************
     298         152 :    SUBROUTINE rtp_fit_mo_coeffs_none(qs_env, admm_env, para_env, matrix_s_aux_fit, matrix_s_mixed, &
     299          76 :                                      mos, mos_aux_fit, rtp, rtp_coeff_aux_fit, geometry_did_change)
     300             : 
     301             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     302             :       TYPE(admm_type), POINTER                           :: admm_env
     303             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     304             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s_aux_fit, matrix_s_mixed
     305             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos, mos_aux_fit
     306             :       TYPE(rt_prop_type), POINTER                        :: rtp
     307             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: rtp_coeff_aux_fit
     308             :       LOGICAL, INTENT(IN)                                :: geometry_did_change
     309             : 
     310             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_fit_mo_coeffs_none'
     311             : 
     312             :       INTEGER                                            :: handle, ispin, nao_aux_fit, nao_orb, &
     313             :                                                             natoms, nmo, nmo_mos, nspins
     314          76 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: occ_num, occ_num_aux
     315          76 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     316             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff, mo_coeff_aux_fit
     317             :       TYPE(dft_control_type), POINTER                    :: dft_control
     318          76 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     319             :       TYPE(section_vals_type), POINTER                   :: input, xc_section
     320             : 
     321          76 :       CALL timeset(routineN, handle)
     322             : 
     323          76 :       NULLIFY (dft_control, qs_kind_set)
     324             : 
     325          76 :       IF (.NOT. (ASSOCIATED(admm_env))) THEN
     326           0 :          CALL get_qs_env(qs_env, input=input, natom=natoms, dft_control=dft_control, qs_kind_set=qs_kind_set)
     327           0 :          CALL get_qs_kind_set(qs_kind_set, nsgf=nao_aux_fit, basis_type="AUX_FIT")
     328           0 :          CALL admm_env_create(admm_env, dft_control%admm_control, mos, para_env, natoms, nao_aux_fit)
     329           0 :          xc_section => section_vals_get_subs_vals(input, "DFT%XC")
     330             :          CALL create_admm_xc_section(x_data=qs_env%x_data, xc_section=xc_section, &
     331           0 :                                      admm_env=admm_env)
     332             :       END IF
     333             : 
     334          76 :       nao_aux_fit = admm_env%nao_aux_fit
     335          76 :       nao_orb = admm_env%nao_orb
     336          76 :       nspins = SIZE(mos)
     337             : 
     338             :       ! *** This part only depends on overlap matrices ==> needs only to be calculated if the geometry changed
     339             : 
     340          76 :       IF (geometry_did_change) THEN
     341          20 :          CALL copy_dbcsr_to_fm(matrix_s_aux_fit(1)%matrix, admm_env%S_inv)
     342          20 :          CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux)
     343          20 :          CALL cp_fm_to_fm(admm_env%S_inv, admm_env%S)
     344             : 
     345          20 :          CALL copy_dbcsr_to_fm(matrix_s_mixed(1)%matrix, admm_env%Q)
     346             : 
     347             :          !! Calculate S'_inverse
     348          20 :          CALL cp_fm_cholesky_decompose(admm_env%S_inv)
     349          20 :          CALL cp_fm_cholesky_invert(admm_env%S_inv)
     350             :          !! Symmetrize the guy
     351          20 :          CALL cp_fm_upper_to_full(admm_env%S_inv, admm_env%work_aux_aux)
     352             :          !! Calculate A=S'^(-1)*P
     353             :          CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
     354             :                             1.0_dp, admm_env%S_inv, admm_env%Q, 0.0_dp, &
     355          20 :                             admm_env%A)
     356             :       END IF
     357             : 
     358             :       ! *** Calculate the mo_coeffs for the fitting basis
     359         152 :       DO ispin = 1, nspins
     360          76 :          nmo = admm_env%nmo(ispin)
     361          76 :          IF (nmo == 0) CYCLE
     362             :          !! Lambda = C^(T)*B*C
     363          76 :          CALL get_rtp(rtp=rtp, mos_new=mos_new)
     364          76 :          CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff, occupation_numbers=occ_num, nmo=nmo_mos)
     365             :          CALL get_mo_set(mos_aux_fit(ispin), mo_coeff=mo_coeff_aux_fit, &
     366          76 :                          occupation_numbers=occ_num_aux)
     367             : 
     368             :          CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
     369             :                             1.0_dp, admm_env%A, mos_new(2*ispin - 1), 0.0_dp, &
     370          76 :                             rtp_coeff_aux_fit(2*ispin - 1))
     371             :          CALL parallel_gemm('N', 'N', nao_aux_fit, nmo, nao_orb, &
     372             :                             1.0_dp, admm_env%A, mos_new(2*ispin), 0.0_dp, &
     373          76 :                             rtp_coeff_aux_fit(2*ispin))
     374             : 
     375         228 :          CALL cp_fm_to_fm(rtp_coeff_aux_fit(2*ispin - 1), mo_coeff_aux_fit)
     376             :       END DO
     377             : 
     378          76 :       CALL timestop(handle)
     379             : 
     380          76 :    END SUBROUTINE rtp_fit_mo_coeffs_none
     381             : 
     382             : ! **************************************************************************************************
     383             : !> \brief ...
     384             : !> \param density_matrix_aux ...
     385             : !> \param density_matrix_aux_im ...
     386             : !> \param rtp_coeff_aux_fit ...
     387             : !> \param ispin ...
     388             : ! **************************************************************************************************
     389         228 :    SUBROUTINE calculate_rtp_admm_density(density_matrix_aux, density_matrix_aux_im, &
     390          76 :                                          rtp_coeff_aux_fit, ispin)
     391             : 
     392             :       TYPE(dbcsr_type), POINTER                          :: density_matrix_aux, density_matrix_aux_im
     393             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: rtp_coeff_aux_fit
     394             :       INTEGER, INTENT(in)                                :: ispin
     395             : 
     396             :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_rtp_admm_density'
     397             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     398             : 
     399             :       INTEGER                                            :: handle, im, ncol, re
     400             :       REAL(KIND=dp)                                      :: alpha
     401             : 
     402          76 :       CALL timeset(routineN, handle)
     403             : 
     404          76 :       re = 2*ispin - 1; im = 2*ispin
     405          76 :       alpha = 3*one - REAL(SIZE(rtp_coeff_aux_fit)/2, dp)
     406          76 :       CALL dbcsr_set(density_matrix_aux, zero)
     407          76 :       CALL cp_fm_get_info(rtp_coeff_aux_fit(re), ncol_global=ncol)
     408             :       CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
     409             :                                  matrix_v=rtp_coeff_aux_fit(re), &
     410             :                                  ncol=ncol, &
     411          76 :                                  alpha=alpha)
     412             : 
     413             :       ! It is actually complex conjugate but i*i=-1 therefore it must be added
     414             :       CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux, &
     415             :                                  matrix_v=rtp_coeff_aux_fit(im), &
     416             :                                  ncol=ncol, &
     417          76 :                                  alpha=alpha)
     418             : 
     419             : !   compute the imaginary part of the dm
     420          76 :       CALL dbcsr_set(density_matrix_aux_im, zero)
     421             :       CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix_aux_im, &
     422             :                                  matrix_v=rtp_coeff_aux_fit(im), &
     423             :                                  matrix_g=rtp_coeff_aux_fit(re), &
     424             :                                  ncol=ncol, &
     425             :                                  alpha=2.0_dp*alpha, &
     426          76 :                                  symmetry_mode=-1)
     427             : 
     428          76 :       CALL timestop(handle)
     429             : 
     430          76 :    END SUBROUTINE calculate_rtp_admm_density
     431             : 
     432             : ! **************************************************************************************************
     433             : !> \brief ...
     434             : !> \param qs_env ...
     435             : ! **************************************************************************************************
     436          76 :    SUBROUTINE rtp_admm_merge_ks_matrix(qs_env)
     437             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     438             : 
     439             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'rtp_admm_merge_ks_matrix'
     440             : 
     441             :       INTEGER                                            :: handle, ispin
     442             :       TYPE(admm_type), POINTER                           :: admm_env
     443          76 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_ks_aux_fit, &
     444          76 :                                                             matrix_ks_aux_fit_im, matrix_ks_im
     445             :       TYPE(dft_control_type), POINTER                    :: dft_control
     446             : 
     447          76 :       NULLIFY (admm_env, dft_control, matrix_ks, matrix_ks_im, matrix_ks_aux_fit, matrix_ks_aux_fit_im)
     448          76 :       CALL timeset(routineN, handle)
     449             : 
     450             :       CALL get_qs_env(qs_env, &
     451             :                       admm_env=admm_env, &
     452             :                       dft_control=dft_control, &
     453             :                       matrix_ks=matrix_ks, &
     454          76 :                       matrix_ks_im=matrix_ks_im)
     455          76 :       CALL get_admm_env(admm_env, matrix_ks_aux_fit=matrix_ks_aux_fit, matrix_ks_aux_fit_im=matrix_ks_aux_fit_im)
     456             : 
     457             :       !note: the GAPW contribution to ks_aux_fit taken care of in qs_ks_methods.F/update_admm_ks_atom
     458             : 
     459         152 :       DO ispin = 1, dft_control%nspins
     460             : 
     461          76 :          SELECT CASE (admm_env%purification_method)
     462             :          CASE (do_admm_purify_none)
     463             :             CALL rt_merge_ks_matrix_none(ispin, admm_env, &
     464          76 :                                          matrix_ks, matrix_ks_aux_fit)
     465             :             CALL rt_merge_ks_matrix_none(ispin, admm_env, &
     466          76 :                                          matrix_ks_im, matrix_ks_aux_fit_im)
     467             :          CASE DEFAULT
     468          76 :             CPWARN("only purification NONE possible with RTP/EMD at the moment")
     469             :          END SELECT
     470             : 
     471             :       END DO !spin loop
     472          76 :       CALL timestop(handle)
     473             : 
     474          76 :    END SUBROUTINE rtp_admm_merge_ks_matrix
     475             : 
     476             : ! **************************************************************************************************
     477             : !> \brief ...
     478             : !> \param ispin ...
     479             : !> \param admm_env ...
     480             : !> \param matrix_ks ...
     481             : !> \param matrix_ks_aux_fit ...
     482             : ! **************************************************************************************************
     483         152 :    SUBROUTINE rt_merge_ks_matrix_none(ispin, admm_env, &
     484             :                                       matrix_ks, matrix_ks_aux_fit)
     485             :       INTEGER, INTENT(IN)                                :: ispin
     486             :       TYPE(admm_type), POINTER                           :: admm_env
     487             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_ks_aux_fit
     488             : 
     489             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'rt_merge_ks_matrix_none'
     490             : 
     491             :       CHARACTER                                          :: matrix_type_fit
     492             :       INTEGER                                            :: handle, nao_aux_fit, nao_orb, nmo
     493             :       INTEGER, SAVE                                      :: counter = 0
     494             :       TYPE(dbcsr_type)                                   :: matrix_ks_nosym
     495             :       TYPE(dbcsr_type), POINTER                          :: matrix_k_tilde
     496             : 
     497         152 :       CALL timeset(routineN, handle)
     498             : 
     499         152 :       counter = counter + 1
     500         152 :       nao_aux_fit = admm_env%nao_aux_fit
     501         152 :       nao_orb = admm_env%nao_orb
     502         152 :       nmo = admm_env%nmo(ispin)
     503             :       CALL dbcsr_create(matrix_ks_nosym, template=matrix_ks_aux_fit(ispin)%matrix, &
     504         152 :                         matrix_type=dbcsr_type_no_symmetry)
     505         152 :       CALL dbcsr_set(matrix_ks_nosym, 0.0_dp)
     506         152 :       CALL dbcsr_desymmetrize(matrix_ks_aux_fit(ispin)%matrix, matrix_ks_nosym)
     507             : 
     508         152 :       CALL copy_dbcsr_to_fm(matrix_ks_nosym, admm_env%K(ispin))
     509             : 
     510             :       !! K*A
     511             :       CALL parallel_gemm('N', 'N', nao_aux_fit, nao_orb, nao_aux_fit, &
     512             :                          1.0_dp, admm_env%K(ispin), admm_env%A, 0.0_dp, &
     513         152 :                          admm_env%work_aux_orb)
     514             :       !! A^T*K*A
     515             :       CALL parallel_gemm('T', 'N', nao_orb, nao_orb, nao_aux_fit, &
     516             :                          1.0_dp, admm_env%A, admm_env%work_aux_orb, 0.0_dp, &
     517         152 :                          admm_env%work_orb_orb)
     518             : 
     519         152 :       CALL dbcsr_get_info(matrix_ks_aux_fit(ispin)%matrix, matrix_type=matrix_type_fit)
     520             : 
     521             :       NULLIFY (matrix_k_tilde)
     522         152 :       ALLOCATE (matrix_k_tilde)
     523             :       CALL dbcsr_create(matrix_k_tilde, template=matrix_ks(ispin)%matrix, &
     524         152 :                         name='MATRIX K_tilde', matrix_type=matrix_type_fit)
     525             : 
     526         152 :       CALL dbcsr_copy(matrix_k_tilde, matrix_ks(ispin)%matrix)
     527         152 :       CALL dbcsr_set(matrix_k_tilde, 0.0_dp)
     528         152 :       CALL copy_fm_to_dbcsr(admm_env%work_orb_orb, matrix_k_tilde, keep_sparsity=.TRUE.)
     529             : 
     530         152 :       CALL dbcsr_add(matrix_ks(ispin)%matrix, matrix_k_tilde, 1.0_dp, 1.0_dp)
     531             : 
     532         152 :       CALL dbcsr_deallocate_matrix(matrix_k_tilde)
     533         152 :       CALL dbcsr_release(matrix_ks_nosym)
     534             : 
     535         152 :       CALL timestop(handle)
     536             : 
     537         152 :    END SUBROUTINE rt_merge_ks_matrix_none
     538             : 
     539             : END MODULE rtp_admm_methods

Generated by: LCOV version 1.15