LCOV - code coverage report
Current view: top level - src - qs_active_space_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 72.7 % 1492 1085
Test Date: 2025-07-25 12:55:17 Functions: 61.5 % 26 16

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Determine active space Hamiltonian
      10              : !> \par History
      11              : !>      04.2016 created [JGH]
      12              : !> \author JGH
      13              : ! **************************************************************************************************
      14              : MODULE qs_active_space_methods
      15              :    USE atomic_kind_types, ONLY: atomic_kind_type
      16              :    USE basis_set_types, ONLY: allocate_sto_basis_set, &
      17              :                               create_gto_from_sto_basis, &
      18              :                               deallocate_sto_basis_set, &
      19              :                               gto_basis_set_type, &
      20              :                               init_orb_basis_set, &
      21              :                               set_sto_basis_set, &
      22              :                               srules, &
      23              :                               sto_basis_set_type
      24              :    USE cell_types, ONLY: cell_type
      25              :    USE cp_blacs_env, ONLY: cp_blacs_env_type, cp_blacs_env_create, cp_blacs_env_release, BLACS_GRID_SQUARE
      26              :    USE cp_control_types, ONLY: dft_control_type, qs_control_type
      27              :    USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t, &
      28              :                                   cp_dbcsr_sm_fm_multiply, &
      29              :                                   dbcsr_allocate_matrix_set, &
      30              :                                   cp_dbcsr_m_by_n_from_template, copy_dbcsr_to_fm
      31              :    USE cp_files, ONLY: close_file, &
      32              :                        file_exists, &
      33              :                        open_file
      34              :    USE cp_fm_struct, ONLY: cp_fm_struct_create, &
      35              :                            cp_fm_struct_release, &
      36              :                            cp_fm_struct_type
      37              :    USE cp_fm_types, ONLY: &
      38              :       cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_init_random, cp_fm_release, &
      39              :       cp_fm_set_all, cp_fm_set_element, cp_fm_to_fm, cp_fm_type
      40              :    USE cp_log_handling, ONLY: cp_get_default_logger, &
      41              :                               cp_logger_get_default_io_unit, &
      42              :                               cp_logger_type
      43              :    USE cp_output_handling, ONLY: &
      44              :       cp_p_file, cp_print_key_finished_output, cp_print_key_should_output, cp_print_key_unit_nr, &
      45              :       debug_print_level, high_print_level, low_print_level, medium_print_level, &
      46              :       silent_print_level
      47              :    USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
      48              :    USE cp_dbcsr_api, ONLY: &
      49              :       dbcsr_copy, dbcsr_csr_create, dbcsr_csr_type, dbcsr_p_type, dbcsr_type, dbcsr_release, &
      50              :       dbcsr_type_no_symmetry, dbcsr_create, dbcsr_set, dbcsr_multiply, dbcsr_iterator_next_block, &
      51              :       dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_blocks_left, &
      52              :       dbcsr_iterator_type, dbcsr_type_symmetric, dbcsr_get_occupation, dbcsr_get_info
      53              :    USE group_dist_types, ONLY: get_group_dist, release_group_dist, group_dist_d1_type
      54              :    USE input_constants, ONLY: &
      55              :       casci_canonical, dmft_model, eri_method_full_gpw, eri_method_gpw_ht, eri_operator_coulomb, &
      56              :       eri_operator_trunc, eri_operator_erf, eri_operator_erfc, eri_operator_gaussian, eri_operator_yukawa, hf_model, &
      57              :       manual_selection, mao_projection, no_solver, qiskit_solver, rsdft_model, wannier_projection
      58              :    USE input_section_types, ONLY: section_vals_get, &
      59              :                                   section_vals_get_subs_vals, &
      60              :                                   section_vals_type, &
      61              :                                   section_vals_val_get
      62              :    USE ISO_C_BINDING, ONLY: c_null_char
      63              :    USE kinds, ONLY: default_path_length, &
      64              :                     default_string_length, &
      65              :                     dp, &
      66              :                     int_8
      67              :    USE machine, ONLY: m_walltime
      68              :    USE mathconstants, ONLY: fourpi
      69              :    USE memory_utilities, ONLY: reallocate
      70              :    USE message_passing, ONLY: mp_comm_type, &
      71              :                               mp_para_env_type, &
      72              :                               mp_para_env_release
      73              :    USE mp2_gpw, ONLY: create_mat_munu, grep_rows_in_subgroups, build_dbcsr_from_rows
      74              :    USE parallel_gemm_api, ONLY: parallel_gemm
      75              :    USE particle_list_types, ONLY: particle_list_type
      76              :    USE particle_types, ONLY: particle_type
      77              :    USE periodic_table, ONLY: ptable
      78              :    USE preconditioner_types, ONLY: preconditioner_type
      79              :    USE pw_env_methods, ONLY: pw_env_create, &
      80              :                              pw_env_rebuild
      81              :    USE pw_env_types, ONLY: pw_env_get, &
      82              :                            pw_env_release, &
      83              :                            pw_env_type
      84              :    USE pw_methods, ONLY: pw_integrate_function, &
      85              :                          pw_multiply, &
      86              :                          pw_multiply_with, &
      87              :                          pw_transfer, &
      88              :                          pw_zero, pw_integral_ab, pw_scale
      89              :    USE pw_poisson_methods, ONLY: pw_poisson_rebuild, &
      90              :                                  pw_poisson_solve
      91              :    USE pw_poisson_types, ONLY: ANALYTIC0D, &
      92              :                                PERIODIC3D, &
      93              :                                greens_fn_type, &
      94              :                                pw_poisson_analytic, &
      95              :                                pw_poisson_periodic, &
      96              :                                pw_poisson_type
      97              :    USE pw_pool_types, ONLY: &
      98              :       pw_pool_type
      99              :    USE pw_types, ONLY: &
     100              :       pw_c1d_gs_type, &
     101              :       pw_r3d_rs_type
     102              :    USE qcschema, ONLY: qcschema_env_create, &
     103              :                        qcschema_env_release, &
     104              :                        qcschema_to_hdf5, &
     105              :                        qcschema_type
     106              :    USE qs_active_space_types, ONLY: active_space_type, &
     107              :                                     create_active_space_type, &
     108              :                                     csr_idx_from_combined, &
     109              :                                     csr_idx_to_combined, &
     110              :                                     eri_type, &
     111              :                                     eri_type_eri_element_func, &
     112              :                                     get_irange_csr
     113              :    USE qs_active_space_utils, ONLY: eri_to_array, &
     114              :                                     subspace_matrix_to_array
     115              :    USE qs_collocate_density, ONLY: calculate_wavefunction
     116              :    USE qs_density_matrices, ONLY: calculate_density_matrix
     117              :    USE qs_energy_types, ONLY: qs_energy_type
     118              :    USE qs_environment_types, ONLY: get_qs_env, &
     119              :                                    qs_environment_type, &
     120              :                                    set_qs_env
     121              :    USE qs_integrate_potential, ONLY: integrate_v_rspace
     122              :    USE qs_kind_types, ONLY: qs_kind_type
     123              :    USE qs_ks_methods, ONLY: qs_ks_update_qs_env
     124              :    USE qs_ks_types, ONLY: qs_ks_did_change, &
     125              :                           qs_ks_env_type
     126              :    USE qs_mo_io, ONLY: write_mo_set_to_output_unit
     127              :    USE qs_mo_methods, ONLY: calculate_subspace_eigenvalues
     128              :    USE qs_mo_types, ONLY: allocate_mo_set, &
     129              :                           get_mo_set, &
     130              :                           init_mo_set, &
     131              :                           mo_set_type
     132              :    USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type, release_neighbor_list_sets
     133              :    USE qs_ot_eigensolver, ONLY: ot_eigensolver
     134              :    USE qs_rho_methods, ONLY: qs_rho_update_rho
     135              :    USE qs_rho_types, ONLY: qs_rho_get, &
     136              :                            qs_rho_type
     137              :    USE qs_subsys_types, ONLY: qs_subsys_get, &
     138              :                               qs_subsys_type
     139              :    USE scf_control_types, ONLY: scf_control_type
     140              : #ifndef __NO_SOCKETS
     141              :    USE sockets_interface, ONLY: accept_socket, &
     142              :                                 close_socket, &
     143              :                                 listen_socket, &
     144              :                                 open_bind_socket, &
     145              :                                 readbuffer, &
     146              :                                 remove_socket_file, &
     147              :                                 writebuffer
     148              : #endif
     149              :    USE task_list_methods, ONLY: generate_qs_task_list
     150              :    USE task_list_types, ONLY: allocate_task_list, &
     151              :                               deallocate_task_list, &
     152              :                               task_list_type
     153              :    USE util, ONLY: get_limit
     154              : #include "./base/base_uses.f90"
     155              : 
     156              :    IMPLICIT NONE
     157              :    PRIVATE
     158              : 
     159              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_methods'
     160              : 
     161              :    PUBLIC :: active_space_main
     162              : 
     163              :    TYPE, EXTENDS(eri_type_eri_element_func) :: eri_fcidump_print
     164              :       INTEGER :: unit_nr = -1, bra_start = -1, ket_start = -1
     165              :    CONTAINS
     166              :       PROCEDURE :: func => eri_fcidump_print_func
     167              :    END TYPE
     168              : 
     169              :    TYPE, EXTENDS(eri_type_eri_element_func) :: eri_fcidump_checksum
     170              :       INTEGER :: bra_start = 0, ket_start = 0
     171              :       REAL(KIND=dp) :: checksum = 0.0_dp
     172              :    CONTAINS
     173              :       PROCEDURE, PASS :: set => eri_fcidump_set
     174              :       PROCEDURE :: func => eri_fcidump_checksum_func
     175              :    END TYPE eri_fcidump_checksum
     176              : 
     177              : CONTAINS
     178              : 
     179              : ! **************************************************************************************************
     180              : !> \brief Sets the starting indices of the bra and ket.
     181              : !> \param this object reference
     182              : !> \param bra_start starting index of the bra
     183              : !> \param ket_start starting index of the ket
     184              : ! **************************************************************************************************
     185           86 :    SUBROUTINE eri_fcidump_set(this, bra_start, ket_start)
     186              :       CLASS(eri_fcidump_checksum) :: this
     187              :       INTEGER, INTENT(IN) :: bra_start, ket_start
     188           86 :       this%bra_start = bra_start
     189           86 :       this%ket_start = ket_start
     190           86 :    END SUBROUTINE eri_fcidump_set
     191              : 
     192              : ! **************************************************************************************************
     193              : !> \brief Main method for determining the active space Hamiltonian
     194              : !> \param qs_env ...
     195              : ! **************************************************************************************************
     196        21701 :    SUBROUTINE active_space_main(qs_env)
     197              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     198              : 
     199              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'active_space_main'
     200              : 
     201              :       CHARACTER(len=10)                                  :: cshell, lnam(5)
     202              :       CHARACTER(len=default_path_length)                 :: qcschema_filename
     203              :       INTEGER :: as_solver, eri_method, eri_operator, eri_print, group_size, handle, i, iatom, &
     204              :          ishell, isp, ispin, iw, j, jm, m, max_orb_ind, mselect, n1, n2, nao, natom, nel, &
     205              :          nelec_active, nelec_inactive, nelec_total, nmo, nmo_active, nmo_available, nmo_inactive, &
     206              :          nmo_inactive_remaining, nmo_occ, nmo_virtual, nn1, nn2, nrow_global, nspins
     207              :       INTEGER, DIMENSION(5)                              :: nshell
     208        21701 :       INTEGER, DIMENSION(:), POINTER                     :: invals
     209              :       LOGICAL                                            :: do_kpoints, ex_operator, ex_perd, &
     210              :                                                             explicit, isolated, stop_after_print, &
     211              :                                                             store_wfn
     212              :       REAL(KIND=dp) :: eri_eps_filter, eri_eps_grid, eri_eps_int, eri_gpw_cutoff, eri_op_param, &
     213              :          eri_rcut, eri_rel_cutoff, fel, focc, maxocc, nze_percentage
     214        21701 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: eigenvalues
     215        21701 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: evals_virtual
     216              :       TYPE(active_space_type), POINTER                   :: active_space_env
     217              :       TYPE(cell_type), POINTER                           :: cell
     218              :       TYPE(cp_blacs_env_type), POINTER                   :: context
     219              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_tmp
     220              :       TYPE(cp_fm_type)                                   :: fm_dummy, mo_virtual
     221              :       TYPE(cp_fm_type), POINTER                          :: fm_target_active, fm_target_inactive, &
     222              :                                                             fmat, mo_coeff, mo_ref, mo_target
     223              :       TYPE(cp_logger_type), POINTER                      :: logger
     224              :       TYPE(dbcsr_csr_type), POINTER                      :: eri_mat
     225        43402 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: ks_matrix, rho_ao, s_matrix
     226              :       TYPE(dbcsr_type), POINTER                          :: denmat
     227              :       TYPE(dft_control_type), POINTER                    :: dft_control
     228        21701 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     229              :       TYPE(mo_set_type), POINTER                         :: mo_set, mo_set_active, mo_set_inactive
     230              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     231              :       TYPE(preconditioner_type), POINTER                 :: local_preconditioner
     232        86804 :       TYPE(qcschema_type)                                :: qcschema_env
     233              :       TYPE(qs_energy_type), POINTER                      :: energy
     234              :       TYPE(qs_rho_type), POINTER                         :: rho
     235              :       TYPE(scf_control_type), POINTER                    :: scf_control
     236              :       TYPE(section_vals_type), POINTER                   :: as_input, input, loc_print, loc_section, &
     237              :                                                             print_orb
     238              : 
     239              :       !--------------------------------------------------------------------------------------------!
     240              : 
     241        21701 :       CALL get_qs_env(qs_env, input=input)
     242        21701 :       as_input => section_vals_get_subs_vals(input, "DFT%ACTIVE_SPACE")
     243        21701 :       CALL section_vals_get(as_input, explicit=explicit)
     244        21701 :       IF (.NOT. explicit) RETURN
     245           62 :       CALL timeset(routineN, handle)
     246              : 
     247           62 :       logger => cp_get_default_logger()
     248           62 :       iw = cp_logger_get_default_io_unit(logger)
     249              : 
     250           62 :       IF (iw > 0) THEN
     251              :          WRITE (iw, '(/,T2,A)') &
     252           31 :             '!-----------------------------------------------------------------------------!'
     253           31 :          WRITE (iw, '(T26,A)') "Generate Active Space Hamiltonian"
     254           31 :          WRITE (iw, '(T27,A)') "Interface to CAS-CI and DMRG-CI"
     255              :          WRITE (iw, '(T2,A)') &
     256           31 :             '!-----------------------------------------------------------------------------!'
     257              :       END IF
     258              : 
     259              :       ! k-points?
     260           62 :       CALL get_qs_env(qs_env, do_kpoints=do_kpoints)
     261           62 :       IF (do_kpoints) THEN
     262            0 :          CALL cp_abort(__LOCATION__, "K-points not allowd for Active Space Interface")
     263              :       END IF
     264              : 
     265           62 :       NULLIFY (active_space_env)
     266           62 :       CALL create_active_space_type(active_space_env)
     267           62 :       active_space_env%energy_total = 0.0_dp
     268           62 :       active_space_env%energy_ref = 0.0_dp
     269           62 :       active_space_env%energy_inactive = 0.0_dp
     270           62 :       active_space_env%energy_active = 0.0_dp
     271              : 
     272              :       ! input options
     273              : 
     274              :       ! figure out what needs to be printed/stored
     275           62 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, as_input, "FCIDUMP"), cp_p_file)) THEN
     276           62 :          active_space_env%fcidump = .TRUE.
     277              :       END IF
     278              : 
     279           62 :       CALL section_vals_val_get(as_input, "QCSCHEMA", c_val=qcschema_filename, explicit=explicit)
     280           62 :       IF (explicit) THEN
     281            0 :          active_space_env%qcschema = .TRUE.
     282            0 :          active_space_env%qcschema_filename = qcschema_filename
     283              :       END IF
     284              : 
     285              :       ! model
     286           62 :       CALL section_vals_val_get(as_input, "MODEL", i_val=active_space_env%model)
     287           62 :       IF (iw > 0) THEN
     288           62 :          SELECT CASE (active_space_env%model)
     289              :          CASE (hf_model)
     290           31 :             WRITE (iw, '(T3,A)') "Hartree-Fock model for interaction Hamiltonian"
     291              :          CASE (rsdft_model)
     292            0 :             WRITE (iw, '(T3,A)') "Range-separated DFT model for interaction Hamiltonian"
     293              :          CASE (dmft_model)
     294            0 :             WRITE (iw, '(T3,A)') "DMFT model for interaction Hamiltonian"
     295              :          CASE DEFAULT
     296           31 :             CPABORT("Unknown Model")
     297              :          END SELECT
     298              :       END IF
     299              : 
     300              :       ! isolated (molecular) system?
     301           62 :       CALL section_vals_val_get(as_input, "ISOLATED_SYSTEM", l_val=isolated)
     302           62 :       active_space_env%molecule = isolated
     303           62 :       IF (iw > 0) THEN
     304           31 :          IF (active_space_env%molecule) THEN
     305           29 :             WRITE (iw, '(T3,A)') "System is treated without periodicity"
     306              :          END IF
     307              :       END IF
     308              : 
     309           62 :       CALL section_vals_val_get(as_input, "ACTIVE_ELECTRONS", i_val=nelec_active)
     310              :       ! actually nelec_spin tells me the number of electrons per spin channel from qs_env
     311              :       ! CALL get_qs_env(qs_env, nelectron_total=nelec_total, nelectron_spin=nelec_spin)
     312           62 :       CALL get_qs_env(qs_env, nelectron_total=nelec_total)
     313              : 
     314           62 :       IF (nelec_active <= 0) CPABORT("Specify a positive number of active electrons.")
     315           62 :       IF (nelec_active > nelec_total) CPABORT("More active electrons than total electrons.")
     316              : 
     317           62 :       nelec_inactive = nelec_total - nelec_active
     318           62 :       IF (MOD(nelec_inactive, 2) /= 0) THEN
     319            0 :          CPABORT("The remaining number of inactive electrons has to be even.")
     320              :       END IF
     321              : 
     322           62 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     323           62 :       nspins = dft_control%nspins
     324           62 :       IF (iw > 0) THEN
     325           31 :          WRITE (iw, '(T3,A,T70,I10)') "Total number of electrons", nelec_total
     326           31 :          WRITE (iw, '(T3,A,T70,I10)') "Number of inactive electrons", nelec_inactive
     327           31 :          WRITE (iw, '(T3,A,T70,I10)') "Number of active electrons", nelec_active
     328              :       END IF
     329              : 
     330           62 :       active_space_env%nelec_active = nelec_active
     331           62 :       active_space_env%nelec_inactive = nelec_inactive
     332           62 :       active_space_env%nelec_total = nelec_total
     333           62 :       active_space_env%nspins = nspins
     334           62 :       active_space_env%multiplicity = dft_control%multiplicity
     335              : 
     336              :       ! define the active/inactive space orbitals
     337           62 :       CALL section_vals_val_get(as_input, "ACTIVE_ORBITALS", explicit=explicit, i_val=nmo_active)
     338           62 :       IF (.NOT. explicit) THEN
     339            0 :          CALL cp_abort(__LOCATION__, "Number of Active Orbitals has to be specified.")
     340              :       END IF
     341           62 :       active_space_env%nmo_active = nmo_active
     342              :       ! this is safe because nelec_inactive is always even
     343           62 :       nmo_inactive = nelec_inactive/2
     344           62 :       active_space_env%nmo_inactive = nmo_inactive
     345              : 
     346           62 :       CALL section_vals_val_get(as_input, "ORBITAL_SELECTION", i_val=mselect)
     347           62 :       IF (iw > 0) THEN
     348            0 :          SELECT CASE (mselect)
     349              :          CASE DEFAULT
     350            0 :             CPABORT("Unknown orbital selection method")
     351              :          CASE (casci_canonical)
     352              :             WRITE (iw, '(/,T3,A)') &
     353           25 :                "Active space orbitals selected using energy ordered canonical orbitals"
     354              :          CASE (wannier_projection)
     355              :             WRITE (iw, '(/,T3,A)') &
     356            0 :                "Active space orbitals selected using projected Wannier orbitals"
     357              :          CASE (mao_projection)
     358              :             WRITE (iw, '(/,T3,A)') &
     359            0 :                "Active space orbitals selected using modified atomic orbitals (MAO)"
     360              :          CASE (manual_selection)
     361              :             WRITE (iw, '(/,T3,A)') &
     362           31 :                "Active space orbitals selected manually"
     363              :          END SELECT
     364              : 
     365           31 :          WRITE (iw, '(T3,A,T70,I10)') "Number of inactive orbitals", nmo_inactive
     366           31 :          WRITE (iw, '(T3,A,T70,I10)') "Number of active orbitals", nmo_active
     367              :       END IF
     368              : 
     369              :       ! get projection spaces
     370           62 :       CALL section_vals_val_get(as_input, "SUBSPACE_ATOM", i_val=iatom, explicit=explicit)
     371           62 :       IF (explicit) THEN
     372            0 :          CALL get_qs_env(qs_env, natom=natom)
     373            0 :          IF (iatom <= 0 .OR. iatom > natom) THEN
     374            0 :             IF (iw > 0) THEN
     375            0 :                WRITE (iw, '(/,T3,A,I3)') "ERROR: SUBSPACE_ATOM number is not valid", iatom
     376              :             END IF
     377            0 :             CPABORT("Select a valid SUBSPACE_ATOM")
     378              :          END IF
     379              :       END IF
     380           62 :       CALL section_vals_val_get(as_input, "SUBSPACE_SHELL", c_val=cshell, explicit=explicit)
     381           62 :       nshell = 0
     382          372 :       lnam = ""
     383           62 :       IF (explicit) THEN
     384            0 :          cshell = ADJUSTL(cshell)
     385            0 :          n1 = 1
     386            0 :          DO i = 1, 5
     387            0 :             ishell = i
     388            0 :             IF (cshell(n1:n1) == " ") THEN
     389           62 :                ishell = ishell - 1
     390              :                EXIT
     391              :             END IF
     392            0 :             READ (cshell(n1:), "(I1,A1)") nshell(i), lnam(i)
     393            0 :             n1 = n1 + 2
     394              :          END DO
     395              :       END IF
     396              : 
     397              :       ! generate orbitals
     398            0 :       SELECT CASE (mselect)
     399              :       CASE DEFAULT
     400            0 :          CPABORT("Unknown orbital selection method")
     401              :       CASE (casci_canonical)
     402           50 :          CALL get_qs_env(qs_env, mos=mos)
     403              : 
     404              :          ! total number of occupied orbitals, i.e. inactive plus active MOs
     405           50 :          nmo_occ = nmo_inactive + nmo_active
     406              : 
     407              :          ! set inactive orbital indices, these are trivially 1...nmo_inactive
     408          160 :          ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
     409          106 :          DO ispin = 1, nspins
     410          146 :             DO i = 1, nmo_inactive
     411           96 :                active_space_env%inactive_orbitals(i, ispin) = i
     412              :             END DO
     413              :          END DO
     414              : 
     415              :          ! set active orbital indices, these are shifted by nmo_inactive
     416          200 :          ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
     417          106 :          DO ispin = 1, nspins
     418          258 :             DO i = 1, nmo_active
     419          208 :                active_space_env%active_orbitals(i, ispin) = nmo_inactive + i
     420              :             END DO
     421              :          END DO
     422              : 
     423              :          ! allocate and initialize inactive and active mo coefficients.
     424              :          ! These are stored in a data structure for the full occupied space:
     425              :          ! for inactive mos, the active subset is set to zero, vice versa for the active mos
     426              :          ! TODO: allocate data structures only for the eaxct number MOs
     427           50 :          maxocc = 2.0_dp
     428           50 :          IF (nspins > 1) maxocc = 1.0_dp
     429          206 :          ALLOCATE (active_space_env%mos_active(nspins))
     430          156 :          ALLOCATE (active_space_env%mos_inactive(nspins))
     431          106 :          DO ispin = 1, nspins
     432           56 :             CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
     433           56 :             CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
     434              :             ! the right number of active electrons per spin channel is initialized further down
     435           56 :             CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, nmo_occ, 0, 0.0_dp, maxocc, 0.0_dp)
     436              :             CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     437           56 :                                      nrow_global=nrow_global, ncol_global=nmo_occ)
     438           56 :             CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO")
     439           56 :             CALL cp_fm_struct_release(fm_struct_tmp)
     440           56 :             IF (nspins == 2) THEN
     441           12 :                nel = nelec_inactive/2
     442              :             ELSE
     443           44 :                nel = nelec_inactive
     444              :             END IF
     445              :             CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, nmo_occ, nel, &
     446           56 :                                  REAL(nel, KIND=dp), maxocc, 0.0_dp)
     447              :             CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     448           56 :                                      nrow_global=nrow_global, ncol_global=nmo_occ)
     449           56 :             CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO")
     450          162 :             CALL cp_fm_struct_release(fm_struct_tmp)
     451              :          END DO
     452              : 
     453              :          ! create canonical orbitals
     454           50 :          IF (dft_control%restricted) THEN
     455            0 :             CPABORT("Unclear how we define MOs in the restricted case ... stopping")
     456              :          ELSE
     457           50 :             IF (dft_control%do_admm) THEN
     458            0 :                CPABORT("ADMM currently not possible for canonical orbital options")
     459              :             END IF
     460              : 
     461          200 :             ALLOCATE (eigenvalues(nmo_occ, nspins))
     462          298 :             eigenvalues = 0.0_dp
     463           50 :             CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control)
     464              : 
     465              :             ! calculate virtual MOs and copy inactive and active orbitals
     466           50 :             IF (iw > 0) THEN
     467           25 :                WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..."
     468              :             END IF
     469          106 :             DO ispin = 1, nspins
     470              :                ! nmo_available is the number of MOs available from the SCF calculation:
     471              :                ! this is at least the number of occupied orbitals in the SCF, plus
     472              :                ! any number of added MOs (virtuals) requested in the SCF section
     473           56 :                CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
     474              : 
     475              :                ! calculate how many extra MOs we still have to compute
     476           56 :                nmo_virtual = nmo_occ - nmo_available
     477           56 :                nmo_virtual = MAX(nmo_virtual, 0)
     478              : 
     479              :                NULLIFY (evals_virtual)
     480          112 :                ALLOCATE (evals_virtual(nmo_virtual))
     481              : 
     482              :                CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, &
     483           56 :                                    nrow_global=nrow_global)
     484              : 
     485              :                CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     486           56 :                                         nrow_global=nrow_global, ncol_global=nmo_virtual)
     487           56 :                CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual")
     488           56 :                CALL cp_fm_struct_release(fm_struct_tmp)
     489           56 :                CALL cp_fm_init_random(mo_virtual, nmo_virtual)
     490              : 
     491           56 :                NULLIFY (local_preconditioner)
     492              : 
     493              :                ! compute missing virtual MOs
     494              :                CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
     495              :                                    matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, &
     496              :                                    eps_gradient=scf_control%eps_lumos, &
     497              :                                    preconditioner=local_preconditioner, &
     498              :                                    iter_max=scf_control%max_iter_lumos, &
     499           56 :                                    size_ortho_space=nmo_available)
     500              : 
     501              :                ! get the eigenvalues
     502           56 :                CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, evals_virtual)
     503              : 
     504              :                ! TODO: double check that we really need this.
     505              :                ! we need to send the copy of MOs to preserve the sign
     506           56 :                CALL cp_fm_create(fm_dummy, mo_ref%matrix_struct)
     507           56 :                CALL cp_fm_to_fm(mo_ref, fm_dummy)
     508              :                CALL calculate_subspace_eigenvalues(fm_dummy, ks_matrix(ispin)%matrix, &
     509           56 :                                                    evals_arg=eigenvalues(:, ispin), do_rotation=.TRUE.)
     510              : 
     511              :                ! copy inactive orbitals
     512           56 :                mo_set => active_space_env%mos_inactive(ispin)
     513           56 :                CALL get_mo_set(mo_set, mo_coeff=mo_target)
     514           96 :                DO i = 1, SIZE(active_space_env%inactive_orbitals, 1)
     515           40 :                   m = active_space_env%inactive_orbitals(i, ispin)
     516           40 :                   CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m)
     517           40 :                   mo_set%eigenvalues(m) = eigenvalues(m, ispin)
     518           96 :                   IF (nspins > 1) THEN
     519           28 :                      mo_set%occupation_numbers(m) = 1.0
     520              :                   ELSE
     521           12 :                      mo_set%occupation_numbers(m) = 2.0
     522              :                   END IF
     523              :                END DO
     524              : 
     525              :                ! copy active orbitals
     526           56 :                mo_set => active_space_env%mos_active(ispin)
     527           56 :                CALL get_mo_set(mo_set, mo_coeff=mo_target)
     528              :                ! for mult > 1, put the polarized electrons in the alpha channel
     529           56 :                IF (nspins == 2) THEN
     530           12 :                   IF (ispin == 1) THEN
     531            6 :                      nel = (nelec_active + active_space_env%multiplicity - 1)/2
     532              :                   ELSE
     533            6 :                      nel = (nelec_active - active_space_env%multiplicity + 1)/2
     534              :                   END IF
     535              :                ELSE
     536           44 :                   nel = nelec_active
     537              :                END IF
     538           56 :                mo_set%nelectron = nel
     539           56 :                mo_set%n_el_f = REAL(nel, KIND=dp)
     540          208 :                DO i = 1, nmo_active
     541          152 :                   m = active_space_env%active_orbitals(i, ispin)
     542          152 :                   IF (m > nmo_available) THEN
     543            0 :                      CALL cp_fm_to_fm(mo_virtual, mo_target, 1, m - nmo_available, m)
     544            0 :                      eigenvalues(m, ispin) = evals_virtual(m - nmo_available)
     545            0 :                      mo_set%occupation_numbers(m) = 0.0
     546              :                   ELSE
     547          152 :                      CALL cp_fm_to_fm(mo_ref, mo_target, 1, m, m)
     548          152 :                      mo_set%occupation_numbers(m) = mos(ispin)%occupation_numbers(m)
     549              :                   END IF
     550          208 :                   mo_set%eigenvalues(m) = eigenvalues(m, ispin)
     551              :                END DO
     552              :                ! Release
     553           56 :                DEALLOCATE (evals_virtual)
     554           56 :                CALL cp_fm_release(fm_dummy)
     555          330 :                CALL cp_fm_release(mo_virtual)
     556              :             END DO
     557              : 
     558           50 :             IF (iw > 0) THEN
     559           53 :                DO ispin = 1, nspins
     560           28 :                   WRITE (iw, '(/,T3,A,I3,T66,A)') "Canonical Orbital Selection for spin", ispin, &
     561           56 :                      "[atomic units]"
     562           36 :                   DO i = 1, nmo_inactive, 4
     563            8 :                      jm = MIN(3, nmo_inactive - i)
     564           56 :                      WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [I]", j=0, jm)
     565              :                   END DO
     566           57 :                   DO i = nmo_inactive + 1, nmo_inactive + nmo_active, 4
     567           29 :                      jm = MIN(3, nmo_inactive + nmo_active - i)
     568          133 :                      WRITE (iw, '(T3,4(F14.6,A5))') (eigenvalues(i + j, ispin), " [A]", j=0, jm)
     569              :                   END DO
     570           28 :                   WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin
     571           82 :                   DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4
     572           29 :                      jm = MIN(3, SIZE(active_space_env%active_orbitals, 1) - i)
     573          133 :                      WRITE (iw, '(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm)
     574              :                   END DO
     575              :                END DO
     576              :             END IF
     577           50 :             DEALLOCATE (eigenvalues)
     578              :          END IF
     579              : 
     580              :       CASE (manual_selection)
     581              :          ! create canonical orbitals
     582           12 :          IF (dft_control%restricted) THEN
     583            0 :             CPABORT("Unclear how we define MOs in the restricted case ... stopping")
     584              :          ELSE
     585           12 :             IF (dft_control%do_admm) THEN
     586            0 :                CPABORT("ADMM currently not possible for manual orbitals selection")
     587              :             END IF
     588              : 
     589           12 :             CALL section_vals_val_get(as_input, "ACTIVE_ORBITAL_INDICES", explicit=explicit, i_vals=invals)
     590           12 :             IF (.NOT. explicit) THEN
     591              :                CALL cp_abort(__LOCATION__, "Manual orbital selection requires to explicitly "// &
     592            0 :                              "set the active orbital indices via ACTIVE_ORBITAL_INDICES")
     593              :             END IF
     594              : 
     595           12 :             IF (nspins == 1) THEN
     596            6 :                CPASSERT(SIZE(invals) == nmo_active)
     597              :             ELSE
     598            6 :                CPASSERT(SIZE(invals) == 2*nmo_active)
     599              :             END IF
     600           36 :             ALLOCATE (active_space_env%inactive_orbitals(nmo_inactive, nspins))
     601           48 :             ALLOCATE (active_space_env%active_orbitals(nmo_active, nspins))
     602              : 
     603           30 :             DO ispin = 1, nspins
     604           66 :                DO i = 1, nmo_active
     605           54 :                   active_space_env%active_orbitals(i, ispin) = invals(i + (ispin - 1)*nmo_active)
     606              :                END DO
     607              :             END DO
     608              : 
     609           12 :             CALL get_qs_env(qs_env, mos=mos)
     610              : 
     611              :             ! include MOs up to the largest index in the list
     612           48 :             max_orb_ind = MAXVAL(invals)
     613           12 :             maxocc = 2.0_dp
     614           12 :             IF (nspins > 1) maxocc = 1.0_dp
     615           54 :             ALLOCATE (active_space_env%mos_active(nspins))
     616           42 :             ALLOCATE (active_space_env%mos_inactive(nspins))
     617           30 :             DO ispin = 1, nspins
     618              :                ! init active orbitals
     619           18 :                CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nao=nao)
     620           18 :                CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, nrow_global=nrow_global)
     621           18 :                CALL allocate_mo_set(active_space_env%mos_active(ispin), nao, max_orb_ind, 0, 0.0_dp, maxocc, 0.0_dp)
     622              :                CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     623           18 :                                         nrow_global=nrow_global, ncol_global=max_orb_ind)
     624           18 :                CALL init_mo_set(active_space_env%mos_active(ispin), fm_struct=fm_struct_tmp, name="Active Space MO")
     625           18 :                CALL cp_fm_struct_release(fm_struct_tmp)
     626              : 
     627              :                ! init inactive orbitals
     628           18 :                IF (nspins == 2) THEN
     629           12 :                   nel = nelec_inactive/2
     630              :                ELSE
     631            6 :                   nel = nelec_inactive
     632              :                END IF
     633           18 :                CALL allocate_mo_set(active_space_env%mos_inactive(ispin), nao, max_orb_ind, nel, REAL(nel, KIND=dp), maxocc, 0.0_dp)
     634              :                CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     635           18 :                                         nrow_global=nrow_global, ncol_global=max_orb_ind)
     636           18 :                CALL init_mo_set(active_space_env%mos_inactive(ispin), fm_struct=fm_struct_tmp, name="Inactive Space MO")
     637              :                ! small hack: set the correct inactive occupations down below
     638           68 :                active_space_env%mos_inactive(ispin)%occupation_numbers = 0.0_dp
     639           48 :                CALL cp_fm_struct_release(fm_struct_tmp)
     640              :             END DO
     641              : 
     642           48 :             ALLOCATE (eigenvalues(max_orb_ind, nspins))
     643           80 :             eigenvalues = 0.0_dp
     644           12 :             CALL get_qs_env(qs_env, matrix_ks=ks_matrix, matrix_s=s_matrix, scf_control=scf_control)
     645              : 
     646              :             ! calculate virtual MOs and copy inactive and active orbitals
     647           12 :             IF (iw > 0) THEN
     648            6 :                WRITE (iw, '(/,T3,A)') "Calculating virtual MOs..."
     649              :             END IF
     650           30 :             DO ispin = 1, nspins
     651           18 :                CALL get_mo_set(mos(ispin), mo_coeff=mo_ref, nmo=nmo_available)
     652           18 :                nmo_virtual = max_orb_ind - nmo_available
     653           18 :                nmo_virtual = MAX(nmo_virtual, 0)
     654              : 
     655              :                NULLIFY (evals_virtual)
     656           36 :                ALLOCATE (evals_virtual(nmo_virtual))
     657              : 
     658              :                CALL cp_fm_get_info(mo_ref, context=context, para_env=para_env, &
     659           18 :                                    nrow_global=nrow_global)
     660              : 
     661              :                CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=context, &
     662           18 :                                         nrow_global=nrow_global, ncol_global=nmo_virtual)
     663           18 :                CALL cp_fm_create(mo_virtual, fm_struct_tmp, name="virtual")
     664           18 :                CALL cp_fm_struct_release(fm_struct_tmp)
     665           18 :                CALL cp_fm_init_random(mo_virtual, nmo_virtual)
     666              : 
     667           18 :                NULLIFY (local_preconditioner)
     668              : 
     669              :                CALL ot_eigensolver(matrix_h=ks_matrix(ispin)%matrix, matrix_s=s_matrix(1)%matrix, &
     670              :                                    matrix_c_fm=mo_virtual, matrix_orthogonal_space_fm=mo_ref, &
     671              :                                    eps_gradient=scf_control%eps_lumos, &
     672              :                                    preconditioner=local_preconditioner, &
     673              :                                    iter_max=scf_control%max_iter_lumos, &
     674           18 :                                    size_ortho_space=nmo_available)
     675              : 
     676              :                CALL calculate_subspace_eigenvalues(mo_virtual, ks_matrix(ispin)%matrix, &
     677           18 :                                                    evals_virtual)
     678              : 
     679              :                ! We need to send the copy of MOs to preserve the sign
     680           18 :                CALL cp_fm_create(fm_dummy, mo_ref%matrix_struct)
     681           18 :                CALL cp_fm_to_fm(mo_ref, fm_dummy)
     682              : 
     683              :                CALL calculate_subspace_eigenvalues(fm_dummy, ks_matrix(ispin)%matrix, &
     684           18 :                                                    evals_arg=eigenvalues(:, ispin), do_rotation=.TRUE.)
     685              : 
     686           18 :                mo_set_active => active_space_env%mos_active(ispin)
     687           18 :                CALL get_mo_set(mo_set_active, mo_coeff=fm_target_active)
     688           18 :                mo_set_inactive => active_space_env%mos_inactive(ispin)
     689           18 :                CALL get_mo_set(mo_set_inactive, mo_coeff=fm_target_inactive)
     690              : 
     691              :                ! copy orbitals
     692           18 :                nmo_inactive_remaining = nmo_inactive
     693           68 :                DO i = 1, max_orb_ind
     694              :                   ! case for i being an active orbital
     695          114 :                   IF (ANY(active_space_env%active_orbitals(:, ispin) == i)) THEN
     696           36 :                      IF (i > nmo_available) THEN
     697            0 :                         CALL cp_fm_to_fm(mo_virtual, fm_target_active, 1, i - nmo_available, i)
     698            0 :                         eigenvalues(i, ispin) = evals_virtual(i - nmo_available)
     699            0 :                         mo_set_active%occupation_numbers(i) = 0.0
     700              :                      ELSE
     701           36 :                         CALL cp_fm_to_fm(fm_dummy, fm_target_active, 1, i, i)
     702           36 :                         mo_set_active%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
     703              :                      END IF
     704           36 :                      mo_set_active%eigenvalues(i) = eigenvalues(i, ispin)
     705              :                      ! if it was not an active orbital, check whether it is an inactive orbital
     706           14 :                   ELSEIF (nmo_inactive_remaining > 0) THEN
     707            0 :                      CALL cp_fm_to_fm(fm_dummy, fm_target_inactive, 1, i, i)
     708              :                      ! store on the fly the mapping of inactive orbitals
     709            0 :                      active_space_env%inactive_orbitals(nmo_inactive - nmo_inactive_remaining + 1, ispin) = i
     710            0 :                      mo_set_inactive%eigenvalues(i) = eigenvalues(i, ispin)
     711            0 :                      mo_set_inactive%occupation_numbers(i) = mos(ispin)%occupation_numbers(i)
     712              :                      ! hack: set homo and lumo manually
     713            0 :                      IF (nmo_inactive_remaining == 1) THEN
     714            0 :                         mo_set_inactive%homo = i
     715            0 :                         mo_set_inactive%lfomo = i + 1
     716              :                      END IF
     717            0 :                      nmo_inactive_remaining = nmo_inactive_remaining - 1
     718              :                   ELSE
     719           14 :                      CYCLE
     720              :                   END IF
     721              :                END DO
     722              : 
     723              :                ! Release
     724           18 :                DEALLOCATE (evals_virtual)
     725           18 :                CALL cp_fm_release(fm_dummy)
     726          102 :                CALL cp_fm_release(mo_virtual)
     727              :             END DO
     728              : 
     729           12 :             IF (iw > 0) THEN
     730           15 :                DO ispin = 1, nspins
     731            9 :                   WRITE (iw, '(/,T3,A,I3,T66,A)') "Orbital Energies and Selection for spin", ispin, "[atomic units]"
     732              : 
     733           18 :                   DO i = 1, max_orb_ind, 4
     734            9 :                      jm = MIN(3, max_orb_ind - i)
     735            9 :                      WRITE (iw, '(T4)', advance="no")
     736           34 :                      DO j = 0, jm
     737           57 :                         IF (ANY(active_space_env%active_orbitals(:, ispin) == i + j)) THEN
     738           18 :                            WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [A]"
     739            7 :                         ELSEIF (ANY(active_space_env%inactive_orbitals(:, ispin) == i + j)) THEN
     740            0 :                            WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [I]"
     741              :                         ELSE
     742            7 :                            WRITE (iw, '(T3,F12.6,A5)', advance="no") eigenvalues(i + j, ispin), " [V]"
     743              :                         END IF
     744              :                      END DO
     745           18 :                      WRITE (iw, *)
     746              :                   END DO
     747            9 :                   WRITE (iw, '(/,T3,A,I3)') "Active Orbital Indices for spin", ispin
     748           24 :                   DO i = 1, SIZE(active_space_env%active_orbitals, 1), 4
     749            9 :                      jm = MIN(3, SIZE(active_space_env%active_orbitals, 1) - i)
     750           36 :                      WRITE (iw, '(T3,4(I4))') (active_space_env%active_orbitals(i + j, ispin), j=0, jm)
     751              :                   END DO
     752              :                END DO
     753              :             END IF
     754           24 :             DEALLOCATE (eigenvalues)
     755              :          END IF
     756              : 
     757              :       CASE (wannier_projection)
     758            0 :          NULLIFY (loc_section, loc_print)
     759            0 :          loc_section => section_vals_get_subs_vals(as_input, "LOCALIZE")
     760            0 :          CPASSERT(ASSOCIATED(loc_section))
     761            0 :          loc_print => section_vals_get_subs_vals(as_input, "LOCALIZE%PRINT")
     762              :          !
     763            0 :          CPABORT("not yet available")
     764              :          !
     765              :       CASE (mao_projection)
     766              :          !
     767           62 :          CPABORT("not yet available")
     768              :          !
     769              :       END SELECT
     770              : 
     771              :       ! Print orbitals on Cube files
     772           62 :       print_orb => section_vals_get_subs_vals(as_input, "PRINT_ORBITAL_CUBES")
     773           62 :       CALL section_vals_get(print_orb, explicit=explicit)
     774           62 :       CALL section_vals_val_get(print_orb, "STOP_AFTER_CUBES", l_val=stop_after_print)
     775           62 :       IF (explicit) THEN
     776              :          !
     777            2 :          CALL print_orbital_cubes(print_orb, qs_env, active_space_env%mos_active)
     778              :          !
     779            2 :          IF (stop_after_print) THEN
     780              : 
     781            0 :             IF (iw > 0) THEN
     782              :                WRITE (iw, '(/,T2,A)') &
     783            0 :                   '!----------------- Early End of Active Space Interface -----------------------!'
     784              :             END IF
     785              : 
     786            0 :             CALL timestop(handle)
     787              : 
     788            0 :             RETURN
     789              :          END IF
     790              :       END IF
     791              : 
     792              :       ! calculate inactive density matrix
     793           62 :       CALL get_qs_env(qs_env, rho=rho)
     794           62 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
     795           62 :       CPASSERT(ASSOCIATED(rho_ao))
     796           62 :       CALL dbcsr_allocate_matrix_set(active_space_env%pmat_inactive, nspins)
     797          136 :       DO ispin = 1, nspins
     798           74 :          ALLOCATE (denmat)
     799           74 :          CALL dbcsr_copy(denmat, rho_ao(ispin)%matrix)
     800           74 :          mo_set => active_space_env%mos_inactive(ispin)
     801           74 :          CALL calculate_density_matrix(mo_set, denmat)
     802          136 :          active_space_env%pmat_inactive(ispin)%matrix => denmat
     803              :       END DO
     804              : 
     805              :       ! generate integrals
     806              :       ! make sure that defaults are set correctly (from basic calculation)
     807              :       ! make sure that periodicity is consistent
     808           62 :       CALL section_vals_val_get(as_input, "ERI%METHOD", i_val=eri_method)
     809           62 :       active_space_env%eri%method = eri_method
     810           62 :       CALL section_vals_val_get(as_input, "ERI%OPERATOR", i_val=eri_operator, explicit=ex_operator)
     811           62 :       active_space_env%eri%operator = eri_operator
     812           62 :       CALL section_vals_val_get(as_input, "ERI%OPERATOR_PARAMETER", r_val=eri_op_param)
     813           62 :       active_space_env%eri%operator_parameter = eri_op_param
     814           62 :       CALL section_vals_val_get(as_input, "ERI%CUTOFF_RADIUS", r_val=eri_rcut)
     815           62 :       active_space_env%eri%cutoff_radius = eri_rcut
     816           62 :       CALL section_vals_val_get(as_input, "ERI%PERIODICITY", i_vals=invals, explicit=ex_perd)
     817           62 :       CALL section_vals_val_get(as_input, "ERI%EPS_INTEGRAL", r_val=eri_eps_int)
     818           62 :       active_space_env%eri%eps_integral = eri_eps_int
     819           62 :       IF (active_space_env%molecule) THEN
     820              :          ! check that we are in a non-periodic setting
     821           58 :          CALL get_qs_env(qs_env, cell=cell)
     822          232 :          IF (SUM(cell%perd) /= 0) THEN
     823            0 :             CPABORT("Active space option ISOLATED_SYSTEM requires non-periodic setting")
     824              :          END IF
     825           58 :          IF (ex_perd) THEN
     826          232 :             IF (SUM(invals) /= 0) THEN
     827            0 :                CPABORT("Active space option ISOLATED_SYSTEM requires non-periodic setting")
     828              :             END IF
     829              :          END IF
     830          232 :          active_space_env%eri%periodicity(1:3) = 0
     831            4 :       ELSE IF (ex_perd) THEN
     832            0 :          IF (SIZE(invals) == 1) THEN
     833            0 :             active_space_env%eri%periodicity(1:3) = invals(1)
     834              :          ELSE
     835            0 :             active_space_env%eri%periodicity(1:3) = invals(1:3)
     836              :          END IF
     837              :       ELSE
     838            4 :          CALL get_qs_env(qs_env, cell=cell)
     839           28 :          active_space_env%eri%periodicity(1:3) = cell%perd(1:3)
     840              :       END IF
     841           62 :       IF (iw > 0) THEN
     842           31 :          WRITE (iw, '(/,T3,A)') "Calculation of Electron Repulsion Integrals"
     843           25 :          SELECT CASE (eri_method)
     844              :          CASE (eri_method_full_gpw)
     845           25 :             WRITE (iw, '(T3,A,T50,A)') "Integration method", "GPW Fourier transform over MOs"
     846              :          CASE (eri_method_gpw_ht)
     847            6 :             WRITE (iw, '(T3,A,T44,A)') "Integration method", "Half transformed integrals from GPW"
     848              :          CASE DEFAULT
     849           31 :             CPABORT("Unknown ERI method")
     850              :          END SELECT
     851           24 :          SELECT CASE (eri_operator)
     852              :          CASE (eri_operator_coulomb)
     853           24 :             WRITE (iw, '(T3,A,T66,A)') "ERI operator", "Coulomb <1/R>"
     854              :          CASE (eri_operator_yukawa)
     855            0 :             WRITE (iw, '(T3,A,T59,A)') "ERI operator", "Yukawa <EXP(-a*R)/R>"
     856            0 :             WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param
     857              :          CASE (eri_operator_erf)
     858            7 :             WRITE (iw, '(T3,A,T53,A)') "ERI operator", "Error function <ERF(a*R)/R>"
     859            7 :             WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param
     860              :          CASE (eri_operator_erfc)
     861            0 :             WRITE (iw, '(T3,A,T45,A)') "ERI operator", "Compl. error function <ERFC(a*R)/R>"
     862            0 :             WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param
     863              :          CASE (eri_operator_gaussian)
     864            0 :             WRITE (iw, '(T3,A,T41,A)') "ERI operator", "Gaussian attenuated <EXP(-a*R^2)/R>"
     865            0 :             WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param
     866              :          CASE (eri_operator_trunc)
     867            0 :             WRITE (iw, '(T3,A,T53,A)') "ERI operator", "Truncated Coulomb <H(a-R)/R>"
     868            0 :             WRITE (iw, '(T3,A,T66,F14.3)') "ERI operator parameter", eri_op_param
     869              :          CASE DEFAULT
     870           31 :             CPABORT("Unknown ERI operator")
     871              :          END SELECT
     872           31 :          WRITE (iw, '(T3,A,T68,E12.4)') "Accuracy of ERI", eri_eps_int
     873           31 :          WRITE (iw, '(T3,A,T71,3I3)') "Periodicity", active_space_env%eri%periodicity(1:3)
     874          124 :          IF (PRODUCT(active_space_env%eri%periodicity(1:3)) == 0) THEN
     875           30 :             IF (eri_rcut > 0.0_dp) WRITE (iw, '(T3,A,T65,F14.6)') "Periodicity (Cutoff)", eri_rcut
     876              :          END IF
     877           31 :          IF (nspins < 2) THEN
     878           25 :             WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI", (nmo_active**4)/8
     879              :          ELSE
     880            6 :             WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|aa)", (nmo_active**4)/8
     881            6 :             WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (bb|bb)", (nmo_active**4)/8
     882            6 :             WRITE (iw, '(T3,A,T68,I12)') "Total Number of ERI (aa|bb)", (nmo_active**4)/4
     883              :          END IF
     884              :       END IF
     885              : 
     886              :       ! allocate container for integrals (CSR matrix)
     887           62 :       CALL get_qs_env(qs_env, para_env=para_env)
     888           62 :       m = (nspins*(nspins + 1))/2
     889          272 :       ALLOCATE (active_space_env%eri%eri(m))
     890          148 :       DO i = 1, m
     891           86 :          CALL get_mo_set(active_space_env%mos_active(1), nmo=nmo)
     892           86 :          ALLOCATE (active_space_env%eri%eri(i)%csr_mat)
     893           86 :          eri_mat => active_space_env%eri%eri(i)%csr_mat
     894           86 :          IF (i == 1) THEN
     895           62 :             n1 = nmo
     896           62 :             n2 = nmo
     897           24 :          ELSEIF (i == 2) THEN
     898           12 :             n1 = nmo
     899           12 :             n2 = nmo
     900              :          ELSE
     901           12 :             n1 = nmo
     902           12 :             n2 = nmo
     903              :          END IF
     904           86 :          nn1 = (n1*(n1 + 1))/2
     905           86 :          nn2 = (n2*(n2 + 1))/2
     906           86 :          CALL dbcsr_csr_create(eri_mat, nn1, nn2, 0_int_8, 0, 0, para_env%get_handle())
     907          234 :          active_space_env%eri%norb = nmo
     908              :       END DO
     909              : 
     910           62 :       SELECT CASE (eri_method)
     911              :       CASE (eri_method_full_gpw, eri_method_gpw_ht)
     912           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%EPS_GRID", r_val=eri_eps_grid)
     913           62 :          active_space_env%eri%eri_gpw%eps_grid = eri_eps_grid
     914           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%EPS_FILTER", r_val=eri_eps_filter)
     915           62 :          active_space_env%eri%eri_gpw%eps_filter = eri_eps_filter
     916           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%CUTOFF", r_val=eri_gpw_cutoff)
     917           62 :          active_space_env%eri%eri_gpw%cutoff = eri_gpw_cutoff
     918           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%REL_CUTOFF", r_val=eri_rel_cutoff)
     919           62 :          active_space_env%eri%eri_gpw%rel_cutoff = eri_rel_cutoff
     920           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%PRINT_LEVEL", i_val=eri_print)
     921           62 :          active_space_env%eri%eri_gpw%print_level = eri_print
     922           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%STORE_WFN", l_val=store_wfn)
     923           62 :          active_space_env%eri%eri_gpw%store_wfn = store_wfn
     924           62 :          CALL section_vals_val_get(as_input, "ERI_GPW%GROUP_SIZE", i_val=group_size)
     925           62 :          active_space_env%eri%eri_gpw%group_size = group_size
     926              :          ! Always redo Poisson solver for now
     927           62 :          active_space_env%eri%eri_gpw%redo_poisson = .TRUE.
     928              :          ! active_space_env%eri%eri_gpw%redo_poisson = (ex_operator .OR. ex_perd)
     929           62 :          IF (iw > 0) THEN
     930           31 :             WRITE (iw, '(/,T2,A,T71,F10.1)') "ERI_GPW| Energy cutoff [Ry]", eri_gpw_cutoff
     931           31 :             WRITE (iw, '(T2,A,T71,F10.1)') "ERI_GPW| Relative energy cutoff [Ry]", eri_rel_cutoff
     932              :          END IF
     933              :          !
     934           62 :          CALL calculate_eri_gpw(active_space_env%mos_active, active_space_env%active_orbitals, active_space_env%eri, qs_env, iw)
     935              :          !
     936              :       CASE DEFAULT
     937           62 :          CPABORT("Unknown ERI method")
     938              :       END SELECT
     939           62 :       IF (iw > 0) THEN
     940           74 :          DO isp = 1, SIZE(active_space_env%eri%eri)
     941           43 :             eri_mat => active_space_env%eri%eri(isp)%csr_mat
     942              :             nze_percentage = 100.0_dp*(REAL(eri_mat%nze_total, KIND=dp) &
     943           43 :                                        /REAL(eri_mat%nrows_total, KIND=dp))/REAL(eri_mat%ncols_total, KIND=dp)
     944           43 :             WRITE (iw, '(/,T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
     945           86 :                "Number of  CSR non-zero elements:", eri_mat%nze_total
     946           43 :             WRITE (iw, '(T2,A,I2,T30,A,T68,F12.4)') "ERI_GPW| Spinmatrix:", isp, &
     947           86 :                "Percentage CSR non-zero elements:", nze_percentage
     948           43 :             WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
     949           86 :                "nrows_total", eri_mat%nrows_total
     950           43 :             WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
     951           86 :                "ncols_total", eri_mat%ncols_total
     952           43 :             WRITE (iw, '(T2,A,I2,T30,A,T68,I12)') "ERI_GPW| Spinmatrix:", isp, &
     953          117 :                "nrows_local", eri_mat%nrows_local
     954              :          END DO
     955              :       END IF
     956              : 
     957              :       ! set the reference active space density matrix
     958           62 :       nspins = active_space_env%nspins
     959          260 :       ALLOCATE (active_space_env%p_active(nspins))
     960          136 :       DO isp = 1, nspins
     961           74 :          mo_set => active_space_env%mos_active(isp)
     962           74 :          CALL get_mo_set(mo_set, mo_coeff=mo_coeff, nmo=nmo)
     963          136 :          CALL create_subspace_matrix(mo_coeff, active_space_env%p_active(isp), nmo)
     964              :       END DO
     965            0 :       SELECT CASE (mselect)
     966              :       CASE DEFAULT
     967            0 :          CPABORT("Unknown orbital selection method")
     968              :       CASE (casci_canonical, manual_selection)
     969           62 :          focc = 2.0_dp
     970           62 :          IF (nspins == 2) focc = 1.0_dp
     971          136 :          DO isp = 1, nspins
     972           74 :             fmat => active_space_env%p_active(isp)
     973           74 :             CALL cp_fm_set_all(fmat, alpha=0.0_dp)
     974           74 :             IF (nspins == 2) THEN
     975           24 :                IF (isp == 1) THEN
     976           12 :                   nel = (active_space_env%nelec_active + active_space_env%multiplicity - 1)/2
     977              :                ELSE
     978           12 :                   nel = (active_space_env%nelec_active - active_space_env%multiplicity + 1)/2
     979              :                END IF
     980              :             ELSE
     981           50 :                nel = active_space_env%nelec_active
     982              :             END IF
     983          324 :             DO i = 1, nmo_active
     984          188 :                m = active_space_env%active_orbitals(i, isp)
     985          188 :                fel = MIN(focc, REAL(nel, KIND=dp))
     986          188 :                CALL cp_fm_set_element(fmat, m, m, fel)
     987          188 :                nel = nel - NINT(fel)
     988          262 :                nel = MAX(nel, 0)
     989              :             END DO
     990              :          END DO
     991              :       CASE (wannier_projection)
     992            0 :          CPABORT("NOT IMPLEMENTED")
     993              :       CASE (mao_projection)
     994           62 :          CPABORT("NOT IMPLEMENTED")
     995              :       END SELECT
     996              : 
     997              :       ! transform KS/Fock, Vxc and Hcore to AS MO basis
     998           62 :       CALL calculate_operators(active_space_env%mos_active, qs_env, active_space_env)
     999              :       ! set the reference energy in the active space
    1000           62 :       CALL get_qs_env(qs_env, energy=energy)
    1001           62 :       active_space_env%energy_ref = energy%total
    1002              :       ! calculate inactive energy and embedding potential
    1003           62 :       CALL subspace_fock_matrix(active_space_env)
    1004              : 
    1005              :       ! associate the active space environment with the qs environment
    1006           62 :       CALL set_qs_env(qs_env, active_space=active_space_env)
    1007              : 
    1008              :       ! Perform the embedding calculation only if qiskit is specified
    1009           62 :       CALL section_vals_val_get(as_input, "AS_SOLVER", i_val=as_solver)
    1010           62 :       SELECT CASE (as_solver)
    1011              :       CASE (no_solver)
    1012           62 :          IF (iw > 0) THEN
    1013           31 :             WRITE (iw, '(/,T3,A)') "No active space solver specified, skipping embedding calculation"
    1014              :          END IF
    1015              :       CASE (qiskit_solver)
    1016            0 :          CALL rsdft_embedding(qs_env, active_space_env, as_input)
    1017              :       CASE DEFAULT
    1018           62 :          CPABORT("Unknown active space solver")
    1019              :       END SELECT
    1020              : 
    1021              :       ! Output a FCIDUMP file if requested
    1022           62 :       IF (active_space_env%fcidump) CALL fcidump(active_space_env, as_input)
    1023              : 
    1024              :       ! Output a QCSchema file if requested
    1025           62 :       IF (active_space_env%qcschema) THEN
    1026            0 :          CALL qcschema_env_create(qcschema_env, qs_env)
    1027            0 :          CALL qcschema_to_hdf5(qcschema_env, active_space_env%qcschema_filename)
    1028            0 :          CALL qcschema_env_release(qcschema_env)
    1029              :       END IF
    1030              : 
    1031           62 :       IF (iw > 0) THEN
    1032              :          WRITE (iw, '(/,T2,A)') &
    1033           31 :             '!-------------------- End of Active Space Interface --------------------------!'
    1034              :       END IF
    1035              : 
    1036           62 :       CALL timestop(handle)
    1037              : 
    1038        87486 :    END SUBROUTINE active_space_main
    1039              : 
    1040              : ! **************************************************************************************************
    1041              : !> \brief computes the one-electron operators in the subspace of the provided orbital set
    1042              : !> \param mos the molecular orbital set within the active subspace
    1043              : !> \param qs_env ...
    1044              : !> \param active_space_env ...
    1045              : !> \par History
    1046              : !>      04.2016 created [JGH]
    1047              : ! **************************************************************************************************
    1048           62 :    SUBROUTINE calculate_operators(mos, qs_env, active_space_env)
    1049              : 
    1050              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1051              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1052              :       TYPE(active_space_type), POINTER                   :: active_space_env
    1053              : 
    1054              :       CHARACTER(len=*), PARAMETER :: routineN = 'calculate_operators'
    1055              : 
    1056              :       INTEGER                                            :: handle, is, nmo, nspins
    1057              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1058           62 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_vxc
    1059           62 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: h_matrix, ks_matrix
    1060              : 
    1061           62 :       CALL timeset(routineN, handle)
    1062              : 
    1063           62 :       nspins = active_space_env%nspins
    1064              : 
    1065              :       ! Kohn-Sham / Fock operator
    1066           62 :       CALL cp_fm_release(active_space_env%ks_sub)
    1067           62 :       CALL get_qs_env(qs_env, matrix_ks_kp=ks_matrix)
    1068           62 :       IF (SIZE(ks_matrix, 2) > 1) THEN
    1069            0 :          CPABORT("No k-points allowed at this point")
    1070              :       END IF
    1071          260 :       ALLOCATE (active_space_env%ks_sub(nspins))
    1072          136 :       DO is = 1, nspins
    1073           74 :          CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
    1074          136 :          CALL subspace_operator(mo_coeff, nmo, ks_matrix(is, 1)%matrix, active_space_env%ks_sub(is))
    1075              :       END DO
    1076              : 
    1077              :       ! Vxc matrix
    1078           62 :       CALL cp_fm_release(active_space_env%vxc_sub)
    1079              : 
    1080           62 :       NULLIFY (matrix_vxc)
    1081           62 :       CALL get_qs_env(qs_env, matrix_vxc=matrix_vxc)
    1082           62 :       IF (ASSOCIATED(matrix_vxc)) THEN
    1083            0 :          ALLOCATE (active_space_env%vxc_sub(nspins))
    1084            0 :          DO is = 1, nspins
    1085            0 :             CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
    1086            0 :             CALL subspace_operator(mo_coeff, nmo, matrix_vxc(is)%matrix, active_space_env%vxc_sub(is))
    1087              :          END DO
    1088              :       END IF
    1089              : 
    1090              :       ! Core Hamiltonian
    1091           62 :       CALL cp_fm_release(active_space_env%h_sub)
    1092              : 
    1093           62 :       NULLIFY (h_matrix)
    1094           62 :       CALL get_qs_env(qs_env=qs_env, matrix_h_kp=h_matrix)
    1095           62 :       IF (SIZE(h_matrix, 2) > 1) THEN
    1096            0 :          CPABORT("No k-points allowed at this point")
    1097              :       END IF
    1098          260 :       ALLOCATE (active_space_env%h_sub(nspins))
    1099          136 :       DO is = 1, nspins
    1100           74 :          CALL get_mo_set(mo_set=mos(is), mo_coeff=mo_coeff, nmo=nmo)
    1101          136 :          CALL subspace_operator(mo_coeff, nmo, h_matrix(1, 1)%matrix, active_space_env%h_sub(is))
    1102              :       END DO
    1103              : 
    1104           62 :       CALL timestop(handle)
    1105              : 
    1106           62 :    END SUBROUTINE calculate_operators
    1107              : 
    1108              : ! **************************************************************************************************
    1109              : !> \brief computes a one-electron operator in the subspace of the provided orbital set
    1110              : !> \param orbitals the orbital coefficient matrix
    1111              : !> \param nmo the number of orbitals
    1112              : !> \param op_matrix operator matrix in AO basis
    1113              : !> \param op_sub operator in orbital basis
    1114              : !> \par History
    1115              : !>      04.2016 created [JGH]
    1116              : ! **************************************************************************************************
    1117          296 :    SUBROUTINE subspace_operator(orbitals, nmo, op_matrix, op_sub)
    1118              : 
    1119              :       TYPE(cp_fm_type), INTENT(IN)                       :: orbitals
    1120              :       INTEGER, INTENT(IN)                                :: nmo
    1121              :       TYPE(dbcsr_type), POINTER                          :: op_matrix
    1122              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: op_sub
    1123              : 
    1124              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'subspace_operator'
    1125              : 
    1126              :       INTEGER                                            :: handle, ncol, nrow
    1127              :       TYPE(cp_fm_type)                                   :: vectors
    1128              : 
    1129          148 :       CALL timeset(routineN, handle)
    1130              : 
    1131          148 :       CALL cp_fm_get_info(matrix=orbitals, ncol_global=ncol, nrow_global=nrow)
    1132          148 :       CPASSERT(nmo <= ncol)
    1133              : 
    1134          148 :       IF (nmo > 0) THEN
    1135              : 
    1136          148 :          CALL cp_fm_create(vectors, orbitals%matrix_struct, "vectors")
    1137              : 
    1138          148 :          CALL create_subspace_matrix(orbitals, op_sub, nmo)
    1139              : 
    1140          148 :          CALL cp_dbcsr_sm_fm_multiply(op_matrix, orbitals, vectors, nmo)
    1141              : 
    1142          148 :          CALL parallel_gemm('T', 'N', nmo, nmo, nrow, 1.0_dp, orbitals, vectors, 0.0_dp, op_sub)
    1143              : 
    1144          148 :          CALL cp_fm_release(vectors)
    1145              : 
    1146              :       END IF
    1147              : 
    1148          148 :       CALL timestop(handle)
    1149              : 
    1150          148 :    END SUBROUTINE subspace_operator
    1151              : 
    1152              : ! **************************************************************************************************
    1153              : !> \brief creates a matrix of subspace size
    1154              : !> \param orbitals the orbital coefficient matrix
    1155              : !> \param op_sub operator in orbital basis
    1156              : !> \param n the number of orbitals
    1157              : !> \par History
    1158              : !>      04.2016 created [JGH]
    1159              : ! **************************************************************************************************
    1160          222 :    SUBROUTINE create_subspace_matrix(orbitals, op_sub, n)
    1161              : 
    1162              :       TYPE(cp_fm_type), INTENT(IN)                       :: orbitals
    1163              :       TYPE(cp_fm_type), INTENT(OUT)                      :: op_sub
    1164              :       INTEGER, INTENT(IN)                                :: n
    1165              : 
    1166              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    1167              : 
    1168          222 :       IF (n > 0) THEN
    1169              : 
    1170          222 :          NULLIFY (fm_struct)
    1171              :          CALL cp_fm_struct_create(fm_struct, nrow_global=n, ncol_global=n, &
    1172              :                                   para_env=orbitals%matrix_struct%para_env, &
    1173          222 :                                   context=orbitals%matrix_struct%context)
    1174          222 :          CALL cp_fm_create(op_sub, fm_struct, name="Subspace operator")
    1175          222 :          CALL cp_fm_struct_release(fm_struct)
    1176              : 
    1177              :       END IF
    1178              : 
    1179          222 :    END SUBROUTINE create_subspace_matrix
    1180              : 
    1181              : ! **************************************************************************************************
    1182              : !> \brief computes a electron repulsion integrals using GPW technology
    1183              : !> \param mos the molecular orbital set within the active subspace
    1184              : !> \param orbitals ...
    1185              : !> \param eri_env ...
    1186              : !> \param qs_env ...
    1187              : !> \param iw ...
    1188              : !> \par History
    1189              : !>      04.2016 created [JGH]
    1190              : ! **************************************************************************************************
    1191           62 :    SUBROUTINE calculate_eri_gpw(mos, orbitals, eri_env, qs_env, iw)
    1192              : 
    1193              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1194              :       INTEGER, DIMENSION(:, :), POINTER                  :: orbitals
    1195              :       TYPE(eri_type)                                     :: eri_env
    1196              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1197              :       INTEGER, INTENT(IN)                                :: iw
    1198              : 
    1199              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'calculate_eri_gpw'
    1200              : 
    1201              :       INTEGER :: col_local, color, handle, i1, i2, i3, i4, i_multigrid, icount2, intcount, &
    1202              :          irange(2), irange_sub(2), isp, isp1, isp2, ispin, iwa1, iwa12, iwa2, iwb1, iwb12, iwb2, &
    1203              :          iwbs, iwbt, iwfn, n_multigrid, ncol_global, ncol_local, nmm, nmo, nmo1, nmo2, &
    1204              :          nrow_global, nrow_local, nspins, number_of_subgroups, nx, row_local, stored_integrals
    1205           62 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: eri_index
    1206           62 :       INTEGER, DIMENSION(:), POINTER                     :: col_indices, row_indices
    1207              :       LOGICAL                                            :: print1, print2, &
    1208              :                                                             skip_load_balance_distributed
    1209              :       REAL(KIND=dp)                                      :: dvol, erint, pair_int, &
    1210              :                                                             progression_factor, rc, rsize
    1211           62 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eri
    1212           62 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1213              :       TYPE(cell_type), POINTER                           :: cell
    1214              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env, blacs_env_sub
    1215              :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct
    1216           62 :       TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:)        :: fm_matrix_pq_rnu, fm_matrix_pq_rs, &
    1217           62 :                                                             fm_mo_coeff_as
    1218              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff, mo_coeff1, mo_coeff2
    1219              :       TYPE(dbcsr_p_type)                                 :: mat_munu
    1220           62 :       TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:)        :: matrix_pq_rnu, mo_coeff_as
    1221              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1222              :       TYPE(mp_comm_type)                                 :: mp_group
    1223              :       TYPE(mp_para_env_type), POINTER                    :: para_env, para_env_sub
    1224              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
    1225           62 :          POINTER                                         :: sab_orb_sub
    1226           62 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1227              :       TYPE(pw_c1d_gs_type)                               :: pot_g, rho_g
    1228              :       TYPE(pw_env_type), POINTER                         :: pw_env_sub
    1229              :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
    1230              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1231              :       TYPE(pw_r3d_rs_type)                               :: rho_r, wfn_r
    1232              :       TYPE(pw_r3d_rs_type), ALLOCATABLE, &
    1233           62 :          DIMENSION(:, :), TARGET                         :: wfn_a
    1234              :       TYPE(pw_r3d_rs_type), POINTER                      :: wfn1, wfn2, wfn3, wfn4
    1235              :       TYPE(qs_control_type), POINTER                     :: qs_control, qs_control_old
    1236           62 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1237              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
    1238              :       TYPE(task_list_type), POINTER                      :: task_list_sub
    1239              : 
    1240           62 :       CALL timeset(routineN, handle)
    1241              : 
    1242              :       ! print levels
    1243          120 :       SELECT CASE (eri_env%eri_gpw%print_level)
    1244              :       CASE (silent_print_level)
    1245           58 :          print1 = .FALSE.
    1246           58 :          print2 = .FALSE.
    1247              :       CASE (low_print_level)
    1248            2 :          print1 = .FALSE.
    1249            2 :          print2 = .FALSE.
    1250              :       CASE (medium_print_level)
    1251            2 :          print1 = .TRUE.
    1252            2 :          print2 = .FALSE.
    1253              :       CASE (high_print_level)
    1254            0 :          print1 = .TRUE.
    1255            0 :          print2 = .TRUE.
    1256              :       CASE (debug_print_level)
    1257            0 :          print1 = .TRUE.
    1258           62 :          print2 = .TRUE.
    1259              :       CASE DEFAULT
    1260              :          ! do nothing
    1261              :       END SELECT
    1262              : 
    1263              :       ! Check the inpuz group
    1264           62 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
    1265           62 :       IF (eri_env%eri_gpw%group_size <= 1) eri_env%eri_gpw%group_size = para_env%num_pe
    1266           62 :       IF (MOD(para_env%num_pe, eri_env%eri_gpw%group_size) /= 0) &
    1267            0 :          CPABORT("Group size must be a divisor of the total number of processes!")
    1268              :       ! Create a new para_env or reuse the old one
    1269           62 :       IF (eri_env%eri_gpw%group_size == para_env%num_pe) THEN
    1270           62 :          para_env_sub => para_env
    1271           62 :          CALL para_env_sub%retain()
    1272           62 :          IF (eri_env%method == eri_method_gpw_ht) THEN
    1273           12 :             blacs_env_sub => blacs_env
    1274           12 :             CALL blacs_env_sub%retain()
    1275              :          END IF
    1276           62 :          number_of_subgroups = 1
    1277           62 :          color = 0
    1278              :       ELSE
    1279            0 :          number_of_subgroups = para_env%num_pe/eri_env%eri_gpw%group_size
    1280            0 :          color = para_env%mepos/eri_env%eri_gpw%group_size
    1281            0 :          ALLOCATE (para_env_sub)
    1282            0 :          CALL para_env_sub%from_split(para_env, color)
    1283            0 :          IF (eri_env%method == eri_method_gpw_ht) THEN
    1284            0 :             NULLIFY (blacs_env_sub)
    1285            0 :             CALL cp_blacs_env_create(blacs_env_sub, para_env_sub, BLACS_GRID_SQUARE, .TRUE.)
    1286              :          END IF
    1287              :       END IF
    1288              : 
    1289              :       ! This should be done differently! Copied from MP2 code
    1290           62 :       CALL get_qs_env(qs_env, dft_control=dft_control)
    1291          186 :       ALLOCATE (qs_control)
    1292           62 :       qs_control_old => dft_control%qs_control
    1293           62 :       qs_control = qs_control_old
    1294           62 :       dft_control%qs_control => qs_control
    1295           62 :       progression_factor = qs_control%progression_factor
    1296           62 :       n_multigrid = SIZE(qs_control%e_cutoff)
    1297           62 :       nspins = SIZE(mos)
    1298              :       ! Allocate new cutoffs (just in private qs_control, not in qs_control_old)
    1299          186 :       ALLOCATE (qs_control%e_cutoff(n_multigrid))
    1300              : 
    1301           62 :       qs_control%cutoff = eri_env%eri_gpw%cutoff*0.5_dp
    1302           62 :       qs_control%e_cutoff(1) = qs_control%cutoff
    1303          248 :       DO i_multigrid = 2, n_multigrid
    1304              :          qs_control%e_cutoff(i_multigrid) = qs_control%e_cutoff(i_multigrid - 1) &
    1305          248 :                                             /progression_factor
    1306              :       END DO
    1307           62 :       qs_control%relative_cutoff = eri_env%eri_gpw%rel_cutoff*0.5_dp
    1308              : 
    1309           62 :       IF (eri_env%method == eri_method_gpw_ht) THEN
    1310              :          ! For now, we will distribute neighbor lists etc. within the global communicator
    1311           12 :          CALL get_qs_env(qs_env, ks_env=ks_env)
    1312              :          CALL create_mat_munu(mat_munu, qs_env, eri_env%eri_gpw%eps_grid, blacs_env_sub, sab_orb_sub=sab_orb_sub, &
    1313           12 :                               do_alloc_blocks_from_nbl=.TRUE., dbcsr_sym_type=dbcsr_type_symmetric)
    1314           12 :          CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
    1315              :       END IF
    1316              : 
    1317              :       ! Generate the appropriate  pw_env
    1318           62 :       NULLIFY (pw_env_sub)
    1319           62 :       CALL pw_env_create(pw_env_sub)
    1320           62 :       CALL pw_env_rebuild(pw_env_sub, qs_env, external_para_env=para_env_sub)
    1321              :       CALL pw_env_get(pw_env_sub, auxbas_pw_pool=auxbas_pw_pool, &
    1322           62 :                       poisson_env=poisson_env)
    1323           62 :       IF (eri_env%eri_gpw%redo_poisson) THEN
    1324              :          ! We need to rebuild the Poisson solver on the fly
    1325          248 :          IF (SUM(eri_env%periodicity) /= 0) THEN
    1326            2 :             poisson_env%parameters%solver = pw_poisson_periodic
    1327              :          ELSE
    1328           60 :             poisson_env%parameters%solver = pw_poisson_analytic
    1329              :          END IF
    1330          248 :          poisson_env%parameters%periodic = eri_env%periodicity
    1331           62 :          CALL pw_poisson_rebuild(poisson_env)
    1332           62 :          IF (eri_env%cutoff_radius > 0.0_dp) THEN
    1333            0 :             poisson_env%green_fft%radius = eri_env%cutoff_radius
    1334              :          ELSE
    1335           62 :             CALL get_qs_env(qs_env, cell=cell)
    1336           62 :             rc = cell%hmat(1, 1)
    1337          248 :             DO iwa1 = 1, 3
    1338          248 :                rc = MIN(rc, 0.5_dp*cell%hmat(iwa1, iwa1))
    1339              :             END DO
    1340           62 :             poisson_env%green_fft%radius = rc
    1341              :          END IF
    1342              : 
    1343           62 :          CALL pw_eri_green_create(poisson_env%green_fft, eri_env)
    1344              : 
    1345           62 :          IF (iw > 0) THEN
    1346           31 :             CALL get_qs_env(qs_env, cell=cell)
    1347          217 :             IF (SUM(cell%perd) /= SUM(eri_env%periodicity)) THEN
    1348            0 :                IF (SUM(eri_env%periodicity) /= 0) THEN
    1349              :                   WRITE (UNIT=iw, FMT="(/,T2,A,T51,A30)") &
    1350            0 :                      "ERI_GPW| Switching Poisson solver to", "PERIODIC"
    1351              :                ELSE
    1352              :                   WRITE (UNIT=iw, FMT="(/,T2,A,T51,A30)") &
    1353            0 :                      "ERI_GPW| Switching Poisson solver to", "ANALYTIC"
    1354              :                END IF
    1355              :             END IF
    1356              :             ! print out the Greens function to check it matches the Poisson solver
    1357           32 :             SELECT CASE (poisson_env%green_fft%method)
    1358              :             CASE (PERIODIC3D)
    1359              :                WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") &
    1360            1 :                   "ERI_GPW| Poisson Greens function", "PERIODIC"
    1361              :             CASE (ANALYTIC0D)
    1362              :                WRITE (UNIT=iw, FMT="(T2,A,T51,A30)") &
    1363           30 :                   "ERI_GPW| Poisson Greens function", "ANALYTIC"
    1364              :             CASE DEFAULT
    1365           31 :                CPABORT("Wrong Greens function setup")
    1366              :             END SELECT
    1367              :          END IF
    1368              :       END IF
    1369              : 
    1370           62 :       IF (eri_env%method == eri_method_gpw_ht) THEN
    1371              :          ! We need a task list
    1372           12 :          NULLIFY (task_list_sub)
    1373           12 :          skip_load_balance_distributed = dft_control%qs_control%skip_load_balance_distributed
    1374           12 :          CALL allocate_task_list(task_list_sub)
    1375              :          CALL generate_qs_task_list(ks_env, task_list_sub, &
    1376              :                                     reorder_rs_grid_ranks=.TRUE., soft_valid=.FALSE., &
    1377              :                                     skip_load_balance_distributed=skip_load_balance_distributed, &
    1378           12 :                                     pw_env_external=pw_env_sub, sab_orb_external=sab_orb_sub)
    1379              : 
    1380              :          ! Create sparse matrices carrying the matrix products, Code borrowed from the MP2 GPW method
    1381              :          ! Create equal distributions for them (no sparsity present)
    1382              :          ! We use the routines from mp2 suggesting that one may replicate the grids later for better performance
    1383              :          ALLOCATE (mo_coeff_as(nspins), matrix_pq_rnu(nspins), &
    1384          156 :                    fm_matrix_pq_rnu(nspins), fm_mo_coeff_as(nspins), fm_matrix_pq_rs(nspins))
    1385           24 :          DO ispin = 1, nspins
    1386           36 :             BLOCK
    1387           12 :                REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: C
    1388           12 :                TYPE(group_dist_d1_type) :: gd_array
    1389              :                TYPE(cp_fm_type), POINTER :: mo_coeff
    1390           12 :                CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
    1391              :                CALL grep_rows_in_subgroups(para_env, para_env_sub, mo_coeff, gd_array, C)
    1392              : 
    1393              :                CALL build_dbcsr_from_rows(para_env_sub, mo_coeff_as(ispin), &
    1394              :                                           C(:, MINVAL(orbitals(:, ispin)):MAXVAL(orbitals(:, ispin))), &
    1395           60 :                                           mat_munu%matrix, gd_array, eri_env%eri_gpw%eps_filter)
    1396           12 :                CALL release_group_dist(gd_array)
    1397           36 :                DEALLOCATE (C)
    1398              :             END BLOCK
    1399              : 
    1400           12 :             CALL dbcsr_create(matrix_pq_rnu(ispin), template=mo_coeff_as(ispin))
    1401           12 :             CALL dbcsr_set(matrix_pq_rnu(ispin), 0.0_dp)
    1402              : 
    1403           12 :             CALL dbcsr_get_info(matrix_pq_rnu(ispin), nfullrows_total=nrow_global, nfullcols_total=ncol_global)
    1404              : 
    1405           12 :             NULLIFY (fm_struct)
    1406              :             CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env_sub, &
    1407           12 :                                      nrow_global=nrow_global, ncol_global=ncol_global)
    1408           12 :             CALL cp_fm_create(fm_matrix_pq_rnu(ispin), fm_struct)
    1409           12 :             CALL cp_fm_create(fm_mo_coeff_as(ispin), fm_struct)
    1410           12 :             CALL cp_fm_struct_release(fm_struct)
    1411              : 
    1412           12 :             CALL copy_dbcsr_to_fm(mo_coeff_as(ispin), fm_mo_coeff_as(ispin))
    1413              : 
    1414           12 :             NULLIFY (fm_struct)
    1415              :             CALL cp_fm_struct_create(fm_struct, context=blacs_env, para_env=para_env_sub, &
    1416           12 :                                      nrow_global=ncol_global, ncol_global=ncol_global)
    1417           12 :             CALL cp_fm_create(fm_matrix_pq_rs(ispin), fm_struct)
    1418           36 :             CALL cp_fm_struct_release(fm_struct)
    1419              :          END DO
    1420              : 
    1421              :          ! Copy the active space of the MOs into DBCSR matrices
    1422              :       END IF
    1423              : 
    1424           62 :       CALL auxbas_pw_pool%create_pw(wfn_r)
    1425           62 :       CALL auxbas_pw_pool%create_pw(rho_g)
    1426              :       CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell, &
    1427           62 :                       particle_set=particle_set, atomic_kind_set=atomic_kind_set)
    1428              : 
    1429           62 :       IF (eri_env%eri_gpw%store_wfn) THEN
    1430              :          ! pre-calculate wavefunctions on reals space grid
    1431           50 :          rsize = 0.0_dp
    1432           50 :          nmo = 0
    1433          112 :          DO ispin = 1, nspins
    1434           62 :             CALL get_mo_set(mo_set=mos(ispin), nmo=nx)
    1435           62 :             nmo = MAX(nmo, nx)
    1436          298 :             rsize = REAL(SIZE(wfn_r%array), KIND=dp)*nx
    1437              :          END DO
    1438           50 :          IF (print1 .AND. iw > 0) THEN
    1439            1 :             rsize = rsize*8._dp/1000000._dp
    1440            1 :             WRITE (iw, "(T4,'ERI_GPW|',' Store active orbitals on real space grid ',T63,F12.3,' MB')") rsize
    1441              :          END IF
    1442          480 :          ALLOCATE (wfn_a(nmo, nspins))
    1443          112 :          DO ispin = 1, nspins
    1444           62 :             CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nmo=nmo)
    1445          276 :             DO i1 = 1, SIZE(orbitals, 1)
    1446          164 :                iwfn = orbitals(i1, ispin)
    1447          164 :                CALL auxbas_pw_pool%create_pw(wfn_a(iwfn, ispin))
    1448              :                CALL calculate_wavefunction(mo_coeff, iwfn, wfn_a(iwfn, ispin), rho_g, atomic_kind_set, &
    1449          164 :                                            qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
    1450          226 :                IF (print2 .AND. iw > 0) THEN
    1451            0 :                   WRITE (iw, "(T4,'ERI_GPW|',' Orbital stored ',I4,'  Spin ',i1)") iwfn, ispin
    1452              :                END IF
    1453              :             END DO
    1454              :          END DO
    1455              :       ELSE
    1456              :          ! Even if we do not store all WFNs, we still need containers for the functions to store
    1457           12 :          ALLOCATE (wfn1, wfn2)
    1458           12 :          CALL auxbas_pw_pool%create_pw(wfn1)
    1459           12 :          CALL auxbas_pw_pool%create_pw(wfn2)
    1460           12 :          IF (eri_env%method /= eri_method_gpw_ht) THEN
    1461            6 :             ALLOCATE (wfn3, wfn4)
    1462            6 :             CALL auxbas_pw_pool%create_pw(wfn3)
    1463            6 :             CALL auxbas_pw_pool%create_pw(wfn4)
    1464              :          END IF
    1465              :       END IF
    1466              : 
    1467              :       ! get some of the grids ready
    1468           62 :       CALL auxbas_pw_pool%create_pw(rho_r)
    1469           62 :       CALL auxbas_pw_pool%create_pw(pot_g)
    1470              : 
    1471              :       ! run the FFT once, to set up buffers and to take into account the memory
    1472           62 :       CALL pw_zero(rho_r)
    1473           62 :       CALL pw_transfer(rho_r, rho_g)
    1474           62 :       dvol = rho_r%pw_grid%dvol
    1475           62 :       CALL mp_group%set_handle(eri_env%eri(1)%csr_mat%mp_group%get_handle())
    1476              : 
    1477              :       ! calculate the integrals
    1478           62 :       stored_integrals = 0
    1479          136 :       DO isp1 = 1, nspins
    1480           74 :          CALL get_mo_set(mo_set=mos(isp1), nmo=nmo1, mo_coeff=mo_coeff1)
    1481           74 :          nmm = (nmo1*(nmo1 + 1))/2
    1482           74 :          IF (eri_env%method == eri_method_gpw_ht) THEN
    1483           36 :             irange = [1, nmm]
    1484              :          ELSE
    1485           62 :             irange = get_irange_csr(nmm, para_env)
    1486              :          END IF
    1487           74 :          irange_sub = get_limit(nmm, number_of_subgroups, color)
    1488          398 :          DO i1 = 1, SIZE(orbitals, 1)
    1489          188 :             iwa1 = orbitals(i1, isp1)
    1490          188 :             IF (eri_env%eri_gpw%store_wfn) THEN
    1491          164 :                wfn1 => wfn_a(iwa1, isp1)
    1492              :             ELSE
    1493              :                CALL calculate_wavefunction(mo_coeff1, iwa1, wfn1, rho_g, atomic_kind_set, &
    1494           24 :                                            qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
    1495              :             END IF
    1496          622 :             DO i2 = i1, SIZE(orbitals, 1)
    1497          360 :                iwa2 = orbitals(i2, isp1)
    1498          360 :                iwa12 = csr_idx_to_combined(iwa1, iwa2, nmo1)
    1499              :                ! Skip calculation directly if the pair is not part of our subgroup
    1500          360 :                IF (iwa12 < irange_sub(1) .OR. iwa12 > irange_sub(2)) CYCLE
    1501          360 :                IF (iwa12 >= irange(1) .AND. iwa12 <= irange(2)) THEN
    1502          198 :                   iwa12 = iwa12 - irange(1) + 1
    1503              :                ELSE
    1504          162 :                   iwa12 = 0
    1505              :                END IF
    1506          360 :                IF (eri_env%eri_gpw%store_wfn) THEN
    1507          324 :                   wfn2 => wfn_a(iwa2, isp1)
    1508              :                ELSE
    1509              :                   CALL calculate_wavefunction(mo_coeff1, iwa2, wfn2, rho_g, atomic_kind_set, &
    1510           36 :                                               qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
    1511              :                END IF
    1512              :                ! calculate charge distribution and potential
    1513          360 :                CALL pw_zero(rho_r)
    1514          360 :                CALL pw_multiply(rho_r, wfn1, wfn2)
    1515          360 :                IF (print2) THEN
    1516            0 :                   erint = pw_integrate_function(rho_r)/dvol
    1517            0 :                   IF (iw > 0) THEN
    1518              :                      WRITE (iw, "(T4,'ERI_GPW| Integral rho_ab ',T32,2I4,' [',I1,']',T58,G20.14)") &
    1519            0 :                         iwa1, iwa2, isp1, erint
    1520              :                   END IF
    1521              :                END IF
    1522          360 :                CALL pw_transfer(rho_r, rho_g)
    1523          360 :                CALL pw_poisson_solve(poisson_env, rho_g, pair_int, pot_g)
    1524              :                ! screening using pair_int
    1525          360 :                IF (pair_int < eri_env%eps_integral) CYCLE
    1526          360 :                CALL pw_transfer(pot_g, rho_r)
    1527              :                !
    1528          908 :                IF (eri_env%method == eri_method_gpw_ht) THEN
    1529           36 :                   CALL pw_scale(rho_r, dvol)
    1530           72 :                   DO isp2 = isp1, nspins
    1531           36 :                      CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2)
    1532           36 :                      nx = (nmo2*(nmo2 + 1))/2
    1533          180 :                      ALLOCATE (eri(nx), eri_index(nx))
    1534           36 :                      CALL dbcsr_set(mat_munu%matrix, 0.0_dp)
    1535              :                      CALL integrate_v_rspace(rho_r, hmat=mat_munu, qs_env=qs_env, &
    1536              :                                              calculate_forces=.FALSE., compute_tau=.FALSE., gapw=.FALSE., &
    1537           36 :                                              pw_env_external=pw_env_sub, task_list_external=task_list_sub)
    1538              : 
    1539              :                      CALL dbcsr_multiply("N", "N", 1.0_dp, mat_munu%matrix, mo_coeff_as(isp2), &
    1540           36 :                                          0.0_dp, matrix_pq_rnu(isp2), filter_eps=eri_env%eri_gpw%eps_filter)
    1541           36 :                      CALL copy_dbcsr_to_fm(matrix_pq_rnu(isp2), fm_matrix_pq_rnu(isp2))
    1542              : 
    1543           36 :                      CALL cp_fm_get_info(fm_matrix_pq_rnu(isp2), ncol_global=ncol_global, nrow_global=nrow_global)
    1544              : 
    1545              :                      CALL parallel_gemm("T", "N", ncol_global, ncol_global, nrow_global, 0.5_dp, &
    1546              :                                         fm_matrix_pq_rnu(isp2), fm_mo_coeff_as(isp2), &
    1547           36 :                                         0.0_dp, fm_matrix_pq_rs(isp2))
    1548              :                      CALL parallel_gemm("T", "N", ncol_global, ncol_global, nrow_global, 0.5_dp, &
    1549              :                                         fm_mo_coeff_as(isp2), fm_matrix_pq_rnu(isp2), &
    1550           36 :                                         1.0_dp, fm_matrix_pq_rs(isp2))
    1551              : 
    1552              :                      CALL cp_fm_get_info(fm_matrix_pq_rs(isp2), ncol_local=ncol_local, nrow_local=nrow_local, &
    1553           36 :                                          col_indices=col_indices, row_indices=row_indices)
    1554              : 
    1555           36 :                      icount2 = 0
    1556          108 :                      DO col_local = 1, ncol_local
    1557           72 :                         iwb2 = orbitals(col_indices(col_local), isp2)
    1558          180 :                         DO row_local = 1, nrow_local
    1559           72 :                            iwb1 = orbitals(row_indices(row_local), isp2)
    1560              : 
    1561          144 :                            IF (iwb1 <= iwb2) THEN
    1562           54 :                               iwb12 = csr_idx_to_combined(iwb1, iwb2, nmo2)
    1563           54 :                               erint = fm_matrix_pq_rs(isp2)%local_data(row_local, col_local)
    1564           54 :                               IF (ABS(erint) > eri_env%eps_integral .AND. (iwa12 <= iwb12 .OR. isp1 /= isp2)) THEN
    1565           30 :                                  icount2 = icount2 + 1
    1566           30 :                                  eri(icount2) = erint
    1567           30 :                                  eri_index(icount2) = iwb12
    1568              :                               END IF
    1569              :                            END IF
    1570              :                         END DO
    1571              :                      END DO
    1572           36 :                      stored_integrals = stored_integrals + icount2
    1573              :                      !
    1574           36 :                      isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
    1575           36 :                      CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
    1576              :                      !
    1577          180 :                      DEALLOCATE (eri, eri_index)
    1578              :                   END DO
    1579          324 :                ELSEIF (eri_env%method == eri_method_full_gpw) THEN
    1580          702 :                   DO isp2 = isp1, nspins
    1581          378 :                      CALL get_mo_set(mo_set=mos(isp2), nmo=nmo2, mo_coeff=mo_coeff2)
    1582          378 :                      nx = (nmo2*(nmo2 + 1))/2
    1583         1890 :                      ALLOCATE (eri(nx), eri_index(nx))
    1584          378 :                      icount2 = 0
    1585          378 :                      iwbs = 1
    1586          378 :                      IF (isp1 == isp2) iwbs = i1
    1587          378 :                      isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
    1588         1334 :                      DO i3 = iwbs, SIZE(orbitals, 1)
    1589          956 :                         iwb1 = orbitals(i3, isp2)
    1590          956 :                         IF (eri_env%eri_gpw%store_wfn) THEN
    1591          926 :                            wfn3 => wfn_a(iwb1, isp2)
    1592              :                         ELSE
    1593              :                            CALL calculate_wavefunction(mo_coeff2, iwb1, wfn3, rho_g, atomic_kind_set, &
    1594           30 :                                                        qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
    1595              :                         END IF
    1596          956 :                         CALL pw_zero(wfn_r)
    1597          956 :                         CALL pw_multiply(wfn_r, rho_r, wfn3)
    1598          956 :                         iwbt = i3
    1599          956 :                         IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
    1600         3020 :                         DO i4 = iwbt, SIZE(orbitals, 1)
    1601         1686 :                            iwb2 = orbitals(i4, isp2)
    1602         1686 :                            IF (eri_env%eri_gpw%store_wfn) THEN
    1603         1650 :                               wfn4 => wfn_a(iwb2, isp2)
    1604              :                            ELSE
    1605              :                               CALL calculate_wavefunction(mo_coeff2, iwb2, wfn4, rho_g, atomic_kind_set, &
    1606           36 :                                                           qs_kind_set, cell, dft_control, particle_set, pw_env_sub)
    1607              :                            END IF
    1608              :                            ! We reduce the amount of communication by collecting the local sums first and sum globally later
    1609         1686 :                            erint = pw_integral_ab(wfn_r, wfn4, local_only=.TRUE.)
    1610         1686 :                            icount2 = icount2 + 1
    1611         1686 :                            eri(icount2) = erint
    1612         2642 :                            eri_index(icount2) = csr_idx_to_combined(iwb1, iwb2, nmo2)
    1613              :                         END DO
    1614              :                      END DO
    1615              :                      ! Now, we sum the integrals globally
    1616          378 :                      CALL wfn_r%pw_grid%para%group%sum(eri)
    1617              :                      ! and we reorder the integrals to prevent storing too small integrals
    1618          378 :                      intcount = 0
    1619          378 :                      icount2 = 0
    1620              :                      iwbs = 1
    1621              :                      IF (isp1 == isp2) iwbs = i1
    1622          378 :                      isp = (isp1 - 1)*isp2 + (isp2 - isp1 + 1)
    1623         1334 :                      DO i3 = iwbs, SIZE(orbitals, 1)
    1624          956 :                         iwb1 = orbitals(i3, isp2)
    1625          956 :                         iwbt = i3
    1626          956 :                         IF (isp1 == isp2 .AND. i1 == i3) iwbt = i2
    1627         3020 :                         DO i4 = iwbt, SIZE(orbitals, 1)
    1628         1686 :                            iwb2 = orbitals(i4, isp2)
    1629         1686 :                            intcount = intcount + 1
    1630         1686 :                            erint = eri(intcount)
    1631         2642 :                            IF (ABS(erint) > eri_env%eps_integral) THEN
    1632         1360 :                               icount2 = icount2 + 1
    1633         1360 :                               eri(icount2) = erint
    1634         1360 :                               eri_index(icount2) = eri_index(intcount)
    1635              :                            END IF
    1636              :                         END DO
    1637              :                      END DO
    1638          378 :                      stored_integrals = stored_integrals + icount2
    1639              :                      !
    1640          378 :                      CALL update_csr_matrix(eri_env%eri(isp)%csr_mat, icount2, eri, eri_index, iwa12)
    1641              :                      !
    1642         1080 :                      DEALLOCATE (eri, eri_index)
    1643              :                   END DO
    1644              :                ELSE
    1645            0 :                   CPABORT("Unknown option")
    1646              :                END IF
    1647              :             END DO
    1648              :          END DO
    1649              :       END DO
    1650              : 
    1651           62 :       IF (print1 .AND. iw > 0) THEN
    1652            1 :          WRITE (iw, "(T4,'ERI_GPW|',' Number of Integrals stored ',T68,I10)") stored_integrals
    1653              :       END IF
    1654              : 
    1655           62 :       IF (eri_env%eri_gpw%store_wfn) THEN
    1656          112 :          DO ispin = 1, nspins
    1657          276 :             DO i1 = 1, SIZE(orbitals, 1)
    1658          164 :                iwfn = orbitals(i1, ispin)
    1659          226 :                CALL wfn_a(iwfn, ispin)%release()
    1660              :             END DO
    1661              :          END DO
    1662           50 :          DEALLOCATE (wfn_a)
    1663              :       ELSE
    1664           12 :          CALL wfn1%release()
    1665           12 :          CALL wfn2%release()
    1666           12 :          DEALLOCATE (wfn1, wfn2)
    1667           12 :          IF (eri_env%method /= eri_method_gpw_ht) THEN
    1668            6 :             CALL wfn3%release()
    1669            6 :             CALL wfn4%release()
    1670            6 :             DEALLOCATE (wfn3, wfn4)
    1671              :          END IF
    1672              :       END IF
    1673           62 :       CALL auxbas_pw_pool%give_back_pw(wfn_r)
    1674           62 :       CALL auxbas_pw_pool%give_back_pw(rho_g)
    1675           62 :       CALL auxbas_pw_pool%give_back_pw(rho_r)
    1676           62 :       CALL auxbas_pw_pool%give_back_pw(pot_g)
    1677           62 :       CALL mp_para_env_release(para_env_sub)
    1678              : 
    1679           62 :       IF (eri_env%method == eri_method_gpw_ht) THEN
    1680           24 :          DO ispin = 1, nspins
    1681           12 :             CALL dbcsr_release(mo_coeff_as(ispin))
    1682           12 :             CALL dbcsr_release(matrix_pq_rnu(ispin))
    1683           12 :             CALL cp_fm_release(fm_mo_coeff_as(ispin))
    1684           12 :             CALL cp_fm_release(fm_matrix_pq_rnu(ispin))
    1685           24 :             CALL cp_fm_release(fm_matrix_pq_rs(ispin))
    1686              :          END DO
    1687            0 :          DEALLOCATE (mo_coeff_as, matrix_pq_rnu, &
    1688           12 :                      fm_mo_coeff_as, fm_matrix_pq_rnu, fm_matrix_pq_rs)
    1689           12 :          CALL deallocate_task_list(task_list_sub)
    1690           12 :          CALL dbcsr_release(mat_munu%matrix)
    1691           12 :          DEALLOCATE (mat_munu%matrix)
    1692           12 :          CALL release_neighbor_list_sets(sab_orb_sub)
    1693           12 :          CALL cp_blacs_env_release(blacs_env_sub)
    1694              :       END IF
    1695           62 :       CALL pw_env_release(pw_env_sub)
    1696              :       ! Return to the old qs_control
    1697           62 :       dft_control%qs_control => qs_control_old
    1698           62 :       DEALLOCATE (qs_control%e_cutoff)
    1699           62 :       DEALLOCATE (qs_control)
    1700              : 
    1701           62 :       CALL timestop(handle)
    1702              : 
    1703          124 :    END SUBROUTINE calculate_eri_gpw
    1704              : 
    1705              : ! **************************************************************************************************
    1706              : !> \brief Sets the Green's function
    1707              : !> \param green ...
    1708              : !> \param eri_env ...
    1709              : !> \par History
    1710              : !>      04.2016 created [JGH]
    1711              : ! **************************************************************************************************
    1712           62 :    SUBROUTINE pw_eri_green_create(green, eri_env)
    1713              : 
    1714              :       TYPE(greens_fn_type), INTENT(INOUT)                :: green
    1715              :       TYPE(eri_type)                                     :: eri_env
    1716              : 
    1717              :       INTEGER                                            :: ig
    1718              :       REAL(KIND=dp)                                      :: a, ea, g2, g3d, ga, gg, rg, rlength
    1719              : 
    1720              :       ! initialize influence function
    1721              :       ASSOCIATE (gf => green%influence_fn, grid => green%influence_fn%pw_grid)
    1722           64 :          SELECT CASE (green%method)
    1723              :          CASE (PERIODIC3D)
    1724              : 
    1725           64 :             SELECT CASE (eri_env%operator)
    1726              :             CASE (eri_operator_coulomb)
    1727       262145 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1728       262143 :                   g2 = grid%gsq(ig)
    1729       262145 :                   gf%array(ig) = fourpi/g2
    1730              :                END DO
    1731            2 :                IF (grid%have_g0) gf%array(1) = 0.0_dp
    1732              :             CASE (eri_operator_yukawa)
    1733            0 :                a = eri_env%operator_parameter**2
    1734            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1735            0 :                   g2 = grid%gsq(ig)
    1736            0 :                   gf%array(ig) = fourpi/(a + g2)
    1737              :                END DO
    1738            0 :                IF (grid%have_g0) gf%array(1) = fourpi/a
    1739              :             CASE (eri_operator_erf)
    1740            0 :                a = eri_env%operator_parameter**2
    1741            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1742            0 :                   g2 = grid%gsq(ig)
    1743            0 :                   ga = -0.25_dp*g2/a
    1744            0 :                   gf%array(ig) = fourpi/g2*EXP(ga)
    1745              :                END DO
    1746            0 :                IF (grid%have_g0) gf%array(1) = 0.0_dp
    1747              :             CASE (eri_operator_erfc)
    1748            0 :                a = eri_env%operator_parameter**2
    1749            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1750            0 :                   g2 = grid%gsq(ig)
    1751            0 :                   ga = -0.25_dp*g2/a
    1752            0 :                   gf%array(ig) = fourpi/g2*(1._dp - EXP(ga))
    1753              :                END DO
    1754            0 :                IF (grid%have_g0) gf%array(1) = 0.25_dp*fourpi/a
    1755              :             CASE (eri_operator_trunc)
    1756            0 :                a = eri_env%operator_parameter
    1757            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1758            0 :                   g2 = grid%gsq(ig)
    1759            0 :                   ga = SQRT(g2)*a
    1760            0 :                   IF (ga >= 0.005_dp) THEN
    1761            0 :                      gf%array(ig) = fourpi/g2*(1.0_dp - COS(ga))
    1762              :                   ELSE
    1763            0 :                      gf%array(ig) = fourpi/g2*ga**2/2.0_dp*(1.0_dp - ga**2/12.0_dp)
    1764              :                   END IF
    1765              :                END DO
    1766            0 :                IF (grid%have_g0) gf%array(1) = 0.0_dp
    1767              :             CASE (eri_operator_gaussian)
    1768            0 :                CPABORT("")
    1769              :             CASE DEFAULT
    1770            2 :                CPABORT("")
    1771              :             END SELECT
    1772              : 
    1773              :          CASE (ANALYTIC0D)
    1774              : 
    1775          106 :             SELECT CASE (eri_env%operator)
    1776              :             CASE (eri_operator_coulomb)
    1777           46 :                rlength = green%radius
    1778     12951437 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1779     12951391 :                   g2 = grid%gsq(ig)
    1780     12951391 :                   gg = SQRT(g2)
    1781     12951391 :                   g3d = fourpi/g2
    1782     12951437 :                   gf%array(ig) = g3d*(1.0_dp - COS(rlength*gg))
    1783              :                END DO
    1784           46 :                IF (grid%have_g0) gf%array(1) = 0.5_dp*fourpi*rlength*rlength
    1785              :             CASE (eri_operator_yukawa)
    1786            0 :                rlength = green%radius
    1787            0 :                a = eri_env%operator_parameter
    1788            0 :                ea = EXP(-a*rlength)
    1789            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1790            0 :                   g2 = grid%gsq(ig)
    1791            0 :                   gg = SQRT(g2)
    1792            0 :                   g3d = fourpi/(a*a + g2)
    1793            0 :                   rg = rlength*gg
    1794            0 :                   gf%array(ig) = g3d*(1.0_dp - ea*(COS(rg) + a/gg*SIN(rg)))
    1795              :                END DO
    1796            0 :                IF (grid%have_g0) gf%array(1) = fourpi/(a*a)*(1.0_dp - ea*(1._dp + a*rlength))
    1797              :             CASE (eri_operator_erf)
    1798           14 :                rlength = green%radius
    1799           14 :                a = eri_env%operator_parameter**2
    1800      1512007 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1801      1511993 :                   g2 = grid%gsq(ig)
    1802      1511993 :                   gg = SQRT(g2)
    1803      1511993 :                   ga = -0.25_dp*g2/a
    1804      1512007 :                   gf%array(ig) = fourpi/g2*EXP(ga)*(1.0_dp - COS(rlength*gg))
    1805              :                END DO
    1806           14 :                IF (grid%have_g0) gf%array(1) = 0.5_dp*fourpi*rlength*rlength
    1807              :             CASE (eri_operator_erfc)
    1808            0 :                rlength = green%radius
    1809            0 :                a = eri_env%operator_parameter**2
    1810            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1811            0 :                   g2 = grid%gsq(ig)
    1812            0 :                   gg = SQRT(g2)
    1813            0 :                   ga = -0.25_dp*g2/a
    1814            0 :                   gf%array(ig) = fourpi/g2*(1._dp - EXP(ga))*(1.0_dp - COS(rlength*gg))
    1815              :                END DO
    1816            0 :                IF (grid%have_g0) gf%array(1) = 0._dp
    1817              :             CASE (eri_operator_trunc)
    1818            0 :                a = eri_env%operator_parameter
    1819            0 :                DO ig = grid%first_gne0, grid%ngpts_cut_local
    1820            0 :                   g2 = grid%gsq(ig)
    1821            0 :                   ga = SQRT(g2)*a
    1822            0 :                   IF (ga >= 0.005_dp) THEN
    1823            0 :                      gf%array(ig) = fourpi/g2*(1.0_dp - COS(ga))
    1824              :                   ELSE
    1825            0 :                      gf%array(ig) = fourpi/g2*ga**2/2.0_dp*(1.0_dp - ga**2/12.0_dp)
    1826              :                   END IF
    1827              :                END DO
    1828            0 :                IF (grid%have_g0) gf%array(1) = 0.0_dp
    1829              :             CASE (eri_operator_gaussian)
    1830            0 :                CPABORT("")
    1831              :             CASE DEFAULT
    1832           60 :                CPABORT("")
    1833              :             END SELECT
    1834              : 
    1835              :          CASE DEFAULT
    1836           62 :             CPABORT("")
    1837              :          END SELECT
    1838              :       END ASSOCIATE
    1839              : 
    1840           62 :    END SUBROUTINE pw_eri_green_create
    1841              : 
    1842              : ! **************************************************************************************************
    1843              : !> \brief Adds data for a new row to the csr matrix
    1844              : !> \param csr_mat ...
    1845              : !> \param nnz ...
    1846              : !> \param rdat ...
    1847              : !> \param rind ...
    1848              : !> \param irow ...
    1849              : !> \par History
    1850              : !>      04.2016 created [JGH]
    1851              : ! **************************************************************************************************
    1852          414 :    SUBROUTINE update_csr_matrix(csr_mat, nnz, rdat, rind, irow)
    1853              : 
    1854              :       TYPE(dbcsr_csr_type), INTENT(INOUT)                :: csr_mat
    1855              :       INTEGER, INTENT(IN)                                :: nnz
    1856              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: rdat
    1857              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: rind
    1858              :       INTEGER, INTENT(IN)                                :: irow
    1859              : 
    1860              :       INTEGER                                            :: k, nrow, nze, nze_new
    1861              : 
    1862          414 :       IF (irow /= 0) THEN
    1863          225 :          nze = csr_mat%nze_local
    1864          225 :          nze_new = nze + nnz
    1865              :          ! values
    1866          225 :          CALL reallocate(csr_mat%nzval_local%r_dp, 1, nze_new)
    1867          935 :          csr_mat%nzval_local%r_dp(nze + 1:nze_new) = rdat(1:nnz)
    1868              :          ! col indices
    1869          225 :          CALL reallocate(csr_mat%colind_local, 1, nze_new)
    1870          935 :          csr_mat%colind_local(nze + 1:nze_new) = rind(1:nnz)
    1871              :          ! rows
    1872          225 :          nrow = csr_mat%nrows_local
    1873          225 :          CALL reallocate(csr_mat%rowptr_local, 1, irow + 1)
    1874          508 :          csr_mat%rowptr_local(nrow + 1:irow) = nze + 1
    1875          225 :          csr_mat%rowptr_local(irow + 1) = nze_new + 1
    1876              :          ! nzerow
    1877          225 :          CALL reallocate(csr_mat%nzerow_local, 1, irow)
    1878          508 :          DO k = nrow + 1, irow
    1879          508 :             csr_mat%nzerow_local(k) = csr_mat%rowptr_local(k + 1) - csr_mat%rowptr_local(k)
    1880              :          END DO
    1881          225 :          csr_mat%nrows_local = irow
    1882          225 :          csr_mat%nze_local = csr_mat%nze_local + nnz
    1883              :       END IF
    1884          414 :       csr_mat%nze_total = csr_mat%nze_total + nnz
    1885          414 :       csr_mat%has_indices = .TRUE.
    1886              : 
    1887          414 :    END SUBROUTINE update_csr_matrix
    1888              : 
    1889              : ! **************************************************************************************************
    1890              : !> \brief Computes and prints the active orbitals on Cube Files
    1891              : !> \param input ...
    1892              : !> \param qs_env the qs_env in which the qs_env lives
    1893              : !> \param mos ...
    1894              : ! **************************************************************************************************
    1895            2 :    SUBROUTINE print_orbital_cubes(input, qs_env, mos)
    1896              :       TYPE(section_vals_type), POINTER                   :: input
    1897              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1898              :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
    1899              : 
    1900              :       CHARACTER(LEN=default_path_length)                 :: filebody, filename, title
    1901              :       INTEGER                                            :: i, imo, isp, nmo, str(3), unit_nr
    1902            2 :       INTEGER, DIMENSION(:), POINTER                     :: alist, blist, istride
    1903              :       LOGICAL                                            :: do_mo, explicit_a, explicit_b
    1904            2 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    1905              :       TYPE(cell_type), POINTER                           :: cell
    1906              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
    1907              :       TYPE(dft_control_type), POINTER                    :: dft_control
    1908              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1909              :       TYPE(particle_list_type), POINTER                  :: particles
    1910            2 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
    1911              :       TYPE(pw_c1d_gs_type)                               :: wf_g
    1912              :       TYPE(pw_env_type), POINTER                         :: pw_env
    1913              :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
    1914              :       TYPE(pw_r3d_rs_type)                               :: wf_r
    1915            2 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1916              :       TYPE(qs_subsys_type), POINTER                      :: subsys
    1917              :       TYPE(section_vals_type), POINTER                   :: dft_section, scf_input
    1918              : 
    1919            2 :       CALL section_vals_val_get(input, "FILENAME", c_val=filebody)
    1920            2 :       CALL section_vals_val_get(input, "STRIDE", i_vals=istride)
    1921            2 :       IF (SIZE(istride) == 1) THEN
    1922            8 :          str(1:3) = istride(1)
    1923            0 :       ELSEIF (SIZE(istride) == 3) THEN
    1924            0 :          str(1:3) = istride(1:3)
    1925              :       ELSE
    1926            0 :          CPABORT("STRIDE arguments inconsistent")
    1927              :       END IF
    1928            2 :       CALL section_vals_val_get(input, "ALIST", i_vals=alist, explicit=explicit_a)
    1929            2 :       CALL section_vals_val_get(input, "BLIST", i_vals=blist, explicit=explicit_b)
    1930              : 
    1931              :       CALL get_qs_env(qs_env=qs_env, &
    1932              :                       dft_control=dft_control, &
    1933              :                       para_env=para_env, &
    1934              :                       subsys=subsys, &
    1935              :                       atomic_kind_set=atomic_kind_set, &
    1936              :                       qs_kind_set=qs_kind_set, &
    1937              :                       cell=cell, &
    1938              :                       particle_set=particle_set, &
    1939              :                       pw_env=pw_env, &
    1940            2 :                       input=scf_input)
    1941              : 
    1942            2 :       CALL qs_subsys_get(subsys, particles=particles)
    1943              :       !
    1944            2 :       CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
    1945            2 :       CALL auxbas_pw_pool%create_pw(wf_r)
    1946            2 :       CALL auxbas_pw_pool%create_pw(wf_g)
    1947              :       !
    1948            2 :       dft_section => section_vals_get_subs_vals(scf_input, "DFT")
    1949              :       !
    1950            4 :       DO isp = 1, SIZE(mos)
    1951            2 :          CALL get_mo_set(mo_set=mos(isp), mo_coeff=mo_coeff, nmo=nmo)
    1952              : 
    1953            2 :          IF (SIZE(mos) > 1) THEN
    1954            0 :             SELECT CASE (isp)
    1955              :             CASE (1)
    1956              :                CALL write_mo_set_to_output_unit(mos(isp), atomic_kind_set, qs_kind_set, particle_set, &
    1957            0 :                                                 dft_section, 4, 0, final_mos=.TRUE., spin="ALPHA")
    1958              :             CASE (2)
    1959              :                CALL write_mo_set_to_output_unit(mos(isp), atomic_kind_set, qs_kind_set, particle_set, &
    1960            0 :                                                 dft_section, 4, 0, final_mos=.TRUE., spin="BETA")
    1961              :             CASE DEFAULT
    1962            0 :                CPABORT("Invalid spin")
    1963              :             END SELECT
    1964              :          ELSE
    1965              :             CALL write_mo_set_to_output_unit(mos(isp), atomic_kind_set, qs_kind_set, particle_set, &
    1966            2 :                                              dft_section, 4, 0, final_mos=.TRUE.)
    1967              :          END IF
    1968              : 
    1969           22 :          DO imo = 1, nmo
    1970           16 :             IF (isp == 1 .AND. explicit_a) THEN
    1971           16 :                IF (alist(1) == -1) THEN
    1972              :                   do_mo = .TRUE.
    1973              :                ELSE
    1974           16 :                   do_mo = .FALSE.
    1975           64 :                   DO i = 1, SIZE(alist)
    1976           64 :                      IF (imo == alist(i)) do_mo = .TRUE.
    1977              :                   END DO
    1978              :                END IF
    1979            0 :             ELSE IF (isp == 2 .AND. explicit_b) THEN
    1980            0 :                IF (blist(1) == -1) THEN
    1981              :                   do_mo = .TRUE.
    1982              :                ELSE
    1983            0 :                   do_mo = .FALSE.
    1984            0 :                   DO i = 1, SIZE(blist)
    1985            0 :                      IF (imo == blist(i)) do_mo = .TRUE.
    1986              :                   END DO
    1987              :                END IF
    1988              :             ELSE
    1989              :                do_mo = .TRUE.
    1990              :             END IF
    1991           16 :             IF (.NOT. do_mo) CYCLE
    1992              :             CALL calculate_wavefunction(mo_coeff, imo, wf_r, wf_g, atomic_kind_set, &
    1993            6 :                                         qs_kind_set, cell, dft_control, particle_set, pw_env)
    1994            6 :             IF (para_env%is_source()) THEN
    1995            3 :                WRITE (filename, '(A,A1,I4.4,A1,I1.1,A)') TRIM(filebody), "_", imo, "_", isp, ".cube"
    1996            3 :                CALL open_file(filename, unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE")
    1997            3 :                WRITE (title, *) "Active Orbital ", imo, " spin ", isp
    1998              :             ELSE
    1999            3 :                unit_nr = -1
    2000              :             END IF
    2001            6 :             CALL cp_pw_to_cube(wf_r, unit_nr, title, particles=particles, stride=istride)
    2002            8 :             IF (para_env%is_source()) THEN
    2003           13 :                CALL close_file(unit_nr)
    2004              :             END IF
    2005              :          END DO
    2006              :       END DO
    2007              : 
    2008            2 :       CALL auxbas_pw_pool%give_back_pw(wf_r)
    2009            2 :       CALL auxbas_pw_pool%give_back_pw(wf_g)
    2010              : 
    2011            2 :    END SUBROUTINE print_orbital_cubes
    2012              : 
    2013              : ! **************************************************************************************************
    2014              : !> \brief Writes a FCIDUMP file
    2015              : !> \param active_space_env ...
    2016              : !> \param as_input ...
    2017              : !> \par History
    2018              : !>      04.2016 created [JGH]
    2019              : ! **************************************************************************************************
    2020           62 :    SUBROUTINE fcidump(active_space_env, as_input)
    2021              : 
    2022              :       TYPE(active_space_type), POINTER                   :: active_space_env
    2023              :       TYPE(section_vals_type), POINTER                   :: as_input
    2024              : 
    2025              :       INTEGER                                            :: i, i1, i2, i3, i4, isym, iw, m1, m2, &
    2026              :                                                             nmo, norb, nspins
    2027              :       REAL(KIND=dp)                                      :: checksum, esub
    2028           62 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: fmat
    2029              :       TYPE(cp_logger_type), POINTER                      :: logger
    2030              :       TYPE(eri_fcidump_checksum)                         :: eri_checksum
    2031              : 
    2032           62 :       checksum = 0.0_dp
    2033              : 
    2034          124 :       logger => cp_get_default_logger()
    2035              :       iw = cp_print_key_unit_nr(logger, as_input, "FCIDUMP", &
    2036           62 :                                 extension=".fcidump", file_status="REPLACE", file_action="WRITE", file_form="FORMATTED")
    2037              :       !
    2038           62 :       nspins = active_space_env%nspins
    2039           62 :       norb = SIZE(active_space_env%active_orbitals, 1)
    2040           62 :       IF (nspins == 1) THEN
    2041              :          ASSOCIATE (ms2 => active_space_env%multiplicity, &
    2042              :                     nelec => active_space_env%nelec_active)
    2043              : 
    2044           50 :             IF (iw > 0) THEN
    2045           25 :                WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
    2046           25 :                isym = 1
    2047           89 :                WRITE (iw, "(A,1000(I1,','))") "  ORBSYM=", (isym, i=1, norb)
    2048           25 :                isym = 0
    2049           25 :                WRITE (iw, "(A,I1,A)") "  ISYM=", isym, ","
    2050           25 :                WRITE (iw, "(A)") " /"
    2051              :             END IF
    2052              :             !
    2053              :             ! Print integrals: ERI
    2054              :             CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, &
    2055           50 :                                                   eri_fcidump_print(iw, 1, 1), 1, 1)
    2056           50 :             CALL eri_checksum%set(1, 1)
    2057           50 :             CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1)
    2058              : 
    2059              :             ! Print integrals: Fij
    2060              :             ! replicate Fock matrix
    2061           50 :             nmo = active_space_env%eri%norb
    2062          200 :             ALLOCATE (fmat(nmo, nmo))
    2063           50 :             CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
    2064           50 :             IF (iw > 0) THEN
    2065           25 :                i3 = 0; i4 = 0
    2066           89 :                DO m1 = 1, SIZE(active_space_env%active_orbitals, 1)
    2067           64 :                   i1 = active_space_env%active_orbitals(m1, 1)
    2068          215 :                   DO m2 = m1, SIZE(active_space_env%active_orbitals, 1)
    2069          126 :                      i2 = active_space_env%active_orbitals(m2, 1)
    2070          126 :                      checksum = checksum + ABS(fmat(i1, i2))
    2071          190 :                      WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
    2072              :                   END DO
    2073              :                END DO
    2074              :             END IF
    2075           50 :             DEALLOCATE (fmat)
    2076              :             ! Print energy
    2077           50 :             esub = active_space_env%energy_inactive
    2078           50 :             i1 = 0; i2 = 0; i3 = 0; i4 = 0
    2079           50 :             checksum = checksum + ABS(esub)
    2080          100 :             IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4
    2081              :          END ASSOCIATE
    2082              : 
    2083              :       ELSE
    2084              :          ASSOCIATE (ms2 => active_space_env%multiplicity, &
    2085              :                     nelec => active_space_env%nelec_active)
    2086              : 
    2087           12 :             IF (iw > 0) THEN
    2088            6 :                WRITE (iw, "(A,A,I4,A,I4,A,I2,A)") "&FCI", " NORB=", norb, ",NELEC=", nelec, ",MS2=", ms2, ","
    2089            6 :                isym = 1
    2090           21 :                WRITE (iw, "(A,1000(I1,','))") "  ORBSYM=", (isym, i=1, norb)
    2091            6 :                isym = 0
    2092            6 :                WRITE (iw, "(A,I1,A)") "  ISYM=", isym, ","
    2093            6 :                WRITE (iw, "(A,I1,A)") "  UHF=", 1, ","
    2094            6 :                WRITE (iw, "(A)") " /"
    2095              :             END IF
    2096              :             !
    2097              :             ! Print integrals: ERI
    2098              :             ! alpha-alpha
    2099              :             CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, &
    2100           12 :                                                   eri_fcidump_print(iw, 1, 1), 1, 1)
    2101           12 :             CALL eri_checksum%set(1, 1)
    2102           12 :             CALL active_space_env%eri%eri_foreach(1, active_space_env%active_orbitals, eri_checksum, 1, 1)
    2103              :             ! alpha-beta
    2104              :             CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, &
    2105           12 :                                                   eri_fcidump_print(iw, 1, norb + 1), 1, 2)
    2106           12 :             CALL eri_checksum%set(1, norb + 1)
    2107           12 :             CALL active_space_env%eri%eri_foreach(2, active_space_env%active_orbitals, eri_checksum, 1, 2)
    2108              :             ! beta-beta
    2109              :             CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, &
    2110           12 :                                                   eri_fcidump_print(iw, norb + 1, norb + 1), 2, 2)
    2111           12 :             CALL eri_checksum%set(norb + 1, norb + 1)
    2112           12 :             CALL active_space_env%eri%eri_foreach(3, active_space_env%active_orbitals, eri_checksum, 2, 2)
    2113              :             ! Print integrals: Fij
    2114              :             ! alpha
    2115           12 :             nmo = active_space_env%eri%norb
    2116           48 :             ALLOCATE (fmat(nmo, nmo))
    2117           12 :             CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(1), fmat)
    2118           12 :             IF (iw > 0) THEN
    2119            6 :                i3 = 0; i4 = 0
    2120           21 :                DO m1 = 1, norb
    2121           15 :                   i1 = active_space_env%active_orbitals(m1, 1)
    2122           48 :                   DO m2 = m1, norb
    2123           27 :                      i2 = active_space_env%active_orbitals(m2, 1)
    2124           27 :                      checksum = checksum + ABS(fmat(i1, i2))
    2125           42 :                      WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1, m2, i3, i4
    2126              :                   END DO
    2127              :                END DO
    2128              :             END IF
    2129           12 :             DEALLOCATE (fmat)
    2130              :             ! beta
    2131           48 :             ALLOCATE (fmat(nmo, nmo))
    2132           12 :             CALL replicate_and_symmetrize_matrix(nmo, active_space_env%fock_sub(2), fmat)
    2133           12 :             IF (iw > 0) THEN
    2134            6 :                i3 = 0; i4 = 0
    2135           21 :                DO m1 = 1, SIZE(active_space_env%active_orbitals, 1)
    2136           15 :                   i1 = active_space_env%active_orbitals(m1, 2)
    2137           48 :                   DO m2 = m1, SIZE(active_space_env%active_orbitals, 1)
    2138           27 :                      i2 = active_space_env%active_orbitals(m2, 2)
    2139           27 :                      checksum = checksum + ABS(fmat(i1, i2))
    2140           42 :                      WRITE (iw, "(ES23.16,4I4)") fmat(i1, i2), m1 + norb, m2 + norb, i3, i4
    2141              :                   END DO
    2142              :                END DO
    2143              :             END IF
    2144           12 :             DEALLOCATE (fmat)
    2145              :             ! Print energy
    2146           12 :             esub = active_space_env%energy_inactive
    2147           12 :             i1 = 0; i2 = 0; i3 = 0; i4 = 0
    2148           12 :             checksum = checksum + ABS(esub)
    2149           24 :             IF (iw > 0) WRITE (iw, "(ES23.16,4I4)") esub, i1, i2, i3, i4
    2150              :          END ASSOCIATE
    2151              :       END IF
    2152              :       !
    2153           62 :       CALL cp_print_key_finished_output(iw, logger, as_input, "FCIDUMP")
    2154              : 
    2155              :       !>>
    2156           62 :       iw = cp_logger_get_default_io_unit(logger)
    2157           62 :       IF (iw > 0) WRITE (iw, '(T4,A,T66,F12.8)') "FCIDUMP| Checksum:", eri_checksum%checksum + checksum
    2158              :       !<<
    2159              : 
    2160          124 :    END SUBROUTINE fcidump
    2161              : 
    2162              : ! **************************************************************************************************
    2163              : !> \brief replicate and symmetrize a matrix
    2164              : !> \param norb the number of orbitals
    2165              : !> \param distributed_matrix ...
    2166              : !> \param replicated_matrix ...
    2167              : ! **************************************************************************************************
    2168          222 :    SUBROUTINE replicate_and_symmetrize_matrix(norb, distributed_matrix, replicated_matrix)
    2169              :       INTEGER, INTENT(IN)                                :: norb
    2170              :       TYPE(cp_fm_type), INTENT(IN)                       :: distributed_matrix
    2171              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: replicated_matrix
    2172              : 
    2173              :       INTEGER                                            :: i1, i2
    2174              :       REAL(dp)                                           :: mval
    2175              : 
    2176         3834 :       replicated_matrix(:, :) = 0.0_dp
    2177          948 :       DO i1 = 1, norb
    2178         2754 :          DO i2 = i1, norb
    2179         1806 :             CALL cp_fm_get_element(distributed_matrix, i1, i2, mval)
    2180         1806 :             replicated_matrix(i1, i2) = mval
    2181         2532 :             replicated_matrix(i2, i1) = mval
    2182              :          END DO
    2183              :       END DO
    2184          222 :    END SUBROUTINE replicate_and_symmetrize_matrix
    2185              : 
    2186              : ! **************************************************************************************************
    2187              : !> \brief Calculates active space Fock matrix and inactive energy
    2188              : !> \param active_space_env ...
    2189              : !> \par History
    2190              : !>      06.2016 created [JGH]
    2191              : ! **************************************************************************************************
    2192           62 :    SUBROUTINE subspace_fock_matrix(active_space_env)
    2193              : 
    2194              :       TYPE(active_space_type), POINTER                   :: active_space_env
    2195              : 
    2196              :       INTEGER                                            :: i1, i2, is, norb, nspins
    2197              :       REAL(KIND=dp)                                      :: eeri, eref, esub, mval
    2198           62 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: ks_a_mat, ks_a_ref, ks_b_mat, ks_b_ref, &
    2199           62 :                                                             ks_mat, ks_ref, p_a_mat, p_b_mat, p_mat
    2200              :       TYPE(cp_fm_type), POINTER                          :: matrix, mo_coef
    2201              :       TYPE(dbcsr_csr_type), POINTER                      :: eri, eri_aa, eri_ab, eri_bb
    2202              : 
    2203           62 :       eref = active_space_env%energy_ref
    2204           62 :       nspins = active_space_env%nspins
    2205              : 
    2206           62 :       IF (nspins == 1) THEN
    2207           50 :          CALL get_mo_set(active_space_env%mos_active(1), nmo=norb, mo_coeff=mo_coef)
    2208              :          !
    2209              :          ! Loop over ERI, calculate subspace HF energy and Fock matrix
    2210              :          !
    2211              :          ! replicate KS, Core, and P matrices
    2212          500 :          ALLOCATE (ks_mat(norb, norb), ks_ref(norb, norb), p_mat(norb, norb))
    2213          726 :          ks_ref = 0.0_dp
    2214              : 
    2215              :          ! ks_mat contains the KS/Fock matrix (of full density) projected onto the AS MO subspace (f_ref in eq. 19)
    2216           50 :          CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_mat)
    2217           50 :          CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_mat)
    2218              : 
    2219              :          ! compute ks_ref = V_H[rho^A] + V_HFX[rho^A]
    2220           50 :          eri => active_space_env%eri%eri(1)%csr_mat
    2221           50 :          CALL build_subspace_fock_matrix(active_space_env%active_orbitals, eri, p_mat, ks_ref, active_space_env%eri%method)
    2222              : 
    2223              :          ! compute eeri = E_H[rho^A] + E_HFX[rho^A] as
    2224              :          ! eeri = 1/2 * (SUM_pq (V_H[rho^A] + V_HFX[rho^A])_pq * D^A_pq)
    2225          726 :          eeri = 0.5_dp*SUM(ks_ref*p_mat)
    2226              : 
    2227              :          ! now calculate the inactive energy acoording to eq. 19, that is
    2228              :          ! esub = E^I = E_ref - f_ref .* D^A + E_H[rho^A] + E_HFX[rho^A]
    2229              :          ! where f^ref = ks_mat, which is the KS/Fock matrix in MO basis, transformed previously
    2230              :          ! and is equal to ks_mat = h^0 + V_core + V_H[rho] + V_HFX[rho]
    2231          726 :          esub = eref - SUM(ks_mat(1:norb, 1:norb)*p_mat(1:norb, 1:norb)) + eeri
    2232              : 
    2233              :          ! reuse ks_mat to store f^I = f^ref - (V_H[rho^A] + V_HFX[rho^A]) according to eq. 20
    2234          726 :          ks_mat(1:norb, 1:norb) = ks_mat(1:norb, 1:norb) - ks_ref(1:norb, 1:norb)
    2235              :          ! this is now the embedding potential for the AS calculation!
    2236              : 
    2237           50 :          active_space_env%energy_inactive = esub
    2238              : 
    2239           50 :          CALL cp_fm_release(active_space_env%fock_sub)
    2240          200 :          ALLOCATE (active_space_env%fock_sub(nspins))
    2241          100 :          DO is = 1, nspins
    2242           50 :             matrix => active_space_env%ks_sub(is)
    2243              :             CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
    2244          100 :                               name="Active Fock operator")
    2245              :          END DO
    2246           50 :          matrix => active_space_env%fock_sub(1)
    2247          246 :          DO i1 = 1, norb
    2248          726 :             DO i2 = 1, norb
    2249          530 :                mval = ks_mat(i1, i2)
    2250          676 :                CALL cp_fm_set_element(matrix, i1, i2, mval)
    2251              :             END DO
    2252              :          END DO
    2253              :       ELSE
    2254              : 
    2255           12 :          CALL get_mo_set(active_space_env%mos_active(1), nmo=norb)
    2256              :          !
    2257              :          ! Loop over ERI, calculate subspace HF energy and Fock matrix
    2258              :          !
    2259              :          ! replicate KS, Core, and P matrices
    2260              :          ALLOCATE (ks_a_mat(norb, norb), ks_b_mat(norb, norb), &
    2261              :               &    ks_a_ref(norb, norb), ks_b_ref(norb, norb), &
    2262          228 :               &     p_a_mat(norb, norb), p_b_mat(norb, norb))
    2263          552 :          ks_a_ref(:, :) = 0.0_dp; ks_b_ref(:, :) = 0.0_dp
    2264              : 
    2265           12 :          CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(1), p_a_mat)
    2266           12 :          CALL replicate_and_symmetrize_matrix(norb, active_space_env%p_active(2), p_b_mat)
    2267           12 :          CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(1), ks_a_mat)
    2268           12 :          CALL replicate_and_symmetrize_matrix(norb, active_space_env%ks_sub(2), ks_b_mat)
    2269              :          !
    2270              :          !
    2271           12 :          eri_aa => active_space_env%eri%eri(1)%csr_mat
    2272           12 :          eri_ab => active_space_env%eri%eri(2)%csr_mat
    2273           12 :          eri_bb => active_space_env%eri%eri(3)%csr_mat
    2274              :          CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_aa, eri_ab, p_a_mat, p_b_mat, ks_a_ref, &
    2275           12 :                                               tr_mixed_eri=.FALSE., eri_method=active_space_env%eri%method)
    2276              :          CALL build_subspace_spin_fock_matrix(active_space_env%active_orbitals, eri_bb, eri_ab, p_b_mat, p_a_mat, ks_b_ref, &
    2277           12 :                                               tr_mixed_eri=.TRUE., eri_method=active_space_env%eri%method)
    2278              :          !
    2279              :          ! calculate energy
    2280           12 :          eeri = 0.0_dp
    2281          552 :          eeri = 0.5_dp*(SUM(ks_a_ref*p_a_mat) + SUM(ks_b_ref*p_b_mat))
    2282          552 :          esub = eref - SUM(ks_a_mat*p_a_mat) - SUM(ks_b_mat*p_b_mat) + eeri
    2283          276 :          ks_a_mat(:, :) = ks_a_mat(:, :) - ks_a_ref(:, :)
    2284          276 :          ks_b_mat(:, :) = ks_b_mat(:, :) - ks_b_ref(:, :)
    2285              :          !
    2286           12 :          active_space_env%energy_inactive = esub
    2287              :          !
    2288           12 :          CALL cp_fm_release(active_space_env%fock_sub)
    2289           60 :          ALLOCATE (active_space_env%fock_sub(nspins))
    2290           36 :          DO is = 1, nspins
    2291           24 :             matrix => active_space_env%ks_sub(is)
    2292              :             CALL cp_fm_create(active_space_env%fock_sub(is), matrix%matrix_struct, &
    2293           36 :                               name="Active Fock operator")
    2294              :          END DO
    2295              : 
    2296           12 :          matrix => active_space_env%fock_sub(1)
    2297           60 :          DO i1 = 1, norb
    2298          276 :             DO i2 = 1, norb
    2299          216 :                mval = ks_a_mat(i1, i2)
    2300          264 :                CALL cp_fm_set_element(matrix, i1, i2, mval)
    2301              :             END DO
    2302              :          END DO
    2303           12 :          matrix => active_space_env%fock_sub(2)
    2304           72 :          DO i1 = 1, norb
    2305          276 :             DO i2 = 1, norb
    2306          216 :                mval = ks_b_mat(i1, i2)
    2307          264 :                CALL cp_fm_set_element(matrix, i1, i2, mval)
    2308              :             END DO
    2309              :          END DO
    2310              : 
    2311              :       END IF
    2312              : 
    2313           62 :    END SUBROUTINE subspace_fock_matrix
    2314              : 
    2315              : ! **************************************************************************************************
    2316              : !> \brief build subspace fockian
    2317              : !> \param active_orbitals the active orbital indices
    2318              : !> \param eri two electon integrals in MO
    2319              : !> \param p_mat density matrix
    2320              : !> \param ks_ref fockian matrix
    2321              : !> \param eri_method ...
    2322              : ! **************************************************************************************************
    2323           50 :    SUBROUTINE build_subspace_fock_matrix(active_orbitals, eri, p_mat, ks_ref, eri_method)
    2324              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: active_orbitals
    2325              :       TYPE(dbcsr_csr_type), INTENT(IN)                   :: eri
    2326              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: p_mat
    2327              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: ks_ref
    2328              :       INTEGER, INTENT(IN)                                :: eri_method
    2329              : 
    2330              :       INTEGER                                            :: i1, i12, i12l, i2, i3, i34, i34l, i4, &
    2331              :                                                             irptr, m1, m2, nindex, nmo_total, norb
    2332              :       INTEGER, DIMENSION(2)                              :: irange
    2333              :       REAL(dp)                                           :: erint
    2334              :       TYPE(mp_comm_type)                                 :: mp_group
    2335              : 
    2336              :       ! Nothing to do
    2337           50 :       norb = SIZE(active_orbitals, 1)
    2338           50 :       nmo_total = SIZE(p_mat, 1)
    2339           50 :       nindex = (nmo_total*(nmo_total + 1))/2
    2340           50 :       CALL mp_group%set_handle(eri%mp_group%get_handle())
    2341           50 :       IF (eri_method == eri_method_gpw_ht) THEN
    2342           36 :          irange = [1, nindex]
    2343              :       ELSE
    2344           38 :          irange = get_irange_csr(nindex, mp_group)
    2345              :       END IF
    2346          178 :       DO m1 = 1, norb
    2347          128 :          i1 = active_orbitals(m1, 1)
    2348          430 :          DO m2 = m1, norb
    2349          252 :             i2 = active_orbitals(m2, 1)
    2350          252 :             i12 = csr_idx_to_combined(i1, i2, nmo_total)
    2351          380 :             IF (i12 >= irange(1) .AND. i12 <= irange(2)) THEN
    2352          144 :                i12l = i12 - irange(1) + 1
    2353          144 :                irptr = eri%rowptr_local(i12l) - 1
    2354          605 :                DO i34l = 1, eri%nzerow_local(i12l)
    2355          461 :                   i34 = eri%colind_local(irptr + i34l)
    2356          461 :                   CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
    2357          461 :                   erint = eri%nzval_local%r_dp(irptr + i34l)
    2358              :                   ! Coulomb
    2359          461 :                   ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
    2360          461 :                   IF (i3 /= i4) THEN
    2361          246 :                      ks_ref(i1, i2) = ks_ref(i1, i2) + erint*p_mat(i3, i4)
    2362              :                   END IF
    2363          461 :                   IF (i12 /= i34) THEN
    2364          335 :                      ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
    2365          335 :                      IF (i1 /= i2) THEN
    2366          220 :                         ks_ref(i3, i4) = ks_ref(i3, i4) + erint*p_mat(i1, i2)
    2367              :                      END IF
    2368              :                   END IF
    2369              :                   ! Exchange
    2370          461 :                   erint = -0.5_dp*erint
    2371          461 :                   ks_ref(i1, i3) = ks_ref(i1, i3) + erint*p_mat(i2, i4)
    2372          461 :                   IF (i1 /= i2) THEN
    2373          282 :                      ks_ref(i2, i3) = ks_ref(i2, i3) + erint*p_mat(i1, i4)
    2374              :                   END IF
    2375          461 :                   IF (i3 /= i4) THEN
    2376          246 :                      ks_ref(i1, i4) = ks_ref(i1, i4) + erint*p_mat(i2, i3)
    2377              :                   END IF
    2378         1066 :                   IF (i1 /= i2 .AND. i3 /= i4) THEN
    2379          193 :                      ks_ref(i2, i4) = ks_ref(i2, i4) + erint*p_mat(i1, i3)
    2380              :                   END IF
    2381              :                END DO
    2382              :             END IF
    2383              :          END DO
    2384              :       END DO
    2385              :       !
    2386          178 :       DO m1 = 1, norb
    2387          128 :          i1 = active_orbitals(m1, 1)
    2388          430 :          DO m2 = m1, norb
    2389          252 :             i2 = active_orbitals(m2, 1)
    2390          380 :             ks_ref(i2, i1) = ks_ref(i1, i2)
    2391              :          END DO
    2392              :       END DO
    2393         1402 :       CALL mp_group%sum(ks_ref)
    2394              : 
    2395           50 :    END SUBROUTINE build_subspace_fock_matrix
    2396              : 
    2397              : ! **************************************************************************************************
    2398              : !> \brief build subspace fockian for unrestricted spins
    2399              : !> \param active_orbitals the active orbital indices
    2400              : !> \param eri_aa two electon integrals in MO with parallel spins
    2401              : !> \param eri_ab two electon integrals in MO with anti-parallel spins
    2402              : !> \param p_a_mat density matrix for up-spin
    2403              : !> \param p_b_mat density matrix for down-spin
    2404              : !> \param ks_a_ref fockian matrix for up-spin
    2405              : !> \param tr_mixed_eri boolean to indicate Coulomb interaction alignment
    2406              : !> \param eri_method ...
    2407              : ! **************************************************************************************************
    2408           24 :    SUBROUTINE build_subspace_spin_fock_matrix(active_orbitals, eri_aa, eri_ab, p_a_mat, p_b_mat, ks_a_ref, tr_mixed_eri, eri_method)
    2409              :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: active_orbitals
    2410              :       TYPE(dbcsr_csr_type), INTENT(IN)                   :: eri_aa, eri_ab
    2411              :       REAL(dp), DIMENSION(:, :), INTENT(IN)              :: p_a_mat, p_b_mat
    2412              :       REAL(dp), DIMENSION(:, :), INTENT(INOUT)           :: ks_a_ref
    2413              :       LOGICAL, INTENT(IN)                                :: tr_mixed_eri
    2414              :       INTEGER, INTENT(IN)                                :: eri_method
    2415              : 
    2416              :       INTEGER                                            :: i1, i12, i12l, i2, i3, i34, i34l, i4, &
    2417              :                                                             irptr, m1, m2, nindex, nmo_total, &
    2418              :                                                             norb, spin1, spin2
    2419              :       INTEGER, DIMENSION(2)                              :: irange
    2420              :       REAL(dp)                                           :: erint
    2421              :       TYPE(mp_comm_type)                                 :: mp_group
    2422              : 
    2423           24 :       norb = SIZE(active_orbitals, 1)
    2424           24 :       nmo_total = SIZE(p_a_mat, 1)
    2425           24 :       nindex = (nmo_total*(nmo_total + 1))/2
    2426           24 :       CALL mp_group%set_handle(eri_aa%mp_group%get_handle())
    2427           24 :       IF (eri_method == eri_method_gpw_ht) THEN
    2428            0 :          irange = [1, nindex]
    2429              :       ELSE
    2430           24 :          irange = get_irange_csr(nindex, mp_group)
    2431              :       END IF
    2432           24 :       IF (tr_mixed_eri) THEN
    2433              :          spin1 = 2
    2434           24 :          spin2 = 1
    2435              :       ELSE
    2436           12 :          spin1 = 1
    2437           12 :          spin2 = 2
    2438              :       END IF
    2439           84 :       DO m1 = 1, norb
    2440           60 :          i1 = active_orbitals(m1, spin1)
    2441          192 :          DO m2 = m1, norb
    2442          108 :             i2 = active_orbitals(m2, spin1)
    2443          108 :             i12 = csr_idx_to_combined(i1, i2, nmo_total)
    2444          168 :             IF (i12 >= irange(1) .AND. i12 <= irange(2)) THEN
    2445           54 :                i12l = i12 - irange(1) + 1
    2446           54 :                irptr = eri_aa%rowptr_local(i12l) - 1
    2447          184 :                DO i34l = 1, eri_aa%nzerow_local(i12l)
    2448          130 :                   i34 = eri_aa%colind_local(irptr + i34l)
    2449          130 :                   CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
    2450          130 :                   erint = eri_aa%nzval_local%r_dp(irptr + i34l)
    2451              :                   ! Coulomb
    2452              :                   !F_ij += (ij|kl)*d_kl
    2453          130 :                   ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_a_mat(i3, i4)
    2454          130 :                   IF (i12 /= i34) THEN
    2455              :                      !F_kl += (ij|kl)*d_ij
    2456           76 :                      ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_a_mat(i1, i2)
    2457              :                   END IF
    2458              :                   ! Exchange
    2459          130 :                   erint = -1.0_dp*erint
    2460              :                   !F_ik -= (ij|kl)*d_jl
    2461          130 :                   ks_a_ref(i1, i3) = ks_a_ref(i1, i3) + erint*p_a_mat(i2, i4)
    2462          130 :                   IF (i1 /= i2) THEN
    2463              :                      !F_jk -= (ij|kl)*d_il
    2464           53 :                      ks_a_ref(i2, i3) = ks_a_ref(i2, i3) + erint*p_a_mat(i1, i4)
    2465              :                   END IF
    2466          130 :                   IF (i3 /= i4) THEN
    2467              :                      !F_il -= (ij|kl)*d_jk
    2468           56 :                      ks_a_ref(i1, i4) = ks_a_ref(i1, i4) + erint*p_a_mat(i2, i3)
    2469              :                   END IF
    2470          314 :                   IF (i1 /= i2 .AND. i3 /= i4) THEN
    2471              :                      !F_jl -= (ij|kl)*d_ik
    2472           33 :                      ks_a_ref(i2, i4) = ks_a_ref(i2, i4) + erint*p_a_mat(i1, i3)
    2473              :                   END IF
    2474              :                END DO
    2475              :             END IF
    2476              :          END DO
    2477              :       END DO
    2478              :       !
    2479              : 
    2480           24 :       CALL mp_group%set_handle(eri_ab%mp_group%get_handle())
    2481           24 :       IF (eri_method == eri_method_gpw_ht) THEN
    2482            0 :          irange = [1, nindex]
    2483              :       ELSE
    2484           24 :          irange = get_irange_csr(nindex, mp_group)
    2485              :       END IF
    2486           84 :       DO m1 = 1, norb
    2487           60 :          i1 = active_orbitals(m1, 1)
    2488          192 :          DO m2 = m1, norb
    2489          108 :             i2 = active_orbitals(m2, 1)
    2490          108 :             i12 = csr_idx_to_combined(i1, i2, nmo_total)
    2491          168 :             IF (i12 >= irange(1) .AND. i12 <= irange(2)) THEN
    2492           54 :                i12l = i12 - irange(1) + 1
    2493           54 :                irptr = eri_ab%rowptr_local(i12l) - 1
    2494          292 :                DO i34l = 1, eri_ab%nzerow_local(i12l)
    2495          238 :                   i34 = eri_ab%colind_local(irptr + i34l)
    2496          238 :                   CALL csr_idx_from_combined(i34, nmo_total, i3, i4)
    2497          238 :                   erint = eri_ab%nzval_local%r_dp(irptr + i34l)
    2498              :                   ! Coulomb
    2499          292 :                   IF (tr_mixed_eri) THEN
    2500              :                      !F_kl += (kl beta|ij alpha )*d_alpha_ij
    2501          119 :                      ks_a_ref(i3, i4) = ks_a_ref(i3, i4) + erint*p_b_mat(i1, i2)
    2502              :                   ELSE
    2503              :                      !F_ij += (ij alpha|kl beta )*d_beta_kl
    2504          119 :                      ks_a_ref(i1, i2) = ks_a_ref(i1, i2) + erint*p_b_mat(i3, i4)
    2505              :                   END IF
    2506              :                END DO
    2507              :             END IF
    2508              :          END DO
    2509              :       END DO
    2510              :       !
    2511           84 :       DO m1 = 1, norb
    2512           60 :          i1 = active_orbitals(m1, spin1)
    2513          192 :          DO m2 = m1, norb
    2514          108 :             i2 = active_orbitals(m2, spin1)
    2515          168 :             ks_a_ref(i2, i1) = ks_a_ref(i1, i2)
    2516              :          END DO
    2517              :       END DO
    2518           24 :       CALL mp_group%set_handle(eri_aa%mp_group%get_handle())
    2519         1080 :       CALL mp_group%sum(ks_a_ref)
    2520              : 
    2521           24 :    END SUBROUTINE build_subspace_spin_fock_matrix
    2522              : 
    2523              : ! **************************************************************************************************
    2524              : !> \brief Creates a local basis
    2525              : !> \param pro_basis_set ...
    2526              : !> \param zval ...
    2527              : !> \param ishell ...
    2528              : !> \param nshell ...
    2529              : !> \param lnam ...
    2530              : !> \par History
    2531              : !>      05.2016 created [JGH]
    2532              : ! **************************************************************************************************
    2533            0 :    SUBROUTINE create_pro_basis(pro_basis_set, zval, ishell, nshell, lnam)
    2534              :       TYPE(gto_basis_set_type), POINTER                  :: pro_basis_set
    2535              :       INTEGER, INTENT(IN)                                :: zval, ishell
    2536              :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nshell
    2537              :       CHARACTER(len=*), DIMENSION(:), INTENT(IN)         :: lnam
    2538              : 
    2539            0 :       CHARACTER(len=6), DIMENSION(:), POINTER            :: sym
    2540              :       INTEGER                                            :: i, l, nj
    2541              :       INTEGER, DIMENSION(4, 7)                           :: ne
    2542            0 :       INTEGER, DIMENSION(:), POINTER                     :: lq, nq
    2543            0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: zet
    2544              :       TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set
    2545              : 
    2546            0 :       CPASSERT(.NOT. ASSOCIATED(pro_basis_set))
    2547            0 :       NULLIFY (sto_basis_set)
    2548              : 
    2549              :       ! electronic configuration
    2550            0 :       ne = 0
    2551            0 :       DO l = 1, 4 !lq(1)+1
    2552            0 :          nj = 2*(l - 1) + 1
    2553            0 :          DO i = l, 7 ! nq(1)
    2554            0 :             ne(l, i) = ptable(zval)%e_conv(l - 1) - 2*nj*(i - l)
    2555            0 :             ne(l, i) = MAX(ne(l, i), 0)
    2556            0 :             ne(l, i) = MIN(ne(l, i), 2*nj)
    2557              :          END DO
    2558              :       END DO
    2559            0 :       ALLOCATE (nq(ishell), lq(ishell), zet(ishell), sym(ishell))
    2560            0 :       DO i = 1, ishell
    2561            0 :          nq(i) = nshell(i)
    2562            0 :          SELECT CASE (lnam(i))
    2563              :          CASE ('S', 's')
    2564            0 :             lq(i) = 0
    2565              :          CASE ('P', 'p')
    2566            0 :             lq(i) = 1
    2567              :          CASE ('D', 'd')
    2568            0 :             lq(i) = 2
    2569              :          CASE ('F', 'f')
    2570            0 :             lq(i) = 3
    2571              :          CASE DEFAULT
    2572            0 :             CPABORT("Wrong l QN")
    2573              :          END SELECT
    2574            0 :          sym(i) = lnam(i)
    2575            0 :          zet(i) = srules(zval, ne, nq(1), lq(1))
    2576              :       END DO
    2577            0 :       CALL allocate_sto_basis_set(sto_basis_set)
    2578            0 :       CALL set_sto_basis_set(sto_basis_set, nshell=1, nq=nq, lq=lq, zet=zet, symbol=sym)
    2579            0 :       CALL create_gto_from_sto_basis(sto_basis_set, pro_basis_set, 6)
    2580            0 :       pro_basis_set%norm_type = 2
    2581            0 :       CALL init_orb_basis_set(pro_basis_set)
    2582            0 :       CALL deallocate_sto_basis_set(sto_basis_set)
    2583              : 
    2584            0 :    END SUBROUTINE create_pro_basis
    2585              : 
    2586              : ! **************************************************************************************************
    2587              : !> \brief Update the density matrix in AO basis with the active density contribution
    2588              : !> \param active_space_env the active space environment
    2589              : !> \param rho_ao the density matrix in AO basis
    2590              : ! **************************************************************************************************
    2591            0 :    SUBROUTINE update_density_ao(active_space_env, rho_ao)
    2592              :       TYPE(active_space_type), POINTER                   :: active_space_env
    2593              :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
    2594              : 
    2595              :       INTEGER                                            :: ispin, nao, nmo, nspins
    2596              :       TYPE(cp_fm_type)                                   :: R, U
    2597              :       TYPE(cp_fm_type), POINTER                          :: C_active, p_active_mo
    2598              :       TYPE(dbcsr_type), POINTER                          :: p_inactive_ao
    2599            0 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos_active
    2600              : 
    2601              :       ! Transform the AS density matrix P_MO to the atomic orbital basis,
    2602              :       ! this is simply C * P_MO * C^T
    2603            0 :       nspins = active_space_env%nspins
    2604            0 :       mos_active => active_space_env%mos_active
    2605            0 :       DO ispin = 1, nspins
    2606              :          ! size of p_inactive_ao is (nao x nao)
    2607            0 :          p_inactive_ao => active_space_env%pmat_inactive(ispin)%matrix
    2608              : 
    2609              :          ! copy p_inactive_ao to rho_ao
    2610            0 :          CALL dbcsr_copy(rho_ao(ispin)%matrix, p_inactive_ao)
    2611              : 
    2612              :          ! size of p_active_mo is (nmo x nmo)
    2613            0 :          p_active_mo => active_space_env%p_active(ispin)
    2614              : 
    2615              :          ! calculate R = p_mo
    2616            0 :          CALL cp_fm_create(R, p_active_mo%matrix_struct)
    2617            0 :          CALL cp_fm_to_fm(p_active_mo, R)
    2618              : 
    2619              :          ! calculate U = C * p_mo
    2620            0 :          CALL get_mo_set(mos_active(ispin), mo_coeff=C_active, nao=nao, nmo=nmo)
    2621            0 :          CALL cp_fm_create(U, C_active%matrix_struct)
    2622            0 :          CALL parallel_gemm("N", "N", nao, nmo, nmo, 1.0_dp, C_active, R, 0.0_dp, U)
    2623              : 
    2624              :          CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=rho_ao(ispin)%matrix, &
    2625            0 :                                     matrix_v=U, matrix_g=C_active, ncol=nmo, alpha=1.0_dp)
    2626              : 
    2627            0 :          CALL cp_fm_release(R)
    2628            0 :          CALL cp_fm_release(U)
    2629              :       END DO
    2630              : 
    2631            0 :    END SUBROUTINE update_density_ao
    2632              : 
    2633              : ! **************************************************************************************************
    2634              : !> \brief Print each value on the master node
    2635              : !> \param this object reference
    2636              : !> \param i i-index
    2637              : !> \param j j-index
    2638              : !> \param k k-index
    2639              : !> \param l l-index
    2640              : !> \param val value of the integral at (i,j,k.l)
    2641              : !> \return always true to dump all integrals
    2642              : ! **************************************************************************************************
    2643         1420 :    LOGICAL FUNCTION eri_fcidump_print_func(this, i, j, k, l, val) RESULT(cont)
    2644              :       CLASS(eri_fcidump_print), INTENT(inout) :: this
    2645              :       INTEGER, INTENT(in)                     :: i, j, k, l
    2646              :       REAL(KIND=dp), INTENT(in)               :: val
    2647              : 
    2648              :       ! write to the actual file only on the master
    2649         1420 :       IF (this%unit_nr > 0) THEN
    2650          710 :          WRITE (this%unit_nr, "(ES23.16,4I4)") val, i + this%bra_start - 1, j + this%bra_start - 1, &
    2651         1420 :               &                                     k + this%ket_start - 1, l + this%ket_start - 1
    2652              :       END IF
    2653              : 
    2654         1420 :       cont = .TRUE.
    2655         1420 :    END FUNCTION eri_fcidump_print_func
    2656              : 
    2657              : ! **************************************************************************************************
    2658              : !> \brief checksum each value on the master node
    2659              : !> \param this object reference
    2660              : !> \param i i-index
    2661              : !> \param j j-index
    2662              : !> \param k k-index
    2663              : !> \param l l-index
    2664              : !> \param val value of the integral at (i,j,k.l)
    2665              : !> \return always true to dump all integrals
    2666              : ! **************************************************************************************************
    2667         1420 :    LOGICAL FUNCTION eri_fcidump_checksum_func(this, i, j, k, l, val) RESULT(cont)
    2668              :       CLASS(eri_fcidump_checksum), INTENT(inout) :: this
    2669              :       INTEGER, INTENT(in)                     :: i, j, k, l
    2670              :       REAL(KIND=dp), INTENT(in)               :: val
    2671              :       MARK_USED(i)
    2672              :       MARK_USED(j)
    2673              :       MARK_USED(k)
    2674              :       MARK_USED(l)
    2675              : 
    2676         1420 :       this%checksum = this%checksum + ABS(val)
    2677              : 
    2678         1420 :       cont = .TRUE.
    2679         1420 :    END FUNCTION eri_fcidump_checksum_func
    2680              : 
    2681              : ! **************************************************************************************************
    2682              : !> \brief Update active space density matrix from a fortran array
    2683              : !> \param p_act_mo density matrix in active space MO basis
    2684              : !> \param active_space_env active space environment
    2685              : !> \param ispin spin index
    2686              : !> \author Vladimir Rybkin
    2687              : ! **************************************************************************************************
    2688            0 :    SUBROUTINE update_active_density(p_act_mo, active_space_env, ispin)
    2689              :       REAL(KIND=dp), DIMENSION(:)                        :: p_act_mo
    2690              :       TYPE(active_space_type), POINTER                   :: active_space_env
    2691              :       INTEGER                                            :: ispin
    2692              : 
    2693              :       INTEGER                                            :: i1, i2, m1, m2, nmo_active
    2694              :       REAL(KIND=dp)                                      :: mval
    2695              :       TYPE(cp_fm_type), POINTER                          :: p_active
    2696              : 
    2697            0 :       p_active => active_space_env%p_active(ispin)
    2698            0 :       nmo_active = active_space_env%nmo_active
    2699              : 
    2700            0 :       DO i1 = 1, nmo_active
    2701            0 :          m1 = active_space_env%active_orbitals(i1, ispin)
    2702            0 :          DO i2 = 1, nmo_active
    2703            0 :             m2 = active_space_env%active_orbitals(i2, ispin)
    2704            0 :             mval = p_act_mo(i2 + (i1 - 1)*nmo_active)
    2705            0 :             CALL cp_fm_set_element(p_active, m1, m2, mval)
    2706              :          END DO
    2707              :       END DO
    2708              : 
    2709            0 :    END SUBROUTINE update_active_density
    2710              : 
    2711              : ! **************************************************************************************************
    2712              : !> \brief ...
    2713              : !> \param qs_env ...
    2714              : !> \param active_space_env ...
    2715              : !> \param as_input ...
    2716              : ! **************************************************************************************************
    2717            0 :    SUBROUTINE rsdft_embedding(qs_env, active_space_env, as_input)
    2718              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    2719              :       TYPE(active_space_type), POINTER                   :: active_space_env
    2720              :       TYPE(section_vals_type), POINTER                   :: as_input
    2721              : 
    2722              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'rsdft_embedding'
    2723              :       INTEGER                                            :: handle
    2724              : 
    2725              : #ifdef __NO_SOCKETS
    2726              :       CALL timeset(routineN, handle)
    2727              :       CPABORT("CP2K was compiled with the __NO_SOCKETS option!")
    2728              :       MARK_USED(qs_env)
    2729              :       MARK_USED(active_space_env)
    2730              :       MARK_USED(as_input)
    2731              : #else
    2732              : 
    2733              :       INTEGER                                            :: iw, client_fd, socket_fd, iter, max_iter
    2734              :       LOGICAL                                            :: converged, do_scf_embedding, ionode
    2735              :       REAL(KIND=dp)                                      :: alpha, delta_E, energy_corr, energy_new, &
    2736              :                                                             energy_old, energy_scf, eps_iter, t1, &
    2737              :                                                             t2
    2738              :       TYPE(cp_logger_type), POINTER                      :: logger
    2739            0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao
    2740            0 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos_active
    2741              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2742              :       TYPE(qs_energy_type), POINTER                      :: energy
    2743              :       TYPE(qs_rho_type), POINTER                         :: rho
    2744              : 
    2745            0 :       CALL timeset(routineN, handle)
    2746              : 
    2747            0 :       t1 = m_walltime()
    2748              : 
    2749            0 :       logger => cp_get_default_logger()
    2750            0 :       iw = cp_logger_get_default_io_unit(logger)
    2751              : 
    2752            0 :       CALL get_qs_env(qs_env, para_env=para_env)
    2753            0 :       ionode = para_env%is_source()
    2754              : 
    2755              :       ! get info from the input
    2756            0 :       CALL section_vals_val_get(as_input, "SCF_EMBEDDING", l_val=do_scf_embedding)
    2757            0 :       active_space_env%do_scf_embedding = do_scf_embedding
    2758            0 :       CALL section_vals_val_get(as_input, "MAX_ITER", i_val=max_iter)
    2759            0 :       CALL section_vals_val_get(as_input, "EPS_ITER", r_val=eps_iter)
    2760            0 :       alpha = 0.0
    2761              : 
    2762              :       ! create the socket and wait for the client to connect
    2763            0 :       CALL initialize_socket(socket_fd, client_fd, as_input, ionode)
    2764            0 :       CALL para_env%sync()
    2765              : 
    2766              :       ! send two-electron integrals to the client
    2767            0 :       CALL send_eri_to_client(client_fd, active_space_env, para_env)
    2768              : 
    2769              :       ! get pointer to density in ao basis
    2770            0 :       CALL get_qs_env(qs_env, rho=rho, energy=energy)
    2771            0 :       CALL qs_rho_get(rho, rho_ao=rho_ao)
    2772              : 
    2773            0 :       IF ((iw > 0)) THEN
    2774              :          WRITE (UNIT=iw, FMT="(/,T3,A,T11,A,T21,A,T34,A,T55,A,T75,A,/,T3,A)") &
    2775            0 :             "Iter", "Update", "Time", "Corr. energy", "Total energy", "Change", REPEAT("-", 78)
    2776              :       END IF
    2777              :       ! CALL cp_add_iter_level(logger%iter_info, "QS_SCF")
    2778              : 
    2779            0 :       iter = 0
    2780            0 :       converged = .FALSE.
    2781              :       ! store the scf energy
    2782            0 :       energy_scf = active_space_env%energy_ref
    2783            0 :       energy_new = energy_scf
    2784            0 :       mos_active => active_space_env%mos_active
    2785              :       ! CALL set_qs_env(qs_env, active_space=active_space_env)
    2786              : 
    2787              :       ! start the self-consistent embedding loop
    2788            0 :       DO WHILE (iter < max_iter)
    2789            0 :          iter = iter + 1
    2790              : 
    2791              :          ! send V_emb and E_ina to the active space solver and update
    2792              :          ! the active space environment with the new active energy and density
    2793            0 :          CALL send_fock_to_client(client_fd, active_space_env, para_env)
    2794              : 
    2795              :          ! update energies
    2796            0 :          energy_old = energy_new
    2797            0 :          energy_new = active_space_env%energy_total
    2798            0 :          energy_corr = energy_new - energy_scf
    2799            0 :          delta_E = energy_new - energy_old
    2800              : 
    2801              :          ! get timer
    2802            0 :          t2 = m_walltime()
    2803              :          ! print out progress
    2804            0 :          IF ((iw > 0)) THEN
    2805              :             WRITE (UNIT=iw, &
    2806              :                    FMT="(T3,I4,T11,A,T21,F4.2,T28,F18.10,T49,F18.10,T70,ES11.2)") &
    2807            0 :                iter, 'P_Mix', t2 - t1, energy_corr, energy_new, delta_E
    2808              :          END IF
    2809              : 
    2810              :          ! update total density in AO basis with the AS contribution
    2811            0 :          CALL update_density_ao(active_space_env, rho_ao) ! rho_ao is updated
    2812              : 
    2813              :          ! calculate F_ks in AO basis (which contains Vxc) with the new density
    2814            0 :          qs_env%requires_matrix_vxc = .TRUE.
    2815            0 :          CALL qs_rho_update_rho(rho, qs_env) ! updates rho_r and rho_g using rho_ao
    2816            0 :          CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.) ! set flags about the change
    2817            0 :          CALL qs_ks_update_qs_env(qs_env) ! this call actually calculates F_ks
    2818              : 
    2819              :          ! update the reference energy
    2820            0 :          active_space_env%energy_ref = energy%total
    2821              : 
    2822              :          ! transform KS/Fock, Vxc and Hcore from AO to MO basis
    2823            0 :          CALL calculate_operators(mos_active, qs_env, active_space_env)
    2824              : 
    2825              :          ! calculate the new inactive energy and embedding potential
    2826            0 :          CALL subspace_fock_matrix(active_space_env)
    2827              : 
    2828              :          ! check if it is a one-shot correction
    2829            0 :          IF (.NOT. active_space_env%do_scf_embedding) THEN
    2830            0 :             IF (iw > 0) THEN
    2831              :                WRITE (UNIT=iw, FMT="(/,T3,A,I5,A/)") &
    2832            0 :                   "*** one-shot embedding correction finished ***"
    2833              :             END IF
    2834              :             converged = .TRUE.
    2835              :             EXIT
    2836              :             ! check for convergence
    2837            0 :          ELSEIF (ABS(delta_E) <= eps_iter) THEN
    2838            0 :             IF (iw > 0) THEN
    2839              :                WRITE (UNIT=iw, FMT="(/,T3,A,I5,A/)") &
    2840            0 :                   "*** rsDFT embedding run converged in ", iter, " iteration(s) ***"
    2841              :             END IF
    2842              :             converged = .TRUE.
    2843              :             EXIT
    2844              :          END IF
    2845              : 
    2846            0 :          t1 = m_walltime()
    2847              :       END DO
    2848              : 
    2849              :       IF (.NOT. converged) THEN
    2850            0 :          IF (iw > 0) THEN
    2851              :             WRITE (UNIT=iw, FMT="(/,T3,A,I5,A/)") &
    2852            0 :                "*** rsDFT embedding did not converged after ", iter, " iteration(s) ***"
    2853              :          END IF
    2854              :       END IF
    2855              : 
    2856              :       ! update qs total energy to the final embedding energy
    2857            0 :       energy%total = active_space_env%energy_total
    2858              : 
    2859            0 :       CALL finalize_socket(socket_fd, client_fd, as_input, ionode)
    2860            0 :       CALL para_env%sync()
    2861              : #endif
    2862              : 
    2863            0 :       CALL timestop(handle)
    2864              : 
    2865            0 :    END SUBROUTINE rsdft_embedding
    2866              : 
    2867              : #ifndef __NO_SOCKETS
    2868              : ! **************************************************************************************************
    2869              : !> \brief Creates the socket, spawns the client and connects to it
    2870              : !> \param socket_fd the socket file descriptor
    2871              : !> \param client_fd the client file descriptor
    2872              : !> \param as_input active space inpute section
    2873              : !> \param ionode logical flag indicating if the process is the master
    2874              : ! **************************************************************************************************
    2875            0 :    SUBROUTINE initialize_socket(socket_fd, client_fd, as_input, ionode)
    2876              :       INTEGER, INTENT(OUT)                               :: socket_fd, client_fd
    2877              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: as_input
    2878              :       LOGICAL, INTENT(IN)                                :: ionode
    2879              : 
    2880              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'initialize_socket'
    2881              :       INTEGER, PARAMETER                                 :: backlog = 10
    2882              : 
    2883              :       CHARACTER(len=default_path_length)                 :: hostname
    2884              :       INTEGER                                            :: handle, iw, port, protocol
    2885              :       LOGICAL                                            :: inet
    2886              :       TYPE(cp_logger_type), POINTER                      :: logger
    2887              : 
    2888            0 :       CALL timeset(routineN, handle)
    2889              : 
    2890            0 :       logger => cp_get_default_logger()
    2891            0 :       iw = cp_logger_get_default_io_unit(logger)
    2892              : 
    2893              :       ! protocol == 0 for UNIX, protocol > 0 for INET
    2894            0 :       CALL section_vals_val_get(as_input, "SOCKET%INET", l_val=inet)
    2895            0 :       IF (inet) THEN
    2896            0 :          protocol = 1
    2897              :       ELSE
    2898            0 :          protocol = 0
    2899              :       END IF
    2900            0 :       CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname)
    2901            0 :       CALL section_vals_val_get(as_input, "SOCKET%PORT", i_val=port)
    2902              : 
    2903            0 :       IF (ionode) THEN
    2904            0 :          CALL open_bind_socket(socket_fd, protocol, port, TRIM(hostname)//C_NULL_CHAR)
    2905            0 :          WRITE (iw, '(/,T3,A,A)') "@SERVER: Created socket with address ", TRIM(hostname)
    2906            0 :          CALL listen_socket(socket_fd, backlog)
    2907              : 
    2908              :          ! wait until a connetion request arrives
    2909            0 :          WRITE (iw, '(T3,A)') "@SERVER: Waiting for requests..."
    2910            0 :          CALL accept_socket(socket_fd, client_fd)
    2911            0 :          WRITE (iw, '(T3,A,I2)') "@SERVER: Accepted socket with fd ", client_fd
    2912              :       END IF
    2913              : 
    2914            0 :       CALL timestop(handle)
    2915              : 
    2916            0 :    END SUBROUTINE initialize_socket
    2917              : 
    2918              : ! **************************************************************************************************
    2919              : !> \brief Closes the connection to the socket and deletes the file
    2920              : !> \param socket_fd the socket file descriptor
    2921              : !> \param client_fd the client file descriptor
    2922              : !> \param as_input active space inpute section
    2923              : !> \param ionode logical flag indicating if the process is the master
    2924              : ! **************************************************************************************************
    2925            0 :    SUBROUTINE finalize_socket(socket_fd, client_fd, as_input, ionode)
    2926              :       INTEGER, INTENT(IN)                                :: socket_fd, client_fd
    2927              :       TYPE(section_vals_type), INTENT(IN), POINTER       :: as_input
    2928              :       LOGICAL, INTENT(IN)                                :: ionode
    2929              : 
    2930              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'finalize_socket'
    2931              :       INTEGER, PARAMETER                                 :: header_len = 12
    2932              : 
    2933              :       CHARACTER(len=default_path_length)                 :: hostname
    2934              :       INTEGER                                            :: handle
    2935              : 
    2936            0 :       CALL timeset(routineN, handle)
    2937              : 
    2938            0 :       CALL section_vals_val_get(as_input, "SOCKET%HOST", c_val=hostname)
    2939              : 
    2940            0 :       IF (ionode) THEN
    2941              :          ! signal the client to quit
    2942            0 :          CALL writebuffer(client_fd, "QUIT        ", header_len)
    2943              :          ! close the connection
    2944            0 :          CALL close_socket(client_fd)
    2945            0 :          CALL close_socket(socket_fd)
    2946              : 
    2947              :          ! delete the socket file
    2948            0 :          IF (file_exists(TRIM(hostname))) THEN
    2949            0 :             CALL remove_socket_file(TRIM(hostname)//C_NULL_CHAR)
    2950              :          END IF
    2951              :       END IF
    2952              : 
    2953            0 :       CALL timestop(handle)
    2954              : 
    2955            0 :    END SUBROUTINE finalize_socket
    2956              : 
    2957              : ! **************************************************************************************************
    2958              : !> \brief Sends the two-electron integrals to the client vie the socket
    2959              : !> \param client_fd the client file descriptor
    2960              : !> \param active_space_env active space environment
    2961              : !> \param para_env parallel environment
    2962              : ! **************************************************************************************************
    2963            0 :    SUBROUTINE send_eri_to_client(client_fd, active_space_env, para_env)
    2964              :       INTEGER, INTENT(IN)                                :: client_fd
    2965              :       TYPE(active_space_type), INTENT(IN), POINTER       :: active_space_env
    2966              :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
    2967              : 
    2968              :       CHARACTER(len=*), PARAMETER :: routineN = 'send_eri_to_client'
    2969              :       INTEGER, PARAMETER                                 :: header_len = 12
    2970              : 
    2971              :       CHARACTER(len=default_string_length)               :: header
    2972              :       INTEGER                                            :: handle, iw
    2973              :       LOGICAL                                            :: ionode
    2974            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: eri_aa, eri_ab, eri_bb
    2975              :       TYPE(cp_logger_type), POINTER                      :: logger
    2976              : 
    2977            0 :       CALL timeset(routineN, handle)
    2978              : 
    2979            0 :       logger => cp_get_default_logger()
    2980            0 :       iw = cp_logger_get_default_io_unit(logger)
    2981            0 :       ionode = para_env%is_source()
    2982              : 
    2983              :       ! TODO: do we really need to allocate the arrays on every process?
    2984            0 :       ALLOCATE (eri_aa(active_space_env%nmo_active**4))
    2985            0 :       CALL eri_to_array(active_space_env%eri, eri_aa, active_space_env%active_orbitals, 1, 1)
    2986            0 :       IF (active_space_env%nspins == 2) THEN
    2987            0 :          ALLOCATE (eri_ab(active_space_env%nmo_active**4))
    2988            0 :          CALL eri_to_array(active_space_env%eri, eri_ab, active_space_env%active_orbitals, 1, 2)
    2989            0 :          ALLOCATE (eri_bb(active_space_env%nmo_active**4))
    2990            0 :          CALL eri_to_array(active_space_env%eri, eri_bb, active_space_env%active_orbitals, 2, 2)
    2991              :       END IF
    2992              : 
    2993              :       ! ask the status of the client
    2994            0 :       IF (ionode) CALL writebuffer(client_fd, "STATUS      ", header_len)
    2995              :       DO
    2996            0 :          header = ""
    2997            0 :          CALL para_env%sync()
    2998            0 :          IF (ionode) THEN
    2999              :             ! IF (iw > 0) WRITE(iw, *) "@SERVER: Waiting for messages..."
    3000            0 :             CALL readbuffer(client_fd, header, header_len)
    3001              :          END IF
    3002            0 :          CALL para_env%bcast(header, para_env%source)
    3003              : 
    3004              :          ! IF (iw > 0) WRITE(iw, *) "@SERVER: Message from client: ", TRIM(header)
    3005              : 
    3006            0 :          IF (TRIM(header) == "READY") THEN
    3007              :             ! if the client is ready, send the data
    3008            0 :             CALL para_env%sync()
    3009            0 :             IF (ionode) THEN
    3010            0 :                CALL writebuffer(client_fd, "TWOBODY     ", header_len)
    3011            0 :                CALL writebuffer(client_fd, active_space_env%nspins)
    3012            0 :                CALL writebuffer(client_fd, active_space_env%nmo_active)
    3013            0 :                CALL writebuffer(client_fd, active_space_env%nelec_active)
    3014            0 :                CALL writebuffer(client_fd, active_space_env%multiplicity)
    3015              :                ! send the alpha component
    3016            0 :                CALL writebuffer(client_fd, eri_aa, SIZE(eri_aa))
    3017              :                ! send the beta part for unrestricted calculations
    3018            0 :                IF (active_space_env%nspins == 2) THEN
    3019            0 :                   CALL writebuffer(client_fd, eri_ab, SIZE(eri_ab))
    3020            0 :                   CALL writebuffer(client_fd, eri_bb, SIZE(eri_bb))
    3021              :                END IF
    3022              :             END IF
    3023            0 :          ELSE IF (TRIM(header) == "RECEIVED") THEN
    3024              :             EXIT
    3025              :          END IF
    3026              :       END DO
    3027              : 
    3028            0 :       DEALLOCATE (eri_aa)
    3029            0 :       IF (active_space_env%nspins == 2) THEN
    3030            0 :          DEALLOCATE (eri_ab)
    3031            0 :          DEALLOCATE (eri_bb)
    3032              :       END IF
    3033              : 
    3034            0 :       CALL para_env%sync()
    3035              : 
    3036            0 :       CALL timestop(handle)
    3037              : 
    3038            0 :    END SUBROUTINE send_eri_to_client
    3039              : 
    3040              : ! **************************************************************************************************
    3041              : !> \brief Sends the one-electron embedding potential and the inactive energy to the client
    3042              : !> \param client_fd the client file descriptor
    3043              : !> \param active_space_env active space environment
    3044              : !> \param para_env parallel environment
    3045              : ! **************************************************************************************************
    3046            0 :    SUBROUTINE send_fock_to_client(client_fd, active_space_env, para_env)
    3047              :       INTEGER, INTENT(IN)                                :: client_fd
    3048              :       TYPE(active_space_type), INTENT(IN), POINTER       :: active_space_env
    3049              :       TYPE(mp_para_env_type), INTENT(IN), POINTER        :: para_env
    3050              : 
    3051              :       CHARACTER(len=*), PARAMETER :: routineN = 'send_fock_to_client'
    3052              :       INTEGER, PARAMETER                                 :: header_len = 12
    3053              : 
    3054              :       CHARACTER(len=default_string_length)               :: header
    3055              :       INTEGER                                            :: handle, iw
    3056              :       LOGICAL                                            :: debug, ionode
    3057            0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: fock_a, fock_b, p_act_mo_a, p_act_mo_b
    3058              :       TYPE(cp_logger_type), POINTER                      :: logger
    3059              : 
    3060            0 :       CALL timeset(routineN, handle)
    3061              : 
    3062              :       ! Set to .TRUE. to activate debug output
    3063            0 :       debug = .FALSE.
    3064              : 
    3065            0 :       logger => cp_get_default_logger()
    3066            0 :       iw = cp_logger_get_default_io_unit(logger)
    3067            0 :       ionode = para_env%is_source()
    3068              : 
    3069            0 :       ALLOCATE (p_act_mo_a(active_space_env%nmo_active**2))
    3070            0 :       ALLOCATE (fock_a(active_space_env%nmo_active**2))
    3071            0 :       IF (active_space_env%nspins == 2) THEN
    3072            0 :          ALLOCATE (p_act_mo_b(active_space_env%nmo_active**2))
    3073            0 :          ALLOCATE (fock_b(active_space_env%nmo_active**2))
    3074              :       END IF
    3075              : 
    3076              :       ! get the fock matrix into Fortran arrays
    3077              :       ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 1))
    3078            0 :          CALL subspace_matrix_to_array(active_space_env%fock_sub(1), fock_a, act_indices, act_indices)
    3079              :       END ASSOCIATE
    3080              : 
    3081            0 :       IF (active_space_env%nspins == 2) THEN
    3082              :          ASSOCIATE (act_indices => active_space_env%active_orbitals(:, 2))
    3083            0 :             CALL subspace_matrix_to_array(active_space_env%fock_sub(2), fock_b, act_indices, act_indices)
    3084              :          END ASSOCIATE
    3085              :       END IF
    3086              : 
    3087              :       ! ask the status of the client
    3088            0 :       IF (ionode) CALL writebuffer(client_fd, "STATUS      ", header_len)
    3089              :       DO
    3090            0 :          header = ""
    3091              : 
    3092            0 :          CALL para_env%sync()
    3093            0 :          IF (ionode) THEN
    3094              :             IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Waiting for messages..."
    3095            0 :             CALL readbuffer(client_fd, header, header_len)
    3096              :          END IF
    3097            0 :          CALL para_env%bcast(header, para_env%source)
    3098              : 
    3099              :          IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Message from client: ", TRIM(header)
    3100              : 
    3101            0 :          IF (TRIM(header) == "READY") THEN
    3102              :             ! if the client is ready, send the data
    3103            0 :             CALL para_env%sync()
    3104            0 :             IF (ionode) THEN
    3105            0 :                CALL writebuffer(client_fd, "ONEBODY     ", header_len)
    3106            0 :                CALL writebuffer(client_fd, active_space_env%energy_inactive)
    3107              :                ! send the alpha component
    3108            0 :                CALL writebuffer(client_fd, fock_a, SIZE(fock_a))
    3109              :                ! send the beta part for unrestricted calculations
    3110            0 :                IF (active_space_env%nspins == 2) THEN
    3111            0 :                   CALL writebuffer(client_fd, fock_b, SIZE(fock_b))
    3112              :                END IF
    3113              :             END IF
    3114              : 
    3115            0 :          ELSE IF (TRIM(header) == "HAVEDATA") THEN
    3116              :             ! qiskit has data to transfer, let them know we want it and wait for it
    3117            0 :             CALL para_env%sync()
    3118            0 :             IF (ionode) THEN
    3119              :                IF (debug .AND. iw > 0) WRITE (iw, *) "@SERVER: Qiskit has data to transfer"
    3120            0 :                CALL writebuffer(client_fd, "GETDENSITY  ", header_len)
    3121              : 
    3122              :                ! read the active energy and density
    3123            0 :                CALL readbuffer(client_fd, active_space_env%energy_active)
    3124            0 :                CALL readbuffer(client_fd, p_act_mo_a, SIZE(p_act_mo_a))
    3125            0 :                IF (active_space_env%nspins == 2) THEN
    3126            0 :                   CALL readbuffer(client_fd, p_act_mo_b, SIZE(p_act_mo_b))
    3127              :                END IF
    3128              :             END IF
    3129              : 
    3130              :             ! broadcast the data to all processors
    3131            0 :             CALL para_env%bcast(active_space_env%energy_active, para_env%source)
    3132            0 :             CALL para_env%bcast(p_act_mo_a, para_env%source)
    3133            0 :             IF (active_space_env%nspins == 2) THEN
    3134            0 :                CALL para_env%bcast(p_act_mo_b, para_env%source)
    3135              :             END IF
    3136              : 
    3137              :             ! update total and reference energies in active space enviornment
    3138            0 :             active_space_env%energy_total = active_space_env%energy_inactive + active_space_env%energy_active
    3139              : 
    3140              :             ! update the active density matrix in the active space environment
    3141            0 :             CALL update_active_density(p_act_mo_a, active_space_env, 1)
    3142            0 :             IF (active_space_env%nspins == 2) THEN
    3143            0 :                CALL update_active_density(p_act_mo_b, active_space_env, 2)
    3144              :             END IF
    3145              : 
    3146              :             ! the non-iterative part is done, we can continue
    3147              :             EXIT
    3148              :          END IF
    3149              : 
    3150              :       END DO
    3151              : 
    3152            0 :       DEALLOCATE (p_act_mo_a)
    3153            0 :       DEALLOCATE (fock_a)
    3154            0 :       IF (active_space_env%nspins == 2) THEN
    3155            0 :          DEALLOCATE (p_act_mo_b)
    3156            0 :          DEALLOCATE (fock_b)
    3157              :       END IF
    3158              : 
    3159            0 :       CALL para_env%sync()
    3160              : 
    3161            0 :       CALL timestop(handle)
    3162              : 
    3163            0 :    END SUBROUTINE send_fock_to_client
    3164              : #endif
    3165              : 
    3166            0 : END MODULE qs_active_space_methods
        

Generated by: LCOV version 2.0-1