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

Generated by: LCOV version 1.15