rt_propagation_forces.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Routines needed for EMD
!> \author Florian Schiffmann (02.09)
! **************************************************************************************************
 
MODULE rt_propagation_forces
   USE admm_types,                      ONLY: admm_type,&
                                              get_admm_env
   USE atomic_kind_types,               ONLY: atomic_kind_type,&
                                              get_atomic_kind_set
   USE cp_control_types,                ONLY: dft_control_type,&
                                              rtp_control_type
   USE cp_dbcsr_api,                    ONLY: &
        dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_get_block_p, &
        dbcsr_get_matrix_type, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
        dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
        dbcsr_p_type, dbcsr_release, dbcsr_type, dbcsr_type_no_symmetry
   USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
   USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
                                              copy_fm_to_dbcsr,&
                                              cp_dbcsr_sm_fm_multiply
   USE cp_fm_types,                     ONLY: cp_fm_create,&
                                              cp_fm_release,&
                                              cp_fm_type
   USE kinds,                           ONLY: dp
   USE mathconstants,                   ONLY: one,&
                                              zero
   USE parallel_gemm_api,               ONLY: parallel_gemm
   USE particle_types,                  ONLY: particle_type
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type
   USE qs_force_types,                  ONLY: add_qs_force,&
                                              qs_force_type
   USE qs_ks_types,                     ONLY: qs_ks_env_type
   USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
   USE qs_overlap,                      ONLY: build_overlap_force
   USE qs_rho_types,                    ONLY: qs_rho_get,&
                                              qs_rho_type
   USE rt_propagation_types,            ONLY: get_rtp,&
                                              rt_prop_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   PUBLIC :: calc_c_mat_force, &
             rt_admm_force
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_forces'
 
CONTAINS
 
