LCOV - code coverage report
Current view: top level - src/emd - rt_propagation_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:15a58fb) Lines: 332 346 96.0 %
Date: 2025-02-18 08:24:35 Functions: 9 9 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines needed for EMD
      10             : !> \author Florian Schiffmann (02.09)
      11             : ! **************************************************************************************************
      12             : 
      13             : MODULE rt_propagation_utils
      14             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      15             :    USE cell_types,                      ONLY: cell_type
      16             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      17             :    USE cp_control_types,                ONLY: dft_control_type,&
      18             :                                               rtp_control_type
      19             :    USE cp_dbcsr_api,                    ONLY: &
      20             :         dbcsr_add, dbcsr_binary_read, dbcsr_checksum, dbcsr_copy, dbcsr_create, &
      21             :         dbcsr_deallocate_matrix, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_filter, &
      22             :         dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      23             :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, &
      24             :         dbcsr_set, dbcsr_type
      25             :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_plus_fm_fm_t,&
      26             :                                               dbcsr_deallocate_matrix_set
      27             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale
      28             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      29             :                                               cp_fm_get_info,&
      30             :                                               cp_fm_release,&
      31             :                                               cp_fm_set_all,&
      32             :                                               cp_fm_to_fm,&
      33             :                                               cp_fm_type
      34             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      35             :                                               cp_logger_get_default_io_unit,&
      36             :                                               cp_logger_get_default_unit_nr,&
      37             :                                               cp_logger_type
      38             :    USE cp_output_handling,              ONLY: cp_p_file,&
      39             :                                               cp_print_key_finished_output,&
      40             :                                               cp_print_key_should_output,&
      41             :                                               cp_print_key_unit_nr
      42             :    USE cp_realspace_grid_cube,          ONLY: cp_pw_to_cube
      43             :    USE input_constants,                 ONLY: use_restart_wfn,&
      44             :                                               use_rt_restart
      45             :    USE input_section_types,             ONLY: section_get_ival,&
      46             :                                               section_get_ivals,&
      47             :                                               section_get_lval,&
      48             :                                               section_vals_get,&
      49             :                                               section_vals_get_subs_vals,&
      50             :                                               section_vals_type,&
      51             :                                               section_vals_val_get
      52             :    USE kinds,                           ONLY: default_path_length,&
      53             :                                               default_string_length,&
      54             :                                               dp
      55             :    USE mathconstants,                   ONLY: zero
      56             :    USE memory_utilities,                ONLY: reallocate
      57             :    USE message_passing,                 ONLY: mp_para_env_type
      58             :    USE orbital_pointers,                ONLY: ncoset
      59             :    USE particle_list_types,             ONLY: particle_list_type
      60             :    USE particle_types,                  ONLY: particle_type
      61             :    USE pw_env_types,                    ONLY: pw_env_get,&
      62             :                                               pw_env_type
      63             :    USE pw_methods,                      ONLY: pw_multiply,&
      64             :                                               pw_zero
      65             :    USE pw_pool_types,                   ONLY: pw_pool_p_type,&
      66             :                                               pw_pool_type
      67             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      68             :                                               pw_r3d_rs_type
      69             :    USE qs_collocate_density,            ONLY: calculate_wavefunction
      70             :    USE qs_density_matrices,             ONLY: calculate_density_matrix
      71             :    USE qs_dftb_matrices,                ONLY: build_dftb_overlap
      72             :    USE qs_environment_types,            ONLY: get_qs_env,&
      73             :                                               qs_environment_type
      74             :    USE qs_kind_types,                   ONLY: qs_kind_type
      75             :    USE qs_ks_types,                     ONLY: qs_ks_did_change,&
      76             :                                               qs_ks_env_type
      77             :    USE qs_mo_io,                        ONLY: read_mo_set_from_restart,&
      78             :                                               read_rt_mos_from_restart,&
      79             :                                               write_mo_set_to_output_unit
      80             :    USE qs_mo_types,                     ONLY: allocate_mo_set,&
      81             :                                               deallocate_mo_set,&
      82             :                                               get_mo_set,&
      83             :                                               init_mo_set,&
      84             :                                               mo_set_type
      85             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      86             :    USE qs_overlap,                      ONLY: build_overlap_matrix
      87             :    USE qs_rho_methods,                  ONLY: qs_rho_update_rho
      88             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      89             :                                               qs_rho_set,&
      90             :                                               qs_rho_type
      91             :    USE qs_scf_wfn_mix,                  ONLY: wfn_mix
      92             :    USE qs_subsys_types,                 ONLY: qs_subsys_get,&
      93             :                                               qs_subsys_type
      94             :    USE rt_propagation_types,            ONLY: get_rtp,&
      95             :                                               rt_prop_type
      96             : #include "../base/base_uses.f90"
      97             : 
      98             :    IMPLICIT NONE
      99             :    PRIVATE
     100             : 
     101             :    PUBLIC :: get_restart_wfn, &
     102             :              calc_S_derivs, &
     103             :              calc_update_rho, &
     104             :              calc_update_rho_sparse, &
     105             :              calculate_P_imaginary, &
     106             :              write_rtp_mos_to_output_unit, &
     107             :              write_rtp_mo_cubes, &
     108             :              warn_section_unused
     109             : 
     110             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'rt_propagation_utils'
     111             : 
     112             : CONTAINS
     113             : 
     114             : ! **************************************************************************************************
     115             : !> \brief Calculates dS/dR respectily the velocity weighted derivatves
     116             : !>        only needed for ehrenfest MD.
     117             : !>
     118             : !> \param qs_env the qs environment
     119             : !> \par History
     120             : !>      02.2009 created [Manuel Guidon]
     121             : !>      02.2014 switched to dbcsr matrices [Samuel Andermatt]
     122             : !> \author Florian Schiffmann
     123             : ! **************************************************************************************************
     124        1222 :    SUBROUTINE calc_S_derivs(qs_env)
     125             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     126             : 
     127             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'calc_S_derivs'
     128             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     129             : 
     130             :       INTEGER                                            :: col_atom, handle, i, j, m, maxder, n, &
     131             :                                                             nder, row_atom
     132             :       INTEGER, DIMENSION(6, 2)                           :: c_map_mat
     133             :       LOGICAL                                            :: return_s_derivatives
     134        1222 :       REAL(dp), DIMENSION(:, :), POINTER                 :: block_values
     135             :       TYPE(dbcsr_iterator_type)                          :: iter
     136        1222 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: C_mat, S_der, s_derivs
     137             :       TYPE(dbcsr_type), POINTER                          :: B_mat, tmp_mat, tmp_mat2
     138             :       TYPE(dft_control_type), POINTER                    :: dft_control
     139             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     140        1222 :          POINTER                                         :: sab_orb
     141        1222 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     142             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     143             :       TYPE(rt_prop_type), POINTER                        :: rtp
     144             : 
     145        1222 :       CALL timeset(routineN, handle)
     146             : 
     147        1222 :       return_s_derivatives = .TRUE.
     148             : 
     149        1222 :       NULLIFY (particle_set)
     150        1222 :       NULLIFY (rtp)
     151        1222 :       NULLIFY (s_derivs)
     152        1222 :       NULLIFY (dft_control)
     153        1222 :       NULLIFY (ks_env)
     154             : 
     155             :       CALL get_qs_env(qs_env=qs_env, &
     156             :                       rtp=rtp, &
     157             :                       particle_set=particle_set, &
     158             :                       sab_orb=sab_orb, &
     159             :                       dft_control=dft_control, &
     160        1222 :                       ks_env=ks_env)
     161             : 
     162        1222 :       CALL get_rtp(rtp=rtp, B_mat=B_mat, C_mat=C_mat, S_der=S_der)
     163             : 
     164        1222 :       nder = 2
     165        1222 :       maxder = ncoset(nder)
     166             : 
     167             :       NULLIFY (tmp_mat)
     168        1222 :       ALLOCATE (tmp_mat)
     169        1222 :       CALL dbcsr_create(tmp_mat, template=S_der(1)%matrix, matrix_type="N")
     170             : 
     171        1222 :       IF (rtp%iter < 2) THEN
     172             :          ! calculate the overlap derivative matrices
     173         342 :          IF (dft_control%qs_control%dftb) THEN
     174          84 :             CALL build_dftb_overlap(qs_env, nder, s_derivs)
     175             :          ELSE
     176             :             CALL build_overlap_matrix(ks_env, nderivative=nder, matrix_s=s_derivs, &
     177         258 :                                       basis_type_a="ORB", basis_type_b="ORB", sab_nl=sab_orb)
     178             :          END IF
     179             : 
     180             :          NULLIFY (tmp_mat2)
     181         342 :          ALLOCATE (tmp_mat2)
     182         342 :          CALL dbcsr_create(tmp_mat2, template=S_der(1)%matrix, matrix_type="S")
     183        3420 :          DO m = 1, 9
     184        3078 :             CALL dbcsr_copy(tmp_mat2, s_derivs(m + 1)%matrix)
     185        3078 :             CALL dbcsr_desymmetrize(tmp_mat2, S_der(m)%matrix)
     186        3078 :             CALL dbcsr_scale(S_der(m)%matrix, -one)
     187        3078 :             CALL dbcsr_filter(S_der(m)%matrix, rtp%filter_eps)
     188             :             !The diagonal should be zero
     189        3078 :             CALL dbcsr_iterator_start(iter, S_der(m)%matrix)
     190       15038 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
     191       11960 :                CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
     192      187460 :                IF (row_atom == col_atom) block_values = 0.0_dp
     193             :             END DO
     194        6498 :             CALL dbcsr_iterator_stop(iter)
     195             :          END DO
     196         342 :          CALL dbcsr_deallocate_matrix_set(s_derivs)
     197         342 :          CALL dbcsr_deallocate_matrix(tmp_mat2)
     198             :       END IF
     199             : 
     200             :       !calculate scalar product v(Rb)*<alpha|d/dRb beta> (B_mat), and store the first derivatives
     201             : 
     202        1222 :       CALL dbcsr_set(B_mat, zero)
     203        4888 :       DO m = 1, 3
     204        3666 :          CALL dbcsr_copy(tmp_mat, S_der(m)%matrix)
     205        3666 :          CALL dbcsr_iterator_start(iter, tmp_mat)
     206       18814 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     207       15148 :             CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
     208      220726 :             IF (row_atom == col_atom) block_values = 0.0_dp
     209      461520 :             block_values = block_values*particle_set(col_atom)%v(m)
     210             :          END DO
     211        3666 :          CALL dbcsr_iterator_stop(iter)
     212        8554 :          CALL dbcsr_add(B_mat, tmp_mat, one, one)
     213             :       END DO
     214        1222 :       CALL dbcsr_filter(B_mat, rtp%filter_eps)
     215             :       !calculate C matrix: v(Rb)*<d/dRa alpha| d/dRb beta>
     216             : 
     217        1222 :       c_map_mat = 0
     218        1222 :       n = 0
     219        4888 :       DO j = 1, 3
     220       12220 :          DO m = j, 3
     221        7332 :             n = n + 1
     222        7332 :             c_map_mat(n, 1) = j
     223        7332 :             IF (m == j) CYCLE
     224       10998 :             c_map_mat(n, 2) = m
     225             :          END DO
     226             :       END DO
     227             : 
     228        4888 :       DO i = 1, 3
     229        4888 :          CALL dbcsr_set(C_mat(i)%matrix, zero)
     230             :       END DO
     231        8554 :       DO m = 1, 6
     232        7332 :          CALL dbcsr_copy(tmp_mat, S_der(m + 3)%matrix)
     233       23218 :          DO j = 1, 2
     234       14664 :             IF (c_map_mat(m, j) == 0) CYCLE
     235       21996 :             CALL dbcsr_add(C_mat(c_map_mat(m, j))%matrix, tmp_mat, one, one)
     236             :          END DO
     237             :       END DO
     238             : 
     239        4888 :       DO m = 1, 3
     240        3666 :          CALL dbcsr_iterator_start(iter, C_mat(m)%matrix)
     241       17746 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     242       14080 :             CALL dbcsr_iterator_next_block(iter, row_atom, col_atom, block_values)
     243      454443 :             block_values = block_values*particle_set(row_atom)%v(m)
     244             :          END DO
     245        3666 :          CALL dbcsr_iterator_stop(iter)
     246        8554 :          CALL dbcsr_filter(C_mat(m)%matrix, rtp%filter_eps)
     247             :       END DO
     248             : 
     249        1222 :       CALL dbcsr_deallocate_matrix(tmp_mat)
     250        1222 :       CALL timestop(handle)
     251        1222 :    END SUBROUTINE
     252             : 
     253             : ! **************************************************************************************************
     254             : !> \brief reads the restart file. At the moment only SCF (means only real)
     255             : !> \param qs_env ...
     256             : !> \author Florian Schiffmann (02.09)
     257             : ! **************************************************************************************************
     258             : 
     259          36 :    SUBROUTINE get_restart_wfn(qs_env)
     260             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     261             : 
     262             :       CHARACTER(LEN=default_path_length)                 :: file_name, project_name
     263             :       INTEGER                                            :: i, id_nr, im, ispin, ncol, nspin, &
     264             :                                                             output_unit, re, unit_nr
     265             :       REAL(KIND=dp)                                      :: alpha, cs_pos
     266          36 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     267             :       TYPE(cp_fm_type)                                   :: mos_occ
     268          36 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_old
     269             :       TYPE(cp_logger_type), POINTER                      :: logger
     270             :       TYPE(dbcsr_distribution_type)                      :: dist
     271          36 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: p_rmpv, rho_new, rho_old
     272             :       TYPE(dft_control_type), POINTER                    :: dft_control
     273          36 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     274             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     275          36 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     276          36 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     277             :       TYPE(qs_rho_type), POINTER                         :: rho_struct
     278             :       TYPE(rt_prop_type), POINTER                        :: rtp
     279             :       TYPE(section_vals_type), POINTER                   :: dft_section, input
     280             : 
     281          36 :       NULLIFY (atomic_kind_set, qs_kind_set, mo_array, particle_set, rho_struct, para_env)
     282             : 
     283             :       CALL get_qs_env(qs_env, &
     284             :                       qs_kind_set=qs_kind_set, &
     285             :                       atomic_kind_set=atomic_kind_set, &
     286             :                       particle_set=particle_set, &
     287             :                       mos=mo_array, &
     288             :                       input=input, &
     289             :                       rtp=rtp, &
     290             :                       dft_control=dft_control, &
     291             :                       rho=rho_struct, &
     292          36 :                       para_env=para_env)
     293          36 :       logger => cp_get_default_logger()
     294          36 :       output_unit = cp_logger_get_default_io_unit(logger)
     295             : 
     296          36 :       IF (logger%para_env%is_source()) THEN
     297          18 :          unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     298             :       ELSE
     299             :          unit_nr = -1
     300             :       END IF
     301             : 
     302          36 :       id_nr = 0
     303          36 :       nspin = SIZE(mo_array)
     304          36 :       CALL qs_rho_get(rho_struct, rho_ao=p_rmpv)
     305          36 :       dft_section => section_vals_get_subs_vals(input, "DFT")
     306          62 :       SELECT CASE (dft_control%rtp_control%initial_wfn)
     307             :       CASE (use_restart_wfn)
     308             :          CALL read_mo_set_from_restart(mo_array, atomic_kind_set, qs_kind_set, particle_set, para_env, &
     309          26 :                                        id_nr=id_nr, multiplicity=dft_control%multiplicity, dft_section=dft_section)
     310          26 :          CALL set_uniform_occupation_mo_array(mo_array, nspin)
     311             : 
     312          26 :          IF (dft_control%rtp_control%apply_wfn_mix_init_restart) &
     313             :             CALL wfn_mix(mo_array, particle_set, dft_section, qs_kind_set, para_env, output_unit, &
     314           4 :                          for_rtp=.TRUE.)
     315             : 
     316          70 :          DO ispin = 1, nspin
     317          70 :             CALL calculate_density_matrix(mo_array(ispin), p_rmpv(ispin)%matrix)
     318             :          END DO
     319          26 :          IF (rtp%linear_scaling) THEN
     320          14 :             CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
     321          34 :             DO ispin = 1, nspin
     322          20 :                re = 2*ispin - 1
     323          20 :                im = 2*ispin
     324          20 :                CALL cp_fm_get_info(mo_array(ispin)%mo_coeff, ncol_global=ncol)
     325             :                CALL cp_fm_create(mos_occ, &
     326             :                                  matrix_struct=mo_array(ispin)%mo_coeff%matrix_struct, &
     327          20 :                                  name="mos_occ")
     328          20 :                CALL cp_fm_to_fm(mo_array(ispin)%mo_coeff, mos_occ)
     329          20 :                IF (mo_array(ispin)%uniform_occupation) THEN
     330          16 :                   alpha = 3.0_dp - REAL(nspin, dp)
     331          78 :                   CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
     332             :                   CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     333             :                                              matrix_v=mos_occ, &
     334             :                                              ncol=ncol, &
     335          16 :                                              alpha=alpha, keep_sparsity=.FALSE.)
     336             :                ELSE
     337           4 :                   alpha = 1.0_dp
     338          88 :                   CALL cp_fm_column_scale(mos_occ, mo_array(ispin)%occupation_numbers/alpha)
     339             :                   CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_old(re)%matrix, &
     340             :                                              matrix_v=mo_array(ispin)%mo_coeff, &
     341             :                                              matrix_g=mos_occ, &
     342             :                                              ncol=ncol, &
     343           4 :                                              alpha=alpha, keep_sparsity=.FALSE.)
     344             :                END IF
     345          20 :                CALL dbcsr_filter(rho_old(re)%matrix, rtp%filter_eps)
     346          20 :                CALL dbcsr_copy(rho_new(re)%matrix, rho_old(re)%matrix)
     347          54 :                CALL cp_fm_release(mos_occ)
     348             :             END DO
     349          14 :             CALL calc_update_rho_sparse(qs_env)
     350             :          ELSE
     351          12 :             CALL get_rtp(rtp=rtp, mos_old=mos_old)
     352          36 :             DO i = 1, SIZE(qs_env%mos)
     353          24 :                CALL cp_fm_to_fm(mo_array(i)%mo_coeff, mos_old(2*i - 1))
     354          36 :                CALL cp_fm_set_all(mos_old(2*i), zero, zero)
     355             :             END DO
     356             :          END IF
     357             :       CASE (use_rt_restart)
     358          36 :          IF (rtp%linear_scaling) THEN
     359           2 :             CALL get_rtp(rtp=rtp, rho_old=rho_old, rho_new=rho_new)
     360           2 :             project_name = logger%iter_info%project_name
     361           4 :             DO ispin = 1, nspin
     362           2 :                re = 2*ispin - 1
     363           2 :                im = 2*ispin
     364           2 :                WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_RE", ispin, "_RESTART.dm"
     365           2 :                CALL dbcsr_get_info(rho_old(re)%matrix, distribution=dist)
     366           2 :                CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=rho_old(re)%matrix)
     367           2 :                cs_pos = dbcsr_checksum(rho_old(re)%matrix, pos=.TRUE.)
     368           2 :                IF (unit_nr > 0) THEN
     369           1 :                   WRITE (unit_nr, '(T2,A,E20.8)') "Read restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     370             :                END IF
     371           2 :                WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_IM", ispin, "_RESTART.dm"
     372           2 :                CALL dbcsr_get_info(rho_old(im)%matrix, distribution=dist)
     373           2 :                CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=rho_old(im)%matrix)
     374           2 :                cs_pos = dbcsr_checksum(rho_old(im)%matrix, pos=.TRUE.)
     375           8 :                IF (unit_nr > 0) THEN
     376           1 :                   WRITE (unit_nr, '(T2,A,E20.8)') "Read restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
     377             :                END IF
     378             :             END DO
     379           6 :             DO i = 1, SIZE(rho_new)
     380           6 :                CALL dbcsr_copy(rho_new(i)%matrix, rho_old(i)%matrix)
     381             :             END DO
     382           2 :             CALL calc_update_rho_sparse(qs_env)
     383             :          ELSE
     384           8 :             CALL get_rtp(rtp=rtp, mos_old=mos_old)
     385             :             CALL read_rt_mos_from_restart(mo_array, mos_old, atomic_kind_set, qs_kind_set, particle_set, para_env, &
     386           8 :                                           id_nr, dft_control%multiplicity, dft_section)
     387           8 :             CALL set_uniform_occupation_mo_array(mo_array, nspin)
     388          16 :             DO ispin = 1, nspin
     389             :                CALL calculate_density_matrix(mo_array(ispin), &
     390          16 :                                              p_rmpv(ispin)%matrix)
     391             :             END DO
     392             :          END IF
     393             :       END SELECT
     394             : 
     395          36 :    END SUBROUTINE get_restart_wfn
     396             : 
     397             : ! **************************************************************************************************
     398             : !> \brief Set mo_array(ispin)%uniform_occupation after a restart
     399             : !> \param mo_array ...
     400             : !> \param nspin ...
     401             : !> \author Guillaume Le Breton (03.23)
     402             : ! **************************************************************************************************
     403             : 
     404          34 :    SUBROUTINE set_uniform_occupation_mo_array(mo_array, nspin)
     405             : 
     406             :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mo_array
     407             :       INTEGER                                            :: nspin
     408             : 
     409             :       INTEGER                                            :: ispin, mo
     410             :       LOGICAL                                            :: is_uniform
     411             : 
     412          86 :       DO ispin = 1, nspin
     413          52 :          is_uniform = .TRUE.
     414         274 :          DO mo = 1, mo_array(ispin)%nmo
     415             :             IF (mo_array(ispin)%occupation_numbers(mo) /= 0.0 .AND. &
     416         222 :                 mo_array(ispin)%occupation_numbers(mo) /= 1.0 .AND. &
     417             :                 mo_array(ispin)%occupation_numbers(mo) /= 2.0) &
     418         100 :                is_uniform = .FALSE.
     419             :          END DO
     420          86 :          mo_array(ispin)%uniform_occupation = is_uniform
     421             :       END DO
     422             : 
     423          34 :    END SUBROUTINE set_uniform_occupation_mo_array
     424             : 
     425             : ! **************************************************************************************************
     426             : !> \brief calculates the density from the complex MOs and passes the density to
     427             : !>        qs_env.
     428             : !> \param qs_env ...
     429             : !> \author Florian Schiffmann (02.09)
     430             : ! **************************************************************************************************
     431             : 
     432        2018 :    SUBROUTINE calc_update_rho(qs_env)
     433             : 
     434             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     435             : 
     436             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calc_update_rho'
     437             :       REAL(KIND=dp), PARAMETER                           :: one = 1.0_dp, zero = 0.0_dp
     438             : 
     439             :       INTEGER                                            :: handle, i, im, ncol, re
     440             :       REAL(KIND=dp)                                      :: alpha
     441             :       TYPE(cp_fm_type)                                   :: mos_occ
     442        2018 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     443        2018 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao, rho_ao_im
     444        2018 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     445             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     446             :       TYPE(qs_rho_type), POINTER                         :: rho
     447             :       TYPE(rt_prop_type), POINTER                        :: rtp
     448             : 
     449        2018 :       CALL timeset(routineN, handle)
     450             : 
     451        2018 :       NULLIFY (rho, ks_env, mos_new, rtp)
     452             :       CALL get_qs_env(qs_env, &
     453             :                       ks_env=ks_env, &
     454             :                       rho=rho, &
     455             :                       rtp=rtp, &
     456        2018 :                       mos=mos)
     457        2018 :       CALL get_rtp(rtp=rtp, mos_new=mos_new)
     458        2018 :       CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
     459        4580 :       DO i = 1, SIZE(mos_new)/2
     460        2562 :          re = 2*i - 1; im = 2*i
     461        2562 :          CALL dbcsr_set(rho_ao(i)%matrix, zero)
     462        2562 :          CALL cp_fm_get_info(mos_new(re), ncol_global=ncol)
     463             :          CALL cp_fm_create(mos_occ, &
     464             :                            matrix_struct=mos(i)%mo_coeff%matrix_struct, &
     465        2562 :                            name="mos_occ")
     466        2562 :          CALL cp_fm_to_fm(mos_new(re), mos_occ)
     467        2562 :          IF (mos(i)%uniform_occupation) THEN
     468        2462 :             alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
     469       10286 :             CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
     470             :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
     471             :                                        matrix_v=mos_occ, &
     472             :                                        ncol=ncol, &
     473        2462 :                                        alpha=alpha)
     474             :          ELSE
     475         100 :             alpha = 1.0_dp
     476         660 :             CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
     477             :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
     478             :                                        matrix_v=mos_new(re), &
     479             :                                        matrix_g=mos_occ, &
     480             :                                        ncol=ncol, &
     481         100 :                                        alpha=alpha)
     482             :          END IF
     483             : 
     484             :          ! It is actually complex conjugate but i*i=-1 therefore it must be added
     485        2562 :          CALL cp_fm_to_fm(mos_new(im), mos_occ)
     486        2562 :          IF (mos(i)%uniform_occupation) THEN
     487        2462 :             alpha = 3*one - REAL(SIZE(mos_new)/2, dp)
     488       10286 :             CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
     489             :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
     490             :                                        matrix_v=mos_occ, &
     491             :                                        ncol=ncol, &
     492        2462 :                                        alpha=alpha)
     493             :          ELSE
     494         100 :             alpha = 1.0_dp
     495         660 :             CALL cp_fm_column_scale(mos_occ, mos(i)%occupation_numbers/alpha)
     496             :             CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(i)%matrix, &
     497             :                                        matrix_v=mos_new(im), &
     498             :                                        matrix_g=mos_occ, &
     499             :                                        ncol=ncol, &
     500         100 :                                        alpha=alpha)
     501             :          END IF
     502        7142 :          CALL cp_fm_release(mos_occ)
     503             :       END DO
     504        2018 :       CALL qs_rho_update_rho(rho, qs_env)
     505             : 
     506        2018 :       IF (rtp%track_imag_density) THEN
     507        1358 :          CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
     508        1358 :          CALL calculate_P_imaginary(qs_env, rtp, rho_ao_im, keep_sparsity=.TRUE.)
     509        1358 :          CALL qs_rho_set(rho, rho_ao_im=rho_ao_im)
     510             :       END IF
     511             : 
     512        2018 :       CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
     513             : 
     514        2018 :       CALL timestop(handle)
     515             : 
     516        2018 :    END SUBROUTINE calc_update_rho
     517             : 
     518             : ! **************************************************************************************************
     519             : !> \brief Copies the density matrix back into the qs_env%rho%rho_ao
     520             : !> \param qs_env ...
     521             : !> \author Samuel Andermatt (3.14)
     522             : ! **************************************************************************************************
     523             : 
     524        1264 :    SUBROUTINE calc_update_rho_sparse(qs_env)
     525             : 
     526             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     527             : 
     528             :       CHARACTER(len=*), PARAMETER :: routineN = 'calc_update_rho_sparse'
     529             :       REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
     530             : 
     531             :       INTEGER                                            :: handle, ispin
     532        1264 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao, rho_ao_im, rho_new
     533             :       TYPE(dft_control_type), POINTER                    :: dft_control
     534             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     535             :       TYPE(qs_rho_type), POINTER                         :: rho
     536             :       TYPE(rt_prop_type), POINTER                        :: rtp
     537             :       TYPE(rtp_control_type), POINTER                    :: rtp_control
     538             : 
     539        1264 :       NULLIFY (rho, ks_env, rtp, dft_control)
     540        1264 :       CALL timeset(routineN, handle)
     541             :       CALL get_qs_env(qs_env, &
     542             :                       ks_env=ks_env, &
     543             :                       rho=rho, &
     544             :                       rtp=rtp, &
     545        1264 :                       dft_control=dft_control)
     546        1264 :       rtp_control => dft_control%rtp_control
     547        1264 :       CALL get_rtp(rtp=rtp, rho_new=rho_new)
     548        1264 :       CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
     549        1264 :       IF (rtp%track_imag_density) CALL qs_rho_get(rho_struct=rho, rho_ao_im=rho_ao_im)
     550        2946 :       DO ispin = 1, SIZE(rho_ao)
     551        1682 :          CALL dbcsr_set(rho_ao(ispin)%matrix, zero)
     552        1682 :          CALL dbcsr_copy(rho_ao(ispin)%matrix, rho_new(ispin*2 - 1)%matrix, keep_sparsity=.TRUE.)
     553        2946 :          IF (rtp%track_imag_density) THEN
     554         488 :             CALL dbcsr_copy(rho_ao_im(ispin)%matrix, rho_new(ispin*2)%matrix, keep_sparsity=.TRUE.)
     555             :          END IF
     556             :       END DO
     557             : 
     558        1264 :       CALL qs_rho_update_rho(rho, qs_env)
     559        1264 :       CALL qs_ks_did_change(ks_env, rho_changed=.TRUE.)
     560             : 
     561        1264 :       CALL timestop(handle)
     562             : 
     563        1264 :    END SUBROUTINE calc_update_rho_sparse
     564             : 
     565             : ! **************************************************************************************************
     566             : !> \brief ...
     567             : !> \param qs_env ...
     568             : !> \param rtp ...
     569             : !> \param matrix_p_im ...
     570             : !> \param keep_sparsity ...
     571             : ! **************************************************************************************************
     572        1370 :    SUBROUTINE calculate_P_imaginary(qs_env, rtp, matrix_p_im, keep_sparsity)
     573             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     574             :       TYPE(rt_prop_type), POINTER                        :: rtp
     575             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p_im
     576             :       LOGICAL, OPTIONAL                                  :: keep_sparsity
     577             : 
     578             :       INTEGER                                            :: i, im, ncol, re
     579             :       LOGICAL                                            :: my_keep_sparsity
     580             :       REAL(KIND=dp)                                      :: alpha
     581        1370 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new, mos_occ
     582        1370 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     583             : 
     584        1370 :       CALL get_rtp(rtp=rtp, mos_new=mos_new)
     585             : 
     586        1370 :       my_keep_sparsity = .FALSE.
     587        1370 :       IF (PRESENT(keep_sparsity)) my_keep_sparsity = keep_sparsity
     588        1370 :       CALL get_qs_env(qs_env, mos=mos)
     589        7482 :       ALLOCATE (mos_occ(SIZE(mos)*2))
     590             : 
     591        3056 :       DO i = 1, SIZE(mos_new)/2
     592        1686 :          re = 2*i - 1; im = 2*i
     593        1686 :          alpha = 3.0_dp - REAL(SIZE(matrix_p_im), dp)
     594             :          CALL cp_fm_create(mos_occ(re), &
     595             :                            matrix_struct=mos(i)%mo_coeff%matrix_struct, &
     596        1686 :                            name="mos_occ")
     597             :          CALL cp_fm_create(mos_occ(im), &
     598             :                            matrix_struct=mos(i)%mo_coeff%matrix_struct, &
     599        1686 :                            name="mos_occ")
     600        1686 :          CALL dbcsr_set(matrix_p_im(i)%matrix, 0.0_dp)
     601        1686 :          CALL cp_fm_get_info(mos_new(re), ncol_global=ncol)
     602        1686 :          CALL cp_fm_to_fm(mos_new(re), mos_occ(re))
     603        7152 :          CALL cp_fm_column_scale(mos_occ(re), mos(i)%occupation_numbers/alpha)
     604        1686 :          CALL cp_fm_to_fm(mos_new(im), mos_occ(im))
     605        7152 :          CALL cp_fm_column_scale(mos_occ(im), mos(i)%occupation_numbers/alpha)
     606             :          CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=matrix_p_im(i)%matrix, &
     607             :                                     matrix_v=mos_occ(im), &
     608             :                                     matrix_g=mos_occ(re), &
     609             :                                     ncol=ncol, &
     610             :                                     keep_sparsity=my_keep_sparsity, &
     611             :                                     alpha=2.0_dp*alpha, &
     612        4742 :                                     symmetry_mode=-1)
     613             :       END DO
     614        1370 :       CALL cp_fm_release(mos_occ)
     615             : 
     616        1370 :    END SUBROUTINE calculate_P_imaginary
     617             : 
     618             : ! **************************************************************************************************
     619             : !> \brief ...
     620             : !> \param qs_env ...
     621             : !> \param rtp ...
     622             : ! **************************************************************************************************
     623         624 :    SUBROUTINE write_rtp_mos_to_output_unit(qs_env, rtp)
     624             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     625             :       TYPE(rt_prop_type), POINTER                        :: rtp
     626             : 
     627             :       CHARACTER(len=*), PARAMETER :: routineN = 'write_rtp_mos_to_output_unit'
     628             : 
     629             :       CHARACTER(LEN=10)                                  :: spin
     630             :       CHARACTER(LEN=2*default_string_length)             :: name
     631             :       INTEGER                                            :: handle, i, ispin, nao, nelectron, nmo, &
     632             :                                                             nspins
     633             :       LOGICAL                                            :: print_eigvecs, print_mo_info
     634             :       REAL(KIND=dp)                                      :: flexible_electron_count, maxocc, n_el_f
     635         312 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     636         312 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     637             :       TYPE(cp_logger_type), POINTER                      :: logger
     638             :       TYPE(mo_set_type)                                  :: mo_set_rtp
     639         312 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     640         312 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     641         312 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     642             :       TYPE(section_vals_type), POINTER                   :: dft_section, input
     643             : 
     644         312 :       NULLIFY (atomic_kind_set, particle_set, qs_kind_set, input, mos, dft_section)
     645             : 
     646         312 :       CALL timeset(routineN, handle)
     647             : 
     648             :       CALL get_qs_env(qs_env, &
     649             :                       atomic_kind_set=atomic_kind_set, &
     650             :                       qs_kind_set=qs_kind_set, &
     651             :                       particle_set=particle_set, &
     652             :                       input=input, &
     653         312 :                       mos=mos)
     654             :       ! Quick return, if no printout of MO information is requested
     655         312 :       dft_section => section_vals_get_subs_vals(input, "DFT")
     656         312 :       CALL section_vals_val_get(dft_section, "PRINT%MO%EIGENVECTORS", l_val=print_eigvecs)
     657             : 
     658         312 :       NULLIFY (logger)
     659         312 :       logger => cp_get_default_logger()
     660             :       print_mo_info = (cp_print_key_should_output(logger%iter_info, &
     661             :                                                   dft_section, "PRINT%MO") /= 0) .OR. &
     662         312 :                       (qs_env%sim_step == 1)
     663             : 
     664         100 :       IF ((.NOT. print_mo_info) .OR. (.NOT. print_eigvecs)) THEN
     665         308 :          CALL timestop(handle)
     666         308 :          RETURN
     667             :       END IF
     668             : 
     669           4 :       CALL get_rtp(rtp=rtp, mos_new=mos_new)
     670             : 
     671           4 :       nspins = SIZE(mos_new)/2
     672             : 
     673           8 :       DO ispin = 1, nspins
     674             :          ! initiate mo_set
     675             :          CALL get_mo_set(mo_set=mos(ispin), nao=nao, nmo=nmo, nelectron=nelectron, &
     676           4 :                          n_el_f=n_el_f, maxocc=maxocc, flexible_electron_count=flexible_electron_count)
     677             : 
     678             :          CALL allocate_mo_set(mo_set_rtp, &
     679             :                               nao=nao, &
     680             :                               nmo=nmo, &
     681             :                               nelectron=nelectron, &
     682             :                               n_el_f=n_el_f, &
     683             :                               maxocc=maxocc, &
     684           4 :                               flexible_electron_count=flexible_electron_count)
     685             : 
     686           4 :          WRITE (name, FMT="(A,I6)") "RTP MO SET, SPIN ", ispin
     687           4 :          CALL init_mo_set(mo_set_rtp, fm_ref=mos_new(2*ispin - 1), name=name)
     688             : 
     689           4 :          IF (nspins > 1) THEN
     690           0 :             IF (ispin == 1) THEN
     691           0 :                spin = "ALPHA SPIN"
     692             :             ELSE
     693           0 :                spin = "BETA SPIN"
     694             :             END IF
     695             :          ELSE
     696           4 :             spin = ""
     697             :          END IF
     698             : 
     699           8 :          mo_set_rtp%occupation_numbers = mos(ispin)%occupation_numbers
     700             : 
     701             :          !loop for real (odd) and imaginary (even) parts
     702          12 :          DO i = 1, 0, -1
     703           8 :             CALL cp_fm_to_fm(mos_new(2*ispin - i), mo_set_rtp%mo_coeff)
     704             :             CALL write_mo_set_to_output_unit(mo_set_rtp, atomic_kind_set, qs_kind_set, particle_set, &
     705          12 :                                           dft_section, 4, 0, rtp=.TRUE., spin=TRIM(spin), cpart=MOD(i, 2), sim_step=qs_env%sim_step)
     706             :          END DO
     707             : 
     708          12 :          CALL deallocate_mo_set(mo_set_rtp)
     709             :       END DO
     710             : 
     711           4 :       CALL timestop(handle)
     712             : 
     713         312 :    END SUBROUTINE write_rtp_mos_to_output_unit
     714             : 
     715             : ! **************************************************************************************************
     716             : !> \brief Write the time dependent amplitude of the MOs in real grid.
     717             : !>        Very close to qs_scf_post_gpw/qs_scf_post_occ_cubes subroutine.
     718             : !> \param qs_env ...
     719             : !> \param rtp ...
     720             : !> \author Guillaume Le Breton (11.22)
     721             : ! **************************************************************************************************
     722         312 :    SUBROUTINE write_rtp_mo_cubes(qs_env, rtp)
     723             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     724             :       TYPE(rt_prop_type), POINTER                        :: rtp
     725             : 
     726             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_rtp_mo_cubes'
     727             : 
     728             :       CHARACTER(LEN=default_path_length)                 :: filename, my_pos_cube, title
     729             :       INTEGER                                            :: handle, homo, i, ir, ispin, ivector, &
     730             :                                                             n_rep, nhomo, nlist, nspins, &
     731             :                                                             rt_time_step, unit_nr
     732         312 :       INTEGER, DIMENSION(:), POINTER                     :: list, list_index
     733             :       LOGICAL                                            :: append_cube, do_kpoints, mpi_io
     734         312 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     735             :       TYPE(cell_type), POINTER                           :: cell
     736             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     737         312 :       TYPE(cp_fm_type), DIMENSION(:), POINTER            :: mos_new
     738             :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     739             :       TYPE(cp_logger_type), POINTER                      :: logger
     740             :       TYPE(dft_control_type), POINTER                    :: dft_control
     741         312 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     742             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     743             :       TYPE(particle_list_type), POINTER                  :: particles
     744         312 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     745             :       TYPE(pw_c1d_gs_type)                               :: wf_g
     746             :       TYPE(pw_env_type), POINTER                         :: pw_env
     747         312 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     748             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     749             :       TYPE(pw_r3d_rs_type)                               :: density_r, wf_r
     750         312 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     751             :       TYPE(qs_subsys_type), POINTER                      :: subsys
     752             :       TYPE(section_vals_type), POINTER                   :: dft_section, input
     753             : 
     754         312 :       CALL timeset(routineN, handle)
     755             : 
     756         312 :       NULLIFY (logger, auxbas_pw_pool, pw_pools, pw_env)
     757             : 
     758             :       ! Get all the info from qs:
     759             :       CALL get_qs_env(qs_env, do_kpoints=do_kpoints, &
     760         312 :                       input=input)
     761             : 
     762             :       ! Kill the run in the case of K points
     763         312 :       IF (do_kpoints) THEN
     764           0 :          CPABORT("K points not handled yet for printing MO_CUBE")
     765             :       END IF
     766             : 
     767         312 :       dft_section => section_vals_get_subs_vals(input, "DFT")
     768         312 :       logger => cp_get_default_logger()
     769             : 
     770             :       ! Quick return if no print required
     771         312 :       IF (.NOT. BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
     772             :                                                  "PRINT%MO_CUBES"), cp_p_file)) THEN
     773         288 :          CALL timestop(handle)
     774         288 :          RETURN
     775             :       END IF
     776             : 
     777             :       CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, &
     778             :                       mos=mos, &
     779             :                       blacs_env=blacs_env, &
     780             :                       qs_kind_set=qs_kind_set, &
     781             :                       pw_env=pw_env, &
     782             :                       subsys=subsys, &
     783             :                       para_env=para_env, &
     784             :                       particle_set=particle_set, &
     785          24 :                       dft_control=dft_control)
     786          24 :       CALL qs_subsys_get(subsys, particles=particles)
     787             : 
     788          24 :       nspins = dft_control%nspins
     789          24 :       rt_time_step = qs_env%sim_step
     790             : 
     791             :       ! Setup the grids needed to compute a wavefunction given a vector
     792             :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
     793          24 :                       pw_pools=pw_pools)
     794          24 :       CALL auxbas_pw_pool%create_pw(wf_r)
     795          24 :       CALL auxbas_pw_pool%create_pw(wf_g)
     796          24 :       CALL auxbas_pw_pool%create_pw(density_r)
     797          24 :       CALL get_rtp(rtp=rtp, mos_new=mos_new)
     798             : 
     799          70 :       DO ispin = 1, nspins
     800          46 :          CALL get_mo_set(mo_set=mos(ispin), homo=homo)
     801             : 
     802          46 :          nhomo = section_get_ival(dft_section, "PRINT%MO_CUBES%NHOMO")
     803          46 :          append_cube = section_get_lval(dft_section, "PRINT%MO_CUBES%APPEND")
     804          46 :          my_pos_cube = "REWIND"
     805          46 :          IF (append_cube) THEN
     806           0 :             my_pos_cube = "APPEND"
     807             :          END IF
     808          46 :          CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", n_rep_val=n_rep)
     809          46 :          IF (n_rep > 0) THEN ! write the cubes of the list
     810           0 :             nlist = 0
     811           0 :             DO ir = 1, n_rep
     812           0 :                NULLIFY (list)
     813             :                CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%HOMO_LIST", i_rep_val=ir, &
     814           0 :                                          i_vals=list)
     815           0 :                IF (ASSOCIATED(list)) THEN
     816           0 :                   CALL reallocate(list_index, 1, nlist + SIZE(list))
     817           0 :                   DO i = 1, SIZE(list)
     818           0 :                      list_index(i + nlist) = list(i)
     819             :                   END DO
     820           0 :                   nlist = nlist + SIZE(list)
     821             :                END IF
     822             :             END DO
     823             :          ELSE
     824             : 
     825          46 :             IF (nhomo == -1) nhomo = homo
     826          46 :             nlist = homo - MAX(1, homo - nhomo + 1) + 1
     827         138 :             ALLOCATE (list_index(nlist))
     828         224 :             DO i = 1, nlist
     829         224 :                list_index(i) = MAX(1, homo - nhomo + 1) + i - 1
     830             :             END DO
     831             :          END IF
     832         224 :          DO i = 1, nlist
     833         178 :             ivector = list_index(i)
     834             :             CALL get_qs_env(qs_env=qs_env, &
     835             :                             atomic_kind_set=atomic_kind_set, &
     836             :                             qs_kind_set=qs_kind_set, &
     837             :                             cell=cell, &
     838             :                             particle_set=particle_set, &
     839         178 :                             pw_env=pw_env)
     840             : 
     841             :             ! density_r contains the density of the MOs
     842         178 :             CALL pw_zero(density_r)
     843         178 :             mo_coeff => mos_new(2*ispin - 1)!Real coeff
     844             :             CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
     845         178 :                                         cell, dft_control, particle_set, pw_env)
     846             :             ! Adding the real part
     847         178 :             CALL pw_multiply(density_r, wf_r, wf_r, 1.0_dp)
     848             : 
     849         178 :             mo_coeff => mos_new(2*ispin) !Im coeff
     850             :             CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
     851         178 :                                         cell, dft_control, particle_set, pw_env)
     852             :             ! Adding the im part
     853         178 :             CALL pw_multiply(density_r, wf_r, wf_r, 1.0_dp)
     854             : 
     855         178 :             WRITE (filename, '(a4,I5.5,a1,I1.1)') "WFN_", ivector, "_", ispin
     856         178 :             mpi_io = .TRUE.
     857             :             unit_nr = cp_print_key_unit_nr(logger, input, "DFT%PRINT%MO_CUBES", extension=".cube", &
     858             :                                            middle_name=TRIM(filename), file_position=my_pos_cube, log_filename=.FALSE., &
     859         178 :                                            mpi_io=mpi_io)
     860         178 :             WRITE (title, *) "DENSITY ", ivector, " spin ", ispin, " i.e. HOMO - ", ivector - homo
     861             :             CALL cp_pw_to_cube(density_r, unit_nr, title, particles=particles, &
     862         178 :                                stride=section_get_ivals(dft_section, "PRINT%MO_CUBES%STRIDE"), mpi_io=mpi_io)
     863         224 :             CALL cp_print_key_finished_output(unit_nr, logger, input, "DFT%PRINT%MO_CUBES", mpi_io=mpi_io)
     864             :          END DO
     865         162 :          IF (ASSOCIATED(list_index)) DEALLOCATE (list_index)
     866             :       END DO
     867             : 
     868             :       ! Deallocate grids needed to compute wavefunctions
     869          24 :       CALL auxbas_pw_pool%give_back_pw(wf_r)
     870          24 :       CALL auxbas_pw_pool%give_back_pw(wf_g)
     871          24 :       CALL auxbas_pw_pool%give_back_pw(density_r)
     872             : 
     873          24 :       CALL timestop(handle)
     874             : 
     875         312 :    END SUBROUTINE write_rtp_mo_cubes
     876             : 
     877             : ! **************************************************************************************************
     878             : !> \brief Warn about unused sections of the print section - only implemented for some of the methods
     879             : !> \param section Parent section
     880             : !> \param subsection_name Name of the subsection which is not used even when explicitly included
     881             : !> \param error_message Message to display as warning
     882             : !> \author Stepan Marek (01.25)
     883             : ! **************************************************************************************************
     884        1704 :    SUBROUTINE warn_section_unused(section, subsection_name, error_message)
     885             :       TYPE(section_vals_type), POINTER                   :: section
     886             :       CHARACTER(len=*)                                   :: subsection_name, error_message
     887             : 
     888             :       LOGICAL                                            :: explicit
     889             :       TYPE(section_vals_type), POINTER                   :: target_section
     890             : 
     891         852 :       target_section => section_vals_get_subs_vals(section, subsection_name)
     892         852 :       CALL section_vals_get(target_section, explicit=explicit)
     893         852 :       IF (explicit) CPWARN(error_message)
     894         852 :    END SUBROUTINE
     895             : 
     896             : END MODULE rt_propagation_utils

Generated by: LCOV version 1.15