LCOV - code coverage report
Current view: top level - src - mixed_cdft_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 931 1038 89.7 %
Date: 2024-03-29 07:50:05 Functions: 12 12 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Utility subroutines for mixed CDFT calculations
      10             : !> \par   History
      11             : !>                 separated from mixed_cdft_methods [01.2017]
      12             : !> \author Nico Holmberg [01.2017]
      13             : ! **************************************************************************************************
      14             : MODULE mixed_cdft_utils
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      16             :    USE cell_types,                      ONLY: cell_type
      17             :    USE cp_array_utils,                  ONLY: cp_1d_i_p_type,&
      18             :                                               cp_1d_r_p_type,&
      19             :                                               cp_2d_r_p_type
      20             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
      21             :                                               cp_blacs_env_type
      22             :    USE cp_control_types,                ONLY: dft_control_type
      23             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      24             :                                               copy_fm_to_dbcsr_bc
      25             :    USE cp_files,                        ONLY: open_file
      26             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      27             :                                               cp_fm_struct_release,&
      28             :                                               cp_fm_struct_type
      29             :    USE cp_fm_types,                     ONLY: cp_fm_copy_general,&
      30             :                                               cp_fm_create,&
      31             :                                               cp_fm_get_info,&
      32             :                                               cp_fm_release,&
      33             :                                               cp_fm_to_fm,&
      34             :                                               cp_fm_type
      35             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      36             :                                               cp_logger_create,&
      37             :                                               cp_logger_set,&
      38             :                                               cp_logger_type,&
      39             :                                               cp_to_string
      40             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      41             :                                               cp_print_key_unit_nr
      42             :    USE cp_realspace_grid_init,          ONLY: init_input_type
      43             :    USE cp_subsys_types,                 ONLY: cp_subsys_get,&
      44             :                                               cp_subsys_type
      45             :    USE cube_utils,                      ONLY: init_cube_info,&
      46             :                                               return_cube_max_iradius
      47             :    USE d3_poly,                         ONLY: init_d3_poly_module
      48             :    USE dbcsr_api,                       ONLY: dbcsr_desymmetrize,&
      49             :                                               dbcsr_get_info,&
      50             :                                               dbcsr_init_p,&
      51             :                                               dbcsr_p_type,&
      52             :                                               dbcsr_release,&
      53             :                                               dbcsr_release_p,&
      54             :                                               dbcsr_type
      55             :    USE force_env_types,                 ONLY: force_env_get,&
      56             :                                               force_env_type,&
      57             :                                               multiple_fe_list
      58             :    USE gaussian_gridlevels,             ONLY: init_gaussian_gridlevel
      59             :    USE global_types,                    ONLY: global_environment_type
      60             :    USE hirshfeld_types,                 ONLY: create_hirshfeld_type,&
      61             :                                               release_hirshfeld_type,&
      62             :                                               set_hirshfeld_info
      63             :    USE input_constants,                 ONLY: becke_cutoff_element,&
      64             :                                               mixed_cdft_parallel,&
      65             :                                               mixed_cdft_parallel_nobuild,&
      66             :                                               mixed_cdft_serial,&
      67             :                                               outer_scf_becke_constraint,&
      68             :                                               outer_scf_hirshfeld_constraint,&
      69             :                                               shape_function_gaussian
      70             :    USE input_section_types,             ONLY: section_vals_duplicate,&
      71             :                                               section_vals_get,&
      72             :                                               section_vals_get_subs_vals,&
      73             :                                               section_vals_release,&
      74             :                                               section_vals_type,&
      75             :                                               section_vals_val_get
      76             :    USE kinds,                           ONLY: default_path_length,&
      77             :                                               dp
      78             :    USE message_passing,                 ONLY: mp_request_type,&
      79             :                                               mp_waitall
      80             :    USE mixed_cdft_types,                ONLY: mixed_cdft_result_type_release,&
      81             :                                               mixed_cdft_result_type_set,&
      82             :                                               mixed_cdft_settings_type,&
      83             :                                               mixed_cdft_type,&
      84             :                                               mixed_cdft_work_type_init
      85             :    USE mixed_environment_types,         ONLY: get_mixed_env,&
      86             :                                               mixed_environment_type
      87             :    USE pw_env_methods,                  ONLY: pw_env_create
      88             :    USE pw_env_types,                    ONLY: pw_env_get,&
      89             :                                               pw_env_type
      90             :    USE pw_grid_types,                   ONLY: HALFSPACE,&
      91             :                                               pw_grid_type
      92             :    USE pw_grids,                        ONLY: do_pw_grid_blocked_false,&
      93             :                                               pw_grid_create,&
      94             :                                               pw_grid_release,&
      95             :                                               pw_grid_setup
      96             :    USE pw_pool_types,                   ONLY: pw_pool_create,&
      97             :                                               pw_pool_p_type,&
      98             :                                               pw_pool_type
      99             :    USE qs_cdft_types,                   ONLY: cdft_control_create,&
     100             :                                               cdft_control_type
     101             :    USE qs_environment_types,            ONLY: get_qs_env,&
     102             :                                               qs_environment_type
     103             :    USE qs_kind_types,                   ONLY: create_qs_kind_set,&
     104             :                                               qs_kind_type
     105             :    USE realspace_grid_types,            ONLY: realspace_grid_desc_p_type,&
     106             :                                               realspace_grid_input_type,&
     107             :                                               realspace_grid_type,&
     108             :                                               rs_grid_create,&
     109             :                                               rs_grid_create_descriptor,&
     110             :                                               rs_grid_print
     111             : #include "./base/base_uses.f90"
     112             : 
     113             :    IMPLICIT NONE
     114             :    PRIVATE
     115             : 
     116             : ! *** Public subroutines ***
     117             : 
     118             :    PUBLIC :: mixed_cdft_parse_settings, mixed_cdft_transfer_settings, &
     119             :              mixed_cdft_init_structures, mixed_cdft_redistribute_arrays, &
     120             :              mixed_cdft_print_couplings, map_permutation_to_states, hfun_zero, &
     121             :              mixed_cdft_release_work, mixed_cdft_read_block_diag, &
     122             :              mixed_cdft_get_blocks, mixed_cdft_diagonalize_blocks, &
     123             :              mixed_cdft_assemble_block_diag
     124             : 
     125             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mixed_cdft_utils'
     126             : 
     127             : CONTAINS
     128             : 
     129             : ! **************************************************************************************************
     130             : !> \brief Parse settings for mixed cdft calculation and check their consistency
     131             : !> \param force_env the force_env that holds the CDFT mixed_env
     132             : !> \param mixed_env the mixed_env that holds the CDFT states
     133             : !> \param mixed_cdft control section for mixed CDFT
     134             : !> \param settings container for settings related to the mixed CDFT calculation
     135             : !> \param natom the total number of atoms
     136             : !> \par History
     137             : !>       01.2017  created [Nico Holmberg]
     138             : ! **************************************************************************************************
     139          72 :    SUBROUTINE mixed_cdft_parse_settings(force_env, mixed_env, mixed_cdft, &
     140             :                                         settings, natom)
     141             :       TYPE(force_env_type), POINTER                      :: force_env
     142             :       TYPE(mixed_environment_type), POINTER              :: mixed_env
     143             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
     144             :       TYPE(mixed_cdft_settings_type)                     :: settings
     145             :       INTEGER                                            :: natom
     146             : 
     147             :       INTEGER                                            :: i, iatom, iforce_eval, igroup, &
     148             :                                                             nforce_eval, nkinds
     149          72 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: constraint_type
     150          72 :       INTEGER, ALLOCATABLE, DIMENSION(:, :, :)           :: array_sizes
     151             :       LOGICAL                                            :: is_match
     152          72 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     153             :       TYPE(cdft_control_type), POINTER                   :: cdft_control
     154          72 :       TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:, :) :: atoms
     155          72 :       TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: coeff
     156             :       TYPE(dft_control_type), POINTER                    :: dft_control
     157             :       TYPE(force_env_type), POINTER                      :: force_env_qs
     158             :       TYPE(pw_env_type), POINTER                         :: pw_env
     159             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     160             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     161             : 
     162          72 :       NULLIFY (dft_control, qs_env, pw_env, auxbas_pw_pool, force_env_qs, &
     163          72 :                cdft_control)
     164             :       ! Allocate storage for temporaries used for checking settings consistency
     165          72 :       settings%max_nkinds = 30
     166          72 :       nforce_eval = SIZE(force_env%sub_force_env)
     167         216 :       ALLOCATE (settings%grid_span(nforce_eval))
     168         216 :       ALLOCATE (settings%npts(3, nforce_eval))
     169         216 :       ALLOCATE (settings%cutoff(nforce_eval))
     170         216 :       ALLOCATE (settings%rel_cutoff(nforce_eval))
     171         216 :       ALLOCATE (settings%spherical(nforce_eval))
     172         216 :       ALLOCATE (settings%rs_dims(2, nforce_eval))
     173         216 :       ALLOCATE (settings%odd(nforce_eval))
     174         288 :       ALLOCATE (settings%atoms(natom, nforce_eval))
     175          72 :       IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
     176          88 :          ALLOCATE (settings%coeffs(natom, nforce_eval))
     177         154 :          settings%coeffs = 0.0_dp
     178             :       END IF
     179             :       ! Some of the checked settings are only defined for certain types of constraints
     180             :       ! We nonetheless use arrays that are large enough to contain settings for all constraints
     181             :       ! This is not completely optimal...
     182         216 :       ALLOCATE (settings%si(6, nforce_eval))
     183         216 :       ALLOCATE (settings%sb(8, nforce_eval))
     184         216 :       ALLOCATE (settings%sr(5, nforce_eval))
     185         216 :       ALLOCATE (settings%cutoffs(settings%max_nkinds, nforce_eval))
     186         216 :       ALLOCATE (settings%radii(settings%max_nkinds, nforce_eval))
     187         240 :       settings%grid_span = 0
     188         744 :       settings%npts = 0
     189         240 :       settings%cutoff = 0.0_dp
     190         240 :       settings%rel_cutoff = 0.0_dp
     191         240 :       settings%spherical = 0
     192          72 :       settings%is_spherical = .FALSE.
     193         576 :       settings%rs_dims = 0
     194         240 :       settings%odd = 0
     195          72 :       settings%is_odd = .FALSE.
     196         614 :       settings%atoms = 0
     197        1248 :       settings%si = 0
     198        1080 :       settings%sr = 0.0_dp
     199        1584 :       settings%sb = .FALSE.
     200        5280 :       settings%cutoffs = 0.0_dp
     201        5280 :       settings%radii = 0.0_dp
     202             :       ! Get information from the sub_force_envs
     203         240 :       DO iforce_eval = 1, nforce_eval
     204         168 :          IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
     205         144 :          force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
     206         144 :          IF (mixed_env%do_mixed_qmmm_cdft) THEN
     207          12 :             qs_env => force_env_qs%qmmm_env%qs_env
     208             :          ELSE
     209         132 :             CALL force_env_get(force_env_qs, qs_env=qs_env)
     210             :          END IF
     211         144 :          CALL get_qs_env(qs_env, pw_env=pw_env, dft_control=dft_control)
     212         144 :          IF (.NOT. dft_control%qs_control%cdft) &
     213             :             CALL cp_abort(__LOCATION__, &
     214             :                           "A mixed CDFT simulation with multiple force_evals was requested, "// &
     215           0 :                           "but CDFT constraints were not active in the QS section of all force_evals!")
     216         144 :          cdft_control => dft_control%qs_control%cdft_control
     217         144 :          CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     218        1440 :          settings%bo = auxbas_pw_pool%pw_grid%bounds_local
     219             :          ! Only the rank 0 process collects info about pw_grid and CDFT
     220         216 :          IF (force_env_qs%para_env%is_source()) THEN
     221             :             ! Grid settings
     222          84 :             settings%grid_span(iforce_eval) = auxbas_pw_pool%pw_grid%grid_span
     223         336 :             settings%npts(:, iforce_eval) = auxbas_pw_pool%pw_grid%npts
     224          84 :             settings%cutoff(iforce_eval) = auxbas_pw_pool%pw_grid%cutoff
     225          84 :             settings%rel_cutoff(iforce_eval) = dft_control%qs_control%relative_cutoff
     226          84 :             IF (auxbas_pw_pool%pw_grid%spherical) settings%spherical(iforce_eval) = 1
     227         252 :             settings%rs_dims(:, iforce_eval) = auxbas_pw_pool%pw_grid%para%rs_dims
     228          84 :             IF (auxbas_pw_pool%pw_grid%grid_span == HALFSPACE) settings%odd(iforce_eval) = 1
     229             :             ! Becke constraint atoms/coeffs
     230          84 :             IF (cdft_control%natoms .GT. SIZE(settings%atoms, 1)) &
     231             :                CALL cp_abort(__LOCATION__, &
     232             :                              "More CDFT constraint atoms than defined in mixed section. "// &
     233           0 :                              "Use default values for MIXED\MAPPING.")
     234         233 :             settings%atoms(1:cdft_control%natoms, iforce_eval) = cdft_control%atoms
     235          84 :             IF (mixed_cdft%run_type == mixed_cdft_parallel) &
     236          66 :                settings%coeffs(1:cdft_control%natoms, iforce_eval) = cdft_control%group(1)%coeff
     237             :             ! Integer type settings
     238          84 :             IF (cdft_control%type == outer_scf_becke_constraint) THEN
     239          76 :                settings%si(1, iforce_eval) = cdft_control%becke_control%cutoff_type
     240          76 :                settings%si(2, iforce_eval) = cdft_control%becke_control%cavity_shape
     241             :             END IF
     242          84 :             settings%si(3, iforce_eval) = dft_control%multiplicity
     243          84 :             settings%si(4, iforce_eval) = SIZE(cdft_control%group)
     244          84 :             settings%si(5, iforce_eval) = cdft_control%type
     245          84 :             IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
     246           8 :                settings%si(6, iforce_eval) = cdft_control%hirshfeld_control%shape_function
     247           8 :                settings%si(6, iforce_eval) = cdft_control%hirshfeld_control%gaussian_shape
     248             :             END IF
     249             :             ! Logicals
     250          84 :             IF (cdft_control%type == outer_scf_becke_constraint) THEN
     251          76 :                settings%sb(1, iforce_eval) = cdft_control%becke_control%cavity_confine
     252          76 :                settings%sb(2, iforce_eval) = cdft_control%becke_control%should_skip
     253          76 :                settings%sb(3, iforce_eval) = cdft_control%becke_control%print_cavity
     254          76 :                settings%sb(4, iforce_eval) = cdft_control%becke_control%in_memory
     255          76 :                settings%sb(5, iforce_eval) = cdft_control%becke_control%adjust
     256          76 :                settings%sb(8, iforce_eval) = cdft_control%becke_control%use_bohr
     257             :             END IF
     258          84 :             IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
     259           8 :                settings%sb(8, iforce_eval) = cdft_control%hirshfeld_control%use_bohr
     260             :             END IF
     261          84 :             settings%sb(6, iforce_eval) = cdft_control%atomic_charges
     262          84 :             settings%sb(7, iforce_eval) = qs_env%has_unit_metric
     263             :             ! Reals
     264          84 :             IF (cdft_control%type == outer_scf_becke_constraint) THEN
     265          76 :                settings%sr(1, iforce_eval) = cdft_control%becke_control%rcavity
     266          76 :                settings%sr(2, iforce_eval) = cdft_control%becke_control%rglobal
     267          76 :                settings%sr(3, iforce_eval) = cdft_control%becke_control%eps_cavity
     268             :             END IF
     269          84 :             IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
     270           8 :                settings%sr(2, iforce_eval) = cdft_control%hirshfeld_control%radius
     271             :             END IF
     272          84 :             settings%sr(4, iforce_eval) = dft_control%qs_control%eps_rho_rspace
     273          84 :             settings%sr(5, iforce_eval) = pw_env%cube_info(pw_env%auxbas_grid)%max_rad_ga
     274          84 :             IF (cdft_control%type == outer_scf_becke_constraint) THEN
     275          76 :                IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
     276          50 :                   nkinds = SIZE(cdft_control%becke_control%cutoffs_tmp)
     277          50 :                   IF (nkinds .GT. settings%max_nkinds) &
     278             :                      CALL cp_abort(__LOCATION__, &
     279             :                                    "More than "//TRIM(cp_to_string(settings%max_nkinds))// &
     280             :                                    " unique elements were defined in BECKE_CONSTRAINT\ELEMENT_CUTOFF. Are you sure"// &
     281           0 :                                    " your input is correct? If yes, please increase max_nkinds and recompile.")
     282         150 :                   settings%cutoffs(1:nkinds, iforce_eval) = cdft_control%becke_control%cutoffs_tmp(:)
     283             :                END IF
     284          76 :                IF (cdft_control%becke_control%adjust) THEN
     285          52 :                   CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
     286          52 :                   IF (.NOT. SIZE(atomic_kind_set) == SIZE(cdft_control%becke_control%radii_tmp)) &
     287             :                      CALL cp_abort(__LOCATION__, &
     288             :                                    "Length of keyword BECKE_CONSTRAINT\ATOMIC_RADII does not "// &
     289           0 :                                    "match number of atomic kinds in the input coordinate file.")
     290          52 :                   nkinds = SIZE(cdft_control%becke_control%radii_tmp)
     291          52 :                   IF (nkinds .GT. settings%max_nkinds) &
     292             :                      CALL cp_abort(__LOCATION__, &
     293             :                                    "More than "//TRIM(cp_to_string(settings%max_nkinds))// &
     294             :                                    " unique elements were defined in BECKE_CONSTRAINT\ATOMIC_RADII. Are you sure"// &
     295           0 :                                    " your input is correct? If yes, please increase max_nkinds and recompile.")
     296         156 :                   settings%radii(1:nkinds, iforce_eval) = cdft_control%becke_control%radii_tmp(:)
     297             :                END IF
     298             :             END IF
     299         108 :             IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
     300           8 :                IF (ASSOCIATED(cdft_control%hirshfeld_control%radii)) THEN
     301           0 :                   CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
     302           0 :                   IF (.NOT. SIZE(atomic_kind_set) == SIZE(cdft_control%hirshfeld_control%radii)) &
     303             :                      CALL cp_abort(__LOCATION__, &
     304             :                                    "Length of keyword HIRSHFELD_CONSTRAINT&RADII does not "// &
     305           0 :                                    "match number of atomic kinds in the input coordinate file.")
     306           0 :                   nkinds = SIZE(cdft_control%hirshfeld_control%radii)
     307           0 :                   IF (nkinds .GT. settings%max_nkinds) &
     308             :                      CALL cp_abort(__LOCATION__, &
     309             :                                    "More than "//TRIM(cp_to_string(settings%max_nkinds))// &
     310             :                                    " unique elements were defined in HIRSHFELD_CONSTRAINT&RADII. Are you sure"// &
     311           0 :                                    " your input is correct? If yes, please increase max_nkinds and recompile.")
     312           0 :                   settings%radii(1:nkinds, iforce_eval) = cdft_control%hirshfeld_control%radii(:)
     313             :                END IF
     314             :             END IF
     315             :          END IF
     316             :       END DO
     317             :       ! Make sure the grids are consistent
     318         408 :       CALL force_env%para_env%sum(settings%grid_span)
     319        1416 :       CALL force_env%para_env%sum(settings%npts)
     320         408 :       CALL force_env%para_env%sum(settings%cutoff)
     321         408 :       CALL force_env%para_env%sum(settings%rel_cutoff)
     322         408 :       CALL force_env%para_env%sum(settings%spherical)
     323        1080 :       CALL force_env%para_env%sum(settings%rs_dims)
     324         408 :       CALL force_env%para_env%sum(settings%odd)
     325          72 :       is_match = .TRUE.
     326         168 :       DO iforce_eval = 2, nforce_eval
     327          96 :          is_match = is_match .AND. (settings%grid_span(1) == settings%grid_span(iforce_eval))
     328          96 :          is_match = is_match .AND. (settings%npts(1, 1) == settings%npts(1, iforce_eval))
     329          96 :          is_match = is_match .AND. (settings%cutoff(1) == settings%cutoff(iforce_eval))
     330          96 :          is_match = is_match .AND. (settings%rel_cutoff(1) == settings%rel_cutoff(iforce_eval))
     331          96 :          is_match = is_match .AND. (settings%spherical(1) == settings%spherical(iforce_eval))
     332          96 :          is_match = is_match .AND. (settings%rs_dims(1, 1) == settings%rs_dims(1, iforce_eval))
     333          96 :          is_match = is_match .AND. (settings%rs_dims(2, 1) == settings%rs_dims(2, iforce_eval))
     334         168 :          is_match = is_match .AND. (settings%odd(1) == settings%odd(iforce_eval))
     335             :       END DO
     336          72 :       IF (.NOT. is_match) &
     337             :          CALL cp_abort(__LOCATION__, &
     338           0 :                        "Mismatch detected in the &MGRID settings of the CDFT force_evals.")
     339          72 :       IF (settings%spherical(1) == 1) settings%is_spherical = .TRUE.
     340          72 :       IF (settings%odd(1) == 1) settings%is_odd = .TRUE.
     341             :       ! Make sure CDFT settings are consistent
     342        1156 :       CALL force_env%para_env%sum(settings%atoms)
     343          72 :       IF (mixed_cdft%run_type == mixed_cdft_parallel) &
     344         286 :          CALL force_env%para_env%sum(settings%coeffs)
     345          72 :       settings%ncdft = 0
     346         224 :       DO i = 1, SIZE(settings%atoms, 1)
     347         374 :          DO iforce_eval = 2, nforce_eval
     348         374 :             IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
     349          44 :                IF (settings%atoms(i, 1) /= settings%atoms(i, iforce_eval)) is_match = .FALSE.
     350          44 :                IF (settings%coeffs(i, 1) /= settings%coeffs(i, iforce_eval)) is_match = .FALSE.
     351             :             END IF
     352             :          END DO
     353         224 :          IF (settings%atoms(i, 1) /= 0) settings%ncdft = settings%ncdft + 1
     354             :       END DO
     355          72 :       IF (.NOT. is_match .AND. mixed_cdft%run_type == mixed_cdft_parallel) &
     356             :          CALL cp_abort(__LOCATION__, &
     357             :                        "Mismatch detected in the &CDFT section of the CDFT force_evals. "// &
     358             :                        "Parallel mode mixed CDFT requires identical constraint definitions in both CDFT states. "// &
     359             :                        "Switch to serial mode or disable keyword PARALLEL_BUILD if you "// &
     360           0 :                        "want to use nonidentical constraint definitions.")
     361        2424 :       CALL force_env%para_env%sum(settings%si)
     362        2088 :       CALL force_env%para_env%sum(settings%sr)
     363         648 :       DO i = 1, SIZE(settings%sb, 1)
     364         576 :          CALL force_env%para_env%sum(settings%sb(i, 1))
     365        1416 :          DO iforce_eval = 2, nforce_eval
     366         768 :             CALL force_env%para_env%sum(settings%sb(i, iforce_eval))
     367        1344 :             IF (settings%sb(i, 1) .NEQV. settings%sb(i, iforce_eval)) is_match = .FALSE.
     368             :          END DO
     369             :       END DO
     370         504 :       DO i = 1, SIZE(settings%si, 1)
     371        1080 :          DO iforce_eval = 2, nforce_eval
     372        1008 :             IF (settings%si(i, 1) /= settings%si(i, iforce_eval)) is_match = .FALSE.
     373             :          END DO
     374             :       END DO
     375         432 :       DO i = 1, SIZE(settings%sr, 1)
     376         912 :          DO iforce_eval = 2, nforce_eval
     377         840 :             IF (settings%sr(i, 1) /= settings%sr(i, iforce_eval)) is_match = .FALSE.
     378             :          END DO
     379             :       END DO
     380          72 :       IF (.NOT. is_match) &
     381             :          CALL cp_abort(__LOCATION__, &
     382           0 :                        "Mismatch detected in the &CDFT settings of the CDFT force_evals.")
     383             :       ! Some CDFT features are currently disabled for mixed calculations: check that these features were not requested
     384          72 :       IF (mixed_cdft%dlb .AND. .NOT. settings%sb(1, 1)) &
     385             :          CALL cp_abort(__LOCATION__, &
     386           0 :                        "Parallel mode mixed CDFT load balancing requires Gaussian cavity confinement.")
     387             :       ! Check for identical constraints in case of run type serial/parallel_nobuild
     388          72 :       IF (mixed_cdft%run_type /= mixed_cdft_parallel) THEN
     389             :          ! Get array sizes
     390         250 :          ALLOCATE (array_sizes(nforce_eval, settings%si(4, 1), 2))
     391         510 :          array_sizes(:, :, :) = 0
     392         174 :          DO iforce_eval = 1, nforce_eval
     393         124 :             IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
     394         122 :             force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
     395         122 :             IF (mixed_env%do_mixed_qmmm_cdft) THEN
     396           8 :                qs_env => force_env_qs%qmmm_env%qs_env
     397             :             ELSE
     398         114 :                CALL force_env_get(force_env_qs, qs_env=qs_env)
     399             :             END IF
     400         122 :             CALL get_qs_env(qs_env, dft_control=dft_control)
     401         122 :             cdft_control => dft_control%qs_control%cdft_control
     402         172 :             IF (force_env_qs%para_env%is_source()) THEN
     403         128 :                DO igroup = 1, SIZE(cdft_control%group)
     404          64 :                   array_sizes(iforce_eval, igroup, 1) = SIZE(cdft_control%group(igroup)%atoms)
     405         126 :                   array_sizes(iforce_eval, igroup, 2) = SIZE(cdft_control%group(igroup)%coeff)
     406             :                END DO
     407             :             END IF
     408             :          END DO
     409             :          ! Sum up array sizes and check consistency
     410          50 :          CALL force_env%para_env%sum(array_sizes)
     411         410 :          IF (ANY(array_sizes(:, :, 1) /= array_sizes(1, 1, 1)) .OR. &
     412             :              ANY(array_sizes(:, :, 2) /= array_sizes(1, 1, 2))) &
     413           0 :             mixed_cdft%identical_constraints = .FALSE.
     414             :          ! Check constraint definitions
     415          50 :          IF (mixed_cdft%identical_constraints) THEN
     416             :             ! Prepare temporary storage
     417         380 :             ALLOCATE (atoms(nforce_eval, settings%si(4, 1)))
     418         380 :             ALLOCATE (coeff(nforce_eval, settings%si(4, 1)))
     419         200 :             ALLOCATE (constraint_type(nforce_eval, settings%si(4, 1)))
     420         230 :             constraint_type(:, :) = 0
     421         174 :             DO iforce_eval = 1, nforce_eval
     422         252 :                DO i = 1, settings%si(4, 1)
     423         128 :                   NULLIFY (atoms(iforce_eval, i)%array)
     424             :                   NULLIFY (coeff(iforce_eval, i)%array)
     425         384 :                   ALLOCATE (atoms(iforce_eval, i)%array(array_sizes(iforce_eval, i, 1)))
     426         384 :                   ALLOCATE (coeff(iforce_eval, i)%array(array_sizes(iforce_eval, i, 1)))
     427         346 :                   atoms(iforce_eval, i)%array(:) = 0
     428         470 :                   coeff(iforce_eval, i)%array(:) = 0
     429             :                END DO
     430             :                ! Get constraint definitions
     431         124 :                IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
     432         122 :                force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
     433         122 :                IF (mixed_env%do_mixed_qmmm_cdft) THEN
     434           8 :                   qs_env => force_env_qs%qmmm_env%qs_env
     435             :                ELSE
     436         114 :                   CALL force_env_get(force_env_qs, qs_env=qs_env)
     437             :                END IF
     438         122 :                CALL get_qs_env(qs_env, dft_control=dft_control)
     439         122 :                cdft_control => dft_control%qs_control%cdft_control
     440         172 :                IF (force_env_qs%para_env%is_source()) THEN
     441         128 :                   DO i = 1, settings%si(4, 1)
     442         173 :                      atoms(iforce_eval, i)%array(:) = cdft_control%group(i)%atoms
     443         173 :                      coeff(iforce_eval, i)%array(:) = cdft_control%group(i)%coeff
     444         126 :                      constraint_type(iforce_eval, i) = cdft_control%group(i)%constraint_type
     445             :                   END DO
     446             :                END IF
     447             :             END DO
     448             :             ! Sum up constraint definitions and check consistency
     449         100 :             DO i = 1, settings%si(4, 1)
     450         180 :                DO iforce_eval = 1, nforce_eval
     451         564 :                   CALL force_env%para_env%sum(atoms(iforce_eval, i)%array)
     452         564 :                   CALL force_env%para_env%sum(coeff(iforce_eval, i)%array)
     453         180 :                   CALL force_env%para_env%sum(constraint_type(iforce_eval, i))
     454             :                END DO
     455         126 :                DO iforce_eval = 2, nforce_eval
     456         194 :                   DO iatom = 1, SIZE(atoms(1, i)%array)
     457         120 :                      IF (atoms(1, i)%array(iatom) /= atoms(iforce_eval, i)%array(iatom)) &
     458           0 :                         mixed_cdft%identical_constraints = .FALSE.
     459         120 :                      IF (coeff(1, i)%array(iatom) /= coeff(iforce_eval, i)%array(iatom)) &
     460           2 :                         mixed_cdft%identical_constraints = .FALSE.
     461         194 :                      IF (.NOT. mixed_cdft%identical_constraints) EXIT
     462             :                   END DO
     463          76 :                   IF (constraint_type(1, i) /= constraint_type(iforce_eval, i)) &
     464           0 :                      mixed_cdft%identical_constraints = .FALSE.
     465         126 :                   IF (.NOT. mixed_cdft%identical_constraints) EXIT
     466             :                END DO
     467         100 :                IF (.NOT. mixed_cdft%identical_constraints) EXIT
     468             :             END DO
     469             :             ! Deallocate temporary storage
     470         174 :             DO iforce_eval = 1, nforce_eval
     471         302 :                DO i = 1, settings%si(4, 1)
     472         128 :                   DEALLOCATE (atoms(iforce_eval, i)%array)
     473         252 :                   DEALLOCATE (coeff(iforce_eval, i)%array)
     474             :                END DO
     475             :             END DO
     476          50 :             DEALLOCATE (atoms)
     477          50 :             DEALLOCATE (coeff)
     478          50 :             DEALLOCATE (constraint_type)
     479             :          END IF
     480          50 :          DEALLOCATE (array_sizes)
     481             :       END IF
     482             :       ! Deallocate some arrays that are no longer needed
     483          72 :       IF (mixed_cdft%identical_constraints .AND. mixed_cdft%run_type /= mixed_cdft_parallel_nobuild) THEN
     484         228 :          DO iforce_eval = 1, nforce_eval
     485         160 :             IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
     486         138 :             force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
     487         138 :             IF (mixed_env%do_mixed_qmmm_cdft) THEN
     488          12 :                qs_env => force_env_qs%qmmm_env%qs_env
     489             :             ELSE
     490         126 :                CALL force_env_get(force_env_qs, qs_env=qs_env)
     491             :             END IF
     492         138 :             CALL get_qs_env(qs_env, dft_control=dft_control)
     493         138 :             cdft_control => dft_control%qs_control%cdft_control
     494         206 :             IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
     495          22 :                IF (.NOT. dft_control%qs_control%gapw) THEN
     496          44 :                   DO i = 1, SIZE(cdft_control%group)
     497          22 :                      DEALLOCATE (cdft_control%group(i)%coeff)
     498          44 :                      DEALLOCATE (cdft_control%group(i)%atoms)
     499             :                   END DO
     500          22 :                   IF (.NOT. cdft_control%atomic_charges) DEALLOCATE (cdft_control%atoms)
     501             :                END IF
     502         116 :             ELSE IF (mixed_cdft%run_type == mixed_cdft_serial) THEN
     503         116 :                IF (iforce_eval == 1) CYCLE
     504         142 :                DO igroup = 1, SIZE(cdft_control%group)
     505         142 :                   IF (.NOT. dft_control%qs_control%gapw) THEN
     506          72 :                      DEALLOCATE (cdft_control%group(igroup)%coeff)
     507          72 :                      DEALLOCATE (cdft_control%group(igroup)%atoms)
     508             :                   END IF
     509             :                END DO
     510          70 :                IF (cdft_control%type == outer_scf_becke_constraint) THEN
     511          64 :                   IF (.NOT. cdft_control%atomic_charges) DEALLOCATE (cdft_control%atoms)
     512          64 :                   IF (cdft_control%becke_control%cavity_confine) &
     513          62 :                      CALL release_hirshfeld_type(cdft_control%becke_control%cavity_env)
     514          64 :                   IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) &
     515          42 :                      DEALLOCATE (cdft_control%becke_control%cutoffs_tmp)
     516          64 :                   IF (cdft_control%becke_control%adjust) &
     517          44 :                      DEALLOCATE (cdft_control%becke_control%radii_tmp)
     518             :                END IF
     519             :             END IF
     520             :          END DO
     521             :       END IF
     522             : 
     523         144 :    END SUBROUTINE mixed_cdft_parse_settings
     524             : 
     525             : ! **************************************************************************************************
     526             : !> \brief Transfer settings to mixed_cdft
     527             : !> \param force_env the force_env that holds the CDFT states
     528             : !> \param mixed_cdft the control section for mixed CDFT calculations
     529             : !> \param settings container for settings related to the mixed CDFT calculation
     530             : !> \par History
     531             : !>       01.2017  created [Nico Holmberg]
     532             : ! **************************************************************************************************
     533          72 :    SUBROUTINE mixed_cdft_transfer_settings(force_env, mixed_cdft, settings)
     534             :       TYPE(force_env_type), POINTER                      :: force_env
     535             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
     536             :       TYPE(mixed_cdft_settings_type)                     :: settings
     537             : 
     538             :       INTEGER                                            :: i, nkinds
     539             :       LOGICAL                                            :: is_match
     540             :       TYPE(cdft_control_type), POINTER                   :: cdft_control
     541             : 
     542          72 :       NULLIFY (cdft_control)
     543          72 :       is_match = .TRUE.
     544             :       ! Transfer global settings
     545          72 :       mixed_cdft%multiplicity = settings%si(3, 1)
     546          72 :       mixed_cdft%has_unit_metric = settings%sb(7, 1)
     547          72 :       mixed_cdft%eps_rho_rspace = settings%sr(4, 1)
     548          72 :       mixed_cdft%nconstraint = settings%si(4, 1)
     549          72 :       settings%radius = settings%sr(5, 1)
     550             :       ! Transfer settings only needed if the constraint should be built in parallel
     551          72 :       IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
     552          22 :          IF (settings%sb(6, 1)) &
     553             :             CALL cp_abort(__LOCATION__, &
     554           0 :                           "Calculation of atomic Becke charges not supported with parallel mode mixed CDFT")
     555          22 :          IF (mixed_cdft%nconstraint /= 1) &
     556             :             CALL cp_abort(__LOCATION__, &
     557           0 :                           "Parallel mode mixed CDFT does not yet support multiple constraints.")
     558             : 
     559          22 :          IF (settings%si(5, 1) /= outer_scf_becke_constraint) &
     560             :             CALL cp_abort(__LOCATION__, &
     561           0 :                           "Parallel mode mixed CDFT does not support Hirshfeld constraints.")
     562             : 
     563          22 :          ALLOCATE (mixed_cdft%cdft_control)
     564          22 :          CALL cdft_control_create(mixed_cdft%cdft_control)
     565          22 :          cdft_control => mixed_cdft%cdft_control
     566          66 :          ALLOCATE (cdft_control%atoms(settings%ncdft))
     567          66 :          cdft_control%atoms = settings%atoms(1:settings%ncdft, 1)
     568          44 :          ALLOCATE (cdft_control%group(1))
     569          66 :          ALLOCATE (cdft_control%group(1)%atoms(settings%ncdft))
     570          66 :          ALLOCATE (cdft_control%group(1)%coeff(settings%ncdft))
     571          22 :          NULLIFY (cdft_control%group(1)%weight)
     572          22 :          NULLIFY (cdft_control%group(1)%gradients)
     573          22 :          NULLIFY (cdft_control%group(1)%integrated)
     574          66 :          cdft_control%group(1)%atoms = cdft_control%atoms
     575          66 :          cdft_control%group(1)%coeff = settings%coeffs(1:settings%ncdft, 1)
     576          22 :          cdft_control%natoms = settings%ncdft
     577          22 :          cdft_control%atomic_charges = settings%sb(6, 1)
     578          22 :          cdft_control%becke_control%cutoff_type = settings%si(1, 1)
     579          22 :          cdft_control%becke_control%cavity_confine = settings%sb(1, 1)
     580          22 :          cdft_control%becke_control%should_skip = settings%sb(2, 1)
     581          22 :          cdft_control%becke_control%print_cavity = settings%sb(3, 1)
     582          22 :          cdft_control%becke_control%in_memory = settings%sb(4, 1)
     583          22 :          cdft_control%becke_control%adjust = settings%sb(5, 1)
     584          22 :          cdft_control%becke_control%cavity_shape = settings%si(2, 1)
     585          22 :          cdft_control%becke_control%use_bohr = settings%sb(8, 1)
     586          22 :          cdft_control%becke_control%rcavity = settings%sr(1, 1)
     587          22 :          cdft_control%becke_control%rglobal = settings%sr(2, 1)
     588          22 :          cdft_control%becke_control%eps_cavity = settings%sr(3, 1)
     589          22 :          nkinds = 0
     590          22 :          IF (cdft_control%becke_control%cutoff_type == becke_cutoff_element) THEN
     591        2250 :             CALL force_env%para_env%sum(settings%cutoffs)
     592         558 :             DO i = 1, SIZE(settings%cutoffs, 1)
     593         540 :                IF (settings%cutoffs(i, 1) /= settings%cutoffs(i, 2)) is_match = .FALSE.
     594         558 :                IF (settings%cutoffs(i, 1) /= 0.0_dp) nkinds = nkinds + 1
     595             :             END DO
     596          18 :             IF (.NOT. is_match) &
     597             :                CALL cp_abort(__LOCATION__, &
     598             :                              "Mismatch detected in the &BECKE_CONSTRAINT "// &
     599           0 :                              "&ELEMENT_CUTOFF settings of the two force_evals.")
     600          54 :             ALLOCATE (cdft_control%becke_control%cutoffs_tmp(nkinds))
     601          54 :             cdft_control%becke_control%cutoffs_tmp = settings%cutoffs(1:nkinds, 1)
     602             :          END IF
     603          22 :          nkinds = 0
     604          22 :          IF (cdft_control%becke_control%adjust) THEN
     605        2250 :             CALL force_env%para_env%sum(settings%radii)
     606         558 :             DO i = 1, SIZE(settings%radii, 1)
     607         540 :                IF (settings%radii(i, 1) /= settings%radii(i, 2)) is_match = .FALSE.
     608         558 :                IF (settings%radii(i, 1) /= 0.0_dp) nkinds = nkinds + 1
     609             :             END DO
     610          18 :             IF (.NOT. is_match) &
     611             :                CALL cp_abort(__LOCATION__, &
     612             :                              "Mismatch detected in the &BECKE_CONSTRAINT "// &
     613           0 :                              "&ATOMIC_RADII settings of the two force_evals.")
     614          54 :             ALLOCATE (cdft_control%becke_control%radii(nkinds))
     615          54 :             cdft_control%becke_control%radii = settings%radii(1:nkinds, 1)
     616             :          END IF
     617             :       END IF
     618             : 
     619          72 :    END SUBROUTINE mixed_cdft_transfer_settings
     620             : 
     621             : ! **************************************************************************************************
     622             : !> \brief Initialize all the structures needed for a mixed CDFT calculation
     623             : !> \param force_env the force_env that holds the CDFT mixed_env
     624             : !> \param force_env_qs the force_env that holds the qs_env, which is CDFT state specific
     625             : !> \param mixed_env the mixed_env that holds the CDFT states
     626             : !> \param mixed_cdft the control section for mixed CDFT calculations
     627             : !> \param settings container for settings related to the mixed CDFT calculation
     628             : !> \par History
     629             : !>       01.2017 created [Nico Holmberg]
     630             : ! **************************************************************************************************
     631          72 :    SUBROUTINE mixed_cdft_init_structures(force_env, force_env_qs, mixed_env, mixed_cdft, settings)
     632             :       TYPE(force_env_type), POINTER                      :: force_env, force_env_qs
     633             :       TYPE(mixed_environment_type), POINTER              :: mixed_env
     634             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
     635             :       TYPE(mixed_cdft_settings_type)                     :: settings
     636             : 
     637             :       CHARACTER(len=default_path_length)                 :: c_val, input_file_path, output_file_path
     638             :       INTEGER                                            :: i, imap, iounit, j, lp, n_force_eval, &
     639             :                                                             ncpu, nforce_eval, ntargets, offset, &
     640             :                                                             unit_nr
     641          72 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: bounds
     642             :       INTEGER, DIMENSION(2, 3)                           :: bo, bo_mixed
     643             :       INTEGER, DIMENSION(3)                              :: higher_grid_layout
     644         144 :       INTEGER, DIMENSION(:), POINTER                     :: i_force_eval, mixed_rs_dims, recvbuffer, &
     645         144 :                                                             recvbuffer2, sendbuffer
     646             :       LOGICAL                                            :: is_match
     647          72 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     648             :       TYPE(cell_type), POINTER                           :: cell_mix
     649             :       TYPE(cp_logger_type), POINTER                      :: logger
     650             :       TYPE(cp_subsys_type), POINTER                      :: subsys_mix
     651             :       TYPE(global_environment_type), POINTER             :: globenv
     652         288 :       TYPE(mp_request_type), DIMENSION(3)                :: req
     653             :       TYPE(pw_env_type), POINTER                         :: pw_env
     654             :       TYPE(pw_grid_type), POINTER                        :: pw_grid
     655          72 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: pw_pools
     656             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
     657             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     658          72 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     659             :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
     660          72 :          POINTER                                         :: rs_descs
     661             :       TYPE(realspace_grid_input_type)                    :: input_settings
     662          72 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: rs_grids
     663             :       TYPE(section_vals_type), POINTER :: force_env_section, force_env_sections, kind_section, &
     664             :          print_section, root_section, rs_grid_section, subsys_section
     665             : 
     666          72 :       NULLIFY (cell_mix, subsys_mix, force_env_section, subsys_section, &
     667          72 :                print_section, root_section, kind_section, force_env_sections, &
     668          72 :                rs_grid_section, auxbas_pw_pool, pw_env, pw_pools, pw_grid, &
     669          72 :                sendbuffer, qs_env, mixed_rs_dims, i_force_eval, recvbuffer, &
     670          72 :                recvbuffer2, globenv, atomic_kind_set, qs_kind_set, rs_descs, &
     671          72 :                rs_grids)
     672             : 
     673         144 :       logger => cp_get_default_logger()
     674          72 :       CALL force_env_get(force_env=force_env, force_env_section=force_env_section)
     675          72 :       print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
     676          72 :       iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
     677          72 :       is_match = .TRUE.
     678          72 :       nforce_eval = SIZE(force_env%sub_force_env)
     679          72 :       ncpu = force_env%para_env%num_pe
     680             :       ! Get infos about the mixed subsys
     681          72 :       IF (.NOT. mixed_env%do_mixed_qmmm_cdft) THEN
     682             :          CALL force_env_get(force_env=force_env, &
     683          64 :                             subsys=subsys_mix)
     684             :       ELSE
     685             :          CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
     686           8 :                          cp_subsys=subsys_mix)
     687             :       END IF
     688             :       ! Init structures only needed when the CDFT states are treated in parallel
     689          72 :       IF (mixed_cdft%run_type == mixed_cdft_parallel) THEN
     690             :          ! Start building the mixed auxbas_pw_pool
     691          22 :          CALL pw_grid_create(pw_grid, force_env%para_env)
     692          22 :          CALL pw_env_create(mixed_cdft%pw_env)
     693             :          ! Decide what kind of layout to use and setup the grid
     694             :          ! Processor mappings currently supported:
     695             :          !  (2np,1)  --> (np,1)
     696             :          !  (nx,2ny) --> (nx,ny)
     697             :          !  (nx,ny)  --> (nx*ny/2,1) (required when xc_smooth is in use and with intermediate proc counts)
     698             :          !
     699             :          ! For cases 2 and 3, dlb redistributes YZ slices from overloaded processors to underloaded processors
     700             :          ! For case 1, XZ slices are redistributed
     701             :          ! TODO: Unify mappings. Now we essentially have separate code for cases 1-2 and 3.
     702             :          !       This leads to very messy code especially with dlb turned on...
     703             :          !       In terms of memory usage, it would be beneficial to replace case 1 with 3
     704             :          !       and implement a similar arbitrary mapping to replace case 2
     705             : 
     706          22 :          mixed_cdft%is_pencil = .FALSE. ! Flag to control the first two mappings
     707          22 :          mixed_cdft%is_special = .FALSE. ! Flag to control the last mapping
     708             :          ! With xc smoothing, the grid is always (ncpu/2,1) distributed
     709             :          ! and correct behavior cannot be guaranteed for ncpu/2 > nx, so we abort...
     710          22 :          IF (ncpu/2 .GT. settings%npts(1, 1)) &
     711           0 :             CPABORT("ncpu/2 => nx: decrease ncpu or disable xc_smoothing")
     712             :          !
     713          22 :          ALLOCATE (mixed_rs_dims(2))
     714          22 :          IF (settings%rs_dims(2, 1) /= 1) mixed_cdft%is_pencil = .TRUE.
     715          22 :          IF (.NOT. mixed_cdft%is_pencil .AND. ncpu .GT. settings%npts(1, 1)) mixed_cdft%is_special = .TRUE.
     716          22 :          IF (mixed_cdft%is_special) THEN
     717           0 :             mixed_rs_dims = (/-1, -1/)
     718          22 :          ELSE IF (mixed_cdft%is_pencil) THEN
     719           0 :             mixed_rs_dims = (/settings%rs_dims(1, 1), 2*settings%rs_dims(2, 1)/)
     720             :          ELSE
     721          66 :             mixed_rs_dims = (/2*settings%rs_dims(1, 1), 1/)
     722             :          END IF
     723          22 :          IF (.NOT. mixed_env%do_mixed_qmmm_cdft) THEN
     724             :             CALL force_env_get(force_env=force_env, &
     725          18 :                                cell=cell_mix)
     726             :          ELSE
     727             :             CALL get_qs_env(force_env_qs%qmmm_env%qs_env, &
     728           4 :                             cell=cell_mix)
     729             :          END IF
     730             :          CALL pw_grid_setup(cell_mix%hmat, pw_grid, grid_span=settings%grid_span(1), &
     731             :                             npts=settings%npts(:, 1), cutoff=settings%cutoff(1), &
     732             :                             spherical=settings%is_spherical, odd=settings%is_odd, &
     733             :                             fft_usage=.TRUE., ncommensurate=0, icommensurate=1, &
     734             :                             blocked=do_pw_grid_blocked_false, rs_dims=mixed_rs_dims, &
     735          22 :                             iounit=iounit)
     736             :          ! Check if the layout was successfully created
     737          22 :          IF (mixed_cdft%is_special) THEN
     738           0 :             IF (.NOT. pw_grid%para%rs_dims(2) /= 1) is_match = .FALSE.
     739          22 :          ELSE IF (mixed_cdft%is_pencil) THEN
     740           0 :             IF (.NOT. pw_grid%para%rs_dims(1) == mixed_rs_dims(1)) is_match = .FALSE.
     741             :          ELSE
     742          22 :             IF (.NOT. pw_grid%para%rs_dims(2) == 1) is_match = .FALSE.
     743             :          END IF
     744             :          IF (.NOT. is_match) &
     745             :             CALL cp_abort(__LOCATION__, &
     746             :                           "Unable to create a suitable grid distribution "// &
     747             :                           "for mixed CDFT calculations. Try decreasing the total number "// &
     748           0 :                           "of processors or disabling xc_smoothing.")
     749          22 :          DEALLOCATE (mixed_rs_dims)
     750             :          ! Create the pool
     751         220 :          bo_mixed = pw_grid%bounds_local
     752          44 :          ALLOCATE (pw_pools(1))
     753          22 :          NULLIFY (pw_pools(1)%pool)
     754          22 :          CALL pw_pool_create(pw_pools(1)%pool, pw_grid=pw_grid)
     755             :          ! Initialize Gaussian cavity confinement
     756          22 :          IF (mixed_cdft%cdft_control%becke_control%cavity_confine) THEN
     757          22 :             CALL create_hirshfeld_type(mixed_cdft%cdft_control%becke_control%cavity_env)
     758             :             CALL set_hirshfeld_info(mixed_cdft%cdft_control%becke_control%cavity_env, &
     759             :                                     shape_function_type=shape_function_gaussian, iterative=.FALSE., &
     760             :                                     radius_type=mixed_cdft%cdft_control%becke_control%cavity_shape, &
     761          22 :                                     use_bohr=mixed_cdft%cdft_control%becke_control%use_bohr)
     762             :          END IF
     763             :          ! Gaussian confinement/wavefunction overlap method needs qs_kind_set
     764             :          ! Gaussian cavity confinement also needs the auxbas_rs_grid
     765          22 :          IF (mixed_cdft%cdft_control%becke_control%cavity_confine .OR. &
     766             :              mixed_cdft%wfn_overlap_method) THEN
     767             :             print_section => section_vals_get_subs_vals(force_env_section, &
     768          22 :                                                         "PRINT%GRID_INFORMATION")
     769          22 :             ALLOCATE (mixed_cdft%pw_env%gridlevel_info)
     770             :             CALL init_gaussian_gridlevel(mixed_cdft%pw_env%gridlevel_info, &
     771             :                                          ngrid_levels=1, cutoff=settings%cutoff, &
     772             :                                          rel_cutoff=settings%rel_cutoff(1), &
     773          22 :                                          print_section=print_section)
     774          44 :             ALLOCATE (rs_descs(1))
     775         374 :             ALLOCATE (rs_grids(1))
     776         638 :             ALLOCATE (mixed_cdft%pw_env%cube_info(1))
     777          22 :             higher_grid_layout = (/-1, -1, -1/)
     778          22 :             CALL init_d3_poly_module()
     779             :             CALL init_cube_info(mixed_cdft%pw_env%cube_info(1), &
     780             :                                 pw_grid%dr(:), pw_grid%dh(:, :), &
     781             :                                 pw_grid%dh_inv(:, :), &
     782          22 :                                 pw_grid%orthorhombic, settings%radius)
     783          22 :             NULLIFY (root_section, force_env_section, force_env_sections, rs_grid_section)
     784          22 :             CALL force_env_get(force_env, root_section=root_section)
     785          22 :             force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
     786          22 :             CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, n_force_eval)
     787             :             CALL section_vals_duplicate(force_env_sections, force_env_section, &
     788          22 :                                         i_force_eval(2), i_force_eval(2))
     789          22 :             rs_grid_section => section_vals_get_subs_vals(force_env_section, "DFT%MGRID%RS_GRID")
     790             :             CALL init_input_type(input_settings, &
     791             :                                  nsmax=2*MAX(1, return_cube_max_iradius(mixed_cdft%pw_env%cube_info(1))) + 1, &
     792             :                                  rs_grid_section=rs_grid_section, ilevel=1, &
     793          22 :                                  higher_grid_layout=higher_grid_layout)
     794          22 :             NULLIFY (rs_descs(1)%rs_desc)
     795          22 :             CALL rs_grid_create_descriptor(rs_descs(1)%rs_desc, pw_grid, input_settings)
     796          22 :             IF (rs_descs(1)%rs_desc%distributed) higher_grid_layout = rs_descs(1)%rs_desc%group_dim
     797          22 :             CALL rs_grid_create(rs_grids(1), rs_descs(1)%rs_desc)
     798          22 :             CALL rs_grid_print(rs_grids(1), iounit)
     799          22 :             mixed_cdft%pw_env%rs_descs => rs_descs
     800          22 :             mixed_cdft%pw_env%rs_grids => rs_grids
     801             :             ! qs_kind_set
     802             :             subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
     803          22 :                                                          i_rep_section=i_force_eval(1))
     804          22 :             kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
     805          22 :             NULLIFY (qs_kind_set)
     806          22 :             CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
     807             :             CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
     808          22 :                                     force_env%para_env, force_env_section)
     809          22 :             mixed_cdft%qs_kind_set => qs_kind_set
     810          22 :             DEALLOCATE (i_force_eval)
     811          22 :             CALL section_vals_release(force_env_section)
     812             :          END IF
     813             :          CALL force_env_get(force_env=force_env, &
     814          22 :                             force_env_section=force_env_section)
     815          22 :          CALL pw_grid_release(pw_grid)
     816          22 :          mixed_cdft%pw_env%auxbas_grid = 1
     817          22 :          NULLIFY (mixed_cdft%pw_env%pw_pools)
     818          22 :          mixed_cdft%pw_env%pw_pools => pw_pools
     819         220 :          bo = settings%bo
     820             :          ! Determine which processors need to exchange data when redistributing the weight/gradient
     821          22 :          IF (.NOT. mixed_cdft%is_special) THEN
     822          22 :             ALLOCATE (mixed_cdft%dest_list(2))
     823          22 :             ALLOCATE (mixed_cdft%source_list(2))
     824          22 :             imap = force_env%para_env%mepos/2
     825          66 :             mixed_cdft%dest_list = (/imap, imap + force_env%para_env%num_pe/2/)
     826             :             imap = MOD(force_env%para_env%mepos, force_env%para_env%num_pe/2) + &
     827          22 :                    MODULO(force_env%para_env%mepos, force_env%para_env%num_pe/2)
     828          66 :             mixed_cdft%source_list = (/imap, imap + 1/)
     829             :             ! Determine bounds of the data that is replicated
     830          22 :             ALLOCATE (mixed_cdft%recv_bo(4))
     831          22 :             ALLOCATE (sendbuffer(2), recvbuffer(2), recvbuffer2(2))
     832          22 :             IF (mixed_cdft%is_pencil) THEN
     833           0 :                sendbuffer = (/bo_mixed(1, 2), bo_mixed(2, 2)/)
     834             :             ELSE
     835          66 :                sendbuffer = (/bo_mixed(1, 1), bo_mixed(2, 1)/)
     836             :             END IF
     837             :             ! Communicate bounds in steps
     838             :             CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dest_list(1), &
     839          22 :                                           request=req(1))
     840             :             CALL force_env%para_env%irecv(msgout=recvbuffer, source=mixed_cdft%source_list(1), &
     841          22 :                                           request=req(2))
     842             :             CALL force_env%para_env%irecv(msgout=recvbuffer2, source=mixed_cdft%source_list(2), &
     843          22 :                                           request=req(3))
     844          22 :             CALL req(1)%wait()
     845             :             CALL force_env%para_env%isend(msgin=sendbuffer, dest=mixed_cdft%dest_list(2), &
     846          22 :                                           request=req(1))
     847          22 :             CALL mp_waitall(req)
     848         110 :             mixed_cdft%recv_bo(1:2) = recvbuffer
     849         110 :             mixed_cdft%recv_bo(3:4) = recvbuffer2
     850          22 :             DEALLOCATE (sendbuffer, recvbuffer, recvbuffer2)
     851             :          ELSE
     852           0 :             IF (mixed_env%do_mixed_qmmm_cdft) THEN
     853           0 :                qs_env => force_env_qs%qmmm_env%qs_env
     854             :             ELSE
     855           0 :                CALL force_env_get(force_env_qs, qs_env=qs_env)
     856             :             END IF
     857           0 :             CALL get_qs_env(qs_env, pw_env=pw_env)
     858           0 :             CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
     859             :             ! work out the pw grid points each proc holds in the two (identical) parallel proc groups
     860             :             ! note we only care about the x dir since we assume the y dir is not subdivided
     861           0 :             ALLOCATE (bounds(0:auxbas_pw_pool%pw_grid%para%group_size - 1, 1:2))
     862           0 :             DO i = 0, auxbas_pw_pool%pw_grid%para%group_size - 1
     863           0 :                bounds(i, 1:2) = auxbas_pw_pool%pw_grid%para%bo(1:2, 1, i, 1)
     864           0 :                bounds(i, 1:2) = bounds(i, 1:2) - auxbas_pw_pool%pw_grid%npts(1)/2 - 1
     865             :             END DO
     866             :             ! work out which procs to send my grid points
     867             :             ! first get the number of target procs per group
     868           0 :             ntargets = 0
     869           0 :             offset = -1
     870           0 :             DO i = 0, auxbas_pw_pool%pw_grid%para%group_size - 1
     871           0 :                IF ((bounds(i, 1) .GE. bo_mixed(1, 1) .AND. bounds(i, 1) .LE. bo_mixed(2, 1)) .OR. &
     872           0 :                    (bounds(i, 2) .GE. bo_mixed(1, 1) .AND. bounds(i, 2) .LE. bo_mixed(2, 1))) THEN
     873           0 :                   ntargets = ntargets + 1
     874           0 :                   IF (offset == -1) offset = i
     875           0 :                ELSE IF (bounds(i, 2) .GT. bo_mixed(2, 1)) THEN
     876             :                   EXIT
     877             :                ELSE
     878           0 :                   CYCLE
     879             :                END IF
     880             :             END DO
     881           0 :             ALLOCATE (mixed_cdft%dest_list(ntargets))
     882           0 :             ALLOCATE (mixed_cdft%dest_list_bo(2, ntargets))
     883             :             ! now determine the actual grid points to send
     884           0 :             j = 1
     885           0 :             DO i = offset, offset + ntargets - 1
     886           0 :                mixed_cdft%dest_list(j) = i
     887             :                mixed_cdft%dest_list_bo(:, j) = (/bo_mixed(1, 1) + (bounds(i, 1) - bo_mixed(1, 1)), &
     888           0 :                                                  bo_mixed(2, 1) + (bounds(i, 2) - bo_mixed(2, 1))/)
     889           0 :                j = j + 1
     890             :             END DO
     891           0 :             ALLOCATE (mixed_cdft%dest_list_save(ntargets), mixed_cdft%dest_bo_save(2, ntargets))
     892             :             ! We need to store backups of these arrays since they might get reallocated during dlb
     893           0 :             mixed_cdft%dest_list_save = mixed_cdft%dest_list
     894           0 :             mixed_cdft%dest_bo_save = mixed_cdft%dest_list_bo
     895             :             ! finally determine which procs will send me grid points
     896             :             ! now we need info about y dir also
     897           0 :             DEALLOCATE (bounds)
     898           0 :             ALLOCATE (bounds(0:pw_pools(1)%pool%pw_grid%para%group_size - 1, 1:4))
     899           0 :             DO i = 0, pw_pools(1)%pool%pw_grid%para%group_size - 1
     900           0 :                bounds(i, 1:2) = pw_pools(1)%pool%pw_grid%para%bo(1:2, 1, i, 1)
     901           0 :                bounds(i, 3:4) = pw_pools(1)%pool%pw_grid%para%bo(1:2, 2, i, 1)
     902           0 :                bounds(i, 1:2) = bounds(i, 1:2) - pw_pools(1)%pool%pw_grid%npts(1)/2 - 1
     903           0 :                bounds(i, 3:4) = bounds(i, 3:4) - pw_pools(1)%pool%pw_grid%npts(2)/2 - 1
     904             :             END DO
     905           0 :             ntargets = 0
     906           0 :             offset = -1
     907           0 :             DO i = 0, pw_pools(1)%pool%pw_grid%para%group_size - 1
     908           0 :                IF ((bo(1, 1) .GE. bounds(i, 1) .AND. bo(1, 1) .LE. bounds(i, 2)) .OR. &
     909           0 :                    (bo(2, 1) .GE. bounds(i, 1) .AND. bo(2, 1) .LE. bounds(i, 2))) THEN
     910           0 :                   ntargets = ntargets + 1
     911           0 :                   IF (offset == -1) offset = i
     912           0 :                ELSE IF (bo(2, 1) .LT. bounds(i, 1)) THEN
     913             :                   EXIT
     914             :                ELSE
     915           0 :                   CYCLE
     916             :                END IF
     917             :             END DO
     918           0 :             ALLOCATE (mixed_cdft%source_list(ntargets))
     919           0 :             ALLOCATE (mixed_cdft%source_list_bo(4, ntargets))
     920           0 :             j = 1
     921           0 :             DO i = offset, offset + ntargets - 1
     922           0 :                mixed_cdft%source_list(j) = i
     923           0 :                IF (bo(1, 1) .GE. bounds(i, 1) .AND. bo(2, 1) .LE. bounds(i, 2)) THEN
     924             :                   mixed_cdft%source_list_bo(:, j) = (/bo(1, 1), bo(2, 1), &
     925           0 :                                                       bounds(i, 3), bounds(i, 4)/)
     926           0 :                ELSE IF (bo(1, 1) .GE. bounds(i, 1) .AND. bo(1, 1) .LE. bounds(i, 2)) THEN
     927             :                   mixed_cdft%source_list_bo(:, j) = (/bo(1, 1), bounds(i, 2), &
     928           0 :                                                       bounds(i, 3), bounds(i, 4)/)
     929             :                ELSE
     930             :                   mixed_cdft%source_list_bo(:, j) = (/bounds(i, 1), bo(2, 1), &
     931           0 :                                                       bounds(i, 3), bounds(i, 4)/)
     932             :                END IF
     933           0 :                j = j + 1
     934             :             END DO
     935           0 :             ALLOCATE (mixed_cdft%source_list_save(ntargets), mixed_cdft%source_bo_save(4, ntargets))
     936             :             ! We need to store backups of these arrays since they might get reallocated during dlb
     937           0 :             mixed_cdft%source_list_save = mixed_cdft%source_list
     938           0 :             mixed_cdft%source_bo_save = mixed_cdft%source_list_bo
     939           0 :             DEALLOCATE (bounds)
     940             :          END IF
     941             :       ELSE
     942             :          ! Create loggers to redirect the output of all CDFT states to different files
     943             :          ! even when the states are treated in serial (the initial print of QS data [basis set etc] for
     944             :          ! all states unfortunately goes to the first log file)
     945          50 :          CALL force_env_get(force_env, root_section=root_section)
     946         224 :          ALLOCATE (mixed_cdft%sub_logger(nforce_eval - 1))
     947         124 :          DO i = 1, nforce_eval - 1
     948          74 :             IF (force_env%para_env%is_source()) THEN
     949             :                CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", &
     950          37 :                                          c_val=input_file_path)
     951          37 :                lp = LEN_TRIM(input_file_path)
     952          37 :                input_file_path(lp + 1:LEN(input_file_path)) = "-r-"//ADJUSTL(cp_to_string(i + 1))
     953          37 :                lp = LEN_TRIM(input_file_path)
     954          37 :                output_file_path = input_file_path(1:lp)//".out"
     955             :                CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
     956             :                               file_action="WRITE", file_position="APPEND", &
     957          37 :                               unit_number=unit_nr)
     958             :             ELSE
     959          37 :                unit_nr = -1
     960             :             END IF
     961             :             CALL cp_logger_create(mixed_cdft%sub_logger(i)%p, &
     962             :                                   para_env=force_env%para_env, &
     963             :                                   default_global_unit_nr=unit_nr, &
     964          74 :                                   close_global_unit_on_dealloc=.FALSE.)
     965             :             ! Try to use better names for the local log if it is not too late
     966             :             CALL section_vals_val_get(root_section, "GLOBAL%OUTPUT_FILE_NAME", &
     967          74 :                                       c_val=c_val)
     968          74 :             IF (c_val /= "") THEN
     969             :                CALL cp_logger_set(mixed_cdft%sub_logger(i)%p, &
     970           0 :                                   local_filename=TRIM(c_val)//"_localLog")
     971             :             END IF
     972          74 :             CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
     973          74 :             IF (c_val /= "") THEN
     974             :                CALL cp_logger_set(mixed_cdft%sub_logger(i)%p, &
     975          74 :                                   local_filename=TRIM(c_val)//"_localLog")
     976             :             END IF
     977          74 :             mixed_cdft%sub_logger(i)%p%iter_info%project_name = c_val
     978             :             CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
     979         124 :                                       i_val=mixed_cdft%sub_logger(i)%p%iter_info%print_level)
     980             :          END DO
     981          50 :          IF (mixed_cdft%wfn_overlap_method) THEN
     982             :             ! qs_kind_set
     983           6 :             NULLIFY (root_section, force_env_section, force_env_sections, rs_grid_section)
     984           6 :             CALL force_env_get(force_env, root_section=root_section)
     985           6 :             force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
     986           6 :             CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, n_force_eval)
     987             :             CALL section_vals_duplicate(force_env_sections, force_env_section, &
     988           6 :                                         i_force_eval(2), i_force_eval(2))
     989             :             subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
     990           6 :                                                          i_rep_section=i_force_eval(1))
     991           6 :             kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
     992           6 :             NULLIFY (qs_kind_set)
     993           6 :             CALL cp_subsys_get(subsys_mix, atomic_kind_set=atomic_kind_set)
     994             :             CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, &
     995           6 :                                     force_env%para_env, force_env_section)
     996           6 :             mixed_cdft%qs_kind_set => qs_kind_set
     997           6 :             DEALLOCATE (i_force_eval)
     998           6 :             CALL section_vals_release(force_env_section)
     999           6 :             mixed_cdft%qs_kind_set => qs_kind_set
    1000             :          END IF
    1001             :          CALL force_env_get(force_env=force_env, &
    1002          50 :                             force_env_section=force_env_section)
    1003             :       END IF
    1004             :       ! Deallocate settings temporaries
    1005          72 :       DEALLOCATE (settings%grid_span)
    1006          72 :       DEALLOCATE (settings%npts)
    1007          72 :       DEALLOCATE (settings%spherical)
    1008          72 :       DEALLOCATE (settings%rs_dims)
    1009          72 :       DEALLOCATE (settings%odd)
    1010          72 :       DEALLOCATE (settings%atoms)
    1011          72 :       IF (mixed_cdft%run_type == mixed_cdft_parallel) &
    1012          22 :          DEALLOCATE (settings%coeffs)
    1013          72 :       DEALLOCATE (settings%cutoffs)
    1014          72 :       DEALLOCATE (settings%radii)
    1015          72 :       DEALLOCATE (settings%si)
    1016          72 :       DEALLOCATE (settings%sr)
    1017          72 :       DEALLOCATE (settings%sb)
    1018          72 :       DEALLOCATE (settings%cutoff)
    1019          72 :       DEALLOCATE (settings%rel_cutoff)
    1020             :       ! Setup mixed blacs_env for redistributing arrays during ET coupling calculation
    1021          72 :       IF (mixed_env%do_mixed_et) THEN
    1022          72 :          NULLIFY (root_section)
    1023          72 :          CALL force_env_get(force_env, globenv=globenv, root_section=root_section)
    1024             :          CALL cp_blacs_env_create(mixed_cdft%blacs_env, force_env%para_env, globenv%blacs_grid_layout, &
    1025          72 :                                   globenv%blacs_repeatable)
    1026             :       END IF
    1027             :       CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
    1028          72 :                                         "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
    1029             : 
    1030         360 :    END SUBROUTINE mixed_cdft_init_structures
    1031             : 
    1032             : ! **************************************************************************************************
    1033             : !> \brief Redistribute arrays needed for an ET coupling calculation from individual CDFT states to
    1034             : !>        the mixed CDFT env, that is, move the arrays to the correct blacs context. For parallel
    1035             : !>        simulations, the array processor distributions also change from N to 2N processors.
    1036             : !> \param force_env the force_env that holds the CDFT states
    1037             : !> \par History
    1038             : !>       01.2017  created [Nico Holmberg]
    1039             : ! **************************************************************************************************
    1040          94 :    SUBROUTINE mixed_cdft_redistribute_arrays(force_env)
    1041             :       TYPE(force_env_type), POINTER                      :: force_env
    1042             : 
    1043             :       INTEGER                                            :: iforce_eval, ispin, ivar, ncol_overlap, &
    1044             :                                                             ncol_wmat, nforce_eval, nrow_overlap, &
    1045             :                                                             nrow_wmat, nspins, nvar
    1046          94 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: ncol_mo, nrow_mo
    1047             :       LOGICAL                                            :: uniform_occupation
    1048          94 :       LOGICAL, ALLOCATABLE, DIMENSION(:)                 :: has_occupation_numbers
    1049          94 :       TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: occno_tmp
    1050             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1051             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_mo, fm_struct_overlap, &
    1052             :                                                             fm_struct_tmp, fm_struct_wmat
    1053             :       TYPE(cp_fm_type)                                   :: matrix_s_tmp, mixed_matrix_s_tmp
    1054          94 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :)     :: matrix_p_tmp, mixed_matrix_p_tmp, &
    1055          94 :                                                             mixed_wmat_tmp, mo_coeff_tmp, wmat_tmp
    1056          94 :       TYPE(cp_fm_type), DIMENSION(:, :), POINTER         :: mixed_mo_coeff
    1057          94 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: density_matrix, w_matrix
    1058             :       TYPE(dbcsr_type)                                   :: desymm_tmp
    1059             :       TYPE(dbcsr_type), POINTER                          :: mixed_matrix_s
    1060             :       TYPE(dft_control_type), POINTER                    :: dft_control
    1061             :       TYPE(force_env_type), POINTER                      :: force_env_qs
    1062             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
    1063             :       TYPE(mixed_environment_type), POINTER              :: mixed_env
    1064             :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1065             : 
    1066          94 :       NULLIFY (mixed_env, mixed_cdft, qs_env, dft_control, fm_struct_mo, &
    1067          94 :                fm_struct_wmat, fm_struct_overlap, fm_struct_tmp, &
    1068          94 :                mixed_mo_coeff, mixed_matrix_s, density_matrix, blacs_env, w_matrix, force_env_qs)
    1069           0 :       CPASSERT(ASSOCIATED(force_env))
    1070          94 :       mixed_env => force_env%mixed_env
    1071          94 :       nforce_eval = SIZE(force_env%sub_force_env)
    1072          94 :       CALL get_mixed_env(mixed_env, cdft_control=mixed_cdft)
    1073          94 :       CPASSERT(ASSOCIATED(mixed_cdft))
    1074          94 :       CALL mixed_cdft_work_type_init(mixed_cdft%matrix)
    1075             :       ! Get nspins and query for non-uniform occupation numbers
    1076         282 :       ALLOCATE (has_occupation_numbers(nforce_eval))
    1077         306 :       has_occupation_numbers = .FALSE.
    1078         306 :       DO iforce_eval = 1, nforce_eval
    1079         212 :          IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
    1080         176 :          force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
    1081         176 :          IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
    1082          24 :             qs_env => force_env_qs%qmmm_env%qs_env
    1083             :          ELSE
    1084         152 :             CALL force_env_get(force_env_qs, qs_env=qs_env)
    1085             :          END IF
    1086         176 :          CALL get_qs_env(qs_env, dft_control=dft_control)
    1087         176 :          CPASSERT(ASSOCIATED(dft_control))
    1088         176 :          nspins = dft_control%nspins
    1089         176 :          IF (force_env_qs%para_env%is_source()) &
    1090         236 :             has_occupation_numbers(iforce_eval) = ALLOCATED(dft_control%qs_control%cdft_control%occupations)
    1091             :       END DO
    1092          94 :       CALL force_env%para_env%sum(has_occupation_numbers(1))
    1093         212 :       DO iforce_eval = 2, nforce_eval
    1094         118 :          CALL force_env%para_env%sum(has_occupation_numbers(iforce_eval))
    1095         118 :          IF (has_occupation_numbers(1) .NEQV. has_occupation_numbers(iforce_eval)) &
    1096             :             CALL cp_abort(__LOCATION__, &
    1097          94 :                           "Mixing of uniform and non-uniform occupations is not allowed.")
    1098             :       END DO
    1099          94 :       uniform_occupation = .NOT. has_occupation_numbers(1)
    1100          94 :       DEALLOCATE (has_occupation_numbers)
    1101             :       ! Get number of weight functions per state as well as the type of each constraint
    1102          94 :       nvar = SIZE(dft_control%qs_control%cdft_control%target)
    1103          94 :       IF (.NOT. ALLOCATED(mixed_cdft%constraint_type)) THEN
    1104         288 :          ALLOCATE (mixed_cdft%constraint_type(nvar, nforce_eval))
    1105         412 :          mixed_cdft%constraint_type(:, :) = 0
    1106          72 :          IF (mixed_cdft%identical_constraints) THEN
    1107         142 :             DO ivar = 1, nvar
    1108             :                mixed_cdft%constraint_type(ivar, :) = &
    1109         310 :                   dft_control%qs_control%cdft_control%group(ivar)%constraint_type
    1110             :             END DO
    1111             :          ELSE
    1112             :             ! Possibly couple spin and charge constraints
    1113           6 :             DO iforce_eval = 1, nforce_eval
    1114           4 :                IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
    1115           4 :                IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
    1116           0 :                   qs_env => force_env%sub_force_env(iforce_eval)%force_env%qmmm_env%qs_env
    1117             :                ELSE
    1118           4 :                   CALL force_env_get(force_env%sub_force_env(iforce_eval)%force_env, qs_env=qs_env)
    1119             :                END IF
    1120           4 :                CALL get_qs_env(qs_env, dft_control=dft_control)
    1121           6 :                IF (force_env%sub_force_env(iforce_eval)%force_env%para_env%is_source()) THEN
    1122           4 :                   DO ivar = 1, nvar
    1123             :                      mixed_cdft%constraint_type(ivar, iforce_eval) = &
    1124           4 :                         dft_control%qs_control%cdft_control%group(ivar)%constraint_type
    1125             :                   END DO
    1126             :                END IF
    1127             :             END DO
    1128           2 :             CALL force_env%para_env%sum(mixed_cdft%constraint_type)
    1129             :          END IF
    1130             :       END IF
    1131             :       ! Transfer data from sub_force_envs to temporaries
    1132         988 :       ALLOCATE (mixed_cdft%matrix%mixed_mo_coeff(nforce_eval, nspins))
    1133          94 :       mixed_mo_coeff => mixed_cdft%matrix%mixed_mo_coeff
    1134         688 :       ALLOCATE (mixed_cdft%matrix%w_matrix(nforce_eval, nvar))
    1135          94 :       w_matrix => mixed_cdft%matrix%w_matrix
    1136          94 :       CALL dbcsr_init_p(mixed_cdft%matrix%mixed_matrix_s)
    1137          94 :       mixed_matrix_s => mixed_cdft%matrix%mixed_matrix_s
    1138          94 :       IF (mixed_cdft%calculate_metric) THEN
    1139         144 :          ALLOCATE (mixed_cdft%matrix%density_matrix(nforce_eval, nspins))
    1140          14 :          density_matrix => mixed_cdft%matrix%density_matrix
    1141             :       END IF
    1142        1582 :       ALLOCATE (mo_coeff_tmp(nforce_eval, nspins), wmat_tmp(nforce_eval, nvar))
    1143         470 :       ALLOCATE (nrow_mo(nspins), ncol_mo(nspins))
    1144         224 :       IF (mixed_cdft%calculate_metric) ALLOCATE (matrix_p_tmp(nforce_eval, nspins))
    1145          94 :       IF (.NOT. uniform_occupation) THEN
    1146         140 :          ALLOCATE (mixed_cdft%occupations(nforce_eval, nspins))
    1147         140 :          ALLOCATE (occno_tmp(nforce_eval, nspins))
    1148             :       END IF
    1149         306 :       DO iforce_eval = 1, nforce_eval
    1150             :          ! Temporary arrays need to be nulled on every process
    1151         636 :          DO ispin = 1, nspins
    1152             :             ! Valgrind 3.12/gfortran 4.8.4 oddly complains here (unconditional jump)
    1153             :             ! if mixed_cdft%calculate_metric = .FALSE. and the need to null the array
    1154             :             ! is queried with IF (mixed_cdft%calculate_metric) &
    1155         424 :             IF (.NOT. uniform_occupation) &
    1156         268 :                NULLIFY (occno_tmp(iforce_eval, ispin)%array)
    1157             :          END DO
    1158         212 :          IF (.NOT. ASSOCIATED(force_env%sub_force_env(iforce_eval)%force_env)) CYCLE
    1159             :          ! From this point onward, we access data local to the sub_force_envs
    1160             :          ! Get qs_env
    1161         176 :          force_env_qs => force_env%sub_force_env(iforce_eval)%force_env
    1162         176 :          IF (force_env%mixed_env%do_mixed_qmmm_cdft) THEN
    1163          24 :             qs_env => force_env_qs%qmmm_env%qs_env
    1164             :          ELSE
    1165         152 :             CALL force_env_get(force_env_qs, qs_env=qs_env)
    1166             :          END IF
    1167         176 :          CALL get_qs_env(qs_env, dft_control=dft_control, blacs_env=blacs_env)
    1168             :          ! Store dimensions of the transferred arrays
    1169             :          CALL dbcsr_get_info(dft_control%qs_control%cdft_control%matrix_s%matrix, &
    1170         176 :                              nfullrows_total=nrow_overlap, nfullcols_total=ncol_overlap)
    1171             :          CALL dbcsr_get_info(dft_control%qs_control%cdft_control%wmat(1)%matrix, &
    1172         176 :                              nfullrows_total=nrow_wmat, nfullcols_total=ncol_wmat)
    1173             :          ! MO Coefficients
    1174         528 :          DO ispin = 1, nspins
    1175             :             CALL cp_fm_get_info(dft_control%qs_control%cdft_control%mo_coeff(ispin), &
    1176         352 :                                 ncol_global=ncol_mo(ispin), nrow_global=nrow_mo(ispin))
    1177             :             CALL cp_fm_create(matrix=mo_coeff_tmp(iforce_eval, ispin), &
    1178             :                               matrix_struct=dft_control%qs_control%cdft_control%mo_coeff(ispin)%matrix_struct, &
    1179             :                               name="MO_COEFF_"//TRIM(ADJUSTL(cp_to_string(iforce_eval)))//"_" &
    1180         352 :                               //TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
    1181             :             CALL cp_fm_to_fm(dft_control%qs_control%cdft_control%mo_coeff(ispin), &
    1182         528 :                              mo_coeff_tmp(iforce_eval, ispin))
    1183             :          END DO
    1184         176 :          CALL cp_fm_release(dft_control%qs_control%cdft_control%mo_coeff)
    1185             :          ! Matrix representation(s) of the weight function(s) (dbcsr -> fm)
    1186             :          CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_wmat, ncol_global=ncol_wmat, context=blacs_env, &
    1187             :                                   para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env, &
    1188         176 :                                   square_blocks=.TRUE.)
    1189         356 :          DO ivar = 1, nvar
    1190         180 :             CALL cp_fm_create(wmat_tmp(iforce_eval, ivar), fm_struct_tmp, name="w_matrix")
    1191         180 :             CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%wmat(ivar)%matrix, desymm_tmp)
    1192         180 :             CALL copy_dbcsr_to_fm(desymm_tmp, wmat_tmp(iforce_eval, ivar))
    1193         180 :             CALL dbcsr_release(desymm_tmp)
    1194         356 :             CALL dbcsr_release_p(dft_control%qs_control%cdft_control%wmat(ivar)%matrix)
    1195             :          END DO
    1196         176 :          DEALLOCATE (dft_control%qs_control%cdft_control%wmat)
    1197         176 :          CALL cp_fm_struct_release(fm_struct_tmp)
    1198             :          ! Overlap matrix is the same for all sub_force_envs, so we just copy the first one (dbcsr -> fm)
    1199         176 :          IF (iforce_eval == 1) THEN
    1200             :             CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_overlap, &
    1201             :                                      ncol_global=ncol_overlap, context=blacs_env, &
    1202          76 :                                      para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env)
    1203          76 :             CALL cp_fm_create(matrix_s_tmp, fm_struct_tmp, name="s_matrix")
    1204          76 :             CALL cp_fm_struct_release(fm_struct_tmp)
    1205          76 :             CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%matrix_s%matrix, desymm_tmp)
    1206          76 :             CALL copy_dbcsr_to_fm(desymm_tmp, matrix_s_tmp)
    1207          76 :             CALL dbcsr_release(desymm_tmp)
    1208             :          END IF
    1209         176 :          CALL dbcsr_release_p(dft_control%qs_control%cdft_control%matrix_s%matrix)
    1210             :          ! Density_matrix (dbcsr -> fm)
    1211         176 :          IF (mixed_cdft%calculate_metric) THEN
    1212          72 :             DO ispin = 1, nspins
    1213             :                ! Size AOxAO
    1214             :                CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=ncol_overlap, &
    1215             :                                         ncol_global=ncol_overlap, context=blacs_env, &
    1216          48 :                                         para_env=force_env%sub_force_env(iforce_eval)%force_env%para_env)
    1217          48 :                CALL cp_fm_create(matrix_p_tmp(iforce_eval, ispin), fm_struct_tmp, name="dm_matrix")
    1218          48 :                CALL cp_fm_struct_release(fm_struct_tmp)
    1219          48 :                CALL dbcsr_desymmetrize(dft_control%qs_control%cdft_control%matrix_p(ispin)%matrix, desymm_tmp)
    1220          48 :                CALL copy_dbcsr_to_fm(desymm_tmp, matrix_p_tmp(iforce_eval, ispin))
    1221          48 :                CALL dbcsr_release(desymm_tmp)
    1222          72 :                CALL dbcsr_release_p(dft_control%qs_control%cdft_control%matrix_p(ispin)%matrix)
    1223             :             END DO
    1224          24 :             DEALLOCATE (dft_control%qs_control%cdft_control%matrix_p)
    1225             :          END IF
    1226             :          ! Occupation numbers
    1227         446 :          IF (.NOT. uniform_occupation) THEN
    1228          84 :             DO ispin = 1, nspins
    1229          56 :                IF (ncol_mo(ispin) /= SIZE(dft_control%qs_control%cdft_control%occupations(ispin)%array)) &
    1230           0 :                   CPABORT("Array dimensions dont match.")
    1231          56 :                IF (force_env_qs%para_env%is_source()) THEN
    1232          84 :                   ALLOCATE (occno_tmp(iforce_eval, ispin)%array(ncol_mo(ispin)))
    1233         126 :                   occno_tmp(iforce_eval, ispin)%array = dft_control%qs_control%cdft_control%occupations(ispin)%array
    1234             :                END IF
    1235          84 :                DEALLOCATE (dft_control%qs_control%cdft_control%occupations(ispin)%array)
    1236             :             END DO
    1237          28 :             DEALLOCATE (dft_control%qs_control%cdft_control%occupations)
    1238             :          END IF
    1239             :       END DO
    1240             :       ! Create needed fm structs
    1241             :       CALL cp_fm_struct_create(fm_struct_wmat, nrow_global=nrow_wmat, ncol_global=ncol_wmat, &
    1242          94 :                                context=mixed_cdft%blacs_env, para_env=force_env%para_env)
    1243             :       CALL cp_fm_struct_create(fm_struct_overlap, nrow_global=nrow_overlap, ncol_global=ncol_overlap, &
    1244          94 :                                context=mixed_cdft%blacs_env, para_env=force_env%para_env)
    1245             :       ! Redistribute arrays with copy_general (this is not optimal for dbcsr matrices but...)
    1246             :       ! We use this method for the serial case (mixed_cdft%run_type == mixed_cdft_serial) as well to move the arrays to the
    1247             :       ! correct blacs_env, which is impossible using a simple copy of the arrays
    1248         688 :       ALLOCATE (mixed_wmat_tmp(nforce_eval, nvar))
    1249          94 :       IF (mixed_cdft%calculate_metric) &
    1250         144 :          ALLOCATE (mixed_matrix_p_tmp(nforce_eval, nspins))
    1251         306 :       DO iforce_eval = 1, nforce_eval
    1252             :          ! MO coefficients
    1253         636 :          DO ispin = 1, nspins
    1254         424 :             NULLIFY (fm_struct_mo)
    1255             :             CALL cp_fm_struct_create(fm_struct_mo, nrow_global=nrow_mo(ispin), ncol_global=ncol_mo(ispin), &
    1256         424 :                                      context=mixed_cdft%blacs_env, para_env=force_env%para_env)
    1257             :             CALL cp_fm_create(matrix=mixed_mo_coeff(iforce_eval, ispin), &
    1258             :                               matrix_struct=fm_struct_mo, &
    1259             :                               name="MO_COEFF_"//TRIM(ADJUSTL(cp_to_string(iforce_eval)))//"_" &
    1260         424 :                               //TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
    1261             :             CALL cp_fm_copy_general(mo_coeff_tmp(iforce_eval, ispin), &
    1262             :                                     mixed_mo_coeff(iforce_eval, ispin), &
    1263         424 :                                     mixed_cdft%blacs_env%para_env)
    1264         424 :             CALL cp_fm_release(mo_coeff_tmp(iforce_eval, ispin))
    1265         636 :             CALL cp_fm_struct_release(fm_struct_mo)
    1266             :          END DO
    1267             :          ! Weight
    1268         428 :          DO ivar = 1, nvar
    1269         216 :             NULLIFY (w_matrix(iforce_eval, ivar)%matrix)
    1270         216 :             CALL dbcsr_init_p(w_matrix(iforce_eval, ivar)%matrix)
    1271             :             CALL cp_fm_create(matrix=mixed_wmat_tmp(iforce_eval, ivar), &
    1272             :                               matrix_struct=fm_struct_wmat, &
    1273         216 :                               name="WEIGHT_"//TRIM(ADJUSTL(cp_to_string(iforce_eval)))//"_MATRIX")
    1274             :             CALL cp_fm_copy_general(wmat_tmp(iforce_eval, ivar), &
    1275             :                                     mixed_wmat_tmp(iforce_eval, ivar), &
    1276         216 :                                     mixed_cdft%blacs_env%para_env)
    1277         216 :             CALL cp_fm_release(wmat_tmp(iforce_eval, ivar))
    1278             :             ! (fm -> dbcsr)
    1279             :             CALL copy_fm_to_dbcsr_bc(mixed_wmat_tmp(iforce_eval, ivar), &
    1280         216 :                                      w_matrix(iforce_eval, ivar)%matrix)
    1281         428 :             CALL cp_fm_release(mixed_wmat_tmp(iforce_eval, ivar))
    1282             :          END DO
    1283             :          ! Density matrix (fm -> dbcsr)
    1284         306 :          IF (mixed_cdft%calculate_metric) THEN
    1285          90 :             DO ispin = 1, nspins
    1286          60 :                NULLIFY (density_matrix(iforce_eval, ispin)%matrix)
    1287          60 :                CALL dbcsr_init_p(density_matrix(iforce_eval, ispin)%matrix)
    1288             :                CALL cp_fm_create(matrix=mixed_matrix_p_tmp(iforce_eval, ispin), &
    1289             :                                  matrix_struct=fm_struct_overlap, &
    1290             :                                  name="DENSITY_"//TRIM(ADJUSTL(cp_to_string(iforce_eval)))//"_"// &
    1291          60 :                                  TRIM(ADJUSTL(cp_to_string(ispin)))//"_MATRIX")
    1292             :                CALL cp_fm_copy_general(matrix_p_tmp(iforce_eval, ispin), &
    1293             :                                        mixed_matrix_p_tmp(iforce_eval, ispin), &
    1294          60 :                                        mixed_cdft%blacs_env%para_env)
    1295          60 :                CALL cp_fm_release(matrix_p_tmp(iforce_eval, ispin))
    1296             :                CALL copy_fm_to_dbcsr_bc(mixed_matrix_p_tmp(iforce_eval, ispin), &
    1297          60 :                                         density_matrix(iforce_eval, ispin)%matrix)
    1298          90 :                CALL cp_fm_release(mixed_matrix_p_tmp(iforce_eval, ispin))
    1299             :             END DO
    1300             :          END IF
    1301             :       END DO
    1302          94 :       CALL cp_fm_struct_release(fm_struct_wmat)
    1303          94 :       DEALLOCATE (mo_coeff_tmp, wmat_tmp, mixed_wmat_tmp)
    1304          94 :       IF (mixed_cdft%calculate_metric) THEN
    1305          14 :          DEALLOCATE (matrix_p_tmp)
    1306          14 :          DEALLOCATE (mixed_matrix_p_tmp)
    1307             :       END IF
    1308             :       ! Overlap (fm -> dbcsr)
    1309             :       CALL cp_fm_create(matrix=mixed_matrix_s_tmp, &
    1310             :                         matrix_struct=fm_struct_overlap, &
    1311          94 :                         name="OVERLAP_MATRIX")
    1312          94 :       CALL cp_fm_struct_release(fm_struct_overlap)
    1313             :       CALL cp_fm_copy_general(matrix_s_tmp, &
    1314             :                               mixed_matrix_s_tmp, &
    1315          94 :                               mixed_cdft%blacs_env%para_env)
    1316          94 :       CALL cp_fm_release(matrix_s_tmp)
    1317          94 :       CALL copy_fm_to_dbcsr_bc(mixed_matrix_s_tmp, mixed_matrix_s)
    1318          94 :       CALL cp_fm_release(mixed_matrix_s_tmp)
    1319             :       ! Occupation numbers
    1320          94 :       IF (.NOT. uniform_occupation) THEN
    1321          42 :          DO iforce_eval = 1, nforce_eval
    1322          98 :             DO ispin = 1, nspins
    1323         168 :                ALLOCATE (mixed_cdft%occupations(iforce_eval, ispin)%array(ncol_mo(ispin)))
    1324         252 :                mixed_cdft%occupations(iforce_eval, ispin)%array = 0.0_dp
    1325          56 :                IF (ASSOCIATED(occno_tmp(iforce_eval, ispin)%array)) THEN
    1326         126 :                   mixed_cdft%occupations(iforce_eval, ispin)%array = occno_tmp(iforce_eval, ispin)%array
    1327          28 :                   DEALLOCATE (occno_tmp(iforce_eval, ispin)%array)
    1328             :                END IF
    1329         476 :                CALL force_env%para_env%sum(mixed_cdft%occupations(iforce_eval, ispin)%array)
    1330             :             END DO
    1331             :          END DO
    1332          14 :          DEALLOCATE (occno_tmp)
    1333             :       END IF
    1334          94 :       DEALLOCATE (ncol_mo, nrow_mo)
    1335             : 
    1336         188 :    END SUBROUTINE mixed_cdft_redistribute_arrays
    1337             : ! **************************************************************************************************
    1338             : !> \brief Routine to print out the electronic coupling(s) between CDFT states.
    1339             : !> \param force_env the force_env that holds the CDFT states
    1340             : !> \par History
    1341             : !>       11.17  created [Nico Holmberg]
    1342             : ! **************************************************************************************************
    1343          94 :    SUBROUTINE mixed_cdft_print_couplings(force_env)
    1344             :       TYPE(force_env_type), POINTER                      :: force_env
    1345             : 
    1346             :       INTEGER                                            :: iounit, ipermutation, istate, ivar, &
    1347             :                                                             jstate, nforce_eval, npermutations, &
    1348             :                                                             nvar
    1349             :       TYPE(cp_logger_type), POINTER                      :: logger
    1350             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
    1351             :       TYPE(section_vals_type), POINTER                   :: force_env_section, print_section
    1352             : 
    1353          94 :       NULLIFY (print_section, mixed_cdft)
    1354             : 
    1355          94 :       logger => cp_get_default_logger()
    1356          94 :       CPASSERT(ASSOCIATED(force_env))
    1357             :       CALL force_env_get(force_env=force_env, &
    1358          94 :                          force_env_section=force_env_section)
    1359          94 :       CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
    1360          94 :       CPASSERT(ASSOCIATED(mixed_cdft))
    1361          94 :       print_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
    1362          94 :       iounit = cp_print_key_unit_nr(logger, print_section, '', extension='.mixedLog')
    1363             :       !
    1364          94 :       CPASSERT(ALLOCATED(mixed_cdft%results%strength))
    1365          94 :       CPASSERT(ALLOCATED(mixed_cdft%results%W_diagonal))
    1366          94 :       CPASSERT(ALLOCATED(mixed_cdft%results%S))
    1367          94 :       CPASSERT(ALLOCATED(mixed_cdft%results%energy))
    1368          94 :       nforce_eval = SIZE(force_env%sub_force_env)
    1369          94 :       nvar = SIZE(mixed_cdft%results%strength, 1)
    1370          94 :       npermutations = nforce_eval*(nforce_eval - 1)/2 ! Size of upper triangular part
    1371          94 :       IF (iounit > 0) THEN
    1372             :          WRITE (iounit, '(/,T3,A,T66)') &
    1373          47 :             '------------------------- CDFT coupling information --------------------------'
    1374             :          WRITE (iounit, '(T3,A,T66,(3X,F12.2))') &
    1375          47 :             'Information at step (fs):', mixed_cdft%sim_step*mixed_cdft%sim_dt
    1376         135 :          DO ipermutation = 1, npermutations
    1377          88 :             CALL map_permutation_to_states(nforce_eval, ipermutation, istate, jstate)
    1378          88 :             WRITE (iounit, '(/,T3,A)') REPEAT('#', 44)
    1379          88 :             WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### CDFT states I =', istate, ' and J = ', jstate, ' ######'
    1380          88 :             WRITE (iounit, '(T3,A)') REPEAT('#', 44)
    1381         177 :             DO ivar = 1, nvar
    1382          89 :                IF (ivar > 1) &
    1383           1 :                   WRITE (iounit, '(A)') ''
    1384          89 :                WRITE (iounit, '(T3,A,T60,(3X,I18))') 'Atomic group:', ivar
    1385             :                WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1386          89 :                   'Strength of constraint I:', mixed_cdft%results%strength(ivar, istate)
    1387             :                WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1388          89 :                   'Strength of constraint J:', mixed_cdft%results%strength(ivar, jstate)
    1389             :                WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1390          89 :                   'Final value of constraint I:', mixed_cdft%results%W_diagonal(ivar, istate)
    1391             :                WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1392         177 :                   'Final value of constraint J:', mixed_cdft%results%W_diagonal(ivar, jstate)
    1393             :             END DO
    1394             :             WRITE (iounit, '(/,T3,A,T60,(3X,F18.12))') &
    1395          88 :                'Overlap between states I and J:', mixed_cdft%results%S(istate, jstate)
    1396             :             WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1397          88 :                'Charge transfer energy (J-I) (Hartree):', (mixed_cdft%results%energy(jstate) - mixed_cdft%results%energy(istate))
    1398          88 :             WRITE (iounit, *)
    1399          88 :             IF (ALLOCATED(mixed_cdft%results%rotation)) THEN
    1400          86 :                IF (ABS(mixed_cdft%results%rotation(ipermutation))*1.0E3_dp .GE. 0.1_dp) THEN
    1401             :                   WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1402          84 :                      'Diabatic electronic coupling (rotation, mHartree):', &
    1403         168 :                      ABS(mixed_cdft%results%rotation(ipermutation)*1.0E3_dp)
    1404             :                ELSE
    1405             :                   WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1406           2 :                      'Diabatic electronic coupling (rotation, microHartree):', &
    1407           4 :                      ABS(mixed_cdft%results%rotation(ipermutation)*1.0E6_dp)
    1408             :                END IF
    1409             :             END IF
    1410          88 :             IF (ALLOCATED(mixed_cdft%results%lowdin)) THEN
    1411          10 :                IF (ABS(mixed_cdft%results%lowdin(ipermutation))*1.0E3_dp .GE. 0.1_dp) THEN
    1412             :                   WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1413           9 :                      'Diabatic electronic coupling (Lowdin, mHartree):', &
    1414          18 :                      ABS(mixed_cdft%results%lowdin(ipermutation)*1.0E3_dp)
    1415             :                ELSE
    1416             :                   WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1417           1 :                      'Diabatic electronic coupling (Lowdin, microHartree):', &
    1418           2 :                      ABS(mixed_cdft%results%lowdin(ipermutation)*1.0E6_dp)
    1419             :                END IF
    1420             :             END IF
    1421          88 :             IF (ALLOCATED(mixed_cdft%results%wfn)) THEN
    1422           6 :                IF (mixed_cdft%results%wfn(ipermutation)*1.0E3_dp .GE. 0.1_dp) THEN
    1423             :                   WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1424           5 :                      'Diabatic electronic coupling (wfn overlap, mHartree):', &
    1425          10 :                      ABS(mixed_cdft%results%wfn(ipermutation)*1.0E3_dp)
    1426             :                ELSE
    1427             :                   WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1428           1 :                      'Diabatic electronic coupling (wfn overlap, microHartree):', &
    1429           2 :                      ABS(mixed_cdft%results%wfn(ipermutation)*1.0E6_dp)
    1430             :                END IF
    1431             :             END IF
    1432          88 :             IF (ALLOCATED(mixed_cdft%results%nonortho)) THEN
    1433             :                WRITE (iounit, '(T3,A,T60,(3X,F18.12))') &
    1434          50 :                   'Diabatic electronic coupling (nonorthogonal, Hartree):', mixed_cdft%results%nonortho(ipermutation)
    1435             :             END IF
    1436         223 :             IF (ALLOCATED(mixed_cdft%results%metric)) THEN
    1437           9 :                WRITE (iounit, *)
    1438           9 :                IF (SIZE(mixed_cdft%results%metric, 2) == 1) THEN
    1439             :                   WRITE (iounit, '(T3,A,T66,(3X,F12.6))') &
    1440           0 :                      'Coupling reliability metric (0 is ideal):', mixed_cdft%results%metric(ipermutation, 1)
    1441             :                ELSE
    1442             :                   WRITE (iounit, '(T3,A,T54,(3X,2F12.6))') &
    1443           9 :                      'Coupling reliability metric (0 is ideal):', &
    1444          18 :                      mixed_cdft%results%metric(ipermutation, 1), mixed_cdft%results%metric(ipermutation, 2)
    1445             :                END IF
    1446             :             END IF
    1447             :          END DO
    1448             :          WRITE (iounit, '(T3,A)') &
    1449          47 :             '------------------------------------------------------------------------------'
    1450             :       END IF
    1451             :       CALL cp_print_key_finished_output(iounit, logger, force_env_section, &
    1452          94 :                                         "MIXED%MIXED_CDFT%PRINT%PROGRAM_RUN_INFO")
    1453             : 
    1454          94 :    END SUBROUTINE mixed_cdft_print_couplings
    1455             : 
    1456             : ! **************************************************************************************************
    1457             : !> \brief Release storage reserved for mixed CDFT matrices
    1458             : !> \param force_env the force_env that holds the CDFT states
    1459             : !> \par History
    1460             : !>       11.17  created [Nico Holmberg]
    1461             : ! **************************************************************************************************
    1462          94 :    SUBROUTINE mixed_cdft_release_work(force_env)
    1463             :       TYPE(force_env_type), POINTER                      :: force_env
    1464             : 
    1465             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
    1466             : 
    1467          94 :       NULLIFY (mixed_cdft)
    1468             : 
    1469          94 :       CPASSERT(ASSOCIATED(force_env))
    1470          94 :       CALL get_mixed_env(force_env%mixed_env, cdft_control=mixed_cdft)
    1471          94 :       CPASSERT(ASSOCIATED(mixed_cdft))
    1472          94 :       CALL mixed_cdft_result_type_release(mixed_cdft%results)
    1473             : 
    1474          94 :    END SUBROUTINE mixed_cdft_release_work
    1475             : 
    1476             : ! **************************************************************************************************
    1477             : !> \brief Given the size of a symmetric matrix and a permutation index, returns indices (i, j) of the
    1478             : !>        off-diagonal element that corresponds to the permutation index. Assumes that the permutation
    1479             : !>        index was computed by going through the upper triangular part of the input matrix row-by-row.
    1480             : !> \param n the size of the symmetric matrix
    1481             : !> \param ipermutation the permutation index
    1482             : !> \param i the row index corresponding to ipermutation
    1483             : !> \param j the column index corresponding to ipermutation
    1484             : ! **************************************************************************************************
    1485        1248 :    SUBROUTINE map_permutation_to_states(n, ipermutation, i, j)
    1486             :       INTEGER, INTENT(IN)                                :: n, ipermutation
    1487             :       INTEGER, INTENT(OUT)                               :: i, j
    1488             : 
    1489             :       INTEGER                                            :: kcol, kpermutation, krow, npermutations
    1490             : 
    1491        1248 :       npermutations = n*(n - 1)/2 ! Size of upper triangular part
    1492        1248 :       IF (ipermutation > npermutations) &
    1493           0 :          CPABORT("Permutation index out of bounds")
    1494        1248 :       kpermutation = 0
    1495        2124 :       DO krow = 1, n
    1496        7566 :          DO kcol = krow + 1, n
    1497        6690 :             kpermutation = kpermutation + 1
    1498        7566 :             IF (kpermutation == ipermutation) THEN
    1499        1248 :                i = krow
    1500        1248 :                j = kcol
    1501        1248 :                RETURN
    1502             :             END IF
    1503             :          END DO
    1504             :       END DO
    1505             : 
    1506             :    END SUBROUTINE map_permutation_to_states
    1507             : 
    1508             : ! **************************************************************************************************
    1509             : !> \brief Determine confinement bounds along confinement dir (hardcoded to be z)
    1510             : !>        and determine the number of nonzero entries
    1511             : !>        Optionally zero entries below a given threshold
    1512             : !> \param fun input 3D potential (real space)
    1513             : !> \param th threshold for screening values
    1514             : !> \param just_zero determines if fun should only be zeroed without returning bounds/work
    1515             : !> \param bounds the confinement bounds: fun is nonzero only between these values along 3rd dimension
    1516             : !> \param work an estimate of the total number of grid points where fun is nonzero
    1517             : ! **************************************************************************************************
    1518          34 :    SUBROUTINE hfun_zero(fun, th, just_zero, bounds, work)
    1519             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT)   :: fun
    1520             :       REAL(KIND=dp), INTENT(IN)                          :: th
    1521             :       LOGICAL                                            :: just_zero
    1522             :       INTEGER, OPTIONAL                                  :: bounds(2), work
    1523             : 
    1524             :       INTEGER                                            :: i1, i2, i3, lb, n1, n2, n3, nzeroed, &
    1525             :                                                             nzeroed_total, ub
    1526             :       LOGICAL                                            :: lb_final, ub_final
    1527             : 
    1528          34 :       n1 = SIZE(fun, 1)
    1529          34 :       n2 = SIZE(fun, 2)
    1530          34 :       n3 = SIZE(fun, 3)
    1531          34 :       nzeroed_total = 0
    1532          34 :       IF (.NOT. just_zero) THEN
    1533          34 :          CPASSERT(PRESENT(bounds))
    1534          34 :          CPASSERT(PRESENT(work))
    1535             :          lb = 1
    1536             :          lb_final = .FALSE.
    1537             :          ub_final = .FALSE.
    1538             :       END IF
    1539        1586 :       DO i3 = 1, n3
    1540        1552 :          IF (.NOT. just_zero) nzeroed = 0
    1541       75920 :          DO i2 = 1, n2
    1542     1956496 :             DO i1 = 1, n1
    1543     1954944 :                IF (fun(i1, i2, i3) < th) THEN
    1544      983300 :                   IF (.NOT. just_zero) THEN
    1545      983300 :                      nzeroed = nzeroed + 1
    1546      983300 :                      nzeroed_total = nzeroed_total + 1
    1547             :                   ELSE
    1548           0 :                      fun(i1, i2, i3) = 0.0_dp
    1549             :                   END IF
    1550             :                END IF
    1551             :             END DO
    1552             :          END DO
    1553        1586 :          IF (.NOT. just_zero) THEN
    1554        1552 :             IF (nzeroed == (n2*n1)) THEN
    1555          80 :                IF (.NOT. lb_final) THEN
    1556             :                   lb = i3
    1557          56 :                ELSE IF (.NOT. ub_final) THEN
    1558           8 :                   ub = i3
    1559           8 :                   ub_final = .TRUE.
    1560             :                END IF
    1561             :             ELSE
    1562             :                IF (.NOT. lb_final) lb_final = .TRUE.
    1563             :                IF (ub_final) ub_final = .FALSE. ! Safeguard against "holes"
    1564             :             END IF
    1565             :          END IF
    1566             :       END DO
    1567          34 :       IF (.NOT. just_zero) THEN
    1568          34 :          IF (.NOT. ub_final) ub = n3
    1569          34 :          bounds(1) = lb
    1570          34 :          bounds(2) = ub
    1571         102 :          bounds = bounds - (n3/2) - 1
    1572          34 :          work = n3*n2*n1 - nzeroed_total
    1573             :       END IF
    1574             : 
    1575          34 :    END SUBROUTINE hfun_zero
    1576             : 
    1577             : ! **************************************************************************************************
    1578             : !> \brief Read input section related to block diagonalization of the mixed CDFT Hamiltonian matrix.
    1579             : !> \param force_env        the force_env that holds the CDFT states
    1580             : !> \param blocks           list of CDFT states defining the matrix blocks
    1581             : !> \param ignore_excited   flag that determines if excited states resulting from the block
    1582             : !>                         diagonalization process should be ignored
    1583             : !> \param nrecursion       integer that determines how many steps of recursive block diagonalization
    1584             : !>                         is performed (1 if disabled)
    1585             : !> \par History
    1586             : !>       01.18  created [Nico Holmberg]
    1587             : ! **************************************************************************************************
    1588           8 :    SUBROUTINE mixed_cdft_read_block_diag(force_env, blocks, ignore_excited, nrecursion)
    1589             :       TYPE(force_env_type), POINTER                      :: force_env
    1590             :       TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:), &
    1591             :          INTENT(OUT)                                     :: blocks
    1592             :       LOGICAL, INTENT(OUT)                               :: ignore_excited
    1593             :       INTEGER, INTENT(OUT)                               :: nrecursion
    1594             : 
    1595             :       INTEGER                                            :: i, j, k, l, nblk, nforce_eval
    1596           8 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
    1597             :       LOGICAL                                            :: do_recursive, explicit, has_duplicates
    1598             :       TYPE(section_vals_type), POINTER                   :: block_section, force_env_section
    1599             : 
    1600             :       EXTERNAL                                           :: dsygv
    1601             : 
    1602           8 :       NULLIFY (force_env_section, block_section)
    1603           0 :       CPASSERT(ASSOCIATED(force_env))
    1604           8 :       nforce_eval = SIZE(force_env%sub_force_env)
    1605             : 
    1606             :       CALL force_env_get(force_env=force_env, &
    1607           8 :                          force_env_section=force_env_section)
    1608           8 :       block_section => section_vals_get_subs_vals(force_env_section, "MIXED%MIXED_CDFT%BLOCK_DIAGONALIZE")
    1609             : 
    1610           8 :       CALL section_vals_get(block_section, explicit=explicit)
    1611           8 :       IF (.NOT. explicit) &
    1612             :          CALL cp_abort(__LOCATION__, &
    1613             :                        "Block diagonalization of CDFT Hamiltonian was requested, but the "// &
    1614           0 :                        "corresponding input section is missing!")
    1615             : 
    1616           8 :       CALL section_vals_val_get(block_section, "BLOCK", n_rep_val=nblk)
    1617          44 :       ALLOCATE (blocks(nblk))
    1618          28 :       DO i = 1, nblk
    1619          20 :          NULLIFY (blocks(i)%array)
    1620          20 :          CALL section_vals_val_get(block_section, "BLOCK", i_rep_val=i, i_vals=tmplist)
    1621          20 :          IF (SIZE(tmplist) < 1) &
    1622           0 :             CPABORT("Each BLOCK must contain at least 1 state.")
    1623          60 :          ALLOCATE (blocks(i)%array(SIZE(tmplist)))
    1624          66 :          blocks(i)%array(:) = tmplist(:)
    1625             :       END DO
    1626           8 :       CALL section_vals_val_get(block_section, "IGNORE_EXCITED", l_val=ignore_excited)
    1627           8 :       CALL section_vals_val_get(block_section, "RECURSIVE_DIAGONALIZATION", l_val=do_recursive)
    1628             :       ! Check that the requested states exist
    1629          28 :       DO i = 1, nblk
    1630          66 :          DO j = 1, SIZE(blocks(i)%array)
    1631          38 :             IF (blocks(i)%array(j) < 1 .OR. blocks(i)%array(j) > nforce_eval) &
    1632          20 :                CPABORT("Requested state does not exist.")
    1633             :          END DO
    1634             :       END DO
    1635             :       ! Check for duplicates
    1636           8 :       has_duplicates = .FALSE.
    1637          28 :       DO i = 1, nblk
    1638             :          ! Within same block
    1639          58 :          DO j = 1, SIZE(blocks(i)%array)
    1640          76 :             DO k = j + 1, SIZE(blocks(i)%array)
    1641          56 :                IF (blocks(i)%array(j) == blocks(i)%array(k)) has_duplicates = .TRUE.
    1642             :             END DO
    1643             :          END DO
    1644             :          ! Within different blocks
    1645          46 :          DO j = i + 1, nblk
    1646          74 :             DO k = 1, SIZE(blocks(i)%array)
    1647         122 :                DO l = 1, SIZE(blocks(j)%array)
    1648         104 :                   IF (blocks(i)%array(k) == blocks(j)%array(l)) has_duplicates = .TRUE.
    1649             :                END DO
    1650             :             END DO
    1651             :          END DO
    1652             :       END DO
    1653           8 :       IF (has_duplicates) CPABORT("Duplicate states are not allowed.")
    1654           8 :       nrecursion = 1
    1655           8 :       IF (do_recursive) THEN
    1656           2 :          IF (MODULO(nblk, 2) /= 0) THEN
    1657             :             CALL cp_warn(__LOCATION__, &
    1658             :                          "Number of blocks not divisible with 2. Recursive diagonalization not possible. "// &
    1659           0 :                          "Calculation proceeds without.")
    1660           0 :             nrecursion = 1
    1661             :          ELSE
    1662           2 :             nrecursion = nblk/2
    1663             :          END IF
    1664           2 :          IF (nrecursion /= 1 .AND. .NOT. ignore_excited) &
    1665             :             CALL cp_abort(__LOCATION__, &
    1666           0 :                           "Keyword IGNORE_EXCITED must be active for recursive diagonalization.")
    1667             :       END IF
    1668             : 
    1669          32 :    END SUBROUTINE mixed_cdft_read_block_diag
    1670             : 
    1671             : ! **************************************************************************************************
    1672             : !> \brief Assembles the matrix blocks from the mixed CDFT Hamiltonian.
    1673             : !> \param mixed_cdft the env that holds the CDFT states
    1674             : !> \param blocks  list of CDFT states defining the matrix blocks
    1675             : !> \param H_block list of Hamiltonian matrix blocks
    1676             : !> \param S_block list of overlap matrix blocks
    1677             : !> \par History
    1678             : !>       01.18  created [Nico Holmberg]
    1679             : ! **************************************************************************************************
    1680          10 :    SUBROUTINE mixed_cdft_get_blocks(mixed_cdft, blocks, H_block, S_block)
    1681             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
    1682             :       TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:)    :: blocks
    1683             :       TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:), &
    1684             :          INTENT(OUT)                                     :: H_block, S_block
    1685             : 
    1686             :       INTEGER                                            :: i, icol, irow, j, k, nblk
    1687             : 
    1688             :       EXTERNAL                                           :: dsygv
    1689             : 
    1690          10 :       CPASSERT(ASSOCIATED(mixed_cdft))
    1691             : 
    1692          10 :       nblk = SIZE(blocks)
    1693          98 :       ALLOCATE (H_block(nblk), S_block(nblk))
    1694          34 :       DO i = 1, nblk
    1695          24 :          NULLIFY (H_block(i)%array)
    1696          24 :          NULLIFY (S_block(i)%array)
    1697          96 :          ALLOCATE (H_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
    1698          96 :          ALLOCATE (S_block(i)%array(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
    1699          24 :          icol = 0
    1700          70 :          DO j = 1, SIZE(blocks(i)%array)
    1701          46 :             irow = 0
    1702          46 :             icol = icol + 1
    1703         160 :             DO k = 1, SIZE(blocks(i)%array)
    1704          90 :                irow = irow + 1
    1705          90 :                H_block(i)%array(irow, icol) = mixed_cdft%results%H(blocks(i)%array(k), blocks(i)%array(j))
    1706         136 :                S_block(i)%array(irow, icol) = mixed_cdft%results%S(blocks(i)%array(k), blocks(i)%array(j))
    1707             :             END DO
    1708             :          END DO
    1709             :          ! Check that none of the interaction energies is repulsive
    1710         160 :          IF (ANY(H_block(i)%array .GE. 0.0_dp)) &
    1711             :             CALL cp_abort(__LOCATION__, &
    1712             :                           "At least one of the interaction energies within block "//TRIM(ADJUSTL(cp_to_string(i)))// &
    1713          10 :                           " is repulsive.")
    1714             :       END DO
    1715             : 
    1716          10 :    END SUBROUTINE mixed_cdft_get_blocks
    1717             : 
    1718             : ! **************************************************************************************************
    1719             : !> \brief Diagonalizes each of the matrix blocks.
    1720             : !> \param blocks  list of CDFT states defining the matrix blocks
    1721             : !> \param H_block list of Hamiltonian matrix blocks
    1722             : !> \param S_block list of overlap matrix blocks
    1723             : !> \param eigenvalues list of eigenvalues for each block
    1724             : !> \par History
    1725             : !>       01.18  created [Nico Holmberg]
    1726             : ! **************************************************************************************************
    1727          10 :    SUBROUTINE mixed_cdft_diagonalize_blocks(blocks, H_block, S_block, eigenvalues)
    1728             :       TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:)    :: blocks
    1729             :       TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:)    :: H_block, S_block
    1730             :       TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:), &
    1731             :          INTENT(OUT)                                     :: eigenvalues
    1732             : 
    1733             :       INTEGER                                            :: i, info, nblk, work_array_size
    1734          10 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: work
    1735          10 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: H_mat_copy, S_mat_copy
    1736             : 
    1737             :       EXTERNAL                                           :: dsygv
    1738             : 
    1739          10 :       nblk = SIZE(blocks)
    1740          54 :       ALLOCATE (eigenvalues(nblk))
    1741          34 :       DO i = 1, nblk
    1742          24 :          NULLIFY (eigenvalues(i)%array)
    1743          72 :          ALLOCATE (eigenvalues(i)%array(SIZE(blocks(i)%array)))
    1744          70 :          eigenvalues(i)%array = 0.0_dp
    1745             :          ! Workspace query
    1746          24 :          ALLOCATE (work(1))
    1747          24 :          info = 0
    1748          96 :          ALLOCATE (H_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
    1749          96 :          ALLOCATE (S_mat_copy(SIZE(blocks(i)%array), SIZE(blocks(i)%array)))
    1750         160 :          H_mat_copy(:, :) = H_block(i)%array(:, :) ! Need explicit copies because dsygv destroys original values
    1751         160 :          S_mat_copy(:, :) = S_block(i)%array(:, :)
    1752             :          CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_mat_copy, SIZE(blocks(i)%array), &
    1753          24 :                     S_mat_copy, SIZE(blocks(i)%array), eigenvalues(i)%array, work, -1, info)
    1754          24 :          work_array_size = NINT(work(1))
    1755          24 :          DEALLOCATE (H_mat_copy, S_mat_copy)
    1756             :          ! Allocate work array
    1757          24 :          DEALLOCATE (work)
    1758          72 :          ALLOCATE (work(work_array_size))
    1759        1588 :          work = 0.0_dp
    1760             :          ! Solve Hc = eSc
    1761          24 :          info = 0
    1762             :          CALL dsygv(1, 'V', 'U', SIZE(blocks(i)%array), H_block(i)%array, SIZE(blocks(i)%array), &
    1763          24 :                     S_block(i)%array, SIZE(blocks(i)%array), eigenvalues(i)%array, work, work_array_size, info)
    1764          24 :          IF (info /= 0) THEN
    1765           0 :             IF (info > SIZE(blocks(i)%array)) THEN
    1766           0 :                CPABORT("Matrix S is not positive definite")
    1767             :             ELSE
    1768           0 :                CPABORT("Diagonalization of H matrix failed.")
    1769             :             END IF
    1770             :          END IF
    1771          34 :          DEALLOCATE (work)
    1772             :       END DO
    1773             : 
    1774          10 :    END SUBROUTINE mixed_cdft_diagonalize_blocks
    1775             : 
    1776             : ! **************************************************************************************************
    1777             : !> \brief Assembles the new block diagonalized mixed CDFT Hamiltonian and overlap matrices.
    1778             : !> \param mixed_cdft the env that holds the CDFT states
    1779             : !> \param blocks  list of CDFT states defining the matrix blocks
    1780             : !> \param H_block list of Hamiltonian matrix blocks
    1781             : !> \param eigenvalues list of eigenvalues for each block
    1782             : !> \param n size of the new Hamiltonian and overlap matrices
    1783             : !> \param iounit the output unit
    1784             : !> \par History
    1785             : !>       01.18  created [Nico Holmberg]
    1786             : ! **************************************************************************************************
    1787          10 :    SUBROUTINE mixed_cdft_assemble_block_diag(mixed_cdft, blocks, H_block, eigenvalues, &
    1788             :                                              n, iounit)
    1789             :       TYPE(mixed_cdft_type), POINTER                     :: mixed_cdft
    1790             :       TYPE(cp_1d_i_p_type), ALLOCATABLE, DIMENSION(:)    :: blocks
    1791             :       TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:)    :: H_block
    1792             :       TYPE(cp_1d_r_p_type), ALLOCATABLE, DIMENSION(:)    :: eigenvalues
    1793             :       INTEGER                                            :: n, iounit
    1794             : 
    1795             :       CHARACTER(LEN=20)                                  :: ilabel, jlabel
    1796             :       CHARACTER(LEN=3)                                   :: tmp
    1797             :       INTEGER                                            :: i, icol, ipermutation, irow, j, k, l, &
    1798             :                                                             nblk, npermutations
    1799             :       LOGICAL                                            :: ignore_excited
    1800          10 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: H_mat, H_offdiag, S_mat, S_offdiag
    1801             : 
    1802             :       EXTERNAL                                           :: dsygv
    1803             : 
    1804          70 :       ALLOCATE (H_mat(n, n), S_mat(n, n))
    1805          10 :       nblk = SIZE(blocks)
    1806          10 :       ignore_excited = (nblk == n)
    1807             :       ! The diagonal contains the eigenvalues of each block
    1808          10 :       IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Eigenvalues of the block diagonalized states"
    1809         126 :       H_mat(:, :) = 0.0_dp
    1810         126 :       S_mat(:, :) = 0.0_dp
    1811          10 :       k = 1
    1812          34 :       DO i = 1, nblk
    1813          24 :          IF (iounit > 0) WRITE (iounit, '(T6,A,I3)') "Block", i
    1814          42 :          DO j = 1, SIZE(eigenvalues(i)%array)
    1815          28 :             H_mat(k, k) = eigenvalues(i)%array(j)
    1816          28 :             S_mat(k, k) = 1.0_dp
    1817          28 :             k = k + 1
    1818          28 :             IF (iounit > 0) THEN
    1819          14 :                IF (j == 1) THEN
    1820          12 :                   WRITE (iounit, '(T9,A,T58,(3X,F20.14))') 'Ground state energy:', eigenvalues(i)%array(j)
    1821             :                ELSE
    1822             :                   WRITE (iounit, '(T9,A,I2,A,T58,(3X,F20.14))') &
    1823           2 :                      'Excited state (', j - 1, ' ) energy:', eigenvalues(i)%array(j)
    1824             :                END IF
    1825             :             END IF
    1826          32 :             IF (ignore_excited .AND. j == 1) EXIT
    1827             :          END DO
    1828             :       END DO
    1829             :       ! Transform the off-diagonal blocks using the eigenvectors of each block
    1830          10 :       npermutations = nblk*(nblk - 1)/2
    1831          10 :       IF (iounit > 0) WRITE (iounit, '(/,T3,A)') "Interactions between block diagonalized states"
    1832          30 :       DO ipermutation = 1, npermutations
    1833          20 :          CALL map_permutation_to_states(nblk, ipermutation, i, j)
    1834             :          ! Get the untransformed off-diagonal block
    1835          80 :          ALLOCATE (H_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array)))
    1836          80 :          ALLOCATE (S_offdiag(SIZE(blocks(i)%array), SIZE(blocks(j)%array)))
    1837          58 :          icol = 0
    1838          58 :          DO k = 1, SIZE(blocks(j)%array)
    1839          38 :             irow = 0
    1840          38 :             icol = icol + 1
    1841         134 :             DO l = 1, SIZE(blocks(i)%array)
    1842          76 :                irow = irow + 1
    1843          76 :                H_offdiag(irow, icol) = mixed_cdft%results%H(blocks(i)%array(l), blocks(j)%array(k))
    1844         114 :                S_offdiag(irow, icol) = mixed_cdft%results%S(blocks(i)%array(l), blocks(j)%array(k))
    1845             :             END DO
    1846             :          END DO
    1847             :          ! Check that none of the interaction energies is repulsive
    1848         134 :          IF (ANY(H_offdiag .GE. 0.0_dp)) &
    1849             :             CALL cp_abort(__LOCATION__, &
    1850             :                           "At least one of the interaction energies between blocks "//TRIM(ADJUSTL(cp_to_string(i)))// &
    1851           0 :                           " and "//TRIM(ADJUSTL(cp_to_string(j)))//" is repulsive.")
    1852             :          ! Now transform: C_i^T * H * C_j
    1853         952 :          H_offdiag(:, :) = MATMUL(H_offdiag, H_block(j)%array)
    1854         972 :          H_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), H_offdiag)
    1855         952 :          S_offdiag(:, :) = MATMUL(S_offdiag, H_block(j)%array)
    1856         972 :          S_offdiag(:, :) = MATMUL(TRANSPOSE(H_block(i)%array), S_offdiag)
    1857             :          ! Make sure the transformation preserves the sign of elements in the S and H matrices
    1858             :          ! The S/H matrices contain only positive/negative values so that any sign flipping occurs in the
    1859             :          ! same elements in both matrices
    1860             :          ! Check for sign flipping using the S matrix
    1861          30 :          IF (ANY(S_offdiag .LT. 0.0_dp)) THEN
    1862          58 :             DO l = 1, SIZE(S_offdiag, 2)
    1863         134 :                DO k = 1, SIZE(S_offdiag, 1)
    1864         114 :                   IF (S_offdiag(k, l) .LT. 0.0_dp) THEN
    1865          36 :                      S_offdiag(k, l) = -1.0_dp*S_offdiag(k, l)
    1866          36 :                      H_offdiag(k, l) = -1.0_dp*H_offdiag(k, l)
    1867             :                   END IF
    1868             :                END DO
    1869             :             END DO
    1870             :          END IF
    1871          20 :          IF (ignore_excited) THEN
    1872          18 :             H_mat(i, j) = H_offdiag(1, 1)
    1873          18 :             H_mat(j, i) = H_mat(i, j)
    1874          18 :             S_mat(i, j) = S_offdiag(1, 1)
    1875          18 :             S_mat(j, i) = S_mat(i, j)
    1876             :          ELSE
    1877           2 :             irow = 1
    1878           2 :             icol = 1
    1879           2 :             DO k = 1, i - 1
    1880           2 :                irow = irow + SIZE(blocks(k)%array)
    1881             :             END DO
    1882           4 :             DO k = 1, j - 1
    1883           4 :                icol = icol + SIZE(blocks(k)%array)
    1884             :             END DO
    1885          14 :             H_mat(irow:irow + SIZE(H_offdiag, 1) - 1, icol:icol + SIZE(H_offdiag, 2) - 1) = H_offdiag(:, :)
    1886          14 :             H_mat(icol:icol + SIZE(H_offdiag, 2) - 1, irow:irow + SIZE(H_offdiag, 1) - 1) = TRANSPOSE(H_offdiag)
    1887          14 :             S_mat(irow:irow + SIZE(H_offdiag, 1) - 1, icol:icol + SIZE(H_offdiag, 2) - 1) = S_offdiag(:, :)
    1888          14 :             S_mat(icol:icol + SIZE(H_offdiag, 2) - 1, irow:irow + SIZE(H_offdiag, 1) - 1) = TRANSPOSE(S_offdiag)
    1889             :          END IF
    1890          20 :          IF (iounit > 0) THEN
    1891          10 :             WRITE (iounit, '(/,T3,A)') REPEAT('#', 39)
    1892          10 :             WRITE (iounit, '(T3,A,I3,A,I3,A)') '###### Blocks I =', i, ' and J = ', j, ' ######'
    1893          10 :             WRITE (iounit, '(T3,A)') REPEAT('#', 39)
    1894          10 :             WRITE (iounit, '(T3,A)') 'Interaction energies'
    1895          21 :             DO irow = 1, SIZE(H_offdiag, 1)
    1896          20 :                ilabel = "(ground state)"
    1897          20 :                IF (irow > 1) THEN
    1898          10 :                   IF (ignore_excited) EXIT
    1899           1 :                   WRITE (tmp, '(I3)') irow - 1
    1900           1 :                   ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
    1901             :                END IF
    1902          34 :                DO icol = 1, SIZE(H_offdiag, 2)
    1903          21 :                   jlabel = "(ground state)"
    1904          21 :                   IF (icol > 1) THEN
    1905          10 :                      IF (ignore_excited) EXIT
    1906           2 :                      WRITE (tmp, '(I3)') icol - 1
    1907           2 :                      jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
    1908             :                   END IF
    1909          24 :                   WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', H_offdiag(irow, icol)
    1910             :                END DO
    1911             :             END DO
    1912          10 :             WRITE (iounit, '(T3,A)') 'Overlaps'
    1913          21 :             DO irow = 1, SIZE(H_offdiag, 1)
    1914          20 :                ilabel = "(ground state)"
    1915          20 :                IF (irow > 1) THEN
    1916          10 :                   IF (ignore_excited) EXIT
    1917           1 :                   ilabel = "(excited state)"
    1918           1 :                   WRITE (tmp, '(I3)') irow - 1
    1919           1 :                   ilabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
    1920             :                END IF
    1921          34 :                DO icol = 1, SIZE(H_offdiag, 2)
    1922          21 :                   jlabel = "(ground state)"
    1923          21 :                   IF (icol > 1) THEN
    1924          10 :                      IF (ignore_excited) EXIT
    1925           2 :                      WRITE (tmp, '(I3)') icol - 1
    1926           2 :                      jlabel = "(excited state "//TRIM(ADJUSTL(tmp))//")"
    1927             :                   END IF
    1928          24 :                   WRITE (iounit, '(T6,A,T58,(3X,F20.14))') TRIM(ilabel)//'-'//TRIM(jlabel)//':', S_offdiag(irow, icol)
    1929             :                END DO
    1930             :             END DO
    1931             :          END IF
    1932          50 :          DEALLOCATE (H_offdiag, S_offdiag)
    1933             :       END DO
    1934          10 :       CALL mixed_cdft_result_type_set(mixed_cdft%results, H=H_mat, S=S_mat)
    1935             :       ! Deallocate work
    1936          10 :       DEALLOCATE (H_mat, S_mat)
    1937             : 
    1938          10 :    END SUBROUTINE mixed_cdft_assemble_block_diag
    1939             : 
    1940             : END MODULE mixed_cdft_utils

Generated by: LCOV version 1.15