! **************************************************************************************************
!> \brief calculates the three additional force contributions needed in EMD
!>        P_imag*C , P_imag*B*S^-1*S_der , P*S^-1*H*S_der
!> \param qs_env ...
!> \par History
!>      02.2014 switched to dbcsr matrices [Samuel Andermatt]
!>      10.2023 merge MO-based and all-atom into one routine [Guillaume Le Breton]
!> \author Florian Schiffmann (02.09)
! **************************************************************************************************
 
   SUBROUTINE calc_c_mat_force(qs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      CHARACTER(LEN=*), PARAMETER                        :: routinen = 'calc_c_mat_force'
      REAL(kind=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
 
      INTEGER                                            :: handle, i, im, ispin, re
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: c_mat, rho_ao, rho_ao_im, rho_new, &
                                                            s_der, sinvb, sinvh, sinvh_imag
      TYPE(dbcsr_type), POINTER                          :: s_inv, tmp
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_rho_type), POINTER                         :: rho
      TYPE(rt_prop_type), POINTER                        :: rtp
      TYPE(rtp_control_type), POINTER                    :: rtp_control
 
      CALL timeset(routinen, handle)
      NULLIFY (rtp, particle_set, atomic_kind_set, dft_control)
 
      CALL get_qs_env(qs_env, &
                      rtp=rtp, &
                      rho=rho, &
                      particle_set=particle_set, &
                      atomic_kind_set=atomic_kind_set, &
                      force=force, &
                      dft_control=dft_control)
 
      rtp_control => dft_control%rtp_control
      CALL get_rtp(rtp=rtp, c_mat=c_mat, s_der=s_der, s_inv=s_inv, &
                   sinvh=sinvh, sinvb=sinvb)
 
      CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
 
      NULLIFY (tmp)
      ALLOCATE (tmp)
      CALL dbcsr_create(tmp, template=sinvb(1)%matrix)
 
      IF (rtp%linear_scaling) THEN
         CALL get_rtp(rtp=rtp, rho_new=rho_new)
      ELSE
         CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
      END IF
 
      ! If SinvH has an imaginary part (the minus sign is already in SinvH_imag)
      IF (rtp%propagate_complex_ks) CALL get_rtp(rtp=rtp, sinvh_imag=sinvh_imag)
 
      DO ispin = 1, SIZE(sinvh)
         re = 2*ispin - 1
         im = 2*ispin
         IF (rtp%linear_scaling) THEN
            CALL dbcsr_multiply("N", "N", one, sinvh(ispin)%matrix, rho_new(re)%matrix, zero, tmp, &
                                filter_eps=rtp%filter_eps)
            IF (rtp%propagate_complex_ks) THEN
               CALL dbcsr_multiply("N", "N", one, sinvh_imag(ispin)%matrix, rho_new(im)%matrix, one, tmp, &
                                   filter_eps=rtp%filter_eps)
            END IF
            CALL dbcsr_multiply("N", "N", one, sinvb(ispin)%matrix, rho_new(im)%matrix, one, tmp, &
                                filter_eps=rtp%filter_eps)
            CALL compute_forces(force, tmp, s_der, rho_new(im)%matrix, c_mat, kind_of, atom_of_kind)
         ELSE
            CALL dbcsr_multiply("N", "N", one, sinvh(ispin)%matrix, rho_ao(ispin)%matrix, zero, tmp)
            IF (rtp%propagate_complex_ks) THEN
               CALL dbcsr_multiply("N", "N", one, sinvh_imag(ispin)%matrix, rho_ao_im(ispin)%matrix, one, tmp)
            END IF
            CALL dbcsr_multiply("N", "N", one, sinvb(ispin)%matrix, rho_ao_im(ispin)%matrix, one, tmp)
            CALL compute_forces(force, tmp, s_der, rho_ao_im(ispin)%matrix, c_mat, kind_of, atom_of_kind)
         END IF
      END DO
 
      ! recall QS forces, at this point have the other sign.
      DO i = 1, SIZE(force)
         force(i)%ehrenfest(:, :) = -force(i)%ehrenfest(:, :)
      END DO
 
      CALL dbcsr_deallocate_matrix(tmp)
 
      CALL timestop(handle)
 
   END SUBROUTINE calc_c_mat_force
 
! **************************************************************************************************
!> \brief ...
!> \param force ...
!> \param tmp ...
!> \param S_der ...
!> \param rho_im ...
!> \param C_mat ...
!> \param kind_of ...
!> \param atom_of_kind ...
! **************************************************************************************************
   SUBROUTINE compute_forces(force, tmp, S_der, rho_im, C_mat, kind_of, atom_of_kind)
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(dbcsr_type), POINTER                          :: tmp
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: s_der
      TYPE(dbcsr_type), POINTER                          :: rho_im
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: c_mat
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of, atom_of_kind
 
      INTEGER                                            :: col_atom, i, ikind, kind_atom, row_atom
      LOGICAL                                            :: found, triangle
      REAL(dp), DIMENSION(:, :), POINTER                 :: block_values, block_values2
      TYPE(dbcsr_iterator_type)                          :: iter
      TYPE(dbcsr_type)                                   :: rho_im_full
 
      ! A block of a (anti)symmetric matrix need not be on the process that holds the same
      ! block of the non-symmetric C_mat, so the stored triangle is taken from a full copy
      triangle = dbcsr_get_matrix_type(rho_im) /= dbcsr_type_no_symmetry
      CALL dbcsr_create(rho_im_full, template=c_mat(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
      CALL dbcsr_desymmetrize(rho_im, rho_im_full)
 
      DO i = 1, 3
         !Calculate the sum over the hadmard product
         !S_der part
 
         CALL dbcsr_iterator_start(iter, tmp)
         DO WHILE (dbcsr_iterator_blocks_left(iter))
            CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
            CALL dbcsr_get_block_p(s_der(i)%matrix, row_atom, col_atom, block_values2, found=found)
            IF (found) THEN
               ikind = kind_of(col_atom)
               kind_atom = atom_of_kind(col_atom)
               !The block_values are in a vector format,
               ! so the dot_product is the sum over all elements of the hamand product, that I need
               force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
                                                      2.0_dp*sum(block_values*block_values2)
            END IF
         END DO
         CALL dbcsr_iterator_stop(iter)
 
         !C_mat part
 
         CALL dbcsr_iterator_start(iter, rho_im_full)
         DO WHILE (dbcsr_iterator_blocks_left(iter))
            CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
            IF (triangle .AND. row_atom > col_atom) cycle
            CALL dbcsr_get_block_p(c_mat(i)%matrix, row_atom, col_atom, block_values2, found=found)
            IF (found) THEN
               ikind = kind_of(col_atom)
               kind_atom = atom_of_kind(col_atom)
               !The block_values are in a vector format, so the dot_product is
               ! the sum over all elements of the hamand product, that I need
               force(ikind)%ehrenfest(i, kind_atom) = force(ikind)%ehrenfest(i, kind_atom) + &
                                                      2.0_dp*sum(block_values*block_values2)
            END IF
         END DO
         CALL dbcsr_iterator_stop(iter)
      END DO
      CALL dbcsr_release(rho_im_full)
 
   END SUBROUTINE compute_forces
 
! **************************************************************************************************
!> \brief Force contribution from the dependence of the ADMM projection on the nuclear positions
!> \param qs_env ...
!> \par History
!>      10.2026 use the density matrix instead of the MO coefficients
! **************************************************************************************************
   SUBROUTINE rt_admm_force(qs_env)
      TYPE(qs_environment_type), POINTER                 :: qs_env
 
      TYPE(admm_type), POINTER                           :: admm_env
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_aux_im, ks_aux_re, matrix_s_aux_fit, &
                                                            matrix_s_aux_fit_vs_orb
      TYPE(rt_prop_type), POINTER                        :: rtp
 
      CALL get_qs_env(qs_env, &
                      admm_env=admm_env, &
                      rtp=rtp)
      CALL get_admm_env(admm_env, matrix_ks_aux_fit=ks_aux_re, &
                        matrix_ks_aux_fit_im=ks_aux_im, &
                        matrix_s_aux_fit=matrix_s_aux_fit, &
                        matrix_s_aux_fit_vs_orb=matrix_s_aux_fit_vs_orb)
 
      ! currently only none option
      CALL rt_admm_forces_none(qs_env, admm_env, rtp, ks_aux_re, ks_aux_im, &
                               matrix_s_aux_fit, matrix_s_aux_fit_vs_orb)
 
   END SUBROUTINE rt_admm_force
 
! **************************************************************************************************
!> \brief ADMM projection force without purification.
!>        With P~ = A*P*A^T and A = S~^-1*Q, the weighted matrices are
!>        W_Q = -2*S~^-1*Re(K~*A*P) and W_S = -W_Q*A^T.
!>        The complex density matrix P is used instead of the MO coefficients, so that
!>        P is consistent with the auxiliary KS matrix K~ and includes the occupations.
!> \param qs_env ...
!> \param admm_env ...
!> \param rtp ...
!> \param KS_aux_re ...
!> \param KS_aux_im ...
!> \param matrix_s_aux_fit ...
!> \param matrix_s_aux_fit_vs_orb ...
! **************************************************************************************************
   SUBROUTINE rt_admm_forces_none(qs_env, admm_env, rtp, KS_aux_re, KS_aux_im, matrix_s_aux_fit, &
                                  matrix_s_aux_fit_vs_orb)
      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(admm_type), POINTER                           :: admm_env
      TYPE(rt_prop_type), POINTER                        :: rtp
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_aux_re, ks_aux_im, matrix_s_aux_fit, &
                                                            matrix_s_aux_fit_vs_orb
 
      INTEGER                                            :: ispin, nao, natom, naux
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :)        :: admm_force
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cp_fm_type)                                   :: tmp_aux_aux, tmp_aux_kap, tmp_aux_nao, &
                                                            tmp_nao_nao
      TYPE(cp_fm_type), DIMENSION(2)                     :: tmp_aux_ap
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao, rho_ao_im, rho_new
      TYPE(dbcsr_type), POINTER                          :: matrix_w_q, matrix_w_s, rho_im, rho_re
      TYPE(neighbor_list_set_p_type), DIMENSION(:), &
         POINTER                                         :: sab_aux_fit_asymm, sab_aux_fit_vs_orb
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_ks_env_type), POINTER                      :: ks_env
      TYPE(qs_rho_type), POINTER                         :: rho
 
      NULLIFY (sab_aux_fit_asymm, sab_aux_fit_vs_orb, ks_env, rho, rho_ao, rho_ao_im, rho_new)
 
      CALL get_qs_env(qs_env, ks_env=ks_env, rho=rho, natom=natom, &
                      atomic_kind_set=atomic_kind_set, force=force)
      CALL get_admm_env(admm_env, sab_aux_fit_asymm=sab_aux_fit_asymm, &
                        sab_aux_fit_vs_orb=sab_aux_fit_vs_orb)
      IF (rtp%linear_scaling) THEN
         CALL get_rtp(rtp=rtp, rho_new=rho_new)
      ELSE
         CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao, rho_ao_im=rho_ao_im)
      END IF
 
      naux = admm_env%nao_aux_fit; nao = admm_env%nao_orb
 
      ALLOCATE (matrix_w_s)
      CALL dbcsr_create(matrix_w_s, template=matrix_s_aux_fit(1)%matrix, &
                        name='W MATRIX AUX S', matrix_type=dbcsr_type_no_symmetry)
      CALL cp_dbcsr_alloc_block_from_nbl(matrix_w_s, sab_aux_fit_asymm)
 
      ALLOCATE (matrix_w_q)
      CALL dbcsr_copy(matrix_w_q, matrix_s_aux_fit_vs_orb(1)%matrix, &
                      "W MATRIX AUX Q")
 
      CALL cp_fm_create(tmp_nao_nao, admm_env%work_orb_orb%matrix_struct, name="too")
      CALL cp_fm_create(tmp_aux_ap(1), admm_env%work_aux_orb%matrix_struct, name="tao")
      CALL cp_fm_create(tmp_aux_ap(2), admm_env%work_aux_orb%matrix_struct, name="tao")
      CALL cp_fm_create(tmp_aux_kap, admm_env%work_aux_orb%matrix_struct, name="tao")
      CALL cp_fm_create(tmp_aux_nao, admm_env%work_aux_orb%matrix_struct, name="tao")
      CALL cp_fm_create(tmp_aux_aux, admm_env%work_aux_aux%matrix_struct, name="taa")
 
      ALLOCATE (admm_force(3, natom))
 
      DO ispin = 1, SIZE(ks_aux_re)
         IF (rtp%linear_scaling) THEN
            rho_re => rho_new(2*ispin - 1)%matrix
            rho_im => rho_new(2*ispin)%matrix
         ELSE
            rho_re => rho_ao(ispin)%matrix
            rho_im => rho_ao_im(ispin)%matrix
         END IF
 
