LCOV - code coverage report
Current view: top level - src - xtb_ehess_force.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:58e3e09) Lines: 237 258 91.9 %
Date: 2024-03-29 07:50:05 Functions: 2 2 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Calculation of forces for Coulomb contributions in response xTB
      10             : !> \author JGH
      11             : ! **************************************************************************************************
      12             : MODULE xtb_ehess_force
      13             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14             :                                               get_atomic_kind_set
      15             :    USE cell_types,                      ONLY: cell_type,&
      16             :                                               get_cell,&
      17             :                                               pbc
      18             :    USE cp_control_types,                ONLY: dft_control_type,&
      19             :                                               xtb_control_type
      20             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      21             :                                               cp_logger_get_default_unit_nr,&
      22             :                                               cp_logger_type
      23             :    USE dbcsr_api,                       ONLY: &
      24             :         dbcsr_add, dbcsr_get_block_p, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      25             :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_type
      26             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      27             :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      28             :                                               ewald_environment_type
      29             :    USE ewald_methods_tb,                ONLY: tb_ewald_overlap,&
      30             :                                               tb_spme_zforce
      31             :    USE ewald_pw_types,                  ONLY: ewald_pw_type
      32             :    USE kinds,                           ONLY: dp
      33             :    USE mathconstants,                   ONLY: oorootpi,&
      34             :                                               pi
      35             :    USE message_passing,                 ONLY: mp_para_env_type
      36             :    USE particle_types,                  ONLY: particle_type
      37             :    USE pw_poisson_types,                ONLY: do_ewald_ewald,&
      38             :                                               do_ewald_none,&
      39             :                                               do_ewald_pme,&
      40             :                                               do_ewald_spme
      41             :    USE qs_energy_types,                 ONLY: qs_energy_type
      42             :    USE qs_environment_types,            ONLY: get_qs_env,&
      43             :                                               qs_environment_type
      44             :    USE qs_force_types,                  ONLY: qs_force_type
      45             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      46             :                                               qs_kind_type
      47             :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      48             :                                               neighbor_list_iterate,&
      49             :                                               neighbor_list_iterator_create,&
      50             :                                               neighbor_list_iterator_p_type,&
      51             :                                               neighbor_list_iterator_release,&
      52             :                                               neighbor_list_set_p_type
      53             :    USE qs_rho_types,                    ONLY: qs_rho_type
      54             :    USE virial_types,                    ONLY: virial_type
      55             :    USE xtb_coulomb,                     ONLY: dgamma_rab_sr,&
      56             :                                               gamma_rab_sr
      57             :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
      58             :                                               xtb_atom_type
      59             : #include "./base/base_uses.f90"
      60             : 
      61             :    IMPLICIT NONE
      62             : 
      63             :    PRIVATE
      64             : 
      65             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_ehess_force'
      66             : 
      67             :    PUBLIC :: calc_xtb_ehess_force
      68             : 
      69             : ! **************************************************************************************************
      70             : 
      71             : CONTAINS
      72             : 
      73             : ! **************************************************************************************************
      74             : !> \brief ...
      75             : !> \param qs_env ...
      76             : !> \param matrix_p0 ...
      77             : !> \param matrix_p1 ...
      78             : !> \param charges0 ...
      79             : !> \param mcharge0 ...
      80             : !> \param charges1 ...
      81             : !> \param mcharge1 ...
      82             : !> \param debug_forces ...
      83             : ! **************************************************************************************************
      84          14 :    SUBROUTINE calc_xtb_ehess_force(qs_env, matrix_p0, matrix_p1, charges0, mcharge0, &
      85          14 :                                    charges1, mcharge1, debug_forces)
      86             : 
      87             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      88             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p0, matrix_p1
      89             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(in)         :: charges0
      90             :       REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: mcharge0
      91             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(in)         :: charges1
      92             :       REAL(KIND=dp), DIMENSION(:), INTENT(in)            :: mcharge1
      93             :       LOGICAL, INTENT(IN)                                :: debug_forces
      94             : 
      95             :       CHARACTER(len=*), PARAMETER :: routineN = 'calc_xtb_ehess_force'
      96             : 
      97             :       INTEGER :: atom_i, atom_j, blk, ewald_type, handle, i, ia, iatom, icol, ikind, iounit, irow, &
      98             :          j, jatom, jkind, la, lb, lmaxa, lmaxb, natom, natorb_a, natorb_b, ni, nimg, nj, nkind, &
      99             :          nmat, za, zb
     100          14 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     101             :       INTEGER, DIMENSION(25)                             :: laoa, laob
     102             :       INTEGER, DIMENSION(3)                              :: cellind, periodic
     103             :       LOGICAL                                            :: calculate_forces, defined, do_ewald, &
     104             :                                                             found, just_energy, use_virial
     105             :       REAL(KIND=dp)                                      :: alpha, deth, dr, etaa, etab, fi, gmij0, &
     106             :                                                             gmij1, kg, rcut, rcuta, rcutb
     107          14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: xgamma
     108          14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gammab, gcij0, gcij1, gmcharge0, &
     109          14 :                                                             gmcharge1
     110          14 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: gchrg0, gchrg1
     111             :       REAL(KIND=dp), DIMENSION(3)                        :: fij, fodeb, rij
     112             :       REAL(KIND=dp), DIMENSION(5)                        :: kappaa, kappab
     113          14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: dsblock, pblock0, pblock1, sblock
     114          14 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     115             :       TYPE(cell_type), POINTER                           :: cell
     116             :       TYPE(cp_logger_type), POINTER                      :: logger
     117             :       TYPE(dbcsr_iterator_type)                          :: iter
     118          14 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     119             :       TYPE(dft_control_type), POINTER                    :: dft_control
     120             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     121             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     122             :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     123             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     124             :       TYPE(neighbor_list_iterator_p_type), &
     125          14 :          DIMENSION(:), POINTER                           :: nl_iterator
     126             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     127          14 :          POINTER                                         :: n_list
     128          14 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     129             :       TYPE(qs_energy_type), POINTER                      :: energy
     130          14 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     131          14 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     132             :       TYPE(qs_rho_type), POINTER                         :: rho
     133             :       TYPE(virial_type), POINTER                         :: virial
     134             :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b, xtb_kind
     135             :       TYPE(xtb_control_type), POINTER                    :: xtb_control
     136             : 
     137          14 :       CALL timeset(routineN, handle)
     138             : 
     139          14 :       logger => cp_get_default_logger()
     140          14 :       IF (logger%para_env%is_source()) THEN
     141           7 :          iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     142             :       ELSE
     143             :          iounit = -1
     144             :       END IF
     145             : 
     146          14 :       CPASSERT(ASSOCIATED(matrix_p1))
     147             : 
     148             :       CALL get_qs_env(qs_env, &
     149             :                       qs_kind_set=qs_kind_set, &
     150             :                       particle_set=particle_set, &
     151             :                       cell=cell, &
     152             :                       rho=rho, &
     153             :                       energy=energy, &
     154             :                       virial=virial, &
     155          14 :                       dft_control=dft_control)
     156             : 
     157          14 :       xtb_control => dft_control%qs_control%xtb_control
     158             : 
     159          14 :       calculate_forces = .TRUE.
     160          14 :       just_energy = .FALSE.
     161          14 :       use_virial = .FALSE.
     162          14 :       nmat = 4
     163          14 :       nimg = dft_control%nimages
     164          14 :       IF (nimg > 1) THEN
     165           0 :          CPABORT('xTB-sTDA forces for k-points not available')
     166             :       END IF
     167             : 
     168          14 :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
     169          56 :       ALLOCATE (gchrg0(natom, 5, nmat))
     170        5030 :       gchrg0 = 0._dp
     171          42 :       ALLOCATE (gmcharge0(natom, nmat))
     172        1006 :       gmcharge0 = 0._dp
     173          56 :       ALLOCATE (gchrg1(natom, 5, nmat))
     174        5030 :       gchrg1 = 0._dp
     175          42 :       ALLOCATE (gmcharge1(natom, nmat))
     176        1006 :       gmcharge1 = 0._dp
     177             : 
     178             :       ! short range contribution (gamma)
     179             :       ! loop over all atom pairs (sab_xtbe)
     180          14 :       kg = xtb_control%kg
     181          14 :       NULLIFY (n_list)
     182          14 :       IF (xtb_control%old_coulomb_damping) THEN
     183           0 :          CALL get_qs_env(qs_env=qs_env, sab_orb=n_list)
     184             :       ELSE
     185          14 :          CALL get_qs_env(qs_env=qs_env, sab_xtbe=n_list)
     186             :       END IF
     187          14 :       CALL neighbor_list_iterator_create(nl_iterator, n_list)
     188       24708 :       DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     189             :          CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     190       24694 :                                 iatom=iatom, jatom=jatom, r=rij, cell=cellind)
     191       24694 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     192       24694 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined, natorb=natorb_a)
     193       24694 :          IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     194       24694 :          CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     195       24694 :          CALL get_xtb_atom_param(xtb_atom_b, defined=defined, natorb=natorb_b)
     196       24694 :          IF (.NOT. defined .OR. natorb_b < 1) CYCLE
     197             :          ! atomic parameters
     198       24694 :          CALL get_xtb_atom_param(xtb_atom_a, eta=etaa, lmax=lmaxa, kappa=kappaa, rcut=rcuta)
     199       24694 :          CALL get_xtb_atom_param(xtb_atom_b, eta=etab, lmax=lmaxb, kappa=kappab, rcut=rcutb)
     200             :          ! gamma matrix
     201       24694 :          ni = lmaxa + 1
     202       24694 :          nj = lmaxb + 1
     203       98776 :          ALLOCATE (gammab(ni, nj))
     204       24694 :          rcut = rcuta + rcutb
     205       98776 :          dr = SQRT(SUM(rij(:)**2))
     206       24694 :          CALL gamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
     207      250847 :          gchrg0(iatom, 1:ni, 1) = gchrg0(iatom, 1:ni, 1) + MATMUL(gammab, charges0(jatom, 1:nj))
     208      250847 :          gchrg1(iatom, 1:ni, 1) = gchrg1(iatom, 1:ni, 1) + MATMUL(gammab, charges1(jatom, 1:nj))
     209       24694 :          IF (iatom /= jatom) THEN
     210      285738 :             gchrg0(jatom, 1:nj, 1) = gchrg0(jatom, 1:nj, 1) + MATMUL(charges0(iatom, 1:ni), gammab)
     211      285738 :             gchrg1(jatom, 1:nj, 1) = gchrg1(jatom, 1:nj, 1) + MATMUL(charges1(iatom, 1:ni), gammab)
     212             :          END IF
     213       24694 :          IF (dr > 1.e-6_dp) THEN
     214       24577 :             CALL dgamma_rab_sr(gammab, dr, ni, kappaa, etaa, nj, kappab, etab, kg, rcut)
     215       98308 :             DO i = 1, 3
     216             :                gchrg0(iatom, 1:ni, i + 1) = gchrg0(iatom, 1:ni, i + 1) &
     217      748626 :                                             + MATMUL(gammab, charges0(jatom, 1:nj))*rij(i)/dr
     218             :                gchrg1(iatom, 1:ni, i + 1) = gchrg1(iatom, 1:ni, i + 1) &
     219      748626 :                                             + MATMUL(gammab, charges1(jatom, 1:nj))*rij(i)/dr
     220       98308 :                IF (iatom /= jatom) THEN
     221             :                   gchrg0(jatom, 1:nj, i + 1) = gchrg0(jatom, 1:nj, i + 1) &
     222      857214 :                                                - MATMUL(charges0(iatom, 1:ni), gammab)*rij(i)/dr
     223             :                   gchrg1(jatom, 1:nj, i + 1) = gchrg1(jatom, 1:nj, i + 1) &
     224      857214 :                                                - MATMUL(charges1(iatom, 1:ni), gammab)*rij(i)/dr
     225             :                END IF
     226             :             END DO
     227             :          END IF
     228       98776 :          DEALLOCATE (gammab)
     229             :       END DO
     230          14 :       CALL neighbor_list_iterator_release(nl_iterator)
     231             : 
     232             :       ! 1/R contribution
     233             : 
     234          14 :       IF (xtb_control%coulomb_lr) THEN
     235          14 :          do_ewald = xtb_control%do_ewald
     236          14 :          IF (do_ewald) THEN
     237             :             ! Ewald sum
     238           6 :             NULLIFY (ewald_env, ewald_pw)
     239             :             CALL get_qs_env(qs_env=qs_env, &
     240           6 :                             ewald_env=ewald_env, ewald_pw=ewald_pw)
     241           6 :             CALL get_cell(cell=cell, periodic=periodic, deth=deth)
     242           6 :             CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type)
     243           6 :             CALL get_qs_env(qs_env=qs_env, sab_tbe=n_list)
     244           6 :             CALL tb_ewald_overlap(gmcharge0, mcharge0, alpha, n_list, virial, use_virial)
     245           6 :             CALL tb_ewald_overlap(gmcharge1, mcharge1, alpha, n_list, virial, use_virial)
     246           0 :             SELECT CASE (ewald_type)
     247             :             CASE DEFAULT
     248           0 :                CPABORT("Invalid Ewald type")
     249             :             CASE (do_ewald_none)
     250           0 :                CPABORT("Not allowed with DFTB")
     251             :             CASE (do_ewald_ewald)
     252           0 :                CPABORT("Standard Ewald not implemented in DFTB")
     253             :             CASE (do_ewald_pme)
     254           0 :                CPABORT("PME not implemented in DFTB")
     255             :             CASE (do_ewald_spme)
     256           6 :                CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge0, mcharge0)
     257          12 :                CALL tb_spme_zforce(ewald_env, ewald_pw, particle_set, cell, gmcharge1, mcharge1)
     258             :             END SELECT
     259             :          ELSE
     260             :             ! direct sum
     261           8 :             CALL get_qs_env(qs_env=qs_env, local_particles=local_particles)
     262          28 :             DO ikind = 1, SIZE(local_particles%n_el)
     263          42 :                DO ia = 1, local_particles%n_el(ikind)
     264          14 :                   iatom = local_particles%list(ikind)%array(ia)
     265          52 :                   DO jatom = 1, iatom - 1
     266          72 :                      rij = particle_set(iatom)%r - particle_set(jatom)%r
     267          72 :                      rij = pbc(rij, cell)
     268          72 :                      dr = SQRT(SUM(rij(:)**2))
     269          32 :                      IF (dr > 1.e-6_dp) THEN
     270          18 :                         gmcharge0(iatom, 1) = gmcharge0(iatom, 1) + mcharge0(jatom)/dr
     271          18 :                         gmcharge0(jatom, 1) = gmcharge0(jatom, 1) + mcharge0(iatom)/dr
     272          18 :                         gmcharge1(iatom, 1) = gmcharge1(iatom, 1) + mcharge1(jatom)/dr
     273          18 :                         gmcharge1(jatom, 1) = gmcharge1(jatom, 1) + mcharge1(iatom)/dr
     274          72 :                         DO i = 2, nmat
     275          54 :                            gmcharge0(iatom, i) = gmcharge0(iatom, i) + rij(i - 1)*mcharge0(jatom)/dr**3
     276          54 :                            gmcharge0(jatom, i) = gmcharge0(jatom, i) - rij(i - 1)*mcharge0(iatom)/dr**3
     277          54 :                            gmcharge1(iatom, i) = gmcharge1(iatom, i) + rij(i - 1)*mcharge1(jatom)/dr**3
     278          72 :                            gmcharge1(jatom, i) = gmcharge1(jatom, i) - rij(i - 1)*mcharge1(iatom)/dr**3
     279             :                         END DO
     280             :                      END IF
     281             :                   END DO
     282             :                END DO
     283             :             END DO
     284             :             CPASSERT(.NOT. use_virial)
     285             :          END IF
     286             :       END IF
     287             : 
     288             :       ! global sum of gamma*p arrays
     289             :       CALL get_qs_env(qs_env=qs_env, &
     290             :                       atomic_kind_set=atomic_kind_set, &
     291          14 :                       force=force, para_env=para_env)
     292          14 :       CALL para_env%sum(gmcharge0(:, 1))
     293          14 :       CALL para_env%sum(gchrg0(:, :, 1))
     294          14 :       CALL para_env%sum(gmcharge1(:, 1))
     295          14 :       CALL para_env%sum(gchrg1(:, :, 1))
     296             : 
     297          14 :       IF (xtb_control%coulomb_lr) THEN
     298          14 :          IF (do_ewald) THEN
     299             :             ! add self charge interaction and background charge contribution
     300         212 :             gmcharge0(:, 1) = gmcharge0(:, 1) - 2._dp*alpha*oorootpi*mcharge0(:)
     301          12 :             IF (ANY(periodic(:) == 1)) THEN
     302         202 :                gmcharge0(:, 1) = gmcharge0(:, 1) - pi/alpha**2/deth
     303             :             END IF
     304         212 :             gmcharge1(:, 1) = gmcharge1(:, 1) - 2._dp*alpha*oorootpi*mcharge1(:)
     305          12 :             IF (ANY(periodic(:) == 1)) THEN
     306         202 :                gmcharge1(:, 1) = gmcharge1(:, 1) - pi/alpha**2/deth
     307             :             END IF
     308             :          END IF
     309             :       END IF
     310             : 
     311             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     312             :                                kind_of=kind_of, &
     313          14 :                                atom_of_kind=atom_of_kind)
     314             : 
     315          14 :       IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     316         248 :       DO iatom = 1, natom
     317         234 :          ikind = kind_of(iatom)
     318         234 :          atom_i = atom_of_kind(iatom)
     319         234 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     320         234 :          CALL get_xtb_atom_param(xtb_kind, lmax=ni)
     321         234 :          ni = ni + 1
     322             :          ! short range
     323         234 :          fij = 0.0_dp
     324         936 :          DO i = 1, 3
     325             :             fij(i) = SUM(charges0(iatom, 1:ni)*gchrg1(iatom, 1:ni, i + 1)) + &
     326        2832 :                      SUM(charges1(iatom, 1:ni)*gchrg0(iatom, 1:ni, i + 1))
     327             :          END DO
     328         234 :          force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
     329         234 :          force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
     330         234 :          force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
     331             :          ! long range
     332         234 :          fij = 0.0_dp
     333         936 :          DO i = 1, 3
     334             :             fij(i) = gmcharge1(iatom, i + 1)*mcharge0(iatom) + &
     335         936 :                      gmcharge0(iatom, i + 1)*mcharge1(iatom)
     336             :          END DO
     337         234 :          force(ikind)%rho_elec(1, atom_i) = force(ikind)%rho_elec(1, atom_i) - fij(1)
     338         234 :          force(ikind)%rho_elec(2, atom_i) = force(ikind)%rho_elec(2, atom_i) - fij(2)
     339         482 :          force(ikind)%rho_elec(3, atom_i) = force(ikind)%rho_elec(3, atom_i) - fij(3)
     340             :       END DO
     341          14 :       IF (debug_forces) THEN
     342           0 :          fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     343           0 :          CALL para_env%sum(fodeb)
     344           0 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dH[Pz]    ", fodeb
     345             :       END IF
     346             : 
     347          14 :       CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
     348             : 
     349          14 :       IF (SIZE(matrix_p0) == 2) THEN
     350             :          CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
     351           0 :                         alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     352             :          CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
     353           0 :                         alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
     354             :       END IF
     355             : 
     356             :       ! no k-points; all matrices have been transformed to periodic bsf
     357          14 :       IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     358          14 :       CALL dbcsr_iterator_start(iter, matrix_s(1, 1)%matrix)
     359        4695 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     360        4681 :          CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
     361        4681 :          ikind = kind_of(irow)
     362        4681 :          jkind = kind_of(icol)
     363             : 
     364             :          ! atomic parameters
     365        4681 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     366        4681 :          CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     367        4681 :          CALL get_xtb_atom_param(xtb_atom_a, z=za, lao=laoa)
     368        4681 :          CALL get_xtb_atom_param(xtb_atom_b, z=zb, lao=laob)
     369             : 
     370        4681 :          ni = SIZE(sblock, 1)
     371        4681 :          nj = SIZE(sblock, 2)
     372       18724 :          ALLOCATE (gcij0(ni, nj))
     373       18724 :          ALLOCATE (gcij1(ni, nj))
     374       19209 :          DO i = 1, ni
     375       52401 :             DO j = 1, nj
     376       33192 :                la = laoa(i) + 1
     377       33192 :                lb = laob(j) + 1
     378       33192 :                gcij0(i, j) = 0.5_dp*(gchrg0(irow, la, 1) + gchrg0(icol, lb, 1))
     379       47720 :                gcij1(i, j) = 0.5_dp*(gchrg1(irow, la, 1) + gchrg1(icol, lb, 1))
     380             :             END DO
     381             :          END DO
     382        4681 :          gmij0 = 0.5_dp*(gmcharge0(irow, 1) + gmcharge0(icol, 1))
     383        4681 :          gmij1 = 0.5_dp*(gmcharge1(irow, 1) + gmcharge1(icol, 1))
     384        4681 :          atom_i = atom_of_kind(irow)
     385        4681 :          atom_j = atom_of_kind(icol)
     386        4681 :          NULLIFY (pblock0)
     387             :          CALL dbcsr_get_block_p(matrix=matrix_p0(1)%matrix, &
     388        4681 :                                 row=irow, col=icol, block=pblock0, found=found)
     389        4681 :          CPASSERT(found)
     390        4681 :          NULLIFY (pblock1)
     391             :          CALL dbcsr_get_block_p(matrix=matrix_p1(1)%matrix, &
     392        4681 :                                 row=irow, col=icol, block=pblock1, found=found)
     393        4681 :          CPASSERT(found)
     394       18724 :          DO i = 1, 3
     395       14043 :             NULLIFY (dsblock)
     396             :             CALL dbcsr_get_block_p(matrix=matrix_s(1 + i, 1)%matrix, &
     397       14043 :                                    row=irow, col=icol, block=dsblock, found=found)
     398       14043 :             CPASSERT(found)
     399             :             ! short range
     400      275571 :             fi = -2.0_dp*SUM(pblock0*dsblock*gcij1) - 2.0_dp*SUM(pblock1*dsblock*gcij0)
     401       14043 :             force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     402       14043 :             force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     403             :             ! long range
     404      275571 :             fi = -2.0_dp*gmij1*SUM(pblock0*dsblock) - 2.0_dp*gmij0*SUM(pblock1*dsblock)
     405       14043 :             force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     406       32767 :             force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     407             :          END DO
     408       23419 :          DEALLOCATE (gcij0, gcij1)
     409             :       END DO
     410          14 :       CALL dbcsr_iterator_stop(iter)
     411          14 :       IF (debug_forces) THEN
     412           0 :          fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     413           0 :          CALL para_env%sum(fodeb)
     414           0 :          IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H[P]*dS  ", fodeb
     415             :       END IF
     416             : 
     417          14 :       IF (xtb_control%tb3_interaction) THEN
     418          14 :          CALL get_qs_env(qs_env, nkind=nkind)
     419          42 :          ALLOCATE (xgamma(nkind))
     420          48 :          DO ikind = 1, nkind
     421          34 :             CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_kind)
     422          48 :             CALL get_xtb_atom_param(xtb_kind, xgamma=xgamma(ikind))
     423             :          END DO
     424             :          ! Diagonal 3rd order correction (DFTB3)
     425          14 :          IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
     426             :          CALL dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
     427          14 :                                            matrix_p0(1)%matrix, matrix_p1(1)%matrix, xgamma)
     428          14 :          IF (debug_forces) THEN
     429           0 :             fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
     430           0 :             CALL para_env%sum(fodeb)
     431           0 :             IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Pz*H3[P]    ", fodeb
     432             :          END IF
     433          14 :          DEALLOCATE (xgamma)
     434             :       END IF
     435             : 
     436          14 :       IF (SIZE(matrix_p0) == 2) THEN
     437             :          CALL dbcsr_add(matrix_p0(1)%matrix, matrix_p0(2)%matrix, &
     438           0 :                         alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
     439             :          CALL dbcsr_add(matrix_p1(1)%matrix, matrix_p1(2)%matrix, &
     440           0 :                         alpha_scalar=1.0_dp, beta_scalar=-1.0_dp)
     441             :       END IF
     442             : 
     443             :       ! QMMM
     444          14 :       IF (qs_env%qmmm .AND. qs_env%qmmm_periodic) THEN
     445           0 :          CPABORT("Not Available")
     446             :       END IF
     447             : 
     448          14 :       DEALLOCATE (gmcharge0, gchrg0, gmcharge1, gchrg1)
     449             : 
     450          14 :       CALL timestop(handle)
     451             : 
     452          28 :    END SUBROUTINE calc_xtb_ehess_force
     453             : 
     454             : ! **************************************************************************************************
     455             : !> \brief ...
     456             : !> \param qs_env ...
     457             : !> \param mcharge0 ...
     458             : !> \param mcharge1 ...
     459             : !> \param matrixp0 ...
     460             : !> \param matrixp1 ...
     461             : !> \param xgamma ...
     462             : ! **************************************************************************************************
     463          14 :    SUBROUTINE dftb3_diagonal_hessian_force(qs_env, mcharge0, mcharge1, &
     464          14 :                                            matrixp0, matrixp1, xgamma)
     465             : 
     466             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     467             :       REAL(dp), DIMENSION(:)                             :: mcharge0, mcharge1
     468             :       TYPE(dbcsr_type), POINTER                          :: matrixp0, matrixp1
     469             :       REAL(dp), DIMENSION(:)                             :: xgamma
     470             : 
     471             :       CHARACTER(len=*), PARAMETER :: routineN = 'dftb3_diagonal_hessian_force'
     472             : 
     473             :       INTEGER                                            :: atom_i, atom_j, blk, handle, i, icol, &
     474             :                                                             ikind, irow, jkind
     475          14 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     476             :       LOGICAL                                            :: found
     477             :       REAL(KIND=dp)                                      :: fi, gmijp, gmijq, ui, uj
     478          14 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: dsblock, p0block, p1block, sblock
     479          14 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     480             :       TYPE(dbcsr_iterator_type)                          :: iter
     481          14 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_s
     482          14 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     483             : 
     484          14 :       CALL timeset(routineN, handle)
     485          14 :       CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s)
     486          14 :       CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set)
     487             :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     488          14 :                                kind_of=kind_of, atom_of_kind=atom_of_kind)
     489          14 :       CALL get_qs_env(qs_env=qs_env, force=force)
     490             :       ! no k-points; all matrices have been transformed to periodic bsf
     491          14 :       CALL dbcsr_iterator_start(iter, matrix_s(1)%matrix)
     492        4695 :       DO WHILE (dbcsr_iterator_blocks_left(iter))
     493        4681 :          CALL dbcsr_iterator_next_block(iter, irow, icol, sblock, blk)
     494        4681 :          ikind = kind_of(irow)
     495        4681 :          atom_i = atom_of_kind(irow)
     496        4681 :          ui = xgamma(ikind)
     497        4681 :          jkind = kind_of(icol)
     498        4681 :          atom_j = atom_of_kind(icol)
     499        4681 :          uj = xgamma(jkind)
     500             :          !
     501        4681 :          gmijp = ui*mcharge0(irow)*mcharge1(irow) + uj*mcharge0(icol)*mcharge1(icol)
     502        4681 :          gmijq = 0.5_dp*ui*mcharge0(irow)**2 + 0.5_dp*uj*mcharge0(icol)**2
     503             :          !
     504        4681 :          NULLIFY (p0block)
     505             :          CALL dbcsr_get_block_p(matrix=matrixp0, &
     506        4681 :                                 row=irow, col=icol, block=p0block, found=found)
     507        4681 :          CPASSERT(found)
     508        4681 :          NULLIFY (p1block)
     509             :          CALL dbcsr_get_block_p(matrix=matrixp1, &
     510        4681 :                                 row=irow, col=icol, block=p1block, found=found)
     511        4681 :          CPASSERT(found)
     512       18738 :          DO i = 1, 3
     513       14043 :             NULLIFY (dsblock)
     514             :             CALL dbcsr_get_block_p(matrix=matrix_s(1 + i)%matrix, &
     515       14043 :                                    row=irow, col=icol, block=dsblock, found=found)
     516       14043 :             CPASSERT(found)
     517      275571 :             fi = gmijp*SUM(p0block*dsblock) + gmijq*SUM(p1block*dsblock)
     518       14043 :             force(ikind)%rho_elec(i, atom_i) = force(ikind)%rho_elec(i, atom_i) + fi
     519       32767 :             force(jkind)%rho_elec(i, atom_j) = force(jkind)%rho_elec(i, atom_j) - fi
     520             :          END DO
     521             :       END DO
     522          14 :       CALL dbcsr_iterator_stop(iter)
     523             : 
     524          14 :       CALL timestop(handle)
     525             : 
     526          28 :    END SUBROUTINE dftb3_diagonal_hessian_force
     527             : 
     528      389834 : END MODULE xtb_ehess_force
     529             : 

Generated by: LCOV version 1.15