LCOV - code coverage report
Current view: top level - src/emd - rt_delta_pulse.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 309 315 98.1 %
Date: 2023-03-30 11:55:16 Functions: 5 5 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines to apply a delta pulse for RTP and EMD
      10             : ! **************************************************************************************************
      11             : 
      12             : MODULE rt_delta_pulse
      13             :    USE bibliography,                    ONLY: Mattiat2019,&
      14             :                                               Mattiat2022,&
      15             :                                               cite_reference
      16             :    USE cell_types,                      ONLY: cell_type
      17             :    USE commutator_rpnl,                 ONLY: build_com_mom_nl,&
      18             :                                               build_com_nl_mag
      19             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      20             :    USE cp_cfm_basic_linalg,             ONLY: cp_cfm_column_scale
      21             :    USE cp_cfm_diag,                     ONLY: cp_cfm_heevd
      22             :    USE cp_cfm_types,                    ONLY: cp_cfm_create,&
      23             :                                               cp_cfm_release,&
      24             :                                               cp_cfm_to_cfm,&
      25             :                                               cp_cfm_type
      26             :    USE cp_control_types,                ONLY: dft_control_type,&
      27             :                                               rtp_control_type
      28             :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      29             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      30             :                                               cp_dbcsr_sm_fm_multiply,&
      31             :                                               dbcsr_allocate_matrix_set,&
      32             :                                               dbcsr_deallocate_matrix_set
      33             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add,&
      34             :                                               cp_fm_triangular_multiply,&
      35             :                                               cp_fm_upper_to_full
      36             :    USE cp_fm_cholesky,                  ONLY: cp_fm_cholesky_decompose,&
      37             :                                               cp_fm_cholesky_invert,&
      38             :                                               cp_fm_cholesky_reduce,&
      39             :                                               cp_fm_cholesky_restore
      40             :    USE cp_fm_diag,                      ONLY: cp_fm_syevd
      41             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      42             :                                               cp_fm_struct_release,&
      43             :                                               cp_fm_struct_type
      44             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      45             :                                               cp_fm_get_info,&
      46             :                                               cp_fm_release,&
      47             :                                               cp_fm_set_all,&
      48             :                                               cp_fm_to_fm,&
      49             :                                               cp_fm_type
      50             :    USE dbcsr_api,                       ONLY: &
      51             :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_deallocate_matrix, dbcsr_get_info, &
      52             :         dbcsr_init_p, dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
      53             :         dbcsr_type_symmetric
      54             :    USE input_section_types,             ONLY: section_get_ival,&
      55             :                                               section_get_lval,&
      56             :                                               section_vals_type,&
      57             :                                               section_vals_val_get
      58             :    USE kinds,                           ONLY: dp
      59             :    USE mathconstants,                   ONLY: one,&
      60             :                                               twopi,&
      61             :                                               zero
      62             :    USE message_passing,                 ONLY: mp_para_env_type
      63             :    USE moments_utils,                   ONLY: get_reference_point
      64             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      65             :    USE particle_types,                  ONLY: particle_type
      66             :    USE qs_dftb_matrices,                ONLY: build_dftb_overlap
      67             :    USE qs_environment_types,            ONLY: get_qs_env,&
      68             :                                               qs_environment_type
      69             :    USE qs_kind_types,                   ONLY: qs_kind_type
      70             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      71             :                                               mo_set_type
      72             :    USE qs_moments,                      ONLY: build_berry_moment_matrix,&
      73             :                                               build_local_magmom_matrix,&
      74             :                                               build_local_moment_matrix
      75             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      76             :    USE rt_propagation_types,            ONLY: get_rtp,&
      77             :                                               rt_prop_create_mos,&
      78             :                                               rt_prop_type
      79             : #include "../base/base_uses.f90"
      80             : 
      81             :    IMPLICIT NONE
      82             : 
      83             :    PRIVATE
      84             : 
      85             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_delta_pulse'
      86             : 
      87             :    PUBLIC :: apply_delta_pulse
      88             : 
      89             : CONTAINS
      90             : 
      91             : ! **************************************************************************************************
      92             : !> \brief Interface to call the delta pulse depending on the type of calculation.
      93             : !> \param qs_env ...
      94             : !> \param rtp ...
      95             : !> \param rtp_control ...
      96             : !> \author Update: Guillaume Le Breton (2023.01)
      97             : ! **************************************************************************************************
      98             : 
      99          54 :    SUBROUTINE apply_delta_pulse(qs_env, rtp, rtp_control)
     100             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     101             :       TYPE(rt_prop_type), POINTER                        :: rtp
     102             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     103             : 
     104             :       LOGICAL                                            :: my_apply_pulse, periodic_cell
     105             :       TYPE(cell_type), POINTER                           :: cell
     106          54 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new, mos_old
     107          54 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     108             :       TYPE(dft_control_type), POINTER                    :: dft_control
     109          54 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     110             : 
     111             :       CALL get_qs_env(qs_env, &
     112             :                       cell=cell, &
     113             :                       dft_control=dft_control, &
     114          54 :                       matrix_s=matrix_s)
     115         174 :       periodic_cell = ANY(cell%perd > 0)
     116          54 :       my_apply_pulse = .TRUE.
     117          54 :       CALL get_qs_env(qs_env, mos=mos)
     118          54 :       IF (rtp%linear_scaling) THEN
     119          40 :          IF (.NOT. ASSOCIATED(mos)) THEN
     120             :             CALL cp_warn(__LOCATION__, "Delta Pulse not implemented for Linear-Scaling based ground "// &
     121             :                          "state calculation. If you want to perform a Linear-Scaling RTP from a "// &
     122             :                          "Linear-Scaling GS calculation you can do the following: (i) LSCF froms "// &
     123             :                          "scratch, (ii) MO-based SCF (for 1 SCF loop for instance) with the LSCF "// &
     124             :                          "result as a restart and (iii) linear scaling RTP + delta kick (for 1 "// &
     125           0 :                          "SCF loop for instance).")
     126             :             my_apply_pulse = .FALSE.
     127             :          ELSE
     128             :             ! create temporary mos_old and mos_new to use delta kick routine designed for MOs-based RTP
     129             :             CALL rt_prop_create_mos(rtp, mos, qs_env%mpools, dft_control, &
     130             :                                     init_mos_old=.TRUE., init_mos_new=.TRUE., &
     131          40 :                                     init_mos_next=.FALSE., init_mos_admn=.FALSE.)
     132             :          END IF
     133             :       END IF
     134             : 
     135             :       IF (my_apply_pulse) THEN
     136          54 :          CALL get_rtp(rtp=rtp, mos_old=mos_old, mos_new=mos_new)
     137          54 :          IF (rtp_control%apply_delta_pulse) THEN
     138          46 :             IF (dft_control%qs_control%dftb) &
     139           0 :                CALL build_dftb_overlap(qs_env, 1, matrix_s)
     140          46 :             IF (rtp_control%periodic) THEN
     141          32 :                CALL apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new)
     142             :             ELSE
     143          14 :                IF (periodic_cell) THEN
     144           4 :                   CPWARN("This application of the delta pulse is not compatible with PBC!")
     145             :                END IF
     146          14 :                CALL apply_delta_pulse_electric(qs_env, mos_old, mos_new)
     147             :             END IF
     148           8 :          ELSE IF (rtp_control%apply_delta_pulse_mag) THEN
     149           8 :             IF (periodic_cell) THEN
     150           0 :                CPWARN("This application of the delta pulse is not compatible with PBC!")
     151             :             END IF
     152           8 :             CALL apply_delta_pulse_mag(qs_env, mos_old, mos_new)
     153             :          ELSE
     154           0 :             CPABORT("Code error: this case should not happen!")
     155             :          END IF
     156             :       END IF
     157             : 
     158          54 :    END SUBROUTINE apply_delta_pulse
     159             : 
     160             : ! **************************************************************************************************
     161             : !> \brief uses perturbation theory to get the proper initial conditions
     162             : !>        The len_rep option is NOT compatible with periodic boundary conditions!
     163             : !> \param qs_env ...
     164             : !> \param mos_old ...
     165             : !> \param mos_new ...
     166             : !> \author Joost & Martin (2011)
     167             : ! **************************************************************************************************
     168             : 
     169          64 :    SUBROUTINE apply_delta_pulse_electric_periodic(qs_env, mos_old, mos_new)
     170             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     171             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_old, mos_new
     172             : 
     173             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric_periodic'
     174             : 
     175             :       INTEGER                                            :: handle, icol, idir, irow, ispin, nao, &
     176             :                                                             ncol_local, nmo, nrow_local, nvirt, &
     177             :                                                             reference
     178          32 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
     179             :       LOGICAL                                            :: com_nl, len_rep, periodic_cell
     180             :       REAL(KIND=dp)                                      :: eps_ppnl, factor
     181             :       REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :), &
     182          32 :          POINTER                                         :: local_data
     183             :       REAL(KIND=dp), DIMENSION(3)                        :: kvec, rcc
     184          32 :       REAL(kind=dp), DIMENSION(:), POINTER               :: eigenvalues, ref_point
     185             :       TYPE(cell_type), POINTER                           :: cell
     186             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct, fm_struct_tmp
     187             :       TYPE(cp_fm_type)                                   :: eigenvectors, mat_ks, mat_tmp, momentum, &
     188             :                                                             S_chol, virtuals
     189          32 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_r, matrix_rv, matrix_s
     190             :       TYPE(dft_control_type), POINTER                    :: dft_control
     191          32 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     192             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     193          32 :          POINTER                                         :: sab_orb, sap_ppnl
     194          32 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     195          32 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     196             :       TYPE(rt_prop_type), POINTER                        :: rtp
     197             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     198             :       TYPE(section_vals_type), POINTER                   :: input
     199             : 
     200          32 :       CALL timeset(routineN, handle)
     201             : 
     202          32 :       NULLIFY (cell, mos, rtp, matrix_s, matrix_ks, input, dft_control, particle_set, fm_struct)
     203             :       ! we need the overlap and ks matrix for a full diagionalization
     204             :       CALL get_qs_env(qs_env, &
     205             :                       cell=cell, &
     206             :                       mos=mos, &
     207             :                       rtp=rtp, &
     208             :                       matrix_s=matrix_s, &
     209             :                       matrix_ks=matrix_ks, &
     210             :                       dft_control=dft_control, &
     211             :                       input=input, &
     212          32 :                       particle_set=particle_set)
     213             : 
     214          32 :       rtp_control => dft_control%rtp_control
     215          98 :       periodic_cell = ANY(cell%perd > 0)
     216             : 
     217             :       ! relevant input parameters
     218          32 :       com_nl = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%COM_NL")
     219          32 :       len_rep = section_get_lval(section_vals=input, keyword_name="DFT%REAL_TIME_PROPAGATION%LEN_REP")
     220             : 
     221             :       ! calculate non-local commutator if necessary
     222          32 :       IF (com_nl) THEN
     223          16 :          CALL cite_reference(Mattiat2019)
     224          16 :          NULLIFY (qs_kind_set, sab_orb, sap_ppnl)
     225             :          CALL get_qs_env(qs_env, &
     226             :                          sap_ppnl=sap_ppnl, &
     227             :                          sab_orb=sab_orb, &
     228          16 :                          qs_kind_set=qs_kind_set)
     229          16 :          eps_ppnl = dft_control%qs_control%eps_ppnl
     230             : 
     231          16 :          NULLIFY (matrix_rv)
     232          16 :          CALL dbcsr_allocate_matrix_set(matrix_rv, 3)
     233          64 :          DO idir = 1, 3
     234          48 :             CALL dbcsr_init_p(matrix_rv(idir)%matrix)
     235             :             CALL dbcsr_create(matrix_rv(idir)%matrix, template=matrix_s(1)%matrix, &
     236          48 :                               matrix_type=dbcsr_type_antisymmetric)
     237          48 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_rv(idir)%matrix, sab_orb)
     238          64 :             CALL dbcsr_set(matrix_rv(idir)%matrix, 0._dp)
     239             :          END DO
     240          16 :          CALL build_com_mom_nl(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, cell, matrix_rv=matrix_rv)
     241             :       END IF
     242             : 
     243             :       ! calculate dipole moment matrix if required, NOT for periodic boundary conditions!
     244          32 :       IF (len_rep) THEN
     245          10 :          CALL cite_reference(Mattiat2022)
     246          10 :          IF (periodic_cell) THEN
     247           2 :             CPWARN("This application of the delta pulse is not compatible with PBC!")
     248             :          END IF
     249             :          ! get reference point
     250             :          reference = section_get_ival(section_vals=input, &
     251          10 :                                       keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
     252          10 :          NULLIFY (ref_point)
     253          10 :          CALL section_vals_val_get(input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
     254          10 :          CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
     255             : 
     256          10 :          NULLIFY (sab_orb)
     257          10 :          CALL get_qs_env(qs_env, sab_orb=sab_orb)
     258             :          ! calculate dipole moment operator
     259          10 :          NULLIFY (matrix_r)
     260          10 :          CALL dbcsr_allocate_matrix_set(matrix_r, 3)
     261          40 :          DO idir = 1, 3
     262          30 :             CALL dbcsr_init_p(matrix_r(idir)%matrix)
     263          30 :             CALL dbcsr_create(matrix_r(idir)%matrix, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_symmetric)
     264          30 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_r(idir)%matrix, sab_orb)
     265          40 :             CALL dbcsr_set(matrix_r(idir)%matrix, 0._dp)
     266             :          END DO
     267          10 :          CALL build_local_moment_matrix(qs_env, matrix_r, 1, rcc)
     268             :       END IF
     269             : 
     270             :       ! the prefactor (strength of the electric field)
     271             :       kvec(:) = cell%h_inv(1, :)*rtp_control%delta_pulse_direction(1) + &
     272             :                 cell%h_inv(2, :)*rtp_control%delta_pulse_direction(2) + &
     273         128 :                 cell%h_inv(3, :)*rtp_control%delta_pulse_direction(3)
     274         128 :       kvec = -kvec*twopi*rtp_control%delta_pulse_scale
     275             : 
     276          32 :       IF (rtp_control%velocity_gauge) THEN
     277           0 :          rtp_control%vec_pot = rtp_control%vec_pot + kvec
     278             :       END IF
     279             : 
     280             :       ! struct for fm matrices
     281          32 :       fm_struct => rtp%ao_ao_fmstruct
     282             : 
     283             :       ! create matrices and get Cholesky decomposition of S
     284          32 :       CALL cp_fm_create(mat_ks, matrix_struct=fm_struct, name="mat_ks")
     285          32 :       CALL cp_fm_create(eigenvectors, matrix_struct=fm_struct, name="eigenvectors")
     286          32 :       CALL cp_fm_create(S_chol, matrix_struct=fm_struct, name="S_chol")
     287          32 :       CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
     288          32 :       CALL cp_fm_cholesky_decompose(S_chol)
     289             : 
     290             :       ! get number of atomic orbitals
     291          32 :       CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
     292             : 
     293          80 :       DO ispin = 1, SIZE(matrix_ks)
     294             :          ! diagonalize KS matrix to get occ and virt mos
     295         144 :          ALLOCATE (eigenvalues(nao))
     296          48 :          CALL cp_fm_create(mat_tmp, matrix_struct=fm_struct, name="mat_tmp")
     297          48 :          CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks)
     298          48 :          CALL cp_fm_cholesky_reduce(mat_ks, S_chol)
     299          48 :          CALL cp_fm_syevd(mat_ks, mat_tmp, eigenvalues)
     300          48 :          CALL cp_fm_cholesky_restore(mat_tmp, nao, S_chol, eigenvectors, "SOLVE")
     301             : 
     302             :          ! virtuals
     303          48 :          CALL get_mo_set(mo_set=mos(ispin), nmo=nmo)
     304          48 :          nvirt = nao - nmo
     305             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=fm_struct%para_env, context=fm_struct%context, &
     306          48 :                                   nrow_global=nao, ncol_global=nvirt)
     307          48 :          CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals")
     308          48 :          CALL cp_fm_struct_release(fm_struct_tmp)
     309          48 :          CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1)
     310             : 
     311             :          ! occupied
     312          48 :          CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
     313             : 
     314             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=fm_struct%para_env, context=fm_struct%context, &
     315          48 :                                   nrow_global=nvirt, ncol_global=nmo)
     316          48 :          CALL cp_fm_create(momentum, matrix_struct=fm_struct_tmp, name="momentum")
     317          48 :          CALL cp_fm_struct_release(fm_struct_tmp)
     318             : 
     319             :          ! the momentum operator (in a given direction)
     320          48 :          CALL cp_fm_set_all(mos_new(2*ispin - 1), 0.0_dp)
     321             : 
     322         192 :          DO idir = 1, 3
     323         144 :             factor = kvec(idir)
     324         192 :             IF (factor .NE. 0.0_dp) THEN
     325          48 :                IF (.NOT. len_rep) THEN
     326             :                   CALL cp_dbcsr_sm_fm_multiply(matrix_s(idir + 1)%matrix, mos_old(2*ispin - 1), &
     327          34 :                                                mos_old(2*ispin), ncol=nmo)
     328             :                ELSE
     329             :                   CALL cp_dbcsr_sm_fm_multiply(matrix_r(idir)%matrix, mos_old(2*ispin - 1), &
     330          14 :                                                mos_old(2*ispin), ncol=nmo)
     331             :                END IF
     332             : 
     333          48 :                CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1), factor, mos_old(2*ispin))
     334          48 :                IF (com_nl) THEN
     335          24 :                   CALL cp_fm_set_all(mos_old(2*ispin), 0.0_dp)
     336             :                   CALL cp_dbcsr_sm_fm_multiply(matrix_rv(idir)%matrix, mos_old(2*ispin - 1), &
     337          24 :                                                mos_old(2*ispin), ncol=nmo)
     338          24 :                   CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1), factor, mos_old(2*ispin))
     339             :                END IF
     340             :             END IF
     341             :          END DO
     342             : 
     343          48 :          CALL parallel_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1), 0.0_dp, momentum)
     344             : 
     345             :          ! the tricky bit ... rescale by the eigenvalue difference
     346          48 :          IF (.NOT. len_rep) THEN
     347             :             CALL cp_fm_get_info(momentum, nrow_local=nrow_local, ncol_local=ncol_local, &
     348          34 :                                 row_indices=row_indices, col_indices=col_indices, local_data=local_data)
     349         232 :             DO icol = 1, ncol_local
     350        4810 :                DO irow = 1, nrow_local
     351        4578 :                   factor = 1/(eigenvalues(col_indices(icol)) - eigenvalues(nmo + row_indices(irow)))
     352        4776 :                   local_data(irow, icol) = factor*local_data(irow, icol)
     353             :                END DO
     354             :             END DO
     355             :          END IF
     356          48 :          CALL cp_fm_release(mat_tmp)
     357          48 :          DEALLOCATE (eigenvalues)
     358             : 
     359             :          ! now obtain the initial condition in mos_old
     360          48 :          CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
     361          48 :          CALL parallel_gemm("N", "N", nao, nmo, nvirt, 1.0_dp, virtuals, momentum, 0.0_dp, mos_old(2*ispin))
     362             : 
     363          48 :          CALL cp_fm_release(virtuals)
     364         128 :          CALL cp_fm_release(momentum)
     365             :       END DO
     366             : 
     367             :       ! release matrices
     368          32 :       CALL cp_fm_release(S_chol)
     369          32 :       CALL cp_fm_release(mat_ks)
     370          32 :       CALL cp_fm_release(eigenvectors)
     371          32 :       IF (com_nl) CALL dbcsr_deallocate_matrix_set(matrix_rv)
     372          32 :       IF (len_rep) CALL dbcsr_deallocate_matrix_set(matrix_r)
     373             : 
     374             :       ! orthonormalize afterwards
     375          32 :       CALL orthonormalize_complex_mos(qs_env, mos_old)
     376             : 
     377          32 :       CALL timestop(handle)
     378             : 
     379          32 :    END SUBROUTINE apply_delta_pulse_electric_periodic
     380             : 
     381             : ! **************************************************************************************************
     382             : !> \brief applies exp(ikr) to the wavefunction.... stored in mos_old...
     383             : !> \param qs_env ...
     384             : !> \param mos_old ...
     385             : !> \param mos_new ...
     386             : !> \author Joost & Martin (2011)
     387             : ! **************************************************************************************************
     388             : 
     389          14 :    SUBROUTINE apply_delta_pulse_electric(qs_env, mos_old, mos_new)
     390             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     391             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_old, mos_new
     392             : 
     393             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_electric'
     394             : 
     395             :       INTEGER                                            :: handle, i, nao, nmo
     396             :       REAL(KIND=dp), DIMENSION(3)                        :: kvec
     397             :       TYPE(cell_type), POINTER                           :: cell
     398             :       TYPE(cp_fm_type)                                   :: S_inv_fm, tmp
     399          14 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     400             :       TYPE(dbcsr_type), POINTER                          :: cosmat, sinmat
     401             :       TYPE(dft_control_type), POINTER                    :: dft_control
     402          14 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     403             :       TYPE(rt_prop_type), POINTER                        :: rtp
     404             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     405             : 
     406          14 :       CALL timeset(routineN, handle)
     407          14 :       NULLIFY (cell, dft_control, matrix_s, mos, rtp, rtp_control)
     408             :       CALL get_qs_env(qs_env, &
     409             :                       cell=cell, &
     410             :                       dft_control=dft_control, &
     411             :                       matrix_s=matrix_s, &
     412             :                       mos=mos, &
     413          14 :                       rtp=rtp)
     414          14 :       rtp_control => dft_control%rtp_control
     415             : 
     416             :       ! direction ... unscaled, this will yield a exp(ikr) that is periodic with the cell
     417             :       kvec(:) = cell%h_inv(1, :)*rtp_control%delta_pulse_direction(1) + &
     418             :                 cell%h_inv(2, :)*rtp_control%delta_pulse_direction(2) + &
     419          56 :                 cell%h_inv(3, :)*rtp_control%delta_pulse_direction(3)
     420          56 :       kvec = -kvec*twopi
     421             :       ! scaling will make the things not periodic with the cell, which would only be good for gas phase systems ?
     422          56 :       kvec(:) = dft_control%rtp_control%delta_pulse_scale*kvec
     423             : 
     424          14 :       IF (rtp_control%velocity_gauge) THEN
     425           0 :          rtp_control%vec_pot = rtp_control%vec_pot + kvec
     426             :       END IF
     427             : 
     428             :       ! calculate exponentials (= Berry moments)
     429          14 :       NULLIFY (cosmat, sinmat)
     430          14 :       ALLOCATE (cosmat, sinmat)
     431          14 :       CALL dbcsr_copy(cosmat, matrix_s(1)%matrix, 'COS MOM')
     432          14 :       CALL dbcsr_copy(sinmat, matrix_s(1)%matrix, 'SIN MOM')
     433          14 :       CALL build_berry_moment_matrix(qs_env, cosmat, sinmat, kvec)
     434             : 
     435             :       ! need inverse of overlap matrix
     436          14 :       CALL cp_fm_create(S_inv_fm, matrix_struct=rtp%ao_ao_fmstruct, name="S_inv_fm")
     437          14 :       CALL cp_fm_create(tmp, matrix_struct=rtp%ao_ao_fmstruct, name="tmp_mat")
     438          14 :       CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_inv_fm)
     439          14 :       CALL cp_fm_cholesky_decompose(S_inv_fm)
     440          14 :       CALL cp_fm_cholesky_invert(S_inv_fm)
     441          14 :       CALL cp_fm_upper_to_full(S_inv_fm, tmp)
     442          14 :       CALL cp_fm_release(tmp)
     443             : 
     444          38 :       DO i = 1, SIZE(mos)
     445             :          ! apply exponentials to mo coefficients
     446          24 :          CALL get_mo_set(mos(i), nao=nao, nmo=nmo)
     447          24 :          CALL cp_dbcsr_sm_fm_multiply(cosmat, mos(i)%mo_coeff, mos_new(2*i - 1), ncol=nmo)
     448          24 :          CALL cp_dbcsr_sm_fm_multiply(sinmat, mos(i)%mo_coeff, mos_new(2*i), ncol=nmo)
     449             : 
     450          24 :          CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i - 1), 0.0_dp, mos_old(2*i - 1))
     451          62 :          CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_inv_fm, mos_new(2*i), 0.0_dp, mos_old(2*i))
     452             :       END DO
     453             : 
     454          14 :       CALL cp_fm_release(S_inv_fm)
     455          14 :       CALL dbcsr_deallocate_matrix(cosmat)
     456          14 :       CALL dbcsr_deallocate_matrix(sinmat)
     457             : 
     458             :       ! orthonormalize afterwards
     459          14 :       CALL orthonormalize_complex_mos(qs_env, mos_old)
     460             : 
     461          14 :       CALL timestop(handle)
     462             : 
     463          14 :    END SUBROUTINE apply_delta_pulse_electric
     464             : 
     465             : ! **************************************************************************************************
     466             : !> \brief apply magnetic delta pulse to linear order
     467             : !> \param qs_env ...
     468             : !> \param mos_old ...
     469             : !> \param mos_new ...
     470             : ! **************************************************************************************************
     471          16 :    SUBROUTINE apply_delta_pulse_mag(qs_env, mos_old, mos_new)
     472             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     473             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_old, mos_new
     474             : 
     475             :       CHARACTER(len=*), PARAMETER :: routineN = 'apply_delta_pulse_mag'
     476             : 
     477             :       INTEGER                                            :: gauge_orig, handle, idir, ispin, nao, &
     478             :                                                             nmo, nrow_global, nvirt
     479             :       REAL(KIND=dp)                                      :: eps_ppnl, factor
     480             :       REAL(KIND=dp), DIMENSION(3)                        :: kvec, rcc
     481           8 :       REAL(kind=dp), DIMENSION(:), POINTER               :: eigenvalues, ref_point
     482             :       TYPE(cell_type), POINTER                           :: cell
     483             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     484             :       TYPE(cp_fm_type)                                   :: eigenvectors, mat_ks, perturbation, &
     485             :                                                             S_chol, virtuals
     486           8 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_mag, matrix_nl, &
     487           8 :                                                             matrix_s
     488             :       TYPE(dft_control_type), POINTER                    :: dft_control
     489           8 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     490             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     491           8 :          POINTER                                         :: sab_all, sab_orb, sap_ppnl
     492           8 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     493           8 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     494             :       TYPE(rt_prop_type), POINTER                        :: rtp
     495             :       TYPE(section_vals_type), POINTER                   :: input
     496             : 
     497           8 :       CALL timeset(routineN, handle)
     498             : 
     499           8 :       CALL cite_reference(Mattiat2022)
     500             : 
     501           8 :       NULLIFY (rtp, dft_control, matrix_ks, matrix_s, input, mos, cell, sab_orb, sab_all, sap_ppnl, &
     502           8 :                qs_kind_set, particle_set)
     503             : 
     504             :       CALL get_qs_env(qs_env, &
     505             :                       rtp=rtp, &
     506             :                       dft_control=dft_control, &
     507             :                       mos=mos, &
     508             :                       matrix_ks=matrix_ks, &
     509             :                       matrix_s=matrix_s, &
     510             :                       input=input, &
     511             :                       cell=cell, &
     512             :                       sab_orb=sab_orb, &
     513             :                       sab_all=sab_all, &
     514           8 :                       sap_ppnl=sap_ppnl)
     515             : 
     516             :       gauge_orig = section_get_ival(section_vals=input, &
     517           8 :                                     keyword_name="DFT%REAL_TIME_PROPAGATION%GAUGE_ORIG")
     518           8 :       NULLIFY (ref_point)
     519           8 :       CALL section_vals_val_get(input, "DFT%REAL_TIME_PROPAGATION%GAUGE_ORIG_MANUAL", r_vals=ref_point)
     520           8 :       CALL get_reference_point(rcc, qs_env=qs_env, reference=gauge_orig, ref_point=ref_point)
     521             : 
     522             :       ! Create fm matrices
     523           8 :       CALL cp_fm_create(S_chol, matrix_struct=rtp%ao_ao_fmstruct, name='Cholesky S')
     524           8 :       CALL cp_fm_create(eigenvectors, matrix_struct=rtp%ao_ao_fmstruct, name="gs evecs fm")
     525           8 :       CALL cp_fm_create(mat_ks, matrix_struct=rtp%ao_ao_fmstruct, name='KS matrix')
     526             : 
     527             :       ! get nrows_global
     528           8 :       CALL cp_fm_get_info(mat_ks, nrow_global=nrow_global)
     529             : 
     530             :       ! cholesky decomposition of overlap matrix
     531           8 :       CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_chol)
     532           8 :       CALL cp_fm_cholesky_decompose(S_chol)
     533             : 
     534             :       ! initiate perturbation matrix
     535           8 :       NULLIFY (matrix_mag)
     536           8 :       CALL dbcsr_allocate_matrix_set(matrix_mag, 3)
     537          32 :       DO idir = 1, 3
     538          24 :          CALL dbcsr_init_p(matrix_mag(idir)%matrix)
     539             :          CALL dbcsr_create(matrix_mag(idir)%matrix, template=matrix_s(1)%matrix, &
     540          24 :                            matrix_type=dbcsr_type_antisymmetric)
     541          24 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_mag(idir)%matrix, sab_orb)
     542          32 :          CALL dbcsr_set(matrix_mag(idir)%matrix, 0._dp)
     543             :       END DO
     544             :       ! construct magnetic dipole moment matrix
     545           8 :       CALL build_local_magmom_matrix(qs_env, matrix_mag, 1, ref_point=rcc)
     546             : 
     547             :       ! work matrix for non-local potential part if necessary
     548           8 :       NULLIFY (matrix_nl)
     549           8 :       IF (ASSOCIATED(sap_ppnl)) THEN
     550           8 :          CALL dbcsr_allocate_matrix_set(matrix_nl, 3)
     551          32 :          DO idir = 1, 3
     552          24 :             CALL dbcsr_init_p(matrix_nl(idir)%matrix)
     553             :             CALL dbcsr_create(matrix_nl(idir)%matrix, template=matrix_s(1)%matrix, &
     554          24 :                               matrix_type=dbcsr_type_antisymmetric)
     555          24 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_nl(idir)%matrix, sab_orb)
     556          32 :             CALL dbcsr_set(matrix_nl(idir)%matrix, 0._dp)
     557             :          END DO
     558             :          ! construct non-local contribution
     559             :          CALL get_qs_env(qs_env, &
     560             :                          qs_kind_set=qs_kind_set, &
     561           8 :                          particle_set=particle_set)
     562           8 :          eps_ppnl = dft_control%qs_control%eps_ppnl
     563             : 
     564           8 :          CALL build_com_nl_mag(qs_kind_set, sab_orb, sap_ppnl, eps_ppnl, particle_set, matrix_nl, rcc, cell)
     565             : 
     566          32 :          DO idir = 1, 3
     567          32 :             CALL dbcsr_add(matrix_mag(idir)%matrix, matrix_nl(idir)%matrix, -one, one)
     568             :          END DO
     569             : 
     570           8 :          CALL dbcsr_deallocate_matrix_set(matrix_nl)
     571             :       END IF
     572             : 
     573          20 :       DO ispin = 1, dft_control%nspins
     574             :          ! allocate eigenvalues
     575             :          NULLIFY (eigenvalues)
     576          36 :          ALLOCATE (eigenvalues(nrow_global))
     577             :          ! diagonalize KS matrix in AO basis using Cholesky decomp. of S
     578          12 :          CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, mat_ks)
     579          12 :          CALL cp_fm_cholesky_reduce(mat_ks, S_chol)
     580          12 :          CALL cp_fm_syevd(mat_ks, eigenvectors, eigenvalues)
     581          12 :          CALL cp_fm_triangular_multiply(S_chol, eigenvectors, invert_tr=.TRUE.)
     582             : 
     583             :          ! virtuals
     584          12 :          CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo)
     585          12 :          nvirt = nao - nmo
     586             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=rtp%ao_ao_fmstruct%para_env, context=rtp%ao_ao_fmstruct%context, &
     587          12 :                                   nrow_global=nrow_global, ncol_global=nvirt)
     588          12 :          CALL cp_fm_create(virtuals, matrix_struct=fm_struct_tmp, name="virtuals")
     589          12 :          CALL cp_fm_struct_release(fm_struct_tmp)
     590          12 :          CALL cp_fm_to_fm(eigenvectors, virtuals, nvirt, nmo + 1, 1)
     591             : 
     592             :          ! occupied
     593          12 :          CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
     594             : 
     595             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=rtp%ao_ao_fmstruct%para_env, context=rtp%ao_ao_fmstruct%context, &
     596          12 :                                   nrow_global=nvirt, ncol_global=nmo)
     597          12 :          CALL cp_fm_create(perturbation, matrix_struct=fm_struct_tmp, name="perturbation")
     598          12 :          CALL cp_fm_struct_release(fm_struct_tmp)
     599             : 
     600             :          ! apply perturbation
     601          12 :          CALL cp_fm_set_all(mos_new(2*ispin - 1), 0.0_dp)
     602             : 
     603             :          ! the prefactor (strength of the magnetic field, should be divided by 2c)
     604             :          kvec(:) = cell%h_inv(1, :)*dft_control%rtp_control%delta_pulse_direction(1) + &
     605             :                    cell%h_inv(2, :)*dft_control%rtp_control%delta_pulse_direction(2) + &
     606          48 :                    cell%h_inv(3, :)*dft_control%rtp_control%delta_pulse_direction(3)
     607          48 :          kvec = -kvec*twopi*dft_control%rtp_control%delta_pulse_scale
     608             : 
     609          48 :          DO idir = 1, 3
     610          36 :             factor = kvec(idir)/2         ! divide by two b.c. of pre-factor for magnetic term
     611          48 :             IF (factor .NE. 0.0_dp) THEN
     612             :                CALL cp_dbcsr_sm_fm_multiply(matrix_mag(idir)%matrix, mos_old(2*ispin - 1), &
     613          12 :                                             mos_old(2*ispin), ncol=nmo)
     614          12 :                CALL cp_fm_scale_and_add(1.0_dp, mos_new(2*ispin - 1), factor, mos_old(2*ispin))
     615             :             END IF
     616             :          END DO
     617             : 
     618          12 :          CALL parallel_gemm('T', 'N', nvirt, nmo, nao, 1.0_dp, virtuals, mos_new(2*ispin - 1), 0.0_dp, perturbation)
     619             : 
     620          12 :          DEALLOCATE (eigenvalues)
     621             : 
     622             :          ! now obtain the initial condition in mos_old
     623          12 :          CALL cp_fm_to_fm(eigenvectors, mos_old(2*ispin - 1), nmo, 1, 1)
     624          12 :          CALL parallel_gemm("N", "N", nao, nmo, nvirt, 1.0_dp, virtuals, perturbation, 0.0_dp, mos_old(2*ispin))
     625             : 
     626          12 :          CALL cp_fm_release(virtuals)
     627          32 :          CALL cp_fm_release(perturbation)
     628             :       END DO
     629             : 
     630             :       ! deallocations
     631           8 :       CALL cp_fm_release(S_chol)
     632           8 :       CALL cp_fm_release(mat_ks)
     633           8 :       CALL cp_fm_release(eigenvectors)
     634           8 :       CALL dbcsr_deallocate_matrix_set(matrix_mag)
     635             : 
     636             :       ! orthonormalize afterwards
     637           8 :       CALL orthonormalize_complex_mos(qs_env, mos_old)
     638             : 
     639           8 :       CALL timestop(handle)
     640             : 
     641           8 :    END SUBROUTINE apply_delta_pulse_mag
     642             : 
     643             : ! **************************************************************************************************
     644             : !> \brief orthonormalize complex mos, e. g. after non-unitary transformations using Löwdin's algorithm
     645             : !> \param qs_env ...
     646             : !> \param coeffs ...
     647             : ! **************************************************************************************************
     648          54 :    SUBROUTINE orthonormalize_complex_mos(qs_env, coeffs)
     649             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     650             :       TYPE(cp_fm_type), DIMENSION(:), INTENT(INOUT), &
     651             :          POINTER                                         :: coeffs
     652             : 
     653          54 :       COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:)        :: eigenvalues_sqrt
     654             :       INTEGER                                            :: im, ispin, j, nao, nmo, nspins, re
     655          54 :       REAL(kind=dp), ALLOCATABLE, DIMENSION(:)           :: eigenvalues
     656             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     657             :       TYPE(cp_cfm_type)                                  :: oo_c, oo_v, oo_vt
     658             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     659             :       TYPE(cp_fm_type)                                   :: oo_1, oo_2, S_fm, tmp
     660         162 :       TYPE(cp_fm_type), DIMENSION(2)                     :: coeffs_tmp
     661          54 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     662             :       TYPE(dft_control_type), POINTER                    :: dft_control
     663          54 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     664             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     665             : 
     666          54 :       NULLIFY (para_env, blacs_env, dft_control, matrix_s, mos)
     667             :       CALL get_qs_env(qs_env, &
     668             :                       blacs_env=blacs_env, &
     669             :                       dft_control=dft_control, &
     670             :                       matrix_s=matrix_s, &
     671             :                       mos=mos, &
     672          54 :                       para_env=para_env)
     673          54 :       nspins = dft_control%nspins
     674          54 :       CALL cp_fm_get_info(coeffs(1), nrow_global=nao)
     675             : 
     676             :       ! get overlap matrix
     677             :       CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nao, ncol_global=nao, &
     678          54 :                                context=blacs_env, para_env=para_env)
     679          54 :       CALL cp_fm_create(S_fm, matrix_struct=fm_struct_tmp, name="overlap fm")
     680          54 :       CALL cp_fm_struct_release(fm_struct_tmp)
     681             :       ! copy overlap matrix
     682          54 :       CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, S_fm)
     683             : 
     684         138 :       DO ispin = 1, nspins
     685          84 :          CALL get_mo_set(mos(ispin), nmo=nmo)
     686             :          CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
     687          84 :                                   nrow_global=nmo, ncol_global=nmo)
     688          84 :          CALL cp_fm_create(oo_1, matrix_struct=fm_struct_tmp, name="oo_1")
     689          84 :          CALL cp_fm_create(oo_2, matrix_struct=fm_struct_tmp, name="oo_2")
     690          84 :          CALL cp_fm_struct_release(fm_struct_tmp)
     691             : 
     692          84 :          CALL cp_fm_create(tmp, matrix_struct=coeffs(2*ispin - 1)%matrix_struct, name="tmp_mat")
     693             :          ! get the complex overlap matrix in MO basis
     694             :          ! x^T S x + y^T S y + i (-y^TS x+x^T S y)
     695          84 :          CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_fm, coeffs(2*ispin - 1), 0.0_dp, tmp)
     696          84 :          CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin - 1), tmp, 0.0_dp, oo_1)
     697          84 :          CALL parallel_gemm("T", "N", nmo, nmo, nao, -1.0_dp, coeffs(2*ispin), tmp, 0.0_dp, oo_2)
     698             : 
     699          84 :          CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, S_fm, coeffs(2*ispin), 0.0_dp, tmp)
     700          84 :          CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin), tmp, 1.0_dp, oo_1)
     701          84 :          CALL parallel_gemm("T", "N", nmo, nmo, nao, 1.0_dp, coeffs(2*ispin - 1), tmp, 1.0_dp, oo_2)
     702          84 :          CALL cp_fm_release(tmp)
     703             : 
     704             :          ! complex Löwdin
     705          84 :          CALL cp_cfm_create(oo_c, oo_1%matrix_struct)
     706          84 :          CALL cp_cfm_create(oo_v, oo_1%matrix_struct)
     707          84 :          CALL cp_cfm_create(oo_vt, oo_1%matrix_struct)
     708        1995 :          oo_c%local_data = CMPLX(oo_1%local_data, oo_2%local_data, KIND=dp)
     709             : 
     710         252 :          ALLOCATE (eigenvalues(nmo))
     711         252 :          ALLOCATE (eigenvalues_sqrt(nmo))
     712          84 :          CALL cp_cfm_heevd(oo_c, oo_v, eigenvalues)
     713         482 :          eigenvalues_sqrt(:) = CMPLX(one/SQRT(eigenvalues(:)), zero, dp)
     714          84 :          CALL cp_cfm_to_cfm(oo_v, oo_vt)
     715          84 :          CALL cp_cfm_column_scale(oo_v, eigenvalues_sqrt)
     716          84 :          DEALLOCATE (eigenvalues)
     717          84 :          DEALLOCATE (eigenvalues_sqrt)
     718             :          CALL parallel_gemm('N', 'C', nmo, nmo, nmo, (1.0_dp, 0.0_dp), &
     719          84 :                             oo_v, oo_vt, (0.0_dp, 0.0_dp), oo_c)
     720        1995 :          oo_1%local_data = REAL(oo_c%local_data, KIND=dp)
     721        1995 :          oo_2%local_data = AIMAG(oo_c%local_data)
     722          84 :          CALL cp_cfm_release(oo_c)
     723          84 :          CALL cp_cfm_release(oo_v)
     724          84 :          CALL cp_cfm_release(oo_vt)
     725             : 
     726             :          ! transform coefficients accordingly
     727         252 :          DO j = 1, 2
     728         252 :             CALL cp_fm_create(coeffs_tmp(j), matrix_struct=coeffs(2*(ispin - 1) + j)%matrix_struct)
     729             :          END DO
     730             : 
     731             :          ! indices for coeffs_tmp
     732          84 :          re = 1
     733          84 :          im = 2
     734          84 :          CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin - 1), oo_1, zero, coeffs_tmp(re))
     735          84 :          CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin - 1), oo_2, zero, coeffs_tmp(im))
     736             : 
     737          84 :          CALL parallel_gemm("N", "N", nao, nmo, nmo, -one, coeffs(2*ispin), oo_2, zero, coeffs(2*ispin - 1))
     738          84 :          CALL cp_fm_scale_and_add(one, coeffs(2*ispin - 1), one, coeffs_tmp(re))
     739             : 
     740          84 :          CALL parallel_gemm("N", "N", nao, nmo, nmo, one, coeffs(2*ispin), oo_1, one, coeffs_tmp(im))
     741          84 :          CALL cp_fm_to_fm(coeffs_tmp(im), coeffs(2*ispin))
     742             : 
     743         252 :          DO j = 1, 2
     744         252 :             CALL cp_fm_release(coeffs_tmp(j))
     745             :          END DO
     746          84 :          CALL cp_fm_release(oo_1)
     747         222 :          CALL cp_fm_release(oo_2)
     748             :       END DO
     749          54 :       CALL cp_fm_release(S_fm)
     750             : 
     751         108 :    END SUBROUTINE orthonormalize_complex_mos
     752             : 
     753             : END MODULE rt_delta_pulse

Generated by: LCOV version 1.15