LCOV - code coverage report
Current view: top level - src/emd - rt_projection_mo_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:04dd428) Lines: 178 188 94.7 %
Date: 2023-09-28 06:57:18 Functions: 5 5 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief 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_trace
      20             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      21             :                                               cp_fm_struct_release,&
      22             :                                               cp_fm_struct_type
      23             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      24             :                                               cp_fm_release,&
      25             :                                               cp_fm_to_fm,&
      26             :                                               cp_fm_type
      27             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      28             :                                               cp_logger_get_default_io_unit,&
      29             :                                               cp_logger_type,&
      30             :                                               cp_to_string
      31             :    USE cp_output_handling,              ONLY: cp_p_file,&
      32             :                                               cp_print_key_finished_output,&
      33             :                                               cp_print_key_generate_filename,&
      34             :                                               cp_print_key_should_output,&
      35             :                                               cp_print_key_unit_nr
      36             :    USE dbcsr_api,                       ONLY: dbcsr_p_type
      37             :    USE input_constants,                 ONLY: proj_mo_ref_scf,&
      38             :                                               proj_mo_ref_xas_tdp
      39             :    USE input_section_types,             ONLY: section_vals_get,&
      40             :                                               section_vals_get_subs_vals,&
      41             :                                               section_vals_type,&
      42             :                                               section_vals_val_get
      43             :    USE kinds,                           ONLY: default_string_length,&
      44             :                                               dp
      45             :    USE message_passing,                 ONLY: mp_para_env_type
      46             :    USE particle_types,                  ONLY: particle_type
      47             :    USE qs_environment_types,            ONLY: get_qs_env,&
      48             :                                               qs_environment_type
      49             :    USE qs_kind_types,                   ONLY: qs_kind_type
      50             :    USE qs_mo_io,                        ONLY: read_mos_restart_low
      51             :    USE qs_mo_types,                     ONLY: deallocate_mo_set,&
      52             :                                               duplicate_mo_set,&
      53             :                                               mo_set_type
      54             : #include "./../base/base_uses.f90"
      55             : 
      56             :    IMPLICIT NONE
      57             :    PRIVATE
      58             : 
      59             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_projection_mo_utils'
      60             : 
      61             :    PUBLIC :: init_mo_projection, compute_and_write_proj_mo
      62             : 
      63             : CONTAINS
      64             : 
      65             : ! **************************************************************************************************
      66             : !> \brief Initialize the mo projection objects for time dependent run
      67             : !> \param qs_env ...
      68             : !> \param rtp_control ...
      69             : !> \author Guillaume Le Breton (04.2023)
      70             : ! **************************************************************************************************
      71           2 :    SUBROUTINE init_mo_projection(qs_env, rtp_control)
      72             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      73             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
      74             : 
      75             :       INTEGER                                            :: i_rep, j_td, n_rep_val, nbr_mo_td_max, &
      76             :                                                             nrep, reftype
      77           2 :       INTEGER, DIMENSION(:), POINTER                     :: tmp_ints
      78             :       TYPE(cp_logger_type), POINTER                      :: logger
      79           2 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
      80             :       TYPE(proj_mo_type), POINTER                        :: proj_mo
      81             :       TYPE(section_vals_type), POINTER                   :: input, print_key, proj_mo_section
      82             : 
      83           2 :       NULLIFY (rtp_control%proj_mo_list, tmp_ints, proj_mo, logger, &
      84           2 :                input, proj_mo_section, print_key, mos)
      85             : 
      86           2 :       IF (.NOT. rtp_control%fixed_ions) &
      87             :          CALL cp_abort(__LOCATION__, &
      88           0 :                        "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO not implemented for EMD yet.")
      89             : 
      90           2 :       CALL get_qs_env(qs_env, input=input, mos=mos)
      91             : 
      92           2 :       proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
      93             : 
      94             :       ! Read the input section and load the reference MOs
      95           2 :       CALL section_vals_get(proj_mo_section, n_repetition=nrep)
      96          26 :       ALLOCATE (rtp_control%proj_mo_list(nrep))
      97             : 
      98          22 :       DO i_rep = 1, nrep
      99          20 :          NULLIFY (rtp_control%proj_mo_list(i_rep)%proj_mo)
     100          20 :          ALLOCATE (rtp_control%proj_mo_list(i_rep)%proj_mo)
     101          20 :          proj_mo => rtp_control%proj_mo_list(i_rep)%proj_mo
     102             : 
     103             :          CALL section_vals_val_get(proj_mo_section, "REFERENCE_TYPE", i_rep_section=i_rep, &
     104          20 :                                    i_val=reftype)
     105             : 
     106             :          CALL section_vals_val_get(proj_mo_section, "REF_MO_FILE_NAME", i_rep_section=i_rep, &
     107          20 :                                    c_val=proj_mo%ref_mo_file_name)
     108             : 
     109          20 :          IF (reftype == proj_mo_ref_scf) THEN
     110             :             ! If no reference .wfn is provided, using the restart SCF file:
     111          18 :             IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
     112          16 :                CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", n_rep_val=n_rep_val)
     113          16 :                IF (n_rep_val > 0) THEN
     114           0 :                   CALL section_vals_val_get(input, "DFT%WFN_RESTART_FILE_NAME", c_val=proj_mo%ref_mo_file_name)
     115             :                ELSE
     116             :                   !try to read from the filename that is generated automatically from the printkey
     117          16 :                   print_key => section_vals_get_subs_vals(input, "DFT%SCF%PRINT%RESTART")
     118          16 :                   logger => cp_get_default_logger()
     119             :                   proj_mo%ref_mo_file_name = cp_print_key_generate_filename(logger, print_key, &
     120          16 :                                                                             extension=".wfn", my_local=.FALSE.)
     121             :                END IF
     122             :             END IF
     123             : 
     124             :             CALL section_vals_val_get(proj_mo_section, "REF_MO_INDEX", i_rep_section=i_rep, &
     125          18 :                                       i_vals=tmp_ints)
     126          82 :             ALLOCATE (proj_mo%ref_mo_index, SOURCE=tmp_ints(:))
     127             :             CALL section_vals_val_get(proj_mo_section, "REF_MO_SPIN", i_rep_section=i_rep, &
     128          18 :                                       i_val=proj_mo%ref_mo_spin)
     129             : 
     130             :             ! Read the SCF mos and store the one required
     131          18 :             CALL read_reference_mo_from_wfn(qs_env, proj_mo)
     132             : 
     133           2 :          ELSE IF (reftype == proj_mo_ref_xas_tdp) THEN
     134           2 :             IF (proj_mo%ref_mo_file_name == "DEFAULT") THEN
     135             :                CALL cp_abort(__LOCATION__, &
     136             :                              "Input error in DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO. "// &
     137             :                              "For REFERENCE_TYPE XAS_TDP one must define the name "// &
     138           0 :                              "of the .wfn file to read the reference MO from. Please define REF_MO_FILE_NAME.")
     139             :             END IF
     140           2 :             ALLOCATE (proj_mo%ref_mo_index(1))
     141             :             ! XAS restart files contain only one excited state
     142           2 :             proj_mo%ref_mo_index(1) = 1
     143           2 :             proj_mo%ref_mo_spin = 1
     144             :             ! Read XAS TDP mos
     145           2 :             CALL read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref=.TRUE.)
     146             : 
     147             :          END IF
     148             : 
     149             :          ! Initialize the other parameters related to the TD mos.
     150             :          CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_REF", i_rep_section=i_rep, &
     151          20 :                                    l_val=proj_mo%sum_on_all_ref)
     152             : 
     153             :          CALL section_vals_val_get(proj_mo_section, "TD_MO_SPIN", i_rep_section=i_rep, &
     154          20 :                                    i_val=proj_mo%td_mo_spin)
     155          20 :          IF (proj_mo%td_mo_spin > SIZE(mos)) &
     156             :             CALL cp_abort(__LOCATION__, &
     157             :                           "You asked to project the time dependent BETA spin while the "// &
     158             :                           "real time DFT run has only one spin defined. "// &
     159           0 :                           "Please set TD_MO_SPIN to 1 or use UKS.")
     160             : 
     161             :          CALL section_vals_val_get(proj_mo_section, "TD_MO_INDEX", i_rep_section=i_rep, &
     162          20 :                                    i_vals=tmp_ints)
     163             : 
     164          90 :          ALLOCATE (proj_mo%td_mo_index, SOURCE=tmp_ints(:))
     165          20 :          nbr_mo_td_max = mos(proj_mo%td_mo_spin)%mo_coeff%matrix_struct%ncol_global
     166          20 :          IF (proj_mo%td_mo_index(1) == -1) THEN
     167           8 :             DEALLOCATE (proj_mo%td_mo_index)
     168          24 :             ALLOCATE (proj_mo%td_mo_index(nbr_mo_td_max))
     169          56 :             DO j_td = 1, nbr_mo_td_max
     170          56 :                proj_mo%td_mo_index(j_td) = j_td
     171             :             END DO
     172             :          ELSE
     173          34 :             DO j_td = 1, SIZE(proj_mo%td_mo_index)
     174          22 :                IF (proj_mo%td_mo_index(j_td) > nbr_mo_td_max) &
     175             :                   CALL cp_abort(__LOCATION__, &
     176             :                                 "The MO number available in the Time Dependent run "// &
     177          12 :                                 "is smaller than the MO number you have required in TD_MO_INDEX.")
     178             :             END DO
     179             :          END IF
     180             : 
     181             :          CALL section_vals_val_get(proj_mo_section, "SUM_ON_ALL_TD", i_rep_section=i_rep, &
     182          42 :                                    l_val=proj_mo%sum_on_all_td)
     183             : 
     184             :       END DO
     185             : 
     186           4 :    END SUBROUTINE init_mo_projection
     187             : 
     188             : ! **************************************************************************************************
     189             : !> \brief Read the MO from .wfn file and store the required mos for TD projections
     190             : !> \param qs_env ...
     191             : !> \param proj_mo ...
     192             : !> \param xas_ref ...
     193             : !> \author Guillaume Le Breton (04.2023)
     194             : ! **************************************************************************************************
     195          20 :    SUBROUTINE read_reference_mo_from_wfn(qs_env, proj_mo, xas_ref)
     196             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     197             :       TYPE(proj_mo_type), POINTER                        :: proj_mo
     198             :       LOGICAL, OPTIONAL                                  :: xas_ref
     199             : 
     200             :       INTEGER                                            :: i_ref, ispin, mo_index, natom, &
     201             :                                                             nbr_mo_max, nbr_ref_mo, nspins, &
     202             :                                                             real_mo_index, restart_unit
     203             :       LOGICAL                                            :: is_file, my_xasref
     204             :       TYPE(cp_fm_struct_type), POINTER                   :: mo_ref_fmstruct
     205             :       TYPE(cp_fm_type)                                   :: mo_coeff_temp
     206          20 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     207             :       TYPE(dft_control_type), POINTER                    :: dft_control
     208          20 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_qs, mo_ref_temp
     209             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     210          20 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     211          20 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     212             : 
     213          20 :       NULLIFY (mo_qs, mo_ref_temp, qs_kind_set, particle_set, para_env, dft_control, &
     214          20 :                mo_ref_fmstruct, matrix_s)
     215             : 
     216          20 :       my_xasref = .FALSE.
     217           2 :       IF (PRESENT(xas_ref)) my_xasref = xas_ref
     218             : 
     219             :       CALL get_qs_env(qs_env, &
     220             :                       qs_kind_set=qs_kind_set, &
     221             :                       particle_set=particle_set, &
     222             :                       dft_control=dft_control, &
     223             :                       matrix_s_kp=matrix_s, &
     224             :                       mos=mo_qs, &
     225          20 :                       para_env=para_env)
     226             : 
     227          20 :       natom = SIZE(particle_set, 1)
     228             : 
     229             :       ! If the restart comes from DFT%XAS_TDP%PRINT%RESTART_WFN, then always 2 spins are saved
     230          20 :       IF (my_xasref) THEN
     231           2 :          nspins = 2
     232           6 :          ALLOCATE (mo_ref_temp(nspins))
     233           6 :          DO ispin = 1, nspins
     234           6 :             CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(1))
     235             :          END DO
     236             :       ELSE
     237          18 :          nspins = SIZE(mo_qs)
     238          90 :          ALLOCATE (mo_ref_temp(nspins))
     239          54 :          DO ispin = 1, nspins
     240          54 :             CALL duplicate_mo_set(mo_ref_temp(ispin), mo_qs(ispin))
     241             :          END DO
     242             :       END IF
     243             : 
     244          20 :       IF (para_env%is_source()) THEN
     245          10 :          INQUIRE (FILE=TRIM(proj_mo%ref_mo_file_name), exist=is_file)
     246          10 :          IF (.NOT. is_file) &
     247             :             CALL cp_abort(__LOCATION__, &
     248           0 :                           "Reference file not found! Name of the file CP2K looked for: "//TRIM(proj_mo%ref_mo_file_name))
     249             : 
     250             :          CALL open_file(file_name=proj_mo%ref_mo_file_name, &
     251             :                         file_action="READ", &
     252             :                         file_form="UNFORMATTED", &
     253             :                         file_status="OLD", &
     254          10 :                         unit_number=restart_unit)
     255             :       END IF
     256             : 
     257             :       CALL read_mos_restart_low(mo_ref_temp, para_env=para_env, qs_kind_set=qs_kind_set, &
     258             :                                 particle_set=particle_set, natom=natom, &
     259          20 :                                 rst_unit=restart_unit)
     260             : 
     261          20 :       IF (para_env%is_source()) CALL close_file(unit_number=restart_unit)
     262             : 
     263          20 :       IF (proj_mo%ref_mo_spin > SIZE(mo_ref_temp)) &
     264             :          CALL cp_abort(__LOCATION__, &
     265             :                        "You asked as reference spin the BETA one while the "// &
     266             :                        "reference .wfn file has only one spin. Use a reference .wfn "// &
     267           0 :                        "with 2 spins separated or set REF_MO_SPIN to 1")
     268             : 
     269             :       ! Store only the mos required
     270          20 :       nbr_mo_max = mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%ncol_global
     271          20 :       IF (proj_mo%ref_mo_index(1) == -1) THEN
     272           6 :          DEALLOCATE (proj_mo%ref_mo_index)
     273          18 :          ALLOCATE (proj_mo%ref_mo_index(nbr_mo_max))
     274          42 :          DO i_ref = 1, nbr_mo_max
     275          42 :             proj_mo%ref_mo_index(i_ref) = i_ref
     276             :          END DO
     277             :       ELSE
     278          38 :          DO i_ref = 1, SIZE(proj_mo%ref_mo_index)
     279          24 :             IF (proj_mo%ref_mo_index(i_ref) > nbr_mo_max) &
     280             :                CALL cp_abort(__LOCATION__, &
     281             :                              "The MO number available in the Reference SCF "// &
     282          14 :                              "is smaller than the MO number you have required in REF_MO_INDEX.")
     283             :          END DO
     284             :       END IF
     285          20 :       nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
     286             : 
     287          20 :       IF (nbr_ref_mo > nbr_mo_max) &
     288             :          CALL cp_abort(__LOCATION__, &
     289           0 :                        "The number of reference mo is larger then the total number of available one in the .wfn file.")
     290             : 
     291             :       ! Store
     292         120 :       ALLOCATE (proj_mo%mo_ref(nbr_ref_mo))
     293             :       CALL cp_fm_struct_create(mo_ref_fmstruct, &
     294             :                                context=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%context, &
     295             :                                nrow_global=mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff%matrix_struct%nrow_global, &
     296          20 :                                ncol_global=1)
     297          20 :       CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_ref')
     298             : 
     299          80 :       DO mo_index = 1, nbr_ref_mo
     300          60 :          real_mo_index = proj_mo%ref_mo_index(mo_index)
     301          60 :          IF (real_mo_index > nbr_mo_max) &
     302             :             CALL cp_abort(__LOCATION__, &
     303           0 :                           "One of reference mo index is larger then the total number of available mo in the .wfn file.")
     304             : 
     305             :          ! fill with the reference mo values
     306          60 :          CALL cp_fm_create(proj_mo%mo_ref(mo_index), mo_ref_fmstruct, 'mo_ref')
     307             :          CALL cp_fm_to_fm(mo_ref_temp(proj_mo%ref_mo_spin)%mo_coeff, mo_coeff_temp, &
     308             :                           ncol=1, &
     309             :                           source_start=real_mo_index, &
     310          60 :                           target_start=1)
     311             : 
     312             :          ! multiply with overlap matrix to save time later on
     313          80 :          CALL cp_dbcsr_sm_fm_multiply(matrix_s(1, 1)%matrix, mo_coeff_temp, proj_mo%mo_ref(mo_index), ncol=1)
     314             :       END DO
     315             : 
     316             :       ! Clean temporary variables
     317          60 :       DO ispin = 1, nspins
     318          60 :          CALL deallocate_mo_set(mo_ref_temp(ispin))
     319             :       END DO
     320          20 :       DEALLOCATE (mo_ref_temp)
     321             : 
     322          20 :       CALL cp_fm_struct_release(mo_ref_fmstruct)
     323          20 :       CALL cp_fm_release(mo_coeff_temp)
     324             : 
     325          20 :    END SUBROUTINE read_reference_mo_from_wfn
     326             : 
     327             : ! **************************************************************************************************
     328             : !> \brief Compute the projection of the current MO coefficients on reference ones
     329             : !>        and write the results.
     330             : !> \param qs_env ...
     331             : !> \param mos_new ...
     332             : !> \param proj_mo ...
     333             : !> \param n_proj ...
     334             : !> \author Guillaume Le Breton
     335             : ! **************************************************************************************************
     336          40 :    SUBROUTINE compute_and_write_proj_mo(qs_env, mos_new, proj_mo, n_proj)
     337             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     338             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     339             :       TYPE(proj_mo_type)                                 :: proj_mo
     340             :       INTEGER                                            :: n_proj
     341             : 
     342             :       INTEGER                                            :: i_ref, nbr_ref_mo, nbr_ref_td
     343          40 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: phase, popu, sum_popu_ref
     344             :       TYPE(cp_logger_type), POINTER                      :: logger
     345             :       TYPE(dft_control_type), POINTER                    :: dft_control
     346             :       TYPE(section_vals_type), POINTER                   :: input, print_mo_section, proj_mo_section
     347             : 
     348          40 :       NULLIFY (dft_control, input, proj_mo_section, print_mo_section, logger)
     349             : 
     350          40 :       logger => cp_get_default_logger()
     351             : 
     352             :       CALL get_qs_env(qs_env, &
     353             :                       dft_control=dft_control, &
     354          40 :                       input=input)
     355             : 
     356             :       ! The general section
     357          40 :       proj_mo_section => section_vals_get_subs_vals(input, "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO")
     358             :       ! The section we are dealing in this particular subroutine call: n_proj.
     359          40 :       print_mo_section => section_vals_get_subs_vals(proj_mo_section, "PRINT", i_rep_section=n_proj)
     360             : 
     361             :       ! STOP if not the required time step
     362          40 :       IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, &
     363             :                                                  print_mo_section, ""), &
     364             :                       cp_p_file)) &
     365             :          RETURN
     366             : 
     367          38 :       IF (.NOT. dft_control%rtp_control%fixed_ions) &
     368             :          CALL cp_abort(__LOCATION__, &
     369           0 :                        "DFT%REAL_TIME_PROPAGATION%PRINT%PROJECTION_MO not implemented for EMD yet.")
     370             : 
     371          38 :       nbr_ref_mo = SIZE(proj_mo%ref_mo_index)
     372          38 :       nbr_ref_td = SIZE(proj_mo%td_mo_index)
     373         114 :       ALLOCATE (popu(nbr_ref_td))
     374         114 :       ALLOCATE (phase(nbr_ref_td))
     375             : 
     376          38 :       IF (proj_mo%sum_on_all_ref) THEN
     377          24 :          ALLOCATE (sum_popu_ref(nbr_ref_td))
     378          56 :          sum_popu_ref(:) = 0.0_dp
     379          56 :          DO i_ref = 1, nbr_ref_mo
     380          48 :             CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
     381         344 :             sum_popu_ref(:) = sum_popu_ref(:) + popu(:)
     382             :          END DO
     383           8 :          IF (proj_mo%sum_on_all_td) THEN
     384          28 :             CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu_tot=SUM(sum_popu_ref), n_proj=n_proj)
     385             :          ELSE
     386           4 :             CALL write_proj_mo(qs_env, print_mo_section, proj_mo, popu=sum_popu_ref, n_proj=n_proj)
     387             :          END IF
     388           8 :          DEALLOCATE (sum_popu_ref)
     389             :       ELSE
     390          98 :          DO i_ref = 1, nbr_ref_mo
     391          68 :             CALL compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
     392             : 
     393          98 :             IF (proj_mo%sum_on_all_td) THEN
     394         196 :                CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu_tot=SUM(popu), n_proj=n_proj)
     395             :             ELSE
     396             : 
     397          40 :                CALL write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref=i_ref, popu=popu, phase=phase, n_proj=n_proj)
     398             :             END IF
     399             :          END DO
     400             :       END IF
     401             : 
     402          38 :       DEALLOCATE (popu)
     403          38 :       DEALLOCATE (phase)
     404             : 
     405          40 :    END SUBROUTINE compute_and_write_proj_mo
     406             : 
     407             : ! **************************************************************************************************
     408             : !> \brief Compute the projection of the current MO coefficients on reference ones
     409             : !> \param popu ...
     410             : !> \param phase ...
     411             : !> \param mos_new ...
     412             : !> \param proj_mo ...
     413             : !> \param i_ref ...
     414             : !> \author Guillaume Le Breton
     415             : ! **************************************************************************************************
     416         116 :    SUBROUTINE compute_proj_mo(popu, phase, mos_new, proj_mo, i_ref)
     417             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: popu, phase
     418             :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     419             :       TYPE(proj_mo_type)                                 :: proj_mo
     420             :       INTEGER                                            :: i_ref
     421             : 
     422             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'compute_proj_mo'
     423             : 
     424             :       INTEGER                                            :: handle, j_td, nbr_ref_td, spin_td
     425             :       REAL(KIND=dp)                                      :: imag_proj, real_proj
     426             :       TYPE(cp_fm_struct_type), POINTER                   :: mo_ref_fmstruct
     427             :       TYPE(cp_fm_type)                                   :: mo_coeff_temp
     428             : 
     429         116 :       CALL timeset(routineN, handle)
     430             : 
     431         116 :       nbr_ref_td = SIZE(popu)
     432         116 :       spin_td = proj_mo%td_mo_spin
     433             : 
     434             :       CALL cp_fm_struct_create(mo_ref_fmstruct, &
     435             :                                context=mos_new(1)%matrix_struct%context, &
     436             :                                nrow_global=mos_new(1)%matrix_struct%nrow_global, &
     437         116 :                                ncol_global=1)
     438         116 :       CALL cp_fm_create(mo_coeff_temp, mo_ref_fmstruct, 'mo_temp')
     439             : 
     440         656 :       DO j_td = 1, nbr_ref_td
     441             :          ! Real part of the projection:
     442             :          real_proj = 0.0_dp
     443             :          CALL cp_fm_to_fm(mos_new(2*spin_td - 1), mo_coeff_temp, &
     444             :                           ncol=1, &
     445             :                           source_start=proj_mo%td_mo_index(j_td), &
     446         540 :                           target_start=1)
     447             : 
     448         540 :          CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), real_proj)
     449             : 
     450             :          ! Imaginary part of the projection
     451             :          imag_proj = 0.0_dp
     452             :          CALL cp_fm_to_fm(mos_new(2*spin_td), mo_coeff_temp, &
     453             :                           ncol=1, &
     454             :                           source_start=proj_mo%td_mo_index(j_td), &
     455         540 :                           target_start=1)
     456             : 
     457         540 :          CALL cp_fm_trace(mo_coeff_temp, proj_mo%mo_ref(i_ref), imag_proj)
     458             : 
     459             :          ! Store the result
     460         540 :          phase(j_td) = ATAN2(imag_proj, real_proj) ! in radians
     461        1196 :          popu(j_td) = real_proj**2 + imag_proj**2
     462             : 
     463             :       END DO
     464             : 
     465         116 :       CALL cp_fm_struct_release(mo_ref_fmstruct)
     466         116 :       CALL cp_fm_release(mo_coeff_temp)
     467             : 
     468         116 :       CALL timestop(handle)
     469             : 
     470         116 :    END SUBROUTINE compute_proj_mo
     471             : 
     472             : ! **************************************************************************************************
     473             : !> \brief Write in one file the projection of (all) the time-dependent MO coefficients
     474             : !>        on one reference ones
     475             : !> \param qs_env ...
     476             : !> \param print_mo_section ...
     477             : !> \param proj_mo ...
     478             : !> \param i_ref ...
     479             : !> \param popu ...
     480             : !> \param phase ...
     481             : !> \param popu_tot ...
     482             : !> \param n_proj ...
     483             : !> \author Guillaume Le Breton
     484             : ! **************************************************************************************************
     485          76 :    SUBROUTINE write_proj_mo(qs_env, print_mo_section, proj_mo, i_ref, popu, phase, popu_tot, n_proj)
     486             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     487             :       TYPE(section_vals_type), POINTER                   :: print_mo_section
     488             :       TYPE(proj_mo_type)                                 :: proj_mo
     489             :       INTEGER, OPTIONAL                                  :: i_ref
     490             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: popu, phase
     491             :       REAL(KIND=dp), OPTIONAL                            :: popu_tot
     492             :       INTEGER, OPTIONAL                                  :: n_proj
     493             : 
     494             :       CHARACTER(LEN=default_string_length)               :: ext, filename
     495             :       INTEGER                                            :: j_td, output_unit, print_unit
     496             :       TYPE(cp_logger_type), POINTER                      :: logger
     497             : 
     498          76 :       NULLIFY (logger)
     499             : 
     500          76 :       logger => cp_get_default_logger()
     501          76 :       output_unit = cp_logger_get_default_io_unit(logger)
     502             : 
     503          76 :       IF (.NOT. (output_unit > 0)) RETURN
     504             : 
     505          38 :       IF (proj_mo%sum_on_all_ref) THEN
     506           4 :          ext = "-"//TRIM(ADJUSTL(cp_to_string(n_proj)))//"-ALL_REF.dat"
     507             :       ELSE
     508             :          ! Filename is update wrt the reference MO number
     509             :          ext = "-"//TRIM(ADJUSTL(cp_to_string(n_proj)))// &
     510             :                "-REF-"// &
     511             :                TRIM(ADJUSTL(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
     512          34 :                ".dat"
     513             :       END IF
     514             : 
     515             :       print_unit = cp_print_key_unit_nr(logger, print_mo_section, "", &
     516          38 :                                         extension=TRIM(ext))
     517             : 
     518          38 :       IF (print_unit /= output_unit) THEN
     519          38 :          INQUIRE (UNIT=print_unit, NAME=filename)
     520             :          WRITE (UNIT=output_unit, FMT="(/,T2,A,2(/,T3,A),/)") &
     521          38 :             "PROJECTION MO", "The projection of the TD MOs is done in the file:", &
     522          76 :             TRIM(filename)
     523             :          WRITE (UNIT=print_unit, FMT="(/,(T2,A,T40,I6))") &
     524          38 :             "Real time propagation step:", qs_env%sim_step
     525             :       ELSE
     526           0 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") "PROJECTION MO"
     527             :       END IF
     528             : 
     529          38 :       IF (proj_mo%sum_on_all_ref) THEN
     530             :          WRITE (print_unit, "(T3,A)") &
     531             :             "Projection on all the required MO number from the reference file "// &
     532           4 :             TRIM(proj_mo%ref_mo_file_name)
     533           4 :          IF (proj_mo%sum_on_all_td) THEN
     534             :             WRITE (print_unit, "(T3, A, E20.12)") &
     535           2 :                "The sum over all the TD MOs population:", popu_tot
     536             :          ELSE
     537             :             WRITE (print_unit, "(T3,A)") &
     538           2 :                "For each TD MOs required is printed: Population "
     539          14 :             DO j_td = 1, SIZE(popu)
     540          14 :                WRITE (print_unit, "(T5,1(E20.12, 1X))") popu(j_td)
     541             :             END DO
     542             :          END IF
     543             :       ELSE
     544             :          WRITE (print_unit, "(T3,A)") &
     545             :             "Projection on the MO number "// &
     546             :             TRIM(ADJUSTL(cp_to_string(proj_mo%ref_mo_index(i_ref))))// &
     547             :             " from the reference file "// &
     548          34 :             TRIM(proj_mo%ref_mo_file_name)
     549             : 
     550          34 :          IF (proj_mo%sum_on_all_td) THEN
     551             :             WRITE (print_unit, "(T3, A, E20.12)") &
     552          14 :                "The sum over all the TD MOs population:", popu_tot
     553             :          ELSE
     554             :             WRITE (print_unit, "(T3,A)") &
     555          20 :                "For each TD MOs required is printed: Population & Phase [rad] "
     556          62 :             DO j_td = 1, SIZE(popu)
     557          62 :                WRITE (print_unit, "(T5,2(E20.12, E16.8, 1X))") popu(j_td), phase(j_td)
     558             :             END DO
     559             :          END IF
     560             :       END IF
     561             : 
     562          38 :       CALL cp_print_key_finished_output(print_unit, logger, print_mo_section, "")
     563             : 
     564             :    END SUBROUTINE write_proj_mo
     565             : 
     566             : END MODULE rt_projection_mo_utils
     567             : 

Generated by: LCOV version 1.15