LCOV - code coverage report
Current view: top level - src - xtb_matrices.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.9 % 746 723
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 4 4

            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 Overlap and Hamiltonian matrices in xTB
      10              : !>        Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov
      11              : !>                   JCTC 13, 1989-2009, (2017)
      12              : !>                   DOI: 10.1021/acs.jctc.7b00118
      13              : !> \author JGH
      14              : ! **************************************************************************************************
      15              : MODULE xtb_matrices
      16              :    USE ai_contraction,                  ONLY: block_add,&
      17              :                                               contraction
      18              :    USE ai_overlap,                      ONLY: overlap_ab
      19              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      20              :                                               get_atomic_kind_set
      21              :    USE atprop_types,                    ONLY: atprop_array_init,&
      22              :                                               atprop_type
      23              :    USE basis_set_types,                 ONLY: gto_basis_set_p_type,&
      24              :                                               gto_basis_set_type
      25              :    USE block_p_types,                   ONLY: block_p_type
      26              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      27              :    USE cp_control_types,                ONLY: dft_control_type,&
      28              :                                               xtb_control_type
      29              :    USE cp_dbcsr_api,                    ONLY: dbcsr_add,&
      30              :                                               dbcsr_create,&
      31              :                                               dbcsr_finalize,&
      32              :                                               dbcsr_get_block_p,&
      33              :                                               dbcsr_p_type
      34              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      35              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set
      36              :    USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_sparse_matrix
      37              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      38              :                                               cp_logger_type
      39              :    USE cp_output_handling,              ONLY: cp_p_file,&
      40              :                                               cp_print_key_finished_output,&
      41              :                                               cp_print_key_should_output,&
      42              :                                               cp_print_key_unit_nr
      43              :    USE eeq_input,                       ONLY: eeq_solver_type
      44              :    USE input_constants,                 ONLY: vdw_pairpot_dftd4
      45              :    USE input_section_types,             ONLY: section_vals_val_get
      46              :    USE kinds,                           ONLY: dp
      47              :    USE kpoint_types,                    ONLY: get_kpoint_info,&
      48              :                                               kpoint_type
      49              :    USE message_passing,                 ONLY: mp_para_env_type
      50              :    USE orbital_pointers,                ONLY: ncoset
      51              :    USE particle_types,                  ONLY: particle_type
      52              :    USE qs_condnum,                      ONLY: overlap_condnum
      53              :    USE qs_dispersion_cnum,              ONLY: cnumber_init,&
      54              :                                               cnumber_release,&
      55              :                                               dcnum_type
      56              :    USE qs_dispersion_pairpot,           ONLY: calculate_dispersion_pairpot
      57              :    USE qs_dispersion_types,             ONLY: qs_dispersion_type
      58              :    USE qs_energy_types,                 ONLY: qs_energy_type
      59              :    USE qs_environment_types,            ONLY: get_qs_env,&
      60              :                                               qs_environment_type,&
      61              :                                               set_qs_env
      62              :    USE qs_force_types,                  ONLY: qs_force_type
      63              :    USE qs_integral_utils,               ONLY: basis_set_list_setup,&
      64              :                                               get_memory_usage
      65              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      66              :                                               qs_kind_type
      67              :    USE qs_ks_types,                     ONLY: get_ks_env,&
      68              :                                               qs_ks_env_type,&
      69              :                                               set_ks_env
      70              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      71              :                                               neighbor_list_iterate,&
      72              :                                               neighbor_list_iterator_create,&
      73              :                                               neighbor_list_iterator_p_type,&
      74              :                                               neighbor_list_iterator_release,&
      75              :                                               neighbor_list_set_p_type
      76              :    USE qs_overlap,                      ONLY: create_sab_matrix
      77              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      78              :                                               qs_rho_type
      79              :    USE virial_methods,                  ONLY: virial_pair_force
      80              :    USE virial_types,                    ONLY: virial_type
      81              :    USE xtb_eeq,                         ONLY: xtb_eeq_calculation,&
      82              :                                               xtb_eeq_forces
      83              :    USE xtb_hcore,                       ONLY: gfn0_huckel,&
      84              :                                               gfn0_kpair,&
      85              :                                               gfn1_huckel,&
      86              :                                               gfn1_kpair
      87              :    USE xtb_potentials,                  ONLY: nonbonded_correction,&
      88              :                                               repulsive_potential,&
      89              :                                               srb_potential,&
      90              :                                               xb_interaction
      91              :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
      92              :                                               xtb_atom_type
      93              : #include "./base/base_uses.f90"
      94              : 
      95              :    IMPLICIT NONE
      96              : 
      97              :    PRIVATE
      98              : 
      99              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_matrices'
     100              : 
     101              :    PUBLIC :: build_xtb_matrices
     102              : 
     103              : CONTAINS
     104              : 
     105              : ! **************************************************************************************************
     106              : !> \brief ...
     107              : !> \param qs_env ...
     108              : !> \param calculate_forces ...
     109              : ! **************************************************************************************************
     110         5224 :    SUBROUTINE build_xtb_matrices(qs_env, calculate_forces)
     111              : 
     112              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     113              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     114              : 
     115              :       INTEGER                                            :: gfn_type
     116              :       TYPE(dft_control_type), POINTER                    :: dft_control
     117              : 
     118         5224 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     119         5224 :       gfn_type = dft_control%qs_control%xtb_control%gfn_type
     120              : 
     121         1646 :       SELECT CASE (gfn_type)
     122              :       CASE (0)
     123         1646 :          CALL build_gfn0_xtb_matrices(qs_env, calculate_forces)
     124              :       CASE (1)
     125         3578 :          CALL build_gfn1_xtb_matrices(qs_env, calculate_forces)
     126              :       CASE (2)
     127            0 :          CPABORT("gfn_type = 2 not yet available")
     128              :       CASE DEFAULT
     129         5224 :          CPABORT("Unknown gfn_type")
     130              :       END SELECT
     131              : 
     132         5224 :    END SUBROUTINE build_xtb_matrices
     133              : 
     134              : ! **************************************************************************************************
     135              : !> \brief ...
     136              : !> \param qs_env ...
     137              : !> \param calculate_forces ...
     138              : ! **************************************************************************************************
     139         1646 :    SUBROUTINE build_gfn0_xtb_matrices(qs_env, calculate_forces)
     140              : 
     141              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     142              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     143              : 
     144              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn0_xtb_matrices'
     145              : 
     146              :       INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
     147              :          j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
     148              :          natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
     149              :          nsetb, sgfa, sgfb, za, zb
     150         3292 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     151              :       INTEGER, DIMENSION(25)                             :: laoa, laob, naoa, naob
     152              :       INTEGER, DIMENSION(3)                              :: cell
     153         1646 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     154         1646 :                                                             npgfb, nsgfa, nsgfb
     155         1646 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     156         1646 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     157              :       LOGICAL                                            :: defined, diagblock, do_nonbonded, found, &
     158              :                                                             use_virial
     159              :       REAL(KIND=dp) :: dfp, dhij, dr, drk, drx, eeq_energy, ef_energy, enonbonded, enscale, erep, &
     160              :          esrb, etaa, etab, f0, f1, f2, fhua, fhub, fhud, foab, fqa, fqb, hij, kf, qlambda, rcova, &
     161              :          rcovab, rcovb, rrab
     162         3292 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: charges, cnumbers, dcharges, qlagrange
     163         3292 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dfblock, dhuckel, dqhuckel, huckel, owork
     164         1646 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: oint, sint
     165         1646 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: kijab
     166              :       REAL(KIND=dp), DIMENSION(3)                        :: fdik, fdika, fdikb, force_ab, rij, rik
     167              :       REAL(KIND=dp), DIMENSION(5)                        :: dpia, dpib, hena, henb, kpolya, kpolyb, &
     168              :                                                             pia, pib
     169         1646 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: eeq_q, set_radius_a, set_radius_b
     170         1646 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: fblock, pblock, rpgfa, rpgfb, sblock, &
     171         1646 :                                                             scon_a, scon_b, wblock, zeta, zetb
     172         1646 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     173              :       TYPE(atprop_type), POINTER                         :: atprop
     174         6584 :       TYPE(block_p_type), DIMENSION(2:4)                 :: dsblocks
     175              :       TYPE(cp_logger_type), POINTER                      :: logger
     176         1646 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p, matrix_s, matrix_w
     177         1646 :       TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:)        :: dcnum
     178              :       TYPE(dft_control_type), POINTER                    :: dft_control
     179              :       TYPE(eeq_solver_type)                              :: eeq_sparam
     180         1646 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     181              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     182              :       TYPE(kpoint_type), POINTER                         :: kpoints
     183              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     184              :       TYPE(neighbor_list_iterator_p_type), &
     185         1646 :          DIMENSION(:), POINTER                           :: nl_iterator
     186              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     187         1646 :          POINTER                                         :: sab_orb, sab_xtb_nonbond
     188         1646 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     189              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     190              :       TYPE(qs_energy_type), POINTER                      :: energy
     191         1646 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     192         1646 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     193              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     194              :       TYPE(qs_rho_type), POINTER                         :: rho
     195              :       TYPE(virial_type), POINTER                         :: virial
     196              :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b
     197              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
     198              : 
     199         1646 :       CALL timeset(routineN, handle)
     200              : 
     201         1646 :       NULLIFY (logger, virial, atprop)
     202         1646 :       logger => cp_get_default_logger()
     203              : 
     204         1646 :       NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
     205         1646 :                qs_kind_set, sab_orb, ks_env)
     206              :       CALL get_qs_env(qs_env=qs_env, &
     207              :                       ks_env=ks_env, &
     208              :                       energy=energy, &
     209              :                       atomic_kind_set=atomic_kind_set, &
     210              :                       qs_kind_set=qs_kind_set, &
     211              :                       matrix_h_kp=matrix_h, &
     212              :                       matrix_s_kp=matrix_s, &
     213              :                       para_env=para_env, &
     214              :                       atprop=atprop, &
     215              :                       dft_control=dft_control, &
     216         1646 :                       sab_orb=sab_orb)
     217              : 
     218         1646 :       nkind = SIZE(atomic_kind_set)
     219         1646 :       xtb_control => dft_control%qs_control%xtb_control
     220         1646 :       do_nonbonded = xtb_control%do_nonbonded
     221         1646 :       nimg = dft_control%nimages
     222         1646 :       nderivatives = 0
     223         1646 :       IF (calculate_forces) nderivatives = 1
     224         1646 :       IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
     225         1646 :       maxder = ncoset(nderivatives)
     226              : 
     227         1646 :       NULLIFY (particle_set)
     228         1646 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
     229         1646 :       natom = SIZE(particle_set)
     230              :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     231         1646 :                                atom_of_kind=atom_of_kind, kind_of=kind_of)
     232              : 
     233         1646 :       IF (calculate_forces) THEN
     234           42 :          NULLIFY (rho, force, matrix_w)
     235              :          CALL get_qs_env(qs_env=qs_env, &
     236              :                          rho=rho, matrix_w_kp=matrix_w, &
     237           42 :                          virial=virial, force=force)
     238           42 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
     239              : 
     240           42 :          IF (SIZE(matrix_p, 1) == 2) THEN
     241            8 :             DO img = 1, nimg
     242              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     243            4 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     244              :                CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
     245            8 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     246              :             END DO
     247              :          END IF
     248           74 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     249              :       END IF
     250              :       ! atomic energy decomposition
     251         1646 :       IF (atprop%energy) THEN
     252            0 :          CALL atprop_array_init(atprop%atecc, natom)
     253              :       END IF
     254              : 
     255         1646 :       NULLIFY (cell_to_index)
     256         1646 :       IF (nimg > 1) THEN
     257          142 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
     258          142 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     259              :       END IF
     260              : 
     261              :       ! set up basis set lists
     262         9022 :       ALLOCATE (basis_set_list(nkind))
     263         1646 :       CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
     264              : 
     265              :       ! allocate overlap matrix
     266         1646 :       CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
     267              :       CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
     268         1646 :                              sab_orb, .TRUE.)
     269         1646 :       CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
     270              : 
     271              :       ! initialize H matrix
     272         1646 :       CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
     273        30552 :       DO img = 1, nimg
     274        28906 :          ALLOCATE (matrix_h(1, img)%matrix)
     275              :          CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
     276        28906 :                            name="HAMILTONIAN MATRIX")
     277        30552 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
     278              :       END DO
     279         1646 :       CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
     280              : 
     281              :       ! Calculate coordination numbers
     282              :       ! needed for effective atomic energy levels
     283              :       ! code taken from D3 dispersion energy
     284         1646 :       CALL cnumber_init(qs_env, cnumbers, dcnum, 2, calculate_forces)
     285              : 
     286         4938 :       ALLOCATE (charges(natom))
     287         9396 :       charges = 0.0_dp
     288         1646 :       CALL xtb_eeq_calculation(qs_env, charges, cnumbers, eeq_sparam, eeq_energy, ef_energy, qlambda)
     289         1646 :       IF (calculate_forces) THEN
     290           84 :          ALLOCATE (dcharges(natom))
     291          242 :          dcharges = qlambda/REAL(para_env%num_pe, KIND=dp)
     292              :       END IF
     293         1646 :       energy%eeq = eeq_energy
     294         1646 :       energy%efield = ef_energy
     295              : 
     296         1646 :       CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
     297              :       ! prepare charges (needed for D4)
     298         1646 :       IF (dispersion_env%pp_type == vdw_pairpot_dftd4) THEN
     299          858 :          dispersion_env%ext_charges = .TRUE.
     300          858 :          IF (ASSOCIATED(dispersion_env%charges)) DEALLOCATE (dispersion_env%charges)
     301         1716 :          ALLOCATE (dispersion_env%charges(natom))
     302         4320 :          dispersion_env%charges = charges
     303          858 :          IF (calculate_forces) THEN
     304           12 :             IF (ASSOCIATED(dispersion_env%dcharges)) DEALLOCATE (dispersion_env%dcharges)
     305           24 :             ALLOCATE (dispersion_env%dcharges(natom))
     306           60 :             dispersion_env%dcharges = 0.0_dp
     307              :          END IF
     308              :       END IF
     309              :       CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
     310         1646 :                                         energy%dispersion, calculate_forces)
     311         1646 :       IF (calculate_forces) THEN
     312           42 :          IF (dispersion_env%pp_type == vdw_pairpot_dftd4 .AND. dispersion_env%ext_charges) THEN
     313           60 :             dcharges(1:natom) = dcharges(1:natom) + dispersion_env%dcharges(1:natom)
     314              :          END IF
     315              :       END IF
     316              : 
     317              :       ! Calculate Huckel parameters
     318         1646 :       CALL gfn0_huckel(qs_env, cnumbers, charges, huckel, dhuckel, dqhuckel, calculate_forces)
     319              : 
     320              :       ! Calculate KAB parameters and electronegativity correction
     321         1646 :       CALL gfn0_kpair(qs_env, kijab)
     322              : 
     323              :       ! loop over all atom pairs with a non-zero overlap (sab_orb)
     324         1646 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     325       422349 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     326              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     327       420703 :                                 iatom=iatom, jatom=jatom, r=rij, cell=cell)
     328       420703 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     329       420703 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
     330       420703 :          IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     331       420703 :          CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     332       420703 :          CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
     333       420703 :          IF (.NOT. defined .OR. natorb_b < 1) CYCLE
     334              : 
     335      1682812 :          dr = SQRT(SUM(rij(:)**2))
     336              : 
     337              :          ! atomic parameters
     338              :          CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
     339       420703 :                                  lmax=lmaxa, nshell=nsa, kpoly=kpolya, hen=hena)
     340              :          CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
     341       420703 :                                  lmax=lmaxb, nshell=nsb, kpoly=kpolyb, hen=henb)
     342              : 
     343       420703 :          IF (nimg == 1) THEN
     344              :             ic = 1
     345              :          ELSE
     346       199870 :             ic = cell_to_index(cell(1), cell(2), cell(3))
     347       199870 :             CPASSERT(ic > 0)
     348              :          END IF
     349              : 
     350       420703 :          icol = MAX(iatom, jatom)
     351       420703 :          irow = MIN(iatom, jatom)
     352       420703 :          NULLIFY (sblock, fblock)
     353              :          CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
     354       420703 :                                 row=irow, col=icol, BLOCK=sblock, found=found)
     355       420703 :          CPASSERT(found)
     356              :          CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
     357       420703 :                                 row=irow, col=icol, BLOCK=fblock, found=found)
     358       420703 :          CPASSERT(found)
     359              : 
     360       420703 :          IF (calculate_forces) THEN
     361        11584 :             NULLIFY (pblock)
     362              :             CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
     363        11584 :                                    row=irow, col=icol, block=pblock, found=found)
     364        11584 :             CPASSERT(ASSOCIATED(pblock))
     365        11584 :             NULLIFY (wblock)
     366              :             CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
     367        11584 :                                    row=irow, col=icol, block=wblock, found=found)
     368        11584 :             CPASSERT(ASSOCIATED(wblock))
     369        46336 :             DO i = 2, 4
     370        34752 :                NULLIFY (dsblocks(i)%block)
     371              :                CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
     372        34752 :                                       row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
     373        46336 :                CPASSERT(found)
     374              :             END DO
     375              :          END IF
     376              : 
     377              :          ! overlap
     378       420703 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     379       420703 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     380       420703 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     381       420703 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     382       420703 :          atom_a = atom_of_kind(iatom)
     383       420703 :          atom_b = atom_of_kind(jatom)
     384              :          ! basis ikind
     385       420703 :          first_sgfa => basis_set_a%first_sgf
     386       420703 :          la_max => basis_set_a%lmax
     387       420703 :          la_min => basis_set_a%lmin
     388       420703 :          npgfa => basis_set_a%npgf
     389       420703 :          nseta = basis_set_a%nset
     390       420703 :          nsgfa => basis_set_a%nsgf_set
     391       420703 :          rpgfa => basis_set_a%pgf_radius
     392       420703 :          set_radius_a => basis_set_a%set_radius
     393       420703 :          scon_a => basis_set_a%scon
     394       420703 :          zeta => basis_set_a%zet
     395              :          ! basis jkind
     396       420703 :          first_sgfb => basis_set_b%first_sgf
     397       420703 :          lb_max => basis_set_b%lmax
     398       420703 :          lb_min => basis_set_b%lmin
     399       420703 :          npgfb => basis_set_b%npgf
     400       420703 :          nsetb = basis_set_b%nset
     401       420703 :          nsgfb => basis_set_b%nsgf_set
     402       420703 :          rpgfb => basis_set_b%pgf_radius
     403       420703 :          set_radius_b => basis_set_b%set_radius
     404       420703 :          scon_b => basis_set_b%scon
     405       420703 :          zetb => basis_set_b%zet
     406              : 
     407       420703 :          ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
     408      3365624 :          ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
     409      2103515 :          ALLOCATE (sint(natorb_a, natorb_b, maxder))
     410     35613059 :          sint = 0.0_dp
     411              : 
     412      1619859 :          DO iset = 1, nseta
     413      1199156 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
     414      1199156 :             n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
     415      1199156 :             sgfa = first_sgfa(1, iset)
     416      5043308 :             DO jset = 1, nsetb
     417      3423449 :                IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
     418      1749537 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
     419      1749537 :                n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
     420      1749537 :                sgfb = first_sgfb(1, jset)
     421      1749537 :                IF (calculate_forces) THEN
     422              :                   CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     423              :                                   lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     424        47176 :                                   rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
     425              :                ELSE
     426              :                   CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
     427              :                                   lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
     428      1702361 :                                   rij, sab=oint(:, :, 1))
     429              :                END IF
     430              :                ! Contraction
     431              :                CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
     432      1749537 :                                 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
     433      1749537 :                CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
     434      2948693 :                IF (calculate_forces) THEN
     435       188704 :                   DO i = 2, 4
     436              :                      CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
     437       141528 :                                       cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
     438       188704 :                      CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
     439              :                   END DO
     440              :                END IF
     441              :             END DO
     442              :          END DO
     443              :          ! forces W matrix
     444       420703 :          IF (calculate_forces) THEN
     445        46336 :             DO i = 1, 3
     446        46336 :                IF (iatom <= jatom) THEN
     447      1466268 :                   force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
     448              :                ELSE
     449      1091826 :                   force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
     450              :                END IF
     451              :             END DO
     452        11584 :             f1 = 2.0_dp
     453        46336 :             force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
     454        46336 :             force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
     455        11584 :             IF (use_virial .AND. dr > 1.e-3_dp) THEN
     456        10966 :                IF (iatom == jatom) f1 = 1.0_dp
     457        10966 :                CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
     458              :             END IF
     459              :          END IF
     460              :          ! update S matrix
     461       420703 :          IF (iatom <= jatom) THEN
     462     18553151 :             sblock(:, :) = sblock(:, :) + sint(:, :, 1)
     463              :          ELSE
     464     14217966 :             sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
     465              :          END IF
     466       420703 :          IF (calculate_forces) THEN
     467        46336 :             DO i = 2, 4
     468        46336 :                IF (iatom <= jatom) THEN
     469      1466268 :                   dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
     470              :                ELSE
     471      1108062 :                   dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
     472              :                END IF
     473              :             END DO
     474              :          END IF
     475              : 
     476              :          ! Calculate Pi = Pia * Pib (Eq. 11)
     477       420703 :          rcovab = rcova + rcovb
     478       420703 :          rrab = SQRT(dr/rcovab)
     479      1619859 :          pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
     480      1615219 :          pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
     481       420703 :          IF (calculate_forces) THEN
     482        11584 :             IF (dr > 1.e-6_dp) THEN
     483        11484 :                drx = 0.5_dp/rrab/rcovab
     484              :             ELSE
     485              :                drx = 0.0_dp
     486              :             END IF
     487        44044 :             dpia(1:nsa) = drx*kpolya(1:nsa)
     488        43952 :             dpib(1:nsb) = drx*kpolyb(1:nsb)
     489              :          END IF
     490              : 
     491              :          ! diagonal block
     492       420703 :          diagblock = .FALSE.
     493       420703 :          IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
     494              :          !
     495              :          ! Eq. 10
     496              :          !
     497              :          IF (diagblock) THEN
     498        23510 :             DO i = 1, natorb_a
     499        19635 :                na = naoa(i)
     500        23510 :                fblock(i, i) = fblock(i, i) + huckel(na, iatom)
     501              :             END DO
     502              :          ELSE
     503      3825345 :             DO j = 1, natorb_b
     504      3408517 :                nb = naob(j)
     505     32481137 :                DO i = 1, natorb_a
     506     28655792 :                   na = naoa(i)
     507     28655792 :                   hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
     508     32064309 :                   IF (iatom <= jatom) THEN
     509     16195035 :                      fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
     510              :                   ELSE
     511     12460757 :                      fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
     512              :                   END IF
     513              :                END DO
     514              :             END DO
     515              :          END IF
     516       420703 :          IF (calculate_forces) THEN
     517        11584 :             f0 = 1.0_dp
     518        11584 :             IF (irow == iatom) f0 = -1.0_dp
     519        11584 :             f2 = 1.0_dp
     520        11584 :             IF (iatom /= jatom) f2 = 2.0_dp
     521              :             ! Derivative wrt coordination number
     522        11584 :             fhua = 0.0_dp
     523        11584 :             fhub = 0.0_dp
     524        11584 :             fhud = 0.0_dp
     525        11584 :             fqa = 0.0_dp
     526        11584 :             fqb = 0.0_dp
     527        11584 :             IF (diagblock) THEN
     528          542 :                DO i = 1, natorb_a
     529          442 :                   la = laoa(i)
     530          442 :                   na = naoa(i)
     531          442 :                   fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
     532          542 :                   fqa = fqa + pblock(i, i)*dqhuckel(na, iatom)
     533              :                END DO
     534          100 :                dcharges(iatom) = dcharges(iatom) + fqa
     535              :             ELSE
     536       102648 :                DO j = 1, natorb_b
     537        91164 :                   lb = laob(j)
     538        91164 :                   nb = naob(j)
     539       849534 :                   DO i = 1, natorb_a
     540       746886 :                      la = laoa(i)
     541       746886 :                      na = naoa(i)
     542       746886 :                      hij = 0.5_dp*pia(na)*pib(nb)
     543       746886 :                      drx = f2*hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)
     544       838050 :                      IF (iatom <= jatom) THEN
     545       425196 :                         fhua = fhua + drx*pblock(i, j)*dhuckel(na, iatom)
     546       425196 :                         fhub = fhub + drx*pblock(i, j)*dhuckel(nb, jatom)
     547       425196 :                         fqa = fqa + drx*pblock(i, j)*dqhuckel(na, iatom)
     548       425196 :                         fqb = fqb + drx*pblock(i, j)*dqhuckel(nb, jatom)
     549              :                      ELSE
     550       321690 :                         fhua = fhua + drx*pblock(j, i)*dhuckel(na, iatom)
     551       321690 :                         fhub = fhub + drx*pblock(j, i)*dhuckel(nb, jatom)
     552       321690 :                         fqa = fqa + drx*pblock(j, i)*dqhuckel(na, iatom)
     553       321690 :                         fqb = fqb + drx*pblock(j, i)*dqhuckel(nb, jatom)
     554              :                      END IF
     555              :                   END DO
     556              :                END DO
     557        11484 :                dcharges(iatom) = dcharges(iatom) + fqa
     558        11484 :                dcharges(jatom) = dcharges(jatom) + fqb
     559              :             END IF
     560              :             ! iatom
     561        11584 :             atom_a = atom_of_kind(iatom)
     562       158260 :             DO i = 1, dcnum(iatom)%neighbors
     563       146676 :                katom = dcnum(iatom)%nlist(i)
     564       146676 :                kkind = kind_of(katom)
     565       146676 :                atom_c = atom_of_kind(katom)
     566       586704 :                rik = dcnum(iatom)%rik(:, i)
     567       586704 :                drk = SQRT(SUM(rik(:)**2))
     568       158260 :                IF (drk > 1.e-3_dp) THEN
     569       586704 :                   fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
     570       586704 :                   force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
     571       586704 :                   force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
     572       586704 :                   fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
     573       586704 :                   force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
     574       586704 :                   force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
     575       146676 :                   IF (use_virial) THEN
     576       583032 :                      fdik = fdika + fdikb
     577       145758 :                      CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
     578              :                   END IF
     579              :                END IF
     580              :             END DO
     581              :             ! jatom
     582        11584 :             atom_b = atom_of_kind(jatom)
     583       157422 :             DO i = 1, dcnum(jatom)%neighbors
     584       145838 :                katom = dcnum(jatom)%nlist(i)
     585       145838 :                kkind = kind_of(katom)
     586       145838 :                atom_c = atom_of_kind(katom)
     587       583352 :                rik = dcnum(jatom)%rik(:, i)
     588       583352 :                drk = SQRT(SUM(rik(:)**2))
     589       157422 :                IF (drk > 1.e-3_dp) THEN
     590       583352 :                   fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
     591       583352 :                   force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
     592       583352 :                   force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
     593       145838 :                   IF (use_virial) THEN
     594       145010 :                      CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
     595              :                   END IF
     596              :                END IF
     597              :             END DO
     598              :             ! force from R dendent Huckel element: Pia*Pib
     599        11584 :             IF (diagblock) THEN
     600          100 :                force_ab = 0._dp
     601              :             ELSE
     602        11484 :                n1 = SIZE(fblock, 1)
     603        11484 :                n2 = SIZE(fblock, 2)
     604        45936 :                ALLOCATE (dfblock(n1, n2))
     605       854946 :                dfblock = 0.0_dp
     606       102648 :                DO j = 1, natorb_b
     607        91164 :                   lb = laob(j)
     608        91164 :                   nb = naob(j)
     609       849534 :                   DO i = 1, natorb_a
     610       746886 :                      la = laoa(i)
     611       746886 :                      na = naoa(i)
     612       746886 :                      dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
     613       838050 :                      IF (iatom <= jatom) THEN
     614       425196 :                         dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
     615              :                      ELSE
     616       321690 :                         dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
     617              :                      END IF
     618              :                   END DO
     619              :                END DO
     620       854946 :                dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
     621        45936 :                DO ir = 1, 3
     622        34452 :                   foab = 2.0_dp*dfp*rij(ir)/dr
     623              :                   ! force from overlap matrix contribution to H
     624       307944 :                   DO j = 1, natorb_b
     625       273492 :                      lb = laob(j)
     626       273492 :                      nb = naob(j)
     627      2548602 :                      DO i = 1, natorb_a
     628      2240658 :                         la = laoa(i)
     629      2240658 :                         na = naoa(i)
     630      2240658 :                         hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
     631      2514150 :                         IF (iatom <= jatom) THEN
     632      1275588 :                            foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
     633              :                         ELSE
     634       965070 :                            foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
     635              :                         END IF
     636              :                      END DO
     637              :                   END DO
     638        45936 :                   force_ab(ir) = foab
     639              :                END DO
     640        11484 :                DEALLOCATE (dfblock)
     641              :             END IF
     642              :          END IF
     643              : 
     644       420703 :          IF (calculate_forces) THEN
     645        11584 :             atom_a = atom_of_kind(iatom)
     646        11584 :             atom_b = atom_of_kind(jatom)
     647        31156 :             IF (irow == iatom) force_ab = -force_ab
     648        46336 :             force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
     649        46336 :             force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
     650        11584 :             IF (use_virial) THEN
     651        11002 :                f1 = 1.0_dp
     652        11002 :                IF (iatom == jatom) f1 = 0.5_dp
     653        11002 :                CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
     654              :             END IF
     655              :          END IF
     656              : 
     657      2525864 :          DEALLOCATE (oint, owork, sint)
     658              : 
     659              :       END DO
     660         1646 :       CALL neighbor_list_iterator_release(nl_iterator)
     661              : 
     662         3292 :       DO i = 1, SIZE(matrix_h, 1)
     663        32198 :          DO img = 1, nimg
     664        28906 :             CALL dbcsr_finalize(matrix_h(i, img)%matrix)
     665        30552 :             CALL dbcsr_finalize(matrix_s(i, img)%matrix)
     666              :          END DO
     667              :       END DO
     668              : 
     669              :       ! EEQ forces (response and direct)
     670         1646 :       IF (calculate_forces) THEN
     671           42 :          CALL para_env%sum(dcharges)
     672           84 :          ALLOCATE (qlagrange(natom))
     673           42 :          CALL xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
     674              :       END IF
     675              : 
     676         1646 :       kf = xtb_control%kf
     677         1646 :       enscale = xtb_control%enscale
     678         1646 :       erep = 0.0_dp
     679         1646 :       CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
     680              : 
     681         1646 :       esrb = 0.0_dp
     682         1646 :       CALL srb_potential(qs_env, esrb, calculate_forces, xtb_control, cnumbers, dcnum)
     683              : 
     684         1646 :       enonbonded = 0.0_dp
     685         1646 :       IF (do_nonbonded) THEN
     686              :          ! nonbonded interactions
     687            0 :          NULLIFY (sab_xtb_nonbond)
     688            0 :          CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
     689              :          CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
     690            0 :                                    atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
     691              :       END IF
     692              : 
     693              :       ! set repulsive energy
     694         1646 :       erep = erep + esrb + enonbonded
     695         1646 :       IF (do_nonbonded) THEN
     696            0 :          CALL para_env%sum(enonbonded)
     697            0 :          energy%xtb_nonbonded = enonbonded
     698              :       END IF
     699         1646 :       CALL para_env%sum(esrb)
     700         1646 :       energy%srb = esrb
     701         1646 :       CALL para_env%sum(erep)
     702         1646 :       energy%repulsive = erep
     703              : 
     704              :       ! save EEQ charges
     705         1646 :       NULLIFY (eeq_q)
     706         1646 :       CALL get_qs_env(qs_env, eeq=eeq_q)
     707         1646 :       IF (ASSOCIATED(eeq_q)) THEN
     708          972 :          CPASSERT(SIZE(eeq_q) == natom)
     709              :       ELSE
     710         1348 :          ALLOCATE (eeq_q(natom))
     711         3432 :          eeq_q(1:natom) = charges(1:natom)
     712              :       END IF
     713         1646 :       CALL set_qs_env(qs_env, eeq=eeq_q)
     714              : 
     715              :       ! deallocate coordination numbers
     716         1646 :       CALL cnumber_release(cnumbers, dcnum, calculate_forces)
     717              : 
     718              :       ! deallocate Huckel parameters
     719         1646 :       DEALLOCATE (huckel)
     720         1646 :       IF (calculate_forces) THEN
     721           42 :          DEALLOCATE (dhuckel, dqhuckel)
     722              :       END IF
     723              :       ! deallocate KAB parameters
     724         1646 :       DEALLOCATE (kijab)
     725              : 
     726              :       ! deallocate charges
     727         1646 :       DEALLOCATE (charges)
     728         1646 :       IF (calculate_forces) THEN
     729           42 :          DEALLOCATE (dcharges, qlagrange)
     730              :       END IF
     731              : 
     732              :       ! AO matrix outputs
     733         1646 :       CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
     734              : 
     735         1646 :       DEALLOCATE (basis_set_list)
     736         1646 :       IF (calculate_forces) THEN
     737           42 :          IF (SIZE(matrix_p, 1) == 2) THEN
     738            8 :             DO img = 1, nimg
     739              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
     740            8 :                               beta_scalar=-1.0_dp)
     741              :             END DO
     742              :          END IF
     743              :       END IF
     744              : 
     745         1646 :       CALL timestop(handle)
     746              : 
     747         4938 :    END SUBROUTINE build_gfn0_xtb_matrices
     748              : 
     749              : ! **************************************************************************************************
     750              : !> \brief ...
     751              : !> \param qs_env ...
     752              : !> \param calculate_forces ...
     753              : ! **************************************************************************************************
     754         3578 :    SUBROUTINE build_gfn1_xtb_matrices(qs_env, calculate_forces)
     755              : 
     756              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     757              :       LOGICAL, INTENT(IN)                                :: calculate_forces
     758              : 
     759              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'build_gfn1_xtb_matrices'
     760              : 
     761              :       INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, ic, icol, ikind, img, ir, irow, iset, &
     762              :          j, jatom, jkind, jset, katom, kkind, la, lb, ldsab, lmaxa, lmaxb, maxder, n1, n2, na, &
     763              :          natom, natorb_a, natorb_b, nb, ncoa, ncob, nderivatives, nimg, nkind, nsa, nsb, nseta, &
     764              :          nsetb, sgfa, sgfb, za, zb
     765         7156 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     766              :       INTEGER, DIMENSION(25)                             :: laoa, laob, naoa, naob
     767              :       INTEGER, DIMENSION(3)                              :: cell
     768         3578 :       INTEGER, DIMENSION(:), POINTER                     :: la_max, la_min, lb_max, lb_min, npgfa, &
     769         3578 :                                                             npgfb, nsgfa, nsgfb
     770         3578 :       INTEGER, DIMENSION(:, :), POINTER                  :: first_sgfa, first_sgfb
     771         3578 :       INTEGER, DIMENSION(:, :, :), POINTER               :: cell_to_index
     772              :       LOGICAL                                            :: defined, diagblock, do_nonbonded, found, &
     773              :                                                             use_virial, xb_inter
     774              :       REAL(KIND=dp)                                      :: dfp, dhij, dr, drk, drx, enonbonded, &
     775              :                                                             enscale, erep, etaa, etab, exb, f0, &
     776              :                                                             f1, fhua, fhub, fhud, foab, hij, kf, &
     777              :                                                             rcova, rcovab, rcovb, rrab
     778         3578 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: cnumbers
     779         7156 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: dfblock, dhuckel, huckel, owork
     780         3578 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: oint, sint
     781         3578 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :)  :: kijab
     782              :       REAL(KIND=dp), DIMENSION(3)                        :: fdik, fdika, fdikb, force_ab, rij, rik
     783              :       REAL(KIND=dp), DIMENSION(5)                        :: dpia, dpib, kpolya, kpolyb, pia, pib
     784         3578 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: set_radius_a, set_radius_b
     785         3578 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: fblock, pblock, rpgfa, rpgfb, sblock, &
     786         3578 :                                                             scon_a, scon_b, wblock, zeta, zetb
     787         3578 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     788              :       TYPE(atprop_type), POINTER                         :: atprop
     789        14312 :       TYPE(block_p_type), DIMENSION(2:4)                 :: dsblocks
     790              :       TYPE(cp_logger_type), POINTER                      :: logger
     791         3578 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_p, matrix_s, matrix_w
     792         3578 :       TYPE(dcnum_type), ALLOCATABLE, DIMENSION(:)        :: dcnum
     793              :       TYPE(dft_control_type), POINTER                    :: dft_control
     794         3578 :       TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER  :: basis_set_list
     795              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_a, basis_set_b
     796              :       TYPE(kpoint_type), POINTER                         :: kpoints
     797              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     798              :       TYPE(neighbor_list_iterator_p_type), &
     799         3578 :          DIMENSION(:), POINTER                           :: nl_iterator
     800              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     801         3578 :          POINTER                                         :: sab_orb, sab_xtb_nonbond
     802         3578 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     803              :       TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
     804              :       TYPE(qs_energy_type), POINTER                      :: energy
     805         3578 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     806         3578 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     807              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     808              :       TYPE(qs_rho_type), POINTER                         :: rho
     809              :       TYPE(virial_type), POINTER                         :: virial
     810              :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b
     811              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
     812              : 
     813         3578 :       CALL timeset(routineN, handle)
     814              : 
     815         3578 :       NULLIFY (logger, virial, atprop)
     816         3578 :       logger => cp_get_default_logger()
     817              : 
     818         3578 :       NULLIFY (matrix_h, matrix_s, matrix_p, matrix_w, atomic_kind_set, &
     819         3578 :                qs_kind_set, sab_orb, ks_env)
     820              : 
     821              :       CALL get_qs_env(qs_env=qs_env, &
     822              :                       ks_env=ks_env, &
     823              :                       energy=energy, &
     824              :                       atomic_kind_set=atomic_kind_set, &
     825              :                       qs_kind_set=qs_kind_set, &
     826              :                       matrix_h_kp=matrix_h, &
     827              :                       matrix_s_kp=matrix_s, &
     828              :                       para_env=para_env, &
     829              :                       atprop=atprop, &
     830              :                       dft_control=dft_control, &
     831         3578 :                       sab_orb=sab_orb)
     832              : 
     833         3578 :       nkind = SIZE(atomic_kind_set)
     834         3578 :       xtb_control => dft_control%qs_control%xtb_control
     835         3578 :       xb_inter = xtb_control%xb_interaction
     836         3578 :       do_nonbonded = xtb_control%do_nonbonded
     837         3578 :       nimg = dft_control%nimages
     838         3578 :       nderivatives = 0
     839         3578 :       IF (calculate_forces) nderivatives = 1
     840         3578 :       IF (dft_control%tddfpt2_control%enabled) nderivatives = 1
     841         3578 :       maxder = ncoset(nderivatives)
     842              : 
     843         3578 :       NULLIFY (particle_set)
     844         3578 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set)
     845         3578 :       natom = SIZE(particle_set)
     846              :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     847         3578 :                                atom_of_kind=atom_of_kind, kind_of=kind_of)
     848              : 
     849         3578 :       IF (calculate_forces) THEN
     850          498 :          NULLIFY (rho, force, matrix_w)
     851              :          CALL get_qs_env(qs_env=qs_env, &
     852              :                          rho=rho, matrix_w_kp=matrix_w, &
     853          498 :                          virial=virial, force=force)
     854          498 :          CALL qs_rho_get(rho, rho_ao_kp=matrix_p)
     855              : 
     856          498 :          IF (SIZE(matrix_p, 1) == 2) THEN
     857          432 :             DO img = 1, nimg
     858              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
     859          414 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     860              :                CALL dbcsr_add(matrix_w(1, img)%matrix, matrix_w(2, img)%matrix, &
     861          432 :                               alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     862              :             END DO
     863              :          END IF
     864          886 :          use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     865              :       END IF
     866              :       ! atomic energy decomposition
     867         3578 :       IF (atprop%energy) THEN
     868           36 :          CALL atprop_array_init(atprop%atecc, natom)
     869              :       END IF
     870              : 
     871         3578 :       NULLIFY (cell_to_index)
     872         3578 :       IF (nimg > 1) THEN
     873          336 :          CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
     874          336 :          CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
     875              :       END IF
     876              : 
     877              :       ! set up basis set lists
     878        19570 :       ALLOCATE (basis_set_list(nkind))
     879         3578 :       CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
     880              : 
     881              :       ! allocate overlap matrix
     882         3578 :       CALL dbcsr_allocate_matrix_set(matrix_s, maxder, nimg)
     883              :       CALL create_sab_matrix(ks_env, matrix_s, "xTB OVERLAP MATRIX", basis_set_list, basis_set_list, &
     884         3578 :                              sab_orb, .TRUE.)
     885         3578 :       CALL set_ks_env(ks_env, matrix_s_kp=matrix_s)
     886              : 
     887              :       ! initialize H matrix
     888         3578 :       CALL dbcsr_allocate_matrix_set(matrix_h, 1, nimg)
     889        50792 :       DO img = 1, nimg
     890        47214 :          ALLOCATE (matrix_h(1, img)%matrix)
     891              :          CALL dbcsr_create(matrix_h(1, img)%matrix, template=matrix_s(1, 1)%matrix, &
     892        47214 :                            name="HAMILTONIAN MATRIX")
     893        50792 :          CALL cp_dbcsr_alloc_block_from_nbl(matrix_h(1, img)%matrix, sab_orb)
     894              :       END DO
     895         3578 :       CALL set_ks_env(ks_env, matrix_h_kp=matrix_h)
     896              : 
     897              :       ! Calculate coordination numbers
     898              :       ! needed for effective atomic energy levels (Eq. 12)
     899              :       ! code taken from D3 dispersion energy
     900         3578 :       CALL cnumber_init(qs_env, cnumbers, dcnum, 1, calculate_forces)
     901              : 
     902              :       ! vdW Potential
     903         3578 :       CALL get_qs_env(qs_env=qs_env, dispersion_env=dispersion_env)
     904              :       CALL calculate_dispersion_pairpot(qs_env, dispersion_env, &
     905         3578 :                                         energy%dispersion, calculate_forces)
     906              : 
     907              :       ! Calculate Huckel parameters
     908         3578 :       CALL gfn1_huckel(qs_env, cnumbers, huckel, dhuckel, calculate_forces)
     909              : 
     910              :       ! Calculate KAB parameters and electronegativity correction
     911         3578 :       CALL gfn1_kpair(qs_env, kijab)
     912              : 
     913              :       ! loop over all atom pairs with a non-zero overlap (sab_orb)
     914         3578 :       CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
     915      1077834 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     916              :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     917      1074256 :                                 iatom=iatom, jatom=jatom, r=rij, cell=cell)
     918      1074256 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     919      1074256 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
     920      1074256 :          IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     921      1074256 :          CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     922      1074256 :          CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
     923      1074256 :          IF (.NOT. defined .OR. natorb_b < 1) CYCLE
     924              : 
     925      4297024 :          dr = SQRT(SUM(rij(:)**2))
     926              : 
     927              :          ! atomic parameters
     928              :          CALL get_xtb_atom_param(xtb_atom_a, z=za, nao=naoa, lao=laoa, rcov=rcova, eta=etaa, &
     929      1074256 :                                  lmax=lmaxa, nshell=nsa, kpoly=kpolya)
     930              :          CALL get_xtb_atom_param(xtb_atom_b, z=zb, nao=naob, lao=laob, rcov=rcovb, eta=etab, &
     931      1074256 :                                  lmax=lmaxb, nshell=nsb, kpoly=kpolyb)
     932              : 
     933      1074256 :          IF (nimg == 1) THEN
     934              :             ic = 1
     935              :          ELSE
     936       241506 :             ic = cell_to_index(cell(1), cell(2), cell(3))
     937       241506 :             CPASSERT(ic > 0)
     938              :          END IF
     939              : 
     940      1074256 :          icol = MAX(iatom, jatom)
     941      1074256 :          irow = MIN(iatom, jatom)
     942      1074256 :          NULLIFY (sblock, fblock)
     943              :          CALL dbcsr_get_block_p(matrix=matrix_s(1, ic)%matrix, &
     944      1074256 :                                 row=irow, col=icol, BLOCK=sblock, found=found)
     945      1074256 :          CPASSERT(found)
     946              :          CALL dbcsr_get_block_p(matrix=matrix_h(1, ic)%matrix, &
     947      1074256 :                                 row=irow, col=icol, BLOCK=fblock, found=found)
     948      1074256 :          CPASSERT(found)
     949              : 
     950      1074256 :          IF (calculate_forces) THEN
     951       242637 :             NULLIFY (pblock)
     952              :             CALL dbcsr_get_block_p(matrix=matrix_p(1, ic)%matrix, &
     953       242637 :                                    row=irow, col=icol, block=pblock, found=found)
     954       242637 :             CPASSERT(found)
     955       242637 :             NULLIFY (wblock)
     956              :             CALL dbcsr_get_block_p(matrix=matrix_w(1, ic)%matrix, &
     957       242637 :                                    row=irow, col=icol, block=wblock, found=found)
     958       242637 :             CPASSERT(found)
     959       970548 :             DO i = 2, 4
     960       727911 :                NULLIFY (dsblocks(i)%block)
     961              :                CALL dbcsr_get_block_p(matrix=matrix_s(i, ic)%matrix, &
     962       727911 :                                       row=irow, col=icol, BLOCK=dsblocks(i)%block, found=found)
     963       970548 :                CPASSERT(found)
     964              :             END DO
     965              :          END IF
     966              : 
     967              :          ! overlap
     968      1074256 :          basis_set_a => basis_set_list(ikind)%gto_basis_set
     969      1074256 :          IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
     970      1074256 :          basis_set_b => basis_set_list(jkind)%gto_basis_set
     971      1074256 :          IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
     972      1074256 :          atom_a = atom_of_kind(iatom)
     973      1074256 :          atom_b = atom_of_kind(jatom)
     974              :          ! basis ikind
     975      1074256 :          first_sgfa => basis_set_a%first_sgf
     976      1074256 :          la_max => basis_set_a%lmax
     977      1074256 :          la_min => basis_set_a%lmin
     978      1074256 :          npgfa => basis_set_a%npgf
     979      1074256 :          nseta = basis_set_a%nset
     980      1074256 :          nsgfa => basis_set_a%nsgf_set
     981      1074256 :          rpgfa => basis_set_a%pgf_radius
     982      1074256 :          set_radius_a => basis_set_a%set_radius
     983      1074256 :          scon_a => basis_set_a%scon
     984      1074256 :          zeta => basis_set_a%zet
     985              :          ! basis jkind
     986      1074256 :          first_sgfb => basis_set_b%first_sgf
     987      1074256 :          lb_max => basis_set_b%lmax
     988      1074256 :          lb_min => basis_set_b%lmin
     989      1074256 :          npgfb => basis_set_b%npgf
     990      1074256 :          nsetb = basis_set_b%nset
     991      1074256 :          nsgfb => basis_set_b%nsgf_set
     992      1074256 :          rpgfb => basis_set_b%pgf_radius
     993      1074256 :          set_radius_b => basis_set_b%set_radius
     994      1074256 :          scon_b => basis_set_b%scon
     995      1074256 :          zetb => basis_set_b%zet
     996              : 
     997      1074256 :          ldsab = get_memory_usage(qs_kind_set, "ORB", "ORB")
     998      8594048 :          ALLOCATE (oint(ldsab, ldsab, maxder), owork(ldsab, ldsab))
     999      5371280 :          ALLOCATE (sint(natorb_a, natorb_b, maxder))
    1000     52484526 :          sint = 0.0_dp
    1001              : 
    1002      3397512 :          DO iset = 1, nseta
    1003      2323256 :             ncoa = npgfa(iset)*ncoset(la_max(iset))
    1004      2323256 :             n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
    1005      2323256 :             sgfa = first_sgfa(1, iset)
    1006      8543482 :             DO jset = 1, nsetb
    1007      5145970 :                IF (set_radius_a(iset) + set_radius_b(jset) < dr) CYCLE
    1008      4132912 :                ncob = npgfb(jset)*ncoset(lb_max(jset))
    1009      4132912 :                n2 = npgfb(jset)*(ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1))
    1010      4132912 :                sgfb = first_sgfb(1, jset)
    1011      4132912 :                IF (calculate_forces) THEN
    1012              :                   CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
    1013              :                                   lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
    1014       953017 :                                   rij, sab=oint(:, :, 1), dab=oint(:, :, 2:4))
    1015              :                ELSE
    1016              :                   CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
    1017              :                                   lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
    1018      3179895 :                                   rij, sab=oint(:, :, 1))
    1019              :                END IF
    1020              :                ! Contraction
    1021              :                CALL contraction(oint(:, :, 1), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
    1022      4132912 :                                 cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
    1023      4132912 :                CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, 1), sgfa, sgfb, trans=.FALSE.)
    1024      6456168 :                IF (calculate_forces) THEN
    1025      3812068 :                   DO i = 2, 4
    1026              :                      CALL contraction(oint(:, :, i), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
    1027      2859051 :                                       cb=scon_b(:, sgfb:), nb=n2, mb=nsgfb(jset), fscale=1.0_dp, trans=.FALSE.)
    1028      3812068 :                      CALL block_add("IN", owork, nsgfa(iset), nsgfb(jset), sint(:, :, i), sgfa, sgfb, trans=.FALSE.)
    1029              :                   END DO
    1030              :                END IF
    1031              :             END DO
    1032              :          END DO
    1033              :          ! forces W matrix
    1034      1074256 :          IF (calculate_forces) THEN
    1035       970548 :             DO i = 1, 3
    1036       970548 :                IF (iatom <= jatom) THEN
    1037     14173335 :                   force_ab(i) = SUM(sint(:, :, i + 1)*wblock(:, :))
    1038              :                ELSE
    1039     10113909 :                   force_ab(i) = SUM(sint(:, :, i + 1)*TRANSPOSE(wblock(:, :)))
    1040              :                END IF
    1041              :             END DO
    1042       242637 :             f1 = 2.0_dp
    1043       970548 :             force(ikind)%overlap(:, atom_a) = force(ikind)%overlap(:, atom_a) - f1*force_ab(:)
    1044       970548 :             force(jkind)%overlap(:, atom_b) = force(jkind)%overlap(:, atom_b) + f1*force_ab(:)
    1045       242637 :             IF (use_virial .AND. dr > 1.e-3_dp) THEN
    1046       148638 :                IF (iatom == jatom) f1 = 1.0_dp
    1047       148638 :                CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
    1048              :             END IF
    1049              :          END IF
    1050              :          ! update S matrix
    1051      1074256 :          IF (iatom <= jatom) THEN
    1052     14410851 :             sblock(:, :) = sblock(:, :) + sint(:, :, 1)
    1053              :          ELSE
    1054     10570739 :             sblock(:, :) = sblock(:, :) + TRANSPOSE(sint(:, :, 1))
    1055              :          END IF
    1056      1074256 :          IF (calculate_forces) THEN
    1057       970548 :             DO i = 2, 4
    1058       970548 :                IF (iatom <= jatom) THEN
    1059     14173335 :                   dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) + sint(:, :, i)
    1060              :                ELSE
    1061      9947919 :                   dsblocks(i)%block(:, :) = dsblocks(i)%block(:, :) - TRANSPOSE(sint(:, :, i))
    1062              :                END IF
    1063              :             END DO
    1064              :          END IF
    1065              : 
    1066              :          ! Calculate Pi = Pia * Pib (Eq. 11)
    1067      1074256 :          rcovab = rcova + rcovb
    1068      1074256 :          rrab = SQRT(dr/rcovab)
    1069      3397512 :          pia(1:nsa) = 1._dp + kpolya(1:nsa)*rrab
    1070      3397517 :          pib(1:nsb) = 1._dp + kpolyb(1:nsb)*rrab
    1071      1074256 :          IF (calculate_forces) THEN
    1072       242637 :             IF (dr > 1.e-6_dp) THEN
    1073       240347 :                drx = 0.5_dp/rrab/rcovab
    1074              :             ELSE
    1075              :                drx = 0.0_dp
    1076              :             END IF
    1077       799289 :             dpia(1:nsa) = drx*kpolya(1:nsa)
    1078       799291 :             dpib(1:nsb) = drx*kpolyb(1:nsb)
    1079              :          END IF
    1080              : 
    1081              :          ! diagonal block
    1082      1074256 :          diagblock = .FALSE.
    1083      1074256 :          IF (iatom == jatom .AND. dr < 0.001_dp) diagblock = .TRUE.
    1084              :          !
    1085              :          ! Eq. 10
    1086              :          !
    1087              :          IF (diagblock) THEN
    1088        64517 :             DO i = 1, natorb_a
    1089        49008 :                na = naoa(i)
    1090        64517 :                fblock(i, i) = fblock(i, i) + huckel(na, iatom)
    1091              :             END DO
    1092              :          ELSE
    1093      4990916 :             DO j = 1, natorb_b
    1094      3932169 :                nb = naob(j)
    1095     24944112 :                DO i = 1, natorb_a
    1096     19953196 :                   na = naoa(i)
    1097     19953196 :                   hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
    1098     23885365 :                   IF (iatom <= jatom) THEN
    1099     11533986 :                      fblock(i, j) = fblock(i, j) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
    1100              :                   ELSE
    1101      8419210 :                      fblock(j, i) = fblock(j, i) + hij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
    1102              :                   END IF
    1103              :                END DO
    1104              :             END DO
    1105              :          END IF
    1106      1074256 :          IF (calculate_forces) THEN
    1107       242637 :             f0 = 1.0_dp
    1108       242637 :             IF (irow == iatom) f0 = -1.0_dp
    1109              :             ! Derivative wrt coordination number
    1110       242637 :             fhua = 0.0_dp
    1111       242637 :             fhub = 0.0_dp
    1112       242637 :             fhud = 0.0_dp
    1113       242637 :             IF (diagblock) THEN
    1114         9920 :                DO i = 1, natorb_a
    1115         7630 :                   la = laoa(i)
    1116         7630 :                   na = naoa(i)
    1117         9920 :                   fhud = fhud + pblock(i, i)*dhuckel(na, iatom)
    1118              :                END DO
    1119              :             ELSE
    1120      1325101 :                DO j = 1, natorb_b
    1121      1084754 :                   lb = laob(j)
    1122      1084754 :                   nb = naob(j)
    1123      8051298 :                   DO i = 1, natorb_a
    1124      6726197 :                      la = laoa(i)
    1125      6726197 :                      na = naoa(i)
    1126      6726197 :                      hij = 0.5_dp*pia(na)*pib(nb)
    1127      7810951 :                      IF (iatom <= jatom) THEN
    1128      3968344 :                         fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(na, iatom)
    1129      3968344 :                         fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(i, j)*dhuckel(nb, jatom)
    1130              :                      ELSE
    1131      2757853 :                         fhua = fhua + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(na, iatom)
    1132      2757853 :                         fhub = fhub + hij*kijab(i, j, ikind, jkind)*sint(i, j, 1)*pblock(j, i)*dhuckel(nb, jatom)
    1133              :                      END IF
    1134              :                   END DO
    1135              :                END DO
    1136       240347 :                IF (iatom /= jatom) THEN
    1137       207279 :                   fhua = 2.0_dp*fhua
    1138       207279 :                   fhub = 2.0_dp*fhub
    1139              :                END IF
    1140              :             END IF
    1141              :             ! iatom
    1142       242637 :             atom_a = atom_of_kind(iatom)
    1143      7931159 :             DO i = 1, dcnum(iatom)%neighbors
    1144      7688522 :                katom = dcnum(iatom)%nlist(i)
    1145      7688522 :                kkind = kind_of(katom)
    1146      7688522 :                atom_c = atom_of_kind(katom)
    1147     30754088 :                rik = dcnum(iatom)%rik(:, i)
    1148     30754088 :                drk = SQRT(SUM(rik(:)**2))
    1149      7931159 :                IF (drk > 1.e-3_dp) THEN
    1150     30754088 :                   fdika(:) = fhua*dcnum(iatom)%dvals(i)*rik(:)/drk
    1151     30754088 :                   force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdika(:)
    1152     30754088 :                   force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdika(:)
    1153     30754088 :                   fdikb(:) = fhud*dcnum(iatom)%dvals(i)*rik(:)/drk
    1154     30754088 :                   force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - fdikb(:)
    1155     30754088 :                   force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdikb(:)
    1156      7688522 :                   IF (use_virial) THEN
    1157     20580152 :                      fdik = fdika + fdikb
    1158      5145038 :                      CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
    1159              :                   END IF
    1160              :                END IF
    1161              :             END DO
    1162              :             ! jatom
    1163       242637 :             atom_b = atom_of_kind(jatom)
    1164      7928354 :             DO i = 1, dcnum(jatom)%neighbors
    1165      7685717 :                katom = dcnum(jatom)%nlist(i)
    1166      7685717 :                kkind = kind_of(katom)
    1167      7685717 :                atom_c = atom_of_kind(katom)
    1168     30742868 :                rik = dcnum(jatom)%rik(:, i)
    1169     30742868 :                drk = SQRT(SUM(rik(:)**2))
    1170      7928354 :                IF (drk > 1.e-3_dp) THEN
    1171     30742868 :                   fdik(:) = fhub*dcnum(jatom)%dvals(i)*rik(:)/drk
    1172     30742868 :                   force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) - fdik(:)
    1173     30742868 :                   force(kkind)%all_potential(:, atom_c) = force(kkind)%all_potential(:, atom_c) + fdik(:)
    1174      7685717 :                   IF (use_virial) THEN
    1175      5144024 :                      CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
    1176              :                   END IF
    1177              :                END IF
    1178              :             END DO
    1179              :             ! force from R dendent Huckel element: Pia*Pib
    1180       242637 :             IF (diagblock) THEN
    1181         2290 :                force_ab = 0._dp
    1182              :             ELSE
    1183       240347 :                n1 = SIZE(fblock, 1)
    1184       240347 :                n2 = SIZE(fblock, 2)
    1185       961388 :                ALLOCATE (dfblock(n1, n2))
    1186      7995968 :                dfblock = 0.0_dp
    1187      1325101 :                DO j = 1, natorb_b
    1188      1084754 :                   lb = laob(j)
    1189      1084754 :                   nb = naob(j)
    1190      8051298 :                   DO i = 1, natorb_a
    1191      6726197 :                      la = laoa(i)
    1192      6726197 :                      na = naoa(i)
    1193      6726197 :                      dhij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*(dpia(na)*pib(nb) + pia(na)*dpib(nb))
    1194      7810951 :                      IF (iatom <= jatom) THEN
    1195      3968344 :                         dfblock(i, j) = dfblock(i, j) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
    1196              :                      ELSE
    1197      2757853 :                         dfblock(j, i) = dfblock(j, i) + dhij*sint(i, j, 1)*kijab(i, j, ikind, jkind)
    1198              :                      END IF
    1199              :                   END DO
    1200              :                END DO
    1201      7995968 :                dfp = f0*SUM(dfblock(:, :)*pblock(:, :))
    1202       961388 :                DO ir = 1, 3
    1203       721041 :                   foab = 2.0_dp*dfp*rij(ir)/dr
    1204              :                   ! force from overlap matrix contribution to H
    1205      3975303 :                   DO j = 1, natorb_b
    1206      3254262 :                      lb = laob(j)
    1207      3254262 :                      nb = naob(j)
    1208     24153894 :                      DO i = 1, natorb_a
    1209     20178591 :                         la = laoa(i)
    1210     20178591 :                         na = naoa(i)
    1211     20178591 :                         hij = 0.5_dp*(huckel(na, iatom) + huckel(nb, jatom))*pia(na)*pib(nb)
    1212     23432853 :                         IF (iatom <= jatom) THEN
    1213     11905032 :                            foab = foab + 2.0_dp*hij*sint(i, j, ir + 1)*pblock(i, j)*kijab(i, j, ikind, jkind)
    1214              :                         ELSE
    1215      8273559 :                            foab = foab - 2.0_dp*hij*sint(i, j, ir + 1)*pblock(j, i)*kijab(i, j, ikind, jkind)
    1216              :                         END IF
    1217              :                      END DO
    1218              :                   END DO
    1219       961388 :                   force_ab(ir) = foab
    1220              :                END DO
    1221       240347 :                DEALLOCATE (dfblock)
    1222              :             END IF
    1223              :          END IF
    1224              : 
    1225      1074256 :          IF (calculate_forces) THEN
    1226       242637 :             atom_a = atom_of_kind(iatom)
    1227       242637 :             atom_b = atom_of_kind(jatom)
    1228       634569 :             IF (irow == iatom) force_ab = -force_ab
    1229       970548 :             force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_ab(:)
    1230       970548 :             force(jkind)%all_potential(:, atom_b) = force(jkind)%all_potential(:, atom_b) + force_ab(:)
    1231       242637 :             IF (use_virial) THEN
    1232       149490 :                f1 = 1.0_dp
    1233       149490 :                IF (iatom == jatom) f1 = 0.5_dp
    1234       149490 :                CALL virial_pair_force(virial%pv_virial, -f1, force_ab, rij)
    1235              :             END IF
    1236              :          END IF
    1237              : 
    1238      6449114 :          DEALLOCATE (oint, owork, sint)
    1239              : 
    1240              :       END DO
    1241         3578 :       CALL neighbor_list_iterator_release(nl_iterator)
    1242              : 
    1243         7156 :       DO i = 1, SIZE(matrix_h, 1)
    1244        54370 :          DO img = 1, nimg
    1245        47214 :             CALL dbcsr_finalize(matrix_h(i, img)%matrix)
    1246        50792 :             CALL dbcsr_finalize(matrix_s(i, img)%matrix)
    1247              :          END DO
    1248              :       END DO
    1249              : 
    1250         3578 :       kf = xtb_control%kf
    1251         3578 :       enscale = xtb_control%enscale
    1252         3578 :       erep = 0.0_dp
    1253         3578 :       CALL repulsive_potential(qs_env, erep, kf, enscale, calculate_forces)
    1254              : 
    1255         3578 :       exb = 0.0_dp
    1256         3578 :       IF (xb_inter) THEN
    1257         3578 :          CALL xb_interaction(qs_env, exb, calculate_forces)
    1258              :       END IF
    1259              : 
    1260         3578 :       enonbonded = 0.0_dp
    1261         3578 :       IF (do_nonbonded) THEN
    1262              :          ! nonbonded interactions
    1263           34 :          NULLIFY (sab_xtb_nonbond)
    1264           34 :          CALL get_qs_env(qs_env=qs_env, sab_xtb_nonbond=sab_xtb_nonbond)
    1265              :          CALL nonbonded_correction(enonbonded, force, qs_env, xtb_control, sab_xtb_nonbond, &
    1266           34 :                                    atomic_kind_set, calculate_forces, use_virial, virial, atprop, atom_of_kind)
    1267              :       END IF
    1268              : 
    1269              :       ! set repulsive energy
    1270         3578 :       erep = erep + exb + enonbonded
    1271         3578 :       IF (xb_inter) THEN
    1272         3578 :          CALL para_env%sum(exb)
    1273         3578 :          energy%xtb_xb_inter = exb
    1274              :       END IF
    1275         3578 :       IF (do_nonbonded) THEN
    1276           34 :          CALL para_env%sum(enonbonded)
    1277           34 :          energy%xtb_nonbonded = enonbonded
    1278              :       END IF
    1279         3578 :       CALL para_env%sum(erep)
    1280         3578 :       energy%repulsive = erep
    1281              : 
    1282              :       ! deallocate coordination numbers
    1283         3578 :       CALL cnumber_release(cnumbers, dcnum, calculate_forces)
    1284              : 
    1285              :       ! deallocate Huckel parameters
    1286         3578 :       DEALLOCATE (huckel)
    1287         3578 :       IF (calculate_forces) THEN
    1288          498 :          DEALLOCATE (dhuckel)
    1289              :       END IF
    1290              :       ! deallocate KAB parameters
    1291         3578 :       DEALLOCATE (kijab)
    1292              : 
    1293              :       ! AO matrix outputs
    1294         3578 :       CALL ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
    1295              : 
    1296         3578 :       DEALLOCATE (basis_set_list)
    1297         3578 :       IF (calculate_forces) THEN
    1298          498 :          IF (SIZE(matrix_p, 1) == 2) THEN
    1299          432 :             DO img = 1, nimg
    1300              :                CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, alpha_scalar=1.0_dp, &
    1301          432 :                               beta_scalar=-1.0_dp)
    1302              :             END DO
    1303              :          END IF
    1304              :       END IF
    1305              : 
    1306         3578 :       CALL timestop(handle)
    1307              : 
    1308        10734 :    END SUBROUTINE build_gfn1_xtb_matrices
    1309              : 
    1310              : ! **************************************************************************************************
    1311              : !> \brief ...
    1312              : !> \param qs_env ...
    1313              : !> \param matrix_h ...
    1314              : !> \param matrix_s ...
    1315              : !> \param calculate_forces ...
    1316              : ! **************************************************************************************************
    1317        10448 :    SUBROUTINE ao_matrix_output(qs_env, matrix_h, matrix_s, calculate_forces)
    1318              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1319              :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_h, matrix_s
    1320              :       LOGICAL, INTENT(IN)                                :: calculate_forces
    1321              : 
    1322              :       INTEGER                                            :: after, i, img, iw, nimg
    1323              :       LOGICAL                                            :: norml1, norml2, omit_headers, use_arnoldi
    1324              :       REAL(KIND=dp), DIMENSION(2)                        :: condnum
    1325              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
    1326              :       TYPE(cp_logger_type), POINTER                      :: logger
    1327              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1328              : 
    1329         5224 :       logger => cp_get_default_logger()
    1330              : 
    1331         5224 :       CALL get_qs_env(qs_env, para_env=para_env)
    1332         5224 :       nimg = SIZE(matrix_h, 2)
    1333         5224 :       CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%OMIT_HEADERS", l_val=omit_headers)
    1334         5224 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
    1335              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN"), cp_p_file)) THEN
    1336              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN", &
    1337            0 :                                    extension=".Log")
    1338            0 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
    1339            0 :          after = MIN(MAX(after, 1), 16)
    1340            0 :          DO img = 1, nimg
    1341              :             CALL cp_dbcsr_write_sparse_matrix(matrix_h(1, img)%matrix, 4, after, qs_env, para_env, &
    1342            0 :                                               output_unit=iw, omit_headers=omit_headers)
    1343              :          END DO
    1344            0 :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/CORE_HAMILTONIAN")
    1345              :       END IF
    1346              : 
    1347         5224 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, &
    1348              :                                            qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP"), cp_p_file)) THEN
    1349              :          iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP", &
    1350            0 :                                    extension=".Log")
    1351            0 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%AO_MATRICES%NDIGITS", i_val=after)
    1352            0 :          after = MIN(MAX(after, 1), 16)
    1353            0 :          DO img = 1, nimg
    1354              :             CALL cp_dbcsr_write_sparse_matrix(matrix_s(1, img)%matrix, 4, after, qs_env, para_env, &
    1355            0 :                                               output_unit=iw, omit_headers=omit_headers)
    1356            0 :             IF (BTEST(cp_print_key_should_output(logger%iter_info, &
    1357            0 :                                                  qs_env%input, "DFT%PRINT%AO_MATRICES/DERIVATIVES"), cp_p_file)) THEN
    1358            0 :                DO i = 2, SIZE(matrix_s, 1)
    1359              :                   CALL cp_dbcsr_write_sparse_matrix(matrix_s(i, img)%matrix, 4, after, qs_env, para_env, &
    1360            0 :                                                     output_unit=iw, omit_headers=omit_headers)
    1361              :                END DO
    1362              :             END IF
    1363              :          END DO
    1364            0 :          CALL cp_print_key_finished_output(iw, logger, qs_env%input, "DFT%PRINT%AO_MATRICES/OVERLAP")
    1365              :       END IF
    1366              : 
    1367              :       ! *** Overlap condition number
    1368         5224 :       IF (.NOT. calculate_forces) THEN
    1369         4684 :          IF (cp_print_key_should_output(logger%iter_info, qs_env%input, &
    1370              :                                         "DFT%PRINT%OVERLAP_CONDITION") .NE. 0) THEN
    1371              :             iw = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%OVERLAP_CONDITION", &
    1372            4 :                                       extension=".Log")
    1373            4 :             CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%1-NORM", l_val=norml1)
    1374            4 :             CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%DIAGONALIZATION", l_val=norml2)
    1375            4 :             CALL section_vals_val_get(qs_env%input, "DFT%PRINT%OVERLAP_CONDITION%ARNOLDI", l_val=use_arnoldi)
    1376            4 :             CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
    1377            4 :             CALL overlap_condnum(matrix_s, condnum, iw, norml1, norml2, use_arnoldi, blacs_env)
    1378              :          END IF
    1379              :       END IF
    1380              : 
    1381         5224 :    END SUBROUTINE ao_matrix_output
    1382              : 
    1383              : END MODULE xtb_matrices
    1384              : 
        

Generated by: LCOV version 2.0-1