! First calculate A*P, real part ends on tmp_aux_ap(1), imaginary part on tmp_aux_ap(2)
         CALL copy_dbcsr_to_fm(rho_re, tmp_nao_nao)
         CALL parallel_gemm('N', 'N', naux, nao, nao, 1.0_dp, admm_env%A, tmp_nao_nao, 0.0_dp, tmp_aux_ap(1))
         CALL copy_dbcsr_to_fm(rho_im, tmp_nao_nao)
         CALL parallel_gemm('N', 'N', naux, nao, nao, 1.0_dp, admm_env%A, tmp_nao_nao, 0.0_dp, tmp_aux_ap(2))
 
! Next step Re(K~*A*P) = K~_re*A*P_re - K~_im*A*P_im
         CALL cp_dbcsr_sm_fm_multiply(ks_aux_re(ispin)%matrix, tmp_aux_ap(1), tmp_aux_kap, nao, 1.0_dp, 0.0_dp)
         CALL cp_dbcsr_sm_fm_multiply(ks_aux_im(ispin)%matrix, tmp_aux_ap(2), tmp_aux_kap, nao, -1.0_dp, 1.0_dp)
 
! Wq=-2*S~^-1*Re(K~*A*P), minus sign due to force = -dE/dR
         CALL parallel_gemm('N', 'N', naux, nao, naux, -2.0_dp, admm_env%S_inv, tmp_aux_kap, 0.0_dp, tmp_aux_nao)
 
