LCOV - code coverage report
Current view: top level - src - qs_core_hamiltonian.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.2 % 279 274
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 6 6

            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 Calculation of the core Hamiltonian integral matrix <a|H|b> over
      10              : !>      Cartesian Gaussian-type functions.
      11              : !>
      12              : !>      <a|H|b> = <a|T|b> + <a|V|b>
      13              : !>
      14              : !>      Kinetic energy:
      15              : !>
      16              : !>      <a|T|b> = <a|-nabla**2/2|b>
      17              : !>                \_______________/
      18              : !>                        |
      19              : !>                     kinetic
      20              : !>
      21              : !>      Nuclear potential energy:
      22              : !>
      23              : !>      a) Allelectron calculation:
      24              : !>
      25              : !>                          erfc(r)
      26              : !>         <a|V|b> = -Z*<a|---------|b>
      27              : !>                             r
      28              : !>
      29              : !>                          1 - erf(r)
      30              : !>                 = -Z*<a|------------|b>
      31              : !>                              r
      32              : !>
      33              : !>                           1           erf(r)
      34              : !>                 = -Z*(<a|---|b> - <a|--------|b>)
      35              : !>                           r             r
      36              : !>
      37              : !>                           1
      38              : !>                 = -Z*(<a|---|b> - N*<ab||c>)
      39              : !>                           r
      40              : !>
      41              : !>                      -Z
      42              : !>                 = <a|---|b> + Z*N*<ab||c>
      43              : !>                       r
      44              : !>                   \_______/       \_____/
      45              : !>                       |              |
      46              : !>                    nuclear        coulomb
      47              : !>
      48              : !>      b) Pseudopotential calculation (Goedecker, Teter and Hutter; GTH):
      49              : !>
      50              : !>         <a|V|b> = <a|(V(local) + V(non-local))|b>
      51              : !>
      52              : !>                 = <a|(V(local)|b> + <a|V(non-local))|b>
      53              : !>
      54              : !>         <a|V(local)|b> = <a|-Z(eff)*erf(SQRT(2)*alpha*r)/r +
      55              : !>                             (C1 + C2*(alpha*r)**2 + C3*(alpha*r)**4 +
      56              : !>                              C4*(alpha*r)**6)*exp(-(alpha*r)**2/2))|b>
      57              : !>
      58              : !>         <a|V(non-local)|b> = <a|p(l,i)>*h(i,j)*<p(l,j)|b>
      59              : !> \par Literature
      60              : !>      S. Goedecker, M. Teter and J. Hutter, Phys. Rev. B 54, 1703 (1996)
      61              : !>      C. Hartwigsen, S. Goedecker and J. Hutter, Phys. Rev. B 58, 3641 (1998)
      62              : !>      M. Krack and M. Parrinello, Phys. Chem. Chem. Phys. 2, 2105 (2000)
      63              : !>      S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)
      64              : !> \par History
      65              : !>      - Joost VandeVondele (April 2003) : added LSD forces
      66              : !>      - Non-redundant calculation of the non-local part of the GTH PP
      67              : !>        (22.05.2003,MK)
      68              : !>      - New parallelization scheme (27.06.2003,MK)
      69              : !>      - OpenMP version (07.12.2003,JGH)
      70              : !>      - Binary search loop for VPPNL operators (09.01.2004,JGH,MK)
      71              : !>      - Refactoring of pseudopotential and nuclear attraction integrals (25.02.2009,JGH)
      72              : !>      - General refactoring (01.10.2010,JGH)
      73              : !>      - Refactoring related to the new kinetic energy and overlap routines (07.2014,JGH)
      74              : !>      - k-point functionality (07.2015,JGH)
      75              : !> \author Matthias Krack (14.09.2000,21.03.02)
      76              : ! **************************************************************************************************
      77              : MODULE qs_core_hamiltonian
      78              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      79              :                                               get_atomic_kind_set
      80              :    USE core_ae,                         ONLY: build_core_ae
      81              :    USE core_ppl,                        ONLY: build_core_ppl
      82              :    USE core_ppnl,                       ONLY: build_core_ppnl
      83              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      84              :    USE cp_control_types,                ONLY: dft_control_type
      85              :    USE cp_dbcsr_api,                    ONLY: &
      86              :         dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_iterator_blocks_left, &
      87              :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      88              :         dbcsr_p_type, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric
      89              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      90              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      91              :                                               dbcsr_deallocate_matrix_set
      92              :    USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_matrix_dist,&
      93              :                                               cp_dbcsr_write_sparse_matrix
      94              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      95              :                                               cp_logger_type
      96              :    USE cp_output_handling,              ONLY: cp_p_file,&
      97              :                                               cp_print_key_finished_output,&
      98              :                                               cp_print_key_should_output,&
      99              :                                               cp_print_key_unit_nr
     100              :    USE input_constants,                 ONLY: do_admm_purify_none,&
     101              :                                               do_ppl_analytic,&
     102              :                                               kg_tnadd_atomic,&
     103              :                                               rel_none,&
     104              :                                               rel_trans_atom
     105              :    USE input_section_types,             ONLY: section_vals_val_get
     106              :    USE kg_environment_types,            ONLY: kg_environment_type
     107              :    USE kg_tnadd_mat,                    ONLY: build_tnadd_mat
     108              :    USE kinds,                           ONLY: default_string_length,&
     109              :                                               dp
     110              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
     111              :                                               kpoint_type
     112              :    USE lri_environment_types,           ONLY: lri_environment_type
     113              :    USE message_passing,                 ONLY: mp_para_env_type
     114              :    USE particle_types,                  ONLY: particle_type
     115              :    USE qs_condnum,                      ONLY: overlap_condnum
     116              :    USE qs_environment_types,            ONLY: get_qs_env,&
     117              :                                               qs_environment_type,&
     118              :                                               set_qs_env
     119              :    USE qs_force_types,                  ONLY: qs_force_type
     120              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
     121              :                                               qs_kind_type
     122              :    USE qs_kinetic,                      ONLY: build_kinetic_matrix
     123              :    USE qs_ks_types,                     ONLY: get_ks_env,&
     124              :                                               qs_ks_env_type,&
     125              :                                               set_ks_env
     126              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
     127              :    USE qs_oce_methods,                  ONLY: build_oce_matrices
     128              :    USE qs_oce_types,                    ONLY: allocate_oce_set,&
     129              :                                               create_oce_set,&
     130              :                                               oce_matrix_type
     131              :    USE qs_overlap,                      ONLY: build_overlap_matrix
     132              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
     133              :                                               qs_rho_type
     134              :    USE virial_types,                    ONLY: virial_type
     135              : #include "./base/base_uses.f90"
     136              : 
     137              :    IMPLICIT NONE
     138              : 
     139              :    PRIVATE
     140              : 
     141              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_hamiltonian'
     142              : 
     143              :    PUBLIC :: build_core_hamiltonian_matrix, core_matrices
     144              :    PUBLIC :: dump_info_core_hamiltonian, qs_matrix_h_allocate_imag_from_real
     145              : 
     146              : CONTAINS
     147              : 
     148              : ! **************************************************************************************************
     149              : !> \brief Cosntruction of the QS Core Hamiltonian Matrix
     150              : !> \param qs_env ...
     151              : !> \param calculate_forces ...
     152              : !> \author Creation (11.03.2002,MK)
     153              : !>      Non-redundant calculation of the non-local part of the GTH PP (22.05.2003,MK)
     154              : !>      New parallelization scheme (27.06.2003,MK)
     155              : ! **************************************************************************************************
     156        16455 :    SUBROUTINE build_core_hamiltonian_matrix(qs_env, calculate_forces)
     157              : 
     158              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     159              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     160              : 
     161              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_core_hamiltonian_matrix'
     162              : 
     163              :       INTEGER                                            :: handle, ic, img, iw, nder, nders, &
     164              :                                                             nimages, nkind
     165        16455 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     166              :       LOGICAL                                            :: h_is_complex, norml1, norml2, ofdft, &
     167              :                                                             use_arnoldi, use_virial
     168              :       REAL(KIND=dp)                                      :: eps_filter, eps_fit
     169              :       REAL(KIND=dp), DIMENSION(2)                        :: condnum
     170        16455 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     171              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     172              :       TYPE(cp_logger_type), POINTER                      :: logger
     173              :       TYPE(dbcsr_distribution_type), POINTER             :: dbcsr_dist
     174        16455 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p, matrix_s, matrix_t, &
     175        16455 :                                                             matrix_w
     176              :       TYPE(dft_control_type), POINTER                    :: dft_control
     177              :       TYPE(kg_environment_type), POINTER                 :: kg_env
     178              :       TYPE(kpoint_type), POINTER                         :: kpoints
     179              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     180        16455 :          POINTER                                         :: sab_orb, sap_oce
     181              :       TYPE(oce_matrix_type), POINTER                     :: oce
     182        16455 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     183        16455 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     184        16455 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     185              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     186              :       TYPE(qs_rho_type), POINTER                         :: rho
     187              :       TYPE(virial_type), POINTER                         :: virial
     188              : 
     189        32910 :       IF (calculate_forces) THEN
     190         5855 :          CALL timeset(routineN//"_forces", handle)
     191              :       ELSE
     192        10600 :          CALL timeset(routineN, handle)
     193              :       END IF
     194              : 
     195        16455 :       NULLIFY (logger)
     196        16455 :       logger => cp_get_default_logger()
     197              : 
     198        16455 :       NULLIFY (dft_control)
     199        16455 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     200              : 
     201              :       ! is this a orbital-free method calculation
     202        16455 :       ofdft = dft_control%qs_control%ofgpw
     203              : 
     204        16455 :       nimages = dft_control%nimages
     205        16455 :       IF (ofdft) THEN
     206            0 :          CPASSERT(nimages == 1)
     207              :       END IF
     208              : 
     209        16455 :       nders = 0
     210        16455 :       IF (calculate_forces) THEN
     211         5855 :          nder = 1
     212              :       ELSE
     213        10600 :          IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
     214              :                                         "DFT%PRINT%AO_MATRICES/DERIVATIVES") /= 0) THEN
     215            4 :             nder = 1
     216              :          ELSE
     217        10596 :             nder = 0
     218              :          END IF
     219              :       END IF
     220              : 
     221        16455 :       IF ((cp_print_key_should_output(logger%iter_info, qs_env%input, &
     222              :                                       "DFT%PRINT%AO_MATRICES/OVERLAP") /= 0 .AND. &
     223              :            BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
     224              :                                             "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file))) THEN
     225            4 :          nders = 1
     226              :       END IF
     227              : 
     228              :       ! the delta pulse in the periodic case needs the momentum operator,
     229              :       ! which is equivalent to the derivative of the overlap matrix
     230        16455 :       IF (ASSOCIATED(dft_control%rtp_control)) THEN
     231         1788 :          IF (dft_control%rtp_control%apply_delta_pulse .AND. &
     232              :              dft_control%rtp_control%periodic) THEN
     233          118 :             nders = 1
     234              :          END IF
     235              :       END IF
     236              : 
     237        16455 :       IF (dft_control%tddfpt2_control%enabled) THEN
     238         1410 :          nders = 1
     239         1410 :          IF (dft_control%do_admm) THEN
     240          318 :             IF (dft_control%admm_control%purification_method /= do_admm_purify_none) &
     241              :                CALL cp_abort(__LOCATION__, &
     242            0 :                              "Only purification method NONE is possible with TDDFT at the moment")
     243              :          END IF
     244              :       END IF
     245              : 
     246              :       ! filter for new matrices
     247        16455 :       eps_filter = dft_control%qs_control%eps_filter_matrix
     248              :       !
     249        16455 :       NULLIFY (ks_env)
     250        16455 :       CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
     251        16455 :       NULLIFY (matrix_s, matrix_t)
     252        16455 :       CALL get_qs_env(qs_env=qs_env, kinetic_kp=matrix_t, matrix_s_kp=matrix_s)
     253        16455 :       NULLIFY (sab_orb)
     254        16455 :       CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
     255        16455 :       NULLIFY (rho, force, matrix_p, matrix_w)
     256        16455 :       IF (calculate_forces) THEN
     257         5855 :          CALL get_qs_env(qs_env=qs_env, force=force, matrix_w_kp=matrix_w)
     258         5855 :          CALL get_qs_env(qs_env=qs_env, rho=rho)
     259         5855 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
     260              :          !     *** If LSD, then combine alpha density and beta density to
     261              :          !     *** total density: alpha <- alpha + beta   and
     262              :          !     *** spin density:   beta <- alpha - beta
     263              :          !     (since all things can be computed based on the sum of these matrices anyway)
     264              :          !     (matrix_p is restored at the end of the run, matrix_w is left in its modified state
     265              :          !     (as it should not be needed afterwards)
     266         5855 :          IF (SIZE(matrix_p, 1) == 2) THEN
     267         2398 :             DO img = 1, nimages
     268              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     269         1580 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     270              :                CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
     271         1580 :                               alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
     272              :                CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
     273         2398 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     274              :             END DO
     275              :          END IF
     276              :          ! S matrix
     277              :          CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
     278              :                                    matrix_name="OVERLAP MATRIX", &
     279              :                                    basis_type_a="ORB", &
     280              :                                    basis_type_b="ORB", &
     281              :                                    sab_nl=sab_orb, calculate_forces=.TRUE., &
     282         5855 :                                    matrixkp_p=matrix_w)
     283              :          ! T matrix
     284         5855 :          IF (.NOT. ofdft) &
     285              :             CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, &
     286              :                                       matrix_name="KINETIC ENERGY MATRIX", &
     287              :                                       basis_type="ORB", &
     288              :                                       sab_nl=sab_orb, calculate_forces=.TRUE., &
     289              :                                       matrixkp_p=matrix_p, &
     290         5855 :                                       eps_filter=eps_filter)
     291              :          ! *** If LSD, then recover alpha density and beta density     ***
     292              :          ! *** from the total density (1) and the spin density (2)     ***
     293              :          ! *** The W matrix is neglected, since it will be destroyed   ***
     294              :          ! *** in the calling force routine after leaving this routine ***
     295         5855 :          IF (SIZE(matrix_p, 1) == 2) THEN
     296         2398 :             DO img = 1, nimages
     297              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     298         1580 :                               alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
     299              :                CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
     300         2398 :                               alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
     301              :             END DO
     302              :          END IF
     303              :       ELSE
     304              :          ! S matrix
     305              :          CALL build_overlap_matrix(ks_env, nderivative=nders, matrixkp_s=matrix_s, &
     306              :                                    matrix_name="OVERLAP MATRIX", &
     307              :                                    basis_type_a="ORB", &
     308              :                                    basis_type_b="ORB", &
     309        10600 :                                    sab_nl=sab_orb)
     310              :          ! T matrix
     311        10600 :          IF (.NOT. ofdft) &
     312              :             CALL build_kinetic_matrix(ks_env, matrixkp_t=matrix_t, &
     313              :                                       matrix_name="KINETIC ENERGY MATRIX", &
     314              :                                       basis_type="ORB", &
     315              :                                       sab_nl=sab_orb, &
     316        10600 :                                       eps_filter=eps_filter)
     317              :       END IF
     318              : 
     319              :       ! (Re-)allocate H matrix based on overlap matrix
     320        16455 :       CALL get_ks_env(ks_env, complex_ks=h_is_complex)
     321        16455 :       CALL qs_matrix_h_allocate(qs_env, matrix_s(1, 1)%matrix, is_complex=h_is_complex)
     322              : 
     323        16455 :       NULLIFY (matrix_h)
     324        16455 :       CALL get_qs_env(qs_env, matrix_h_kp=matrix_h)
     325              : 
     326        16455 :       IF (.NOT. ofdft) THEN
     327        66234 :          DO img = 1, nimages
     328              :             CALL dbcsr_copy(matrix_h(1, img)%matrix, matrix_t(1, img)%matrix, &
     329        66234 :                             keep_sparsity=.TRUE., name="CORE HAMILTONIAN MATRIX")
     330              :          END DO
     331              :       END IF
     332              : 
     333        16455 :       NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
     334              :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
     335        16455 :                       particle_set=particle_set)
     336              : 
     337        16455 :       IF (.NOT. ofdft) THEN
     338              :          ! relativistic atomic correction to kinetic energy
     339        16455 :          IF (qs_env%rel_control%rel_method /= rel_none) THEN
     340           58 :             IF (qs_env%rel_control%rel_transformation == rel_trans_atom) THEN
     341           58 :                IF (nimages == 1) THEN
     342           58 :                   ic = 1
     343              :                ELSE
     344            4 :                   CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
     345            4 :                   CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     346            4 :                   ic = cell_to_index(0, 0, 0)
     347              :                END IF
     348           58 :                CALL build_atomic_relmat(matrix_h(1, 1)%matrix, atomic_kind_set, qs_kind_set)
     349              :             ELSE
     350            0 :                CPABORT("Relativistic corrections of this type are currently not implemented")
     351              :             END IF
     352              :          END IF
     353              :       END IF
     354              : 
     355              :       ! *** core and pseudopotentials
     356        16455 :       CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
     357              : 
     358              :       ! *** GAPW one-center-expansion (oce) matrices
     359        16455 :       NULLIFY (sap_oce)
     360        16455 :       CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)
     361        16455 :       NULLIFY (oce)
     362        16455 :       IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
     363         2128 :          CALL get_qs_env(qs_env=qs_env, oce=oce)
     364         2128 :          CALL create_oce_set(oce)
     365         2128 :          nkind = SIZE(atomic_kind_set)
     366         2128 :          CALL allocate_oce_set(oce, nkind)
     367         2128 :          eps_fit = dft_control%qs_control%gapw_control%eps_fit
     368         2128 :          IF (ASSOCIATED(sap_oce)) &
     369              :             CALL build_oce_matrices(oce%intac, calculate_forces, nder, qs_kind_set, particle_set, &
     370         2094 :                                     sap_oce, eps_fit)
     371              :       END IF
     372              : 
     373              :       ! *** KG atomic potentials for nonadditive kinetic energy
     374        16455 :       IF (dft_control%qs_control%do_kg) THEN
     375          182 :          IF (qs_env%kg_env%tnadd_method == kg_tnadd_atomic) THEN
     376           30 :             CALL get_qs_env(qs_env=qs_env, kg_env=kg_env, virial=virial, dbcsr_dist=dbcsr_dist)
     377           30 :             use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     378              :             CALL build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
     379           30 :                                  qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
     380              :          END IF
     381              :       END IF
     382              : 
     383              :       ! *** Put the core Hamiltonian matrix in the QS environment ***
     384        16455 :       CALL set_qs_env(qs_env, oce=oce)
     385        16455 :       CALL set_ks_env(ks_env, matrix_s_kp=matrix_s, kinetic_kp=matrix_t, matrix_h_kp=matrix_h)
     386              : 
     387              :       ! *** Print matrices if requested
     388        16455 :       CALL dump_info_core_hamiltonian(qs_env, calculate_forces)
     389              : 
     390              :       ! *** Overlap condition number
     391        16455 :       IF (.NOT. calculate_forces) THEN
     392        10600 :          IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
     393              :                                         "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN
     394              :             iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
     395           36 :                                       extension=".Log")
     396           36 :             CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
     397           36 :             CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
     398           36 :             CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
     399           36 :             CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
     400           36 :             CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
     401              :          END IF
     402              :       END IF
     403              : 
     404        16455 :       CALL timestop(handle)
     405              : 
     406        16455 :    END SUBROUTINE build_core_hamiltonian_matrix
     407              : 
     408              : ! **************************************************************************************************
     409              : !> \brief ...
     410              : !> \param qs_env ...
     411              : !> \param matrix_h ...
     412              : !> \param matrix_p ...
     413              : !> \param calculate_forces ...
     414              : !> \param nder ...
     415              : !> \param atcore ...
     416              : ! **************************************************************************************************
     417        16521 :    SUBROUTINE core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder, atcore)
     418              : 
     419              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     420              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p
     421              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     422              :       INTEGER, INTENT(IN)                                :: nder
     423              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL              :: atcore
     424              : 
     425              :       INTEGER                                            :: natom, nimages
     426        16521 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     427              :       LOGICAL                                            :: all_present, my_gt_nl, ppl_present, &
     428              :                                                             ppnl_present, use_virial
     429              :       REAL(KIND=dp)                                      :: eps_ppnl
     430        16521 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     431              :       TYPE(dft_control_type), POINTER                    :: dft_control
     432              :       TYPE(kpoint_type), POINTER                         :: kpoints
     433              :       TYPE(lri_environment_type), POINTER                :: lri_env
     434              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     435        16521 :          POINTER                                         :: sab_orb, sac_ae, sac_ppl, sap_ppnl
     436        16521 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     437        16521 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     438        16521 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     439              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     440              :       TYPE(virial_type), POINTER                         :: virial
     441              : 
     442        16521 :       NULLIFY (dft_control)
     443        16521 :       CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, dft_control=dft_control, natom=natom)
     444        16521 :       nimages = dft_control%nimages
     445        16521 :       IF (PRESENT(atcore)) THEN
     446           66 :          CPASSERT(SIZE(atcore) >= natom)
     447              :       END IF
     448              : 
     449              :       ! check whether a gauge transformed version of the non-local potential part has to be used
     450        16521 :       my_gt_nl = .FALSE.
     451        16521 :       IF (qs_env%run_rtp) THEN
     452         1238 :          CPASSERT(ASSOCIATED(dft_control%rtp_control))
     453         1238 :          IF (dft_control%rtp_control%velocity_gauge) THEN
     454           54 :             my_gt_nl = dft_control%rtp_control%nl_gauge_transform
     455              :          END IF
     456              :       END IF
     457              : 
     458              :       ! prepare for k-points
     459        16521 :       NULLIFY (cell_to_index)
     460        16521 :       IF (nimages > 1) THEN
     461          310 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
     462          310 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     463          310 :          dft_control%qs_control%do_ppl_method = do_ppl_analytic
     464              :       END IF
     465              :       ! force analytic ppl calculation for GAPW methods
     466        16521 :       IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
     467         2148 :          dft_control%qs_control%do_ppl_method = do_ppl_analytic
     468              :       END IF
     469              : 
     470              :       ! force
     471        16521 :       NULLIFY (force)
     472        16521 :       IF (calculate_forces) CALL get_qs_env(qs_env=qs_env, force=force)
     473              :       ! virial
     474        16521 :       CALL get_qs_env(qs_env=qs_env, virial=virial)
     475        16521 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     476              : 
     477        16521 :       NULLIFY (qs_kind_set, atomic_kind_set, particle_set)
     478              :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, atomic_kind_set=atomic_kind_set, &
     479        16521 :                       particle_set=particle_set)
     480              : 
     481        16521 :       NULLIFY (sab_orb, sac_ae, sac_ppl, sap_ppnl)
     482              :       CALL get_qs_env(qs_env=qs_env, &
     483              :                       sab_orb=sab_orb, &
     484              :                       sac_ae=sac_ae, &
     485              :                       sac_ppl=sac_ppl, &
     486        16521 :                       sap_ppnl=sap_ppnl)
     487              : 
     488              :       ! *** compute the nuclear attraction contribution to the core hamiltonian ***
     489        16521 :       all_present = ASSOCIATED(sac_ae)
     490        16521 :       IF (all_present) THEN
     491              :          CALL build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
     492              :                             qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
     493         2056 :                             nimages, cell_to_index, "ORB", atcore=atcore)
     494              :       END IF
     495              : 
     496              :       ! *** compute the ppl contribution to the core hamiltonian ***
     497        16521 :       ppl_present = ASSOCIATED(sac_ppl)
     498        16521 :       IF (ppl_present) THEN
     499        15563 :          IF (dft_control%qs_control%do_ppl_method == do_ppl_analytic) THEN
     500        15539 :             IF (dft_control%qs_control%lrigpw) THEN
     501           80 :                CALL get_qs_env(qs_env, lri_env=lri_env)
     502           80 :                IF (lri_env%ppl_ri) THEN
     503            4 :                   IF (lri_env%exact_1c_terms) THEN
     504            0 :                      CPABORT("not implemented")
     505              :                   END IF
     506              :                ELSE
     507              :                   CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
     508              :                                       qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
     509          152 :                                       nimages, cell_to_index, "ORB", atcore=atcore)
     510              :                END IF
     511              :             ELSE
     512              :                CALL build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
     513              :                                    qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
     514        30856 :                                    nimages, cell_to_index, "ORB", atcore=atcore)
     515              :             END IF
     516              :          END IF
     517              :       END IF
     518              : 
     519              :       ! *** compute the ppnl contribution to the core hamiltonian ***
     520        16521 :       eps_ppnl = dft_control%qs_control%eps_ppnl
     521        16521 :       ppnl_present = ASSOCIATED(sap_ppnl)
     522        16521 :       IF (ppnl_present) THEN
     523        12684 :          IF (.NOT. my_gt_nl) THEN
     524              :             CALL build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, &
     525              :                                  qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, &
     526        25234 :                                  nimages, cell_to_index, "ORB", atcore=atcore)
     527              :          END IF
     528              :       END IF
     529              : 
     530        16521 :    END SUBROUTINE core_matrices
     531              : 
     532              : ! **************************************************************************************************
     533              : !> \brief Adds atomic blocks of relativistic correction for the kinetic energy
     534              : !> \param matrix_h ...
     535              : !> \param atomic_kind_set ...
     536              : !> \param qs_kind_set ...
     537              : ! **************************************************************************************************
     538           58 :    SUBROUTINE build_atomic_relmat(matrix_h, atomic_kind_set, qs_kind_set)
     539              :       TYPE(dbcsr_type), POINTER                          :: matrix_h
     540              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     541              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     542              : 
     543              :       INTEGER                                            :: iatom, ikind, jatom
     544           58 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
     545           58 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: hblock, reltmat
     546              :       TYPE(dbcsr_iterator_type)                          :: iter
     547              : 
     548           58 :       CALL get_atomic_kind_set(atomic_kind_set, kind_of=kind_of)
     549              : 
     550           58 :       CALL dbcsr_iterator_start(iter, matrix_h)
     551          221 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     552          163 :          CALL dbcsr_iterator_next_block(iter, iatom, jatom, hblock)
     553          221 :          IF (iatom == jatom) THEN
     554           83 :             ikind = kind_of(iatom)
     555           83 :             CALL get_qs_kind(qs_kind_set(ikind), reltmat=reltmat)
     556       192766 :             IF (ASSOCIATED(reltmat)) hblock = hblock + reltmat
     557              :          END IF
     558              :       END DO
     559           58 :       CALL dbcsr_iterator_stop(iter)
     560              : 
     561          116 :    END SUBROUTINE build_atomic_relmat
     562              : 
     563              : ! **************************************************************************************************
     564              : !> \brief Possibly prints matrices after the construction of the Core
     565              : !>     Hamiltonian Matrix
     566              : !> \param qs_env ...
     567              : !> \param calculate_forces ...
     568              : ! **************************************************************************************************
     569        32910 :    SUBROUTINE dump_info_core_hamiltonian(qs_env, calculate_forces)
     570              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     571              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     572              : 
     573              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'dump_info_core_hamiltonian'
     574              : 
     575              :       INTEGER                                            :: after, handle, i, ic, iw, output_unit
     576              :       LOGICAL                                            :: omit_headers
     577              :       TYPE(cp_logger_type), POINTER                      :: logger
     578        16455 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_v
     579        16455 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrixkp_h, matrixkp_s, matrixkp_t
     580              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     581              : 
     582        16455 :       CALL timeset(routineN, handle)
     583              : 
     584        16455 :       NULLIFY (logger, matrix_v, para_env)
     585        16455 :       logger => cp_get_default_logger()
     586        16455 :       CALL get_qs_env(qs_env, para_env=para_env)
     587              : 
     588              :       ! Print the distribution of the overlap matrix blocks
     589              :       ! this duplicates causes duplicate printing at the force calc
     590        16455 :       IF (.NOT. calculate_forces) THEN
     591        10600 :          IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     592              :                                               qs_env%input, "PRINT%DISTRIBUTION"), cp_p_file)) THEN
     593              :             output_unit = cp_print_key_unit_nr(logger, qs_env%input, "PRINT%DISTRIBUTION", &
     594           92 :                                                extension=".distribution")
     595           92 :             CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
     596           92 :             CALL cp_dbcsr_write_matrix_dist(matrixkp_s(1, 1)%matrix, output_unit, para_env)
     597           92 :             CALL cp_print_key_finished_output(output_unit, logger, qs_env%input, "PRINT%DISTRIBUTION")
     598              :          END IF
     599              :       END IF
     600              : 
     601        16455 :       CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
     602              :       ! Print the overlap integral matrix, if requested
     603        16455 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     604              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
     605              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
     606            4 :                                    extension=".Log")
     607            4 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
     608            4 :          after = MIN(MAX(after, 1), 16)
     609            4 :          CALL get_qs_env(qs_env, matrix_s_kp=matrixkp_s)
     610            4 :          IF (ASSOCIATED(matrixkp_s)) THEN
     611            8 :             DO ic = 1, SIZE(matrixkp_s, 2)
     612              :                CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(1, ic)%matrix, 4, after, qs_env, para_env, &
     613            8 :                                                  output_unit=iw, omit_headers=omit_headers)
     614              :             END DO
     615            4 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, qs_env%input, &
     616              :                                                  "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
     617            8 :                DO ic = 1, SIZE(matrixkp_s, 2)
     618           20 :                   DO i = 2, SIZE(matrixkp_s, 1)
     619              :                      CALL cp_dbcsr_write_sparse_matrix(matrixkp_s(i, ic)%matrix, 4, after, qs_env, para_env, &
     620           16 :                                                        output_unit=iw, omit_headers=omit_headers)
     621              :                   END DO
     622              :                END DO
     623              :             END IF
     624              :          END IF
     625              :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
     626            4 :                                            "DFT%PRINT%AO_MATRICES/OVERLAP")
     627              :       END IF
     628              : 
     629              :       ! Print the kinetic energy integral matrix, if requested
     630        16455 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     631              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY"), cp_p_file)) THEN
     632              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY", &
     633           48 :                                    extension=".Log")
     634           48 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
     635           48 :          after = MIN(MAX(after, 1), 16)
     636           48 :          CALL get_qs_env(qs_env, kinetic_kp=matrixkp_t)
     637           48 :          IF (ASSOCIATED(matrixkp_t)) THEN
     638           96 :             DO ic = 1, SIZE(matrixkp_t, 2)
     639              :                CALL cp_dbcsr_write_sparse_matrix(matrixkp_t(1, ic)%matrix, 4, after, qs_env, para_env, &
     640           96 :                                                  output_unit=iw, omit_headers=omit_headers)
     641              :             END DO
     642              :          END IF
     643              :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
     644           48 :                                            "DFT%PRINT%AO_MATRICES/KINETIC_ENERGY")
     645              :       END IF
     646              : 
     647              :       ! Print the potential energy matrix, if requested
     648        16455 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     649              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY"), cp_p_file)) THEN
     650              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY", &
     651           48 :                                    extension=".Log")
     652           48 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
     653           48 :          after = MIN(MAX(after, 1), 16)
     654           48 :          CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h, kinetic_kp=matrixkp_t)
     655           48 :          IF (ASSOCIATED(matrixkp_h)) THEN
     656           48 :             IF (SIZE(matrixkp_h, 2) == 1) THEN
     657           48 :                CALL dbcsr_allocate_matrix_set(matrix_v, 1)
     658           48 :                ALLOCATE (matrix_v(1)%matrix)
     659           48 :                CALL dbcsr_copy(matrix_v(1)%matrix, matrixkp_h(1, 1)%matrix, name="POTENTIAL ENERGY MATRIX")
     660              :                CALL dbcsr_add(matrix_v(1)%matrix, matrixkp_t(1, 1)%matrix, &
     661           48 :                               alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
     662              :                CALL cp_dbcsr_write_sparse_matrix(matrix_v(1)%matrix, 4, after, qs_env, &
     663           48 :                                                  para_env, output_unit=iw, omit_headers=omit_headers)
     664           48 :                CALL dbcsr_deallocate_matrix_set(matrix_v)
     665              :             ELSE
     666            0 :                CPWARN("Printing of potential energy matrix not implemented for k-points")
     667              :             END IF
     668              :          END IF
     669              :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
     670           48 :                                            "DFT%PRINT%AO_MATRICES/POTENTIAL_ENERGY")
     671              :       END IF
     672              : 
     673              :       ! Print the core Hamiltonian matrix, if requested
     674        16455 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
     675              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
     676              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
     677           48 :                                    extension=".Log")
     678           48 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
     679           48 :          after = MIN(MAX(after, 1), 16)
     680           48 :          CALL get_qs_env(qs_env, matrix_h_kp=matrixkp_h)
     681           48 :          IF (ASSOCIATED(matrixkp_h)) THEN
     682           96 :             DO ic = 1, SIZE(matrixkp_h, 2)
     683              :                CALL cp_dbcsr_write_sparse_matrix(matrixkp_h(1, ic)%matrix, 4, after, qs_env, para_env, &
     684           96 :                                                  output_unit=iw, omit_headers=omit_headers)
     685              :             END DO
     686              :          END IF
     687              :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, &
     688           48 :                                            "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
     689              :       END IF
     690              : 
     691        16455 :       CALL timestop(handle)
     692              : 
     693        16455 :    END SUBROUTINE dump_info_core_hamiltonian
     694              : 
     695              : ! **************************************************************************************************
     696              : !> \brief (Re-)allocate matrix_h based on the template (typically the overlap matrix)
     697              : !> \param qs_env ...
     698              : !> \param template ...
     699              : !> \param is_complex ...
     700              : ! **************************************************************************************************
     701        16455 :    SUBROUTINE qs_matrix_h_allocate(qs_env, template, is_complex)
     702              :       TYPE(qs_environment_type)                          :: qs_env
     703              :       TYPE(dbcsr_type), INTENT(in)                       :: template
     704              :       LOGICAL, INTENT(in)                                :: is_complex
     705              : 
     706              :       CHARACTER(LEN=default_string_length)               :: headline
     707              :       INTEGER                                            :: img, nimages
     708        16455 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_h_im
     709              :       TYPE(dft_control_type), POINTER                    :: dft_control
     710              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     711        16455 :          POINTER                                         :: sab_orb
     712              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     713              : 
     714        16455 :       NULLIFY (matrix_h, matrix_h_im, sab_orb, dft_control, ks_env)
     715              :       CALL get_qs_env(qs_env=qs_env, &
     716              :                       matrix_h_kp=matrix_h, &
     717              :                       matrix_h_im_kp=matrix_h_im, &
     718              :                       sab_orb=sab_orb, &
     719              :                       dft_control=dft_control, &
     720        16455 :                       ks_env=ks_env)
     721              : 
     722        16455 :       nimages = dft_control%nimages
     723        16455 :       CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimages)
     724        16455 :       headline = "CORE HAMILTONIAN MATRIX"
     725        66234 :       DO img = 1, nimages
     726        49779 :          ALLOCATE (matrix_h(1, img)%matrix)
     727        49779 :          CALL dbcsr_create(matrix_h(1, img)%matrix, name=TRIM(headline), template=template)
     728        49779 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
     729        66234 :          CALL dbcsr_set(matrix_h(1, img)%matrix, 0.0_dp)
     730              :       END DO
     731        16455 :       CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
     732              : 
     733        16455 :       IF (is_complex) THEN
     734          352 :          headline = "IMAGINARY PART OF CORE HAMILTONIAN MATRIX"
     735          352 :          CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
     736          704 :          DO img = 1, nimages
     737          352 :             ALLOCATE (matrix_h_im(1, img)%matrix)
     738              :             CALL dbcsr_create(matrix_h_im(1, img)%matrix, name=TRIM(headline), template=template, &
     739          352 :                               matrix_type=dbcsr_type_antisymmetric)
     740          352 :             CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, img)%matrix, sab_orb)
     741          704 :             CALL dbcsr_set(matrix_h_im(1, img)%matrix, 0.0_dp)
     742              :          END DO
     743          352 :          CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
     744              :       END IF
     745              : 
     746        16455 :    END SUBROUTINE qs_matrix_h_allocate
     747              : 
     748              : ! **************************************************************************************************
     749              : !> \brief (Re-)allocates matrix_h_im from matrix_h
     750              : !> \param qs_env ...
     751              : ! **************************************************************************************************
     752            8 :    SUBROUTINE qs_matrix_h_allocate_imag_from_real(qs_env)
     753              :       TYPE(qs_environment_type)                          :: qs_env
     754              : 
     755              :       CHARACTER(LEN=default_string_length)               :: headline
     756              :       INTEGER                                            :: image, nimages
     757            8 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_h_im
     758              :       TYPE(dbcsr_type), POINTER                          :: template
     759              :       TYPE(dft_control_type), POINTER                    :: dft_control
     760              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     761            8 :          POINTER                                         :: sab_orb
     762              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     763              : 
     764            8 :       NULLIFY (matrix_h_im, matrix_h, dft_control, template, sab_orb, ks_env)
     765              : 
     766              :       CALL get_qs_env(qs_env, &
     767              :                       matrix_h_im_kp=matrix_h_im, &
     768              :                       matrix_h_kp=matrix_h, &
     769              :                       dft_control=dft_control, &
     770              :                       sab_orb=sab_orb, &
     771            8 :                       ks_env=ks_env)
     772              : 
     773            8 :       nimages = dft_control%nimages
     774              : 
     775            8 :       CPASSERT(nimages .EQ. SIZE(matrix_h, 2))
     776              : 
     777            8 :       CALL dbcsr_allocate_matrix_set(matrix_h_im, 1, nimages)
     778              : 
     779           16 :       DO image = 1, nimages
     780            8 :          headline = "IMAGINARY CORE HAMILTONIAN MATRIX"
     781            8 :          ALLOCATE (matrix_h_im(1, image)%matrix)
     782            8 :          template => matrix_h(1, image)%matrix ! base on real part, but anti-symmetric
     783              :          CALL dbcsr_create(matrix=matrix_h_im(1, image)%matrix, template=template, &
     784            8 :                            name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
     785            8 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_h_im(1, image)%matrix, sab_orb)
     786           16 :          CALL dbcsr_set(matrix_h_im(1, image)%matrix, 0.0_dp)
     787              :       END DO
     788            8 :       CALL set_ks_env(ks_env, matrix_h_im_kp=matrix_h_im)
     789              : 
     790            8 :    END SUBROUTINE qs_matrix_h_allocate_imag_from_real
     791              : 
     792              : END MODULE qs_core_hamiltonian
        

Generated by: LCOV version 2.0-1