LCOV - code coverage report
Current view: top level - src/emd - rt_projection_mo_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 229 239 95.8 %
Date: 2024-04-20 06:29:22 Functions: 6 6 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 Function related to MO projection in RTP calculations
      10             : !> \author Guillaume Le Breton 04.2023
      11             : ! **************************************************************************************************
      12             : MODULE rt_projection_mo_utils
      13             :    USE cp_control_types,                ONLY: dft_control_type,&
      14             :                                               proj_mo_type,&
      15             :                                               rtp_control_type
      16             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply
      17             :    USE cp_files,                        ONLY: close_file,&
      18             :                                               open_file
      19             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_scale_and_add,&
      20             :                                               cp_fm_trace
      21             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      22             :                                               cp_fm_struct_release,&
      23             :                                               cp_fm_struct_type
      24             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      25             :                                               cp_fm_release,&
      26             :                                               cp_fm_to_fm,&
      27             :                                               cp_fm_type
      28             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      29             :                                               cp_logger_get_default_io_unit,&
      30             :                                               cp_logger_type,&
      31             :                                               cp_to_string
      32             :    USE cp_output_handling,              ONLY: cp_p_file,&
      33             :                                               cp_print_key_finished_output,&
      34             :                                               cp_print_key_generate_filename,&
      35             :                                               cp_print_key_should_output,&
      36             :                                               cp_print_key_unit_nr
      37             :    USE dbcsr_api,                       ONLY: dbcsr_p_type
      38             :    USE input_constants,                 ONLY: proj_mo_ref_scf,&
      39             :                                               proj_mo_ref_xas_tdp
      40             :    USE input_section_types,             ONLY: section_vals_get,&
      41             :                                               section_vals_get_subs_vals,&
      42             :                                               section_vals_type,&
      43             :                                               section_vals_val_get
      44             :    USE kinds,                           ONLY: default_string_length,&
      45             :                                               dp
      46             :    USE message_passing,                 ONLY: mp_para_env_type
      47             :    USE particle_types,                  ONLY: particle_type
      48             :    USE qs_environment_types,            ONLY: get_qs_env,&
      49             :                                               qs_environment_type
      50             :    USE qs_kind_types,                   ONLY: qs_kind_type
      51             :    USE qs_mo_io,                        ONLY: read_mos_restart_low
      52             :    USE qs_mo_types,                     ONLY: deallocate_mo_set,&
      53             :                                               mo_set_type
      54             :    USE rt_propagation_types,            ONLY: get_rtp,&
      55             :                                               rt_prop_type
      56             : #include "./../base/base_uses.f90"
      57             : 
      58             :    IMPLICIT NONE
      59             :    PRIVATE
      60             : 
      61             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_projection_mo_utils'
      62             : 
      63             :    PUBLIC :: init_mo_projection, compute_and_write_proj_mo
      64             : 
      65             : CONTAINS
      66             : 
      67             : ! **************************************************************************************************
      68             : !> \brief Initialize the mo projection objects for time dependent run
      69             : !> \param qs_env ...
      70             : !> \param rtp_control ...
      71             : !> \author Guillaume Le Breton (04.2023)
      72             : ! **************************************************************************************************
      73           4 :    SUBROUTINE init_mo_projection(qs_env, rtp_control)
      74             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      75             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
      76             : 
      77             :       INTEGER                                            :: i_rep, j_td, n_rep_val, nbr_mo_td_max, &
      78             :                                                             nrep, reftype
      79           4 :       INTEGER, DIMENSION(:), POINTER                     :: tmp_ints
      80             :       TYPE(cp_logger_type), POINTER                      :: logger
      81           4 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
      82             :       TYPE(proj_mo_type), POINTER                        :: proj_mo
      83             :       TYPE(section_vals_type), POINTER                   :: input, print_key, proj_mo_section
      84             : 
      85           4 :       NULLIFY (rtp_control%proj_mo_list, tmp_ints, proj_mo, logger, &
      86           4 :                input, proj_mo_section, print_key, mos)
      87             : 
      88           4 :       CALL get_qs_env(qs_env, input=input, mos=mos)
      89             : 
      90           4 :       proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
      91             : 
      92             :       ! Read the input section and load the reference MOs
      93           4 :       CALL section_vals_get(proj_mo_section, n_repetition=nrep)
      94          64 :       ALLOCATE (rtp_control%proj_mo_list(nrep))
      95             : 
      96          56 :       DO i_rep = 1, nrep
      97          52 :          NULLIFY (rtp_control%proj_mo_list(i_rep)%proj_mo)
      98          52 :          ALLOCATE (rtp_control%proj_mo_list(i_rep)%proj_mo)
      99          52 :          proj_mo => rtp_control%proj_mo_list(i_rep)%proj_mo
     100             : 
     101             :          CALL section_vals_val_get(proj_mo_section, "REFERENCE_TYPE", i_rep_section=i_rep, &
     102          52 :                                    i_val=reftype)
     103             : 
     104             :          CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, &
     105          52 :                                    c_val=proj_mo%ref_mo_file_name)
     106             : 
     107             :          CALL section_vals_val_get(proj_mo_section, "REF_ADD_LUMO", i_rep_section=i_rep, &
     108          52 :                                    i_val=proj_mo%ref_nlumo)
     109             : 
     110             :          ! Relevent only in EMD
     111          52 :          IF (.NOT. rtp_control%fixed_ions) &
     112             :             CALL section_vals_val_get(proj_mo_section, "PROPAGATE_REF", i_rep_section=i_rep, &
     113          28 :                                       l_val=proj_mo%propagate_ref)
     114             : 
     115          52 :          IF (reftype == proj_mo_ref_scf) THEN
     116             :             ! If no reference .wfn is provided, using the restart SCF file:
     117          46 :             IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
     118          38 :                CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
     119          38 :                IF (n_rep_val > 0) THEN
     120           0 :                   CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", c_val=proj_mo%ref_mo_file_name)
     121             :                ELSE
     122             :                   !try to read from the filename that is generated automatically from the printkey
     123          38 :                   print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
     124          38 :                   logger => cp_get_default_logger()
     125             :                   proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
     126          38 :                                                                             extension=".wfn", my_local=.FALSE.)
     127             :                END IF
     128             :             END IF
     129             : 
     130             :             CALL section_vals_val_get(proj_mo_section, "REF_MO_INDEX", i_rep_section=i_rep, &
     131          46 :                                       i_vals=tmp_ints)
     132         194 :             ALLOCATE (proj_mo%ref_mo_index, SOURCE=tmp_ints(:))
     133             :             CALL section_vals_val_get(proj_mo_section, "REF_MO_SPIN", i_rep_section=i_rep, &
     134          46 :                                       i_val=proj_mo%ref_mo_spin)
     135             : 
     136             :             ! Read the SCF mos and store the one required
     137          46 :             CALL read_reference_mo_from_wfn(qs_env, proj_mo)
     138             : 
     139           6 :          ELSE IF (reftype == proj_mo_ref_xas_tdp) THEN
     140           6 :             IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
     141             :                CALL cp_abort(__LOCATION__, &
     142             :                              "Input error in DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO. "// &
     143             :                              "For REFERENCE_TYPE XAS_TDP one must define the name "// &
     144           0 :                              "of the .wfn file to read the reference MO from. Please define REF_MO_FILE_NAME.")
     145             :             END IF
     146           6 :             ALLOCATE (proj_mo%ref_mo_index(1))
     147             :             ! XAS restart files contain only one excited state
     148           6 :             proj_mo%ref_mo_index(1) = 1
     149           6 :             proj_mo%ref_mo_spin = 1
     150             :             ! Read XAS TDP mos
     151           6 :             CALL read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref=.TRUE.)
     152             : 
     153             :          END IF
     154             : 
     155             :          ! Initialize the other parameters related to the TD mos.
     156             :          CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_REF", i_rep_section=i_rep, &
     157          52 :                                    l_val=proj_mo%sum_on_all_ref)
     158             : 
     159             :          CALL section_vals_val_get(proj_mo_section, "TD_MO_SPIN", i_rep_section=i_rep, &
     160          52 :                                    i_val=proj_mo%td_mo_spin)
     161          52 :          IF (proj_mo%td_mo_spin > SIZE(mos)) &
     162             :             CALL cp_abort(__LOCATION__, &
     163             :                           "You asked to project the time dependent BETA spin while the "// &
     164             :                           "real time DFT run has only one spin defined. "// &
     165           0 :                           "Please set TD_MO_SPIN to 1 or use UKS.")
     166             : 
     167             :          CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, &
     168          52 :                                    i_vals=tmp_ints)
     169             : 
     170          52 :          nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
     171             : 
     172         218 :          ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:))
     173          52 :          IF (proj_mo%td_mo_index(1) == -1) THEN
     174          24 :             DEALLOCATE (proj_mo%td_mo_index)
     175          72 :             ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
     176          72 :             ALLOCATE (proj_mo%td_mo_occ(nbr_mo_td_max))
     177         168 :             DO j_td = 1, nbr_mo_td_max
     178         144 :                proj_mo%td_mo_index(j_td) = j_td
     179         168 :                proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
     180             :             END DO
     181             :          ELSE
     182          84 :             ALLOCATE (proj_mo%td_mo_occ(SIZE(proj_mo%td_mo_index)))
     183          66 :             proj_mo%td_mo_occ(:) = 0.0_dp
     184          66 :             DO j_td = 1, SIZE(proj_mo%td_mo_index)
     185          38 :                IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) &
     186             :                   CALL cp_abort(__LOCATION__, &
     187             :                                 "The MO number available in the Time Dependent run "// &
     188           0 :                                 "is smaller than the MO number you have required in TD_MO_INDEX.")
     189          66 :                proj_mo%td_mo_occ(j_td) = mos(proj_mo%td_mo_spin)%occupation_numbers(proj_mo%td_mo_index(j_td))
     190             :             END DO
     191             :          END IF
     192             : 
     193             :          CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_TD", i_rep_section=i_rep, &
     194         108 :                                    l_val=proj_mo%sum_on_all_td)
     195             : 
     196             :       END DO
     197             : 
     198           8 :    END SUBROUTINE init_mo_projection
     199             : 
     200             : ! **************************************************************************************************
     201             : !> \brief Read the MO from .wfn file and store the required mos for TD projections
     202             : !> \param qs_env ...
     203             : !> \param proj_mo ...
     204             : !> \param xas_ref ...
     205             : !> \author Guillaume Le Breton (04.2023)
     206             : ! **************************************************************************************************
     207          52 :    SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref)
     208             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     209             :       TYPE(proj_mo_type), POINTER                        :: proj_mo
     210             :       LOGICAL, OPTIONAL                                  :: xas_ref
     211             : 
     212             :       INTEGER                                            :: i_ref, ispin, mo_index, natom, &
     213             :                                                             nbr_mo_max, nbr_ref_mo, nspins, &
     214             :                                                             real_mo_index, restart_unit
     215             :       LOGICAL                                            :: is_file, my_xasref
     216             :       TYPE(cp_fm_struct_type), POINTER                   :: mo_ref_fmstruct
     217             :       TYPE(cp_fm_type)                                   :: mo_coeff_temp
     218          52 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     219             :       TYPE(dft_control_type), POINTER                    :: dft_control
     220          52 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_qs, mo_ref_temp
     221             :       TYPE(mo_set_type), POINTER                         :: mo_set
     222             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     223          52 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     224          52 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     225             : 
     226          52 :       NULLIFY (mo_qs, mo_ref_temp, mo_set, qs_kind_set, particle_set, para_env, dft_control, &
     227          52 :                mo_ref_fmstruct, matrix_s)
     228             : 
     229          52 :       my_xasref = .FALSE.
     230           6 :       IF (PRESENT(xas_ref)) my_xasref = xas_ref
     231             : 
     232             :       CALL get_qs_env(qs_env, &
     233             :                       qs_kind_set=qs_kind_set, &
     234             :                       particle_set=particle_set, &
     235             :                       dft_control=dft_control, &
     236             :                       matrix_s_kp=matrix_s, &
     237             :                       mos=mo_qs, &
     238          52 :                       para_env=para_env)
     239             : 
     240          52 :       natom = SIZE(particle_set, 1)
     241             : 
     242          52 :       nspins = SIZE(mo_qs)
     243             :       ! If the restart comes from DFT%XAS_TDP%PRINT%RESTART_WFN, then always 2 spins are saved
     244          52 :       IF (my_xasref .AND. nspins < 2) THEN
     245           0 :          nspins = 2
     246             :       END IF
     247         260 :       ALLOCATE (mo_ref_temp(nspins))
     248             : 
     249         156 :       DO ispin = 1, nspins
     250         104 :          IF (my_xasref) THEN
     251          12 :             mo_set => mo_qs(1)
     252             :          ELSE
     253          92 :             mo_set => mo_qs(ispin)
     254             :          END IF
     255         104 :          mo_ref_temp(ispin)%nmo = mo_set%nmo + proj_mo%ref_nlumo
     256         104 :          NULLIFY (mo_ref_fmstruct)
     257             :          CALL cp_fm_struct_create(mo_ref_fmstruct, nrow_global=mo_set%nao, &
     258         104 :                                ncol_global=mo_ref_temp(ispin)%nmo, para_env=para_env, context=mo_set%mo_coeff%matrix_struct%context)
     259         104 :          NULLIFY (mo_ref_temp(ispin)%mo_coeff)
     260         104 :          ALLOCATE (mo_ref_temp(ispin)%mo_coeff)
     261         104 :          CALL cp_fm_create(mo_ref_temp(ispin)%mo_coeff, mo_ref_fmstruct)
     262         104 :          CALL cp_fm_struct_release(mo_ref_fmstruct)
     263             : 
     264         104 :          mo_ref_temp(ispin)%nao = mo_set%nao
     265         104 :          mo_ref_temp(ispin)%homo = mo_set%homo
     266         104 :          mo_ref_temp(ispin)%nelectron = mo_set%nelectron
     267         312 :          ALLOCATE (mo_ref_temp(ispin)%eigenvalues(mo_ref_temp(ispin)%nmo))
     268         312 :          ALLOCATE (mo_ref_temp(ispin)%occupation_numbers(mo_ref_temp(ispin)%nmo))
     269         156 :          NULLIFY (mo_set)
     270             :       END DO
     271             : 
     272             : !         DO ispin = 1, nspins
     273             : !            CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(1))
     274             : !         END DO
     275             : !      ELSE
     276             : !         DO ispin = 1, nspins
     277             : !            CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(ispin))
     278             : !         END DO
     279             : !      END IF
     280             : 
     281          52 :       IF (para_env%is_source()) THEN
     282          26 :          INQUIRE (FILE=TRIM(proj_mo%ref_mo_file_name), exist=is_file)
     283          26 :          IF (.NOT. is_file) &
     284             :             CALL cp_abort(__LOCATION__, &
     285           0 :                           "Reference file not found! Name of the file CP2K looked for: "//TRIM(proj_mo%ref_mo_file_name))
     286             : 
     287             :          CALL open_file(file_name=proj_mo%ref_mo_file_name, &
     288             :                         file_action="READ", &
     289             :                         file_form="UNFORMATTED", &
     290             :                         file_status="OLD", &
     291          26 :                         unit_number=restart_unit)
     292             :       END IF
     293             : 
     294             :       CALL read_mos_restart_low(mo_ref_temp, para_env=para_env, qs_kind_set=qs_kind_set, &
     295             :                                 particle_set=particle_set, natom=natom, &
     296          52 :                                 rst_unit=restart_unit)
     297             : 
     298          52 :       IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
     299             : 
     300          52 :       IF (proj_mo%ref_mo_spin > SIZE(mo_ref_temp)) &
     301             :          CALL cp_abort(__LOCATION__, &
     302             :                        "You asked as reference spin the BETA one while the "// &
     303             :                        "reference .wfn file has only one spin. Use a reference .wfn "// &
     304           0 :                        "with 2 spins separated or set REF_MO_SPIN to 1")
     305             : 
     306             :       ! Store only the mos required
     307          52 :       nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
     308          52 :       IF (proj_mo%ref_mo_index(1) == -1) THEN
     309          18 :          DEALLOCATE (proj_mo%ref_mo_index)
     310          54 :          ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
     311         126 :          DO i_ref = 1, nbr_mo_max
     312         126 :             proj_mo%ref_mo_index(i_ref) = i_ref
     313             :          END DO
     314             :       ELSE
     315          78 :          DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
     316          44 :             IF (proj_mo%ref_mo_index(i_ref) > nbr_mo_max) &
     317             :                CALL cp_abort(__LOCATION__, &
     318             :                              "The MO number available in the Reference SCF "// &
     319          34 :                              "is smaller than the MO number you have required in REF_MO_INDEX.")
     320             :          END DO
     321             :       END IF
     322          52 :       nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
     323             : 
     324          52 :       IF (nbr_ref_mo > nbr_mo_max) &
     325             :          CALL cp_abort(__LOCATION__, &
     326           0 :                        "The number of reference mo is larger then the total number of available one in the .wfn file.")
     327             : 
     328             :       ! Store
     329         308 :       ALLOCATE (proj_mo%mo_ref(nbr_ref_mo))
     330             :       CALL cp_fm_struct_create(mo_ref_fmstruct, &
     331             :                                context=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%context, &
     332             :                                nrow_global=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%nrow_global, &
     333          52 :                                ncol_global=1)
     334             : 
     335          52 :       IF (dft_control%rtp_control%fixed_ions) &
     336          24 :          CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
     337             : 
     338         204 :       DO mo_index = 1, nbr_ref_mo
     339         152 :          real_mo_index = proj_mo%ref_mo_index(mo_index)
     340         152 :          IF (real_mo_index > nbr_mo_max) &
     341             :             CALL cp_abort(__LOCATION__, &
     342           0 :                           "One of reference mo index is larger then the total number of available mo in the .wfn file.")
     343             : 
     344             :          ! fill with the reference mo values
     345         152 :          CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
     346         204 :          IF (dft_control%rtp_control%fixed_ions) THEN
     347             :             ! multiply with overlap matrix to save time later on: proj_mo%mo_ref is SxMO_ref
     348             :             CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, mo_coeff_temp, &
     349             :                              ncol=1, &
     350             :                              source_start=real_mo_index, &
     351          64 :                              target_start=1)
     352          64 :             CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff_temp, proj_mo%mo_ref(mo_index), ncol=1)
     353             :          ELSE
     354             :             ! the AO will change with times: proj_mo%mo_ref are really the MOs coeffs
     355             :             CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, proj_mo%mo_ref(mo_index), &
     356             :                              ncol=1, &
     357             :                              source_start=real_mo_index, &
     358          88 :                              target_start=1)
     359             :          END IF
     360             :       END DO
     361             : 
     362             :       ! Clean temporary variables
     363         156 :       DO ispin = 1, nspins
     364         156 :          CALL deallocate_mo_set(mo_ref_temp(ispin))
     365             :       END DO
     366          52 :       DEALLOCATE (mo_ref_temp)
     367             : 
     368          52 :       CALL cp_fm_struct_release(mo_ref_fmstruct)
     369          52 :       IF (dft_control%rtp_control%fixed_ions) &
     370          24 :          CALL cp_fm_release(mo_coeff_temp)
     371             : 
     372          52 :    END SUBROUTINE read_reference_mo_from_wfn
     373             : 
     374             : ! **************************************************************************************************
     375             : !> \brief Compute the projection of the current MO coefficients on reference ones
     376             : !>        and write the results.
     377             : !> \param qs_env ...
     378             : !> \param mos_new ...
     379             : !> \param proj_mo ...
     380             : !> \param n_proj ...
     381             : !> \author Guillaume Le Breton
     382             : ! **************************************************************************************************
     383         104 :    SUBROUTINE compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
     384             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     385             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     386             :       TYPE(proj_mo_type)                                 :: proj_mo
     387             :       INTEGER                                            :: n_proj
     388             : 
     389             :       INTEGER                                            :: i_ref, nbr_ref_mo, nbr_ref_td
     390         104 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: phase, popu, sum_popu_ref
     391             :       TYPE(cp_fm_struct_type), POINTER                   :: mo_ref_fmstruct
     392             :       TYPE(cp_fm_type)                                   :: S_mo_ref
     393             :       TYPE(cp_logger_type), POINTER                      :: logger
     394         104 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     395             :       TYPE(dft_control_type), POINTER                    :: dft_control
     396             :       TYPE(section_vals_type), POINTER                   :: input, print_mo_section, proj_mo_section
     397             : 
     398         104 :       NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
     399             : 
     400         208 :       logger => cp_get_default_logger()
     401             : 
     402             :       CALL get_qs_env(qs_env, &
     403             :                       dft_control=dft_control, &
     404         104 :                       input=input)
     405             : 
     406             :       ! The general section
     407         104 :       proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
     408             :       ! The section we are dealing in this particular subroutine call: n_proj.
     409         104 :       print_mo_section => section_vals_get_subs_vals(proj_mo_section, "PRINT", i_rep_section=n_proj)
     410             : 
     411             :       ! Propagate the reference MO if required at each time step
     412         104 :       IF (proj_mo%propagate_ref) CALL propagate_ref_mo(qs_env, proj_mo)
     413             : 
     414             :       ! Does not compute the projection if not the required time step
     415         104 :       IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
     416             :                                                  print_mo_section, ""), &
     417             :                       cp_p_file)) &
     418             :          RETURN
     419             : 
     420         102 :       IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
     421             :          CALL get_qs_env(qs_env, &
     422          56 :                          matrix_s_kp=matrix_s)
     423             :          CALL cp_fm_struct_create(mo_ref_fmstruct, &
     424             :                                   context=proj_mo%mo_ref(1)%matrix_struct%context, &
     425             :                                   nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
     426          56 :                                   ncol_global=1)
     427          56 :          CALL cp_fm_create(S_mo_ref, mo_ref_fmstruct, 'S_mo_ref')
     428             :       END IF
     429             : 
     430         102 :       nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
     431         102 :       nbr_ref_td = SIZE(proj_mo%td_mo_index)
     432         306 :       ALLOCATE (popu(nbr_ref_td))
     433         306 :       ALLOCATE (phase(nbr_ref_td))
     434             : 
     435         102 :       IF (proj_mo%sum_on_all_ref) THEN
     436          72 :          ALLOCATE (sum_popu_ref(nbr_ref_td))
     437         168 :          sum_popu_ref(:) = 0.0_dp
     438         168 :          DO i_ref = 1, nbr_ref_mo
     439             :             ! Compute SxMO_ref for the upcoming projection later on
     440         144 :             IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
     441          96 :                CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), S_mo_ref, ncol=1)
     442          96 :                CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref=S_mo_ref)
     443             :             ELSE
     444          48 :                CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
     445             :             END IF
     446        1032 :             sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
     447             :          END DO
     448          24 :          IF (proj_mo%sum_on_all_td) THEN
     449          84 :             CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=SUM(sum_popu_ref), n_proj=n_proj)
     450             :          ELSE
     451          12 :             CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
     452             :          END IF
     453          24 :          DEALLOCATE (sum_popu_ref)
     454             :       ELSE
     455         234 :          DO i_ref = 1, nbr_ref_mo
     456         156 :             IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
     457          80 :                CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, proj_mo%mo_ref(i_ref), S_mo_ref, ncol=1)
     458          80 :                CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref=S_mo_ref)
     459             :             ELSE
     460          76 :                CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
     461             :             END IF
     462         234 :             IF (proj_mo%sum_on_all_td) THEN
     463         588 :                CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu_tot=SUM(popu), n_proj=n_proj)
     464             :             ELSE
     465             : 
     466          72 :                CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu=popu, phase=phase, n_proj=n_proj)
     467             :             END IF
     468             :          END DO
     469             :       END IF
     470             : 
     471         102 :       IF (.NOT. dft_control%rtp_control%fixed_ions) THEN
     472          56 :          CALL cp_fm_struct_release(mo_ref_fmstruct)
     473          56 :          CALL cp_fm_release(S_mo_ref)
     474             :       END IF
     475         102 :       DEALLOCATE (popu)
     476         102 :       DEALLOCATE (phase)
     477             : 
     478         208 :    END SUBROUTINE compute_and_write_proj_mo
     479             : 
     480             : ! **************************************************************************************************
     481             : !> \brief Compute the projection of the current MO coefficients on reference ones
     482             : !> \param popu ...
     483             : !> \param phase ...
     484             : !> \param mos_new ...
     485             : !> \param proj_mo ...
     486             : !> \param i_ref ...
     487             : !> \param S_mo_ref ...
     488             : !> \author Guillaume Le Breton
     489             : ! **************************************************************************************************
     490         300 :    SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref, S_mo_ref)
     491             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: popu, phase
     492             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     493             :       TYPE(proj_mo_type)                                 :: proj_mo
     494             :       INTEGER                                            :: i_ref
     495             :       TYPE(cp_fm_type), OPTIONAL                         :: S_mo_ref
     496             : 
     497             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_proj_mo'
     498             : 
     499             :       INTEGER                                            :: handle, j_td, nbr_ref_td, spin_td
     500             :       LOGICAL                                            :: is_emd
     501             :       REAL(KIND=dp)                                      :: imag_proj, real_proj
     502             :       TYPE(cp_fm_struct_type), POINTER                   :: mo_ref_fmstruct
     503             :       TYPE(cp_fm_type)                                   :: mo_coeff_temp
     504             : 
     505         300 :       CALL timeset(routineN, handle)
     506             : 
     507         300 :       is_emd = .FALSE.
     508         300 :       IF (PRESENT(S_mo_ref)) is_emd = .TRUE.
     509             : 
     510         300 :       nbr_ref_td = SIZE(popu)
     511         300 :       spin_td = proj_mo%td_mo_spin
     512             : 
     513             :       CALL cp_fm_struct_create(mo_ref_fmstruct, &
     514             :                                context=mos_new(1)%matrix_struct%context, &
     515             :                                nrow_global=mos_new(1)%matrix_struct%nrow_global, &
     516         300 :                                ncol_global=1)
     517         300 :       CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
     518             : 
     519        1784 :       DO j_td = 1, nbr_ref_td
     520             :          ! Real part of the projection:
     521             :          real_proj = 0.0_dp
     522             :          CALL cp_fm_to_fm(mos_new(2*spin_td - 1), mo_coeff_temp, &
     523             :                           ncol=1, &
     524             :                           source_start=proj_mo%td_mo_index(j_td), &
     525        1484 :                           target_start=1)
     526        1484 :          IF (is_emd) THEN
     527             :             ! The reference MO have to be propagated in the new basis, so the projection
     528         936 :             CALL cp_fm_trace(mo_coeff_temp, S_mo_ref, real_proj)
     529             :          ELSE
     530             :             ! The reference MO is time independent. proj_mo%mo_ref(i_ref) is in fact SxMO_ref already
     531         548 :             CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), real_proj)
     532             :          END IF
     533             : 
     534             :          ! Imaginary part of the projection
     535             :          imag_proj = 0.0_dp
     536             :          CALL cp_fm_to_fm(mos_new(2*spin_td), mo_coeff_temp, &
     537             :                           ncol=1, &
     538             :                           source_start=proj_mo%td_mo_index(j_td), &
     539        1484 :                           target_start=1)
     540             : 
     541        1484 :          IF (is_emd) THEN
     542         936 :             CALL cp_fm_trace(mo_coeff_temp, S_mo_ref, imag_proj)
     543             :          ELSE
     544         548 :             CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
     545             :          END IF
     546             : 
     547             :          ! Store the result
     548        1484 :          phase(j_td) = ATAN2(imag_proj, real_proj) ! in radians
     549        1784 :          popu(j_td) = proj_mo%td_mo_occ(j_td)*(real_proj**2 + imag_proj**2)
     550             :       END DO
     551             : 
     552         300 :       CALL cp_fm_struct_release(mo_ref_fmstruct)
     553         300 :       CALL cp_fm_release(mo_coeff_temp)
     554             : 
     555         300 :       CALL timestop(handle)
     556             : 
     557         300 :    END SUBROUTINE compute_proj_mo
     558             : 
     559             : ! **************************************************************************************************
     560             : !> \brief Write in one file the projection of (all) the time-dependent MO coefficients
     561             : !>        on one reference ones
     562             : !> \param qs_env ...
     563             : !> \param print_mo_section ...
     564             : !> \param proj_mo ...
     565             : !> \param i_ref ...
     566             : !> \param popu ...
     567             : !> \param phase ...
     568             : !> \param popu_tot ...
     569             : !> \param n_proj ...
     570             : !> \author Guillaume Le Breton
     571             : ! **************************************************************************************************
     572         180 :    SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, popu_tot, n_proj)
     573             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     574             :       TYPE(section_vals_type), POINTER                   :: print_mo_section
     575             :       TYPE(proj_mo_type)                                 :: proj_mo
     576             :       INTEGER, OPTIONAL                                  :: i_ref
     577             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: popu, phase
     578             :       REAL(KIND=dp), OPTIONAL                            :: popu_tot
     579             :       INTEGER, OPTIONAL                                  :: n_proj
     580             : 
     581             :       CHARACTER(LEN=default_string_length)               :: ext, filename
     582             :       INTEGER                                            :: j_td, output_unit, print_unit
     583             :       TYPE(cp_logger_type), POINTER                      :: logger
     584             : 
     585         180 :       NULLIFY (logger)
     586             : 
     587         180 :       logger => cp_get_default_logger()
     588         180 :       output_unit = cp_logger_get_default_io_unit(logger)
     589             : 
     590         180 :       IF (.NOT. (output_unit > 0)) RETURN
     591             : 
     592          90 :       IF (proj_mo%sum_on_all_ref) THEN
     593          12 :          ext = "-"//TRIM(ADJUSTL(cp_to_string(n_proj)))//"-ALL_REF.dat"
     594             :       ELSE
     595             :          ! Filename is update wrt the reference MO number
     596             :          ext = "-"//TRIM(ADJUSTL(cp_to_string(n_proj)))// &
     597             :                "-REF-"// &
     598             :                TRIM(ADJUSTL(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
     599          78 :                ".dat"
     600             :       END IF
     601             : 
     602             :       print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
     603          90 :                                         extension=TRIM(ext))
     604             : 
     605          90 :       IF (print_unit /= output_unit) THEN
     606          90 :          INQUIRE (UNIT=print_unit, NAME=filename)
     607             : !         WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
     608             : !            "PROJECTION MO", "The projection of the TD MOs is done in the file:", &
     609             : !            TRIM(filename)
     610             :          WRITE (UNIT=print_unit, FMT="(/,(T2,A,T40,I6))") &
     611          90 :             "Real time propagation step:", qs_env%sim_step
     612             :       ELSE
     613           0 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") "PROJECTION MO"
     614             :       END IF
     615             : 
     616          90 :       IF (proj_mo%sum_on_all_ref) THEN
     617             :          WRITE (print_unit, "(T3,A)") &
     618             :             "Projection on all the required MO number from the reference file "// &
     619          12 :             TRIM(proj_mo%ref_mo_file_name)
     620          12 :          IF (proj_mo%sum_on_all_td) THEN
     621             :             WRITE (print_unit, "(T3, A, E20.12)") &
     622           6 :                "The sum over all the TD MOs population:", popu_tot
     623             :          ELSE
     624             :             WRITE (print_unit, "(T3,A)") &
     625           6 :                "For each TD MOs required is printed: Population "
     626          42 :             DO j_td = 1, SIZE(popu)
     627          42 :                WRITE (print_unit, "(T5,1(E20.12, 1X))") popu(j_td)
     628             :             END DO
     629             :          END IF
     630             :       ELSE
     631             :          WRITE (print_unit, "(T3,A)") &
     632             :             "Projection on the MO number "// &
     633             :             TRIM(ADJUSTL(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
     634             :             " from the reference file "// &
     635          78 :             TRIM(proj_mo%ref_mo_file_name)
     636             : 
     637          78 :          IF (proj_mo%sum_on_all_td) THEN
     638             :             WRITE (print_unit, "(T3, A, E20.12)") &
     639          42 :                "The sum over all the TD MOs population:", popu_tot
     640             :          ELSE
     641             :             WRITE (print_unit, "(T3,A)") &
     642          36 :                "For each TD MOs required is printed: Population & Phase [rad] "
     643          94 :             DO j_td = 1, SIZE(popu)
     644          94 :                WRITE (print_unit, "(T5,2(E20.12, E16.8, 1X))") popu(j_td), phase(j_td)
     645             :             END DO
     646             :          END IF
     647             :       END IF
     648             : 
     649          90 :       CALL cp_print_key_finished_output(print_unit, logger, print_mo_section, "")
     650             : 
     651             :    END SUBROUTINE write_proj_mo
     652             : 
     653             : ! **************************************************************************************************
     654             : !> \brief Propagate the reference MO in case of EMD: since the nuclei moves, the MO coeff can be
     655             : !>        propagate to represent the same MO (because the AO move with the nuclei).
     656             : !>        To do so, we use the same formula as for the electrons of the system, but without the
     657             : !>        Hamiltonian:
     658             : !>        dc^j_alpha/dt = - sum_{beta, gamma} S^{-1}_{alpha, beta} B_{beta,gamma} c^j_gamma
     659             : !> \param qs_env ...
     660             : !> \param proj_mo ...
     661             : !> \author Guillaume Le Breton
     662             : ! **************************************************************************************************
     663          56 :    SUBROUTINE propagate_ref_mo(qs_env, proj_mo)
     664             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     665             :       TYPE(proj_mo_type)                                 :: proj_mo
     666             : 
     667             :       INTEGER                                            :: i_ref
     668             :       REAL(Kind=dp)                                      :: dt
     669             :       TYPE(cp_fm_struct_type), POINTER                   :: mo_ref_fmstruct
     670             :       TYPE(cp_fm_type)                                   :: d_mo
     671          28 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: SinvB
     672             :       TYPE(rt_prop_type), POINTER                        :: rtp
     673             : 
     674          28 :       CALL get_qs_env(qs_env, rtp=rtp)
     675          28 :       CALL get_rtp(rtp=rtp, SinvB=SinvB, dt=dt)
     676             : 
     677             :       CALL cp_fm_struct_create(mo_ref_fmstruct, &
     678             :                                context=proj_mo%mo_ref(1)%matrix_struct%context, &
     679             :                                nrow_global=proj_mo%mo_ref(1)%matrix_struct%nrow_global, &
     680          28 :                                ncol_global=1)
     681          28 :       CALL cp_fm_create(d_mo, mo_ref_fmstruct, 'd_mo')
     682             : 
     683         116 :       DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
     684             :          ! MO(t+dt) = MO(t) - dtxS_inv.B(t).MO(t)
     685          88 :          CALL cp_dbcsr_sm_fm_multiply(SinvB(1)%matrix, proj_mo%mo_ref(i_ref), d_mo, ncol=1, alpha=-dt)
     686         116 :          CALL cp_fm_scale_and_add(1.0_dp, proj_mo%mo_ref(i_ref), 1.0_dp, d_mo)
     687             :       END DO
     688             : 
     689          28 :       CALL cp_fm_struct_release(mo_ref_fmstruct)
     690          28 :       CALL cp_fm_release(d_mo)
     691             : 
     692          28 :    END SUBROUTINE propagate_ref_mo
     693             : 
     694             : END MODULE rt_projection_mo_utils
     695             : 

Generated by: LCOV version 1.15