! Ws=-Wq*A^T
         CALL parallel_gemm('N', 'T', naux, naux, nao, -1.0_dp, tmp_aux_nao, admm_env%A, 0.0_dp, tmp_aux_aux)
 
         ! *** copy to sparse matrix
         CALL copy_fm_to_dbcsr(tmp_aux_nao, matrix_w_q, keep_sparsity=.true.)
 
         ! *** copy to sparse matrix
         CALL copy_fm_to_dbcsr(tmp_aux_aux, matrix_w_s, keep_sparsity=.true.)
 
         admm_force = 0.0_dp
         CALL build_overlap_force(ks_env, admm_force, &
                                  basis_type_a="AUX_FIT", basis_type_b="AUX_FIT", &
                                  sab_nl=sab_aux_fit_asymm, matrix_p=matrix_w_s)
         CALL build_overlap_force(ks_env, admm_force, &
                                  basis_type_a="AUX_FIT", basis_type_b="ORB", &
                                  sab_nl=sab_aux_fit_vs_orb, matrix_p=matrix_w_q)
         ! add forces
         CALL add_qs_force(admm_force, force, "overlap_admm", atomic_kind_set)
 
      END DO
 
      DEALLOCATE (admm_force)
 
      ! *** Deallocated weighted density matrices
      CALL dbcsr_deallocate_matrix(matrix_w_s)
      CALL dbcsr_deallocate_matrix(matrix_w_q)
 
      CALL cp_fm_release(tmp_nao_nao)
      CALL cp_fm_release(tmp_aux_ap(1))
      CALL cp_fm_release(tmp_aux_ap(2))
      CALL cp_fm_release(tmp_aux_kap)
      CALL cp_fm_release(tmp_aux_nao)
      CALL cp_fm_release(tmp_aux_aux)
 
   END SUBROUTINE rt_admm_forces_none
 
END MODULE rt_propagation_forces