LCOV - code coverage report
Current view: top level - src - xtb_eeq.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 94.3 % 264 249
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 3 3

            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 charge equilibration in xTB
      10              : !> \author JGH
      11              : ! **************************************************************************************************
      12              : MODULE xtb_eeq
      13              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      14              :                                               get_atomic_kind_set
      15              :    USE atprop_types,                    ONLY: atprop_array_init,&
      16              :                                               atprop_type
      17              :    USE cell_types,                      ONLY: cell_type,&
      18              :                                               pbc
      19              :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      20              :    USE cp_control_types,                ONLY: dft_control_type,&
      21              :                                               xtb_control_type
      22              :    USE cp_log_handling,                 ONLY: cp_logger_get_default_unit_nr
      23              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      24              :    USE eeq_input,                       ONLY: eeq_solver_type
      25              :    USE eeq_method,                      ONLY: eeq_efield_energy,&
      26              :                                               eeq_efield_force_loc,&
      27              :                                               eeq_efield_force_periodic,&
      28              :                                               eeq_efield_pot,&
      29              :                                               eeq_solver
      30              :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      31              :                                               ewald_environment_type
      32              :    USE ewald_pw_types,                  ONLY: ewald_pw_type
      33              :    USE kinds,                           ONLY: dp
      34              :    USE mathconstants,                   ONLY: oorootpi
      35              :    USE message_passing,                 ONLY: mp_para_env_type
      36              :    USE particle_types,                  ONLY: particle_type
      37              :    USE qs_dispersion_cnum,              ONLY: dcnum_type
      38              :    USE qs_environment_types,            ONLY: get_qs_env,&
      39              :                                               qs_environment_type
      40              :    USE qs_force_types,                  ONLY: qs_force_type
      41              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      42              :                                               qs_kind_type
      43              :    USE qs_neighbor_list_types,          ONLY: get_iterator_info,&
      44              :                                               neighbor_list_iterate,&
      45              :                                               neighbor_list_iterator_create,&
      46              :                                               neighbor_list_iterator_p_type,&
      47              :                                               neighbor_list_iterator_release,&
      48              :                                               neighbor_list_set_p_type
      49              :    USE spme,                            ONLY: spme_forces,&
      50              :                                               spme_virial
      51              :    USE virial_methods,                  ONLY: virial_pair_force
      52              :    USE virial_types,                    ONLY: virial_type
      53              :    USE xtb_types,                       ONLY: get_xtb_atom_param,&
      54              :                                               xtb_atom_type
      55              : #include "./base/base_uses.f90"
      56              : 
      57              :    IMPLICIT NONE
      58              : 
      59              :    PRIVATE
      60              : 
      61              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xtb_eeq'
      62              : 
      63              :    PUBLIC :: xtb_eeq_calculation, xtb_eeq_forces, xtb_eeq_lagrange
      64              : 
      65              : CONTAINS
      66              : 
      67              : ! **************************************************************************************************
      68              : !> \brief ...
      69              : !> \param qs_env ...
      70              : !> \param charges ...
      71              : !> \param cnumbers ...
      72              : !> \param eeq_sparam ...
      73              : !> \param eeq_energy ...
      74              : !> \param ef_energy ...
      75              : !> \param lambda ...
      76              : ! **************************************************************************************************
      77         1702 :    SUBROUTINE xtb_eeq_calculation(qs_env, charges, cnumbers, &
      78              :                                   eeq_sparam, eeq_energy, ef_energy, lambda)
      79              : 
      80              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      81              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: charges
      82              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: cnumbers
      83              :       TYPE(eeq_solver_type), INTENT(IN)                  :: eeq_sparam
      84              :       REAL(KIND=dp), INTENT(INOUT)                       :: eeq_energy, ef_energy, lambda
      85              : 
      86              :       CHARACTER(len=*), PARAMETER :: routineN = 'xtb_eeq_calculation'
      87              : 
      88              :       INTEGER                                            :: enshift_type, handle, iatom, ikind, &
      89              :                                                             iunit, jkind, natom, nkind
      90         1702 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
      91              :       LOGICAL                                            :: defined, do_ewald
      92              :       REAL(KIND=dp)                                      :: ala, alb, cn, esg, gama, kappa, scn, &
      93              :                                                             sgamma, totalcharge, xi
      94         1702 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: chia, efr, gam
      95              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gab
      96         1702 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      97              :       TYPE(atprop_type), POINTER                         :: atprop
      98              :       TYPE(cell_type), POINTER                           :: cell
      99              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     100              :       TYPE(dft_control_type), POINTER                    :: dft_control
     101              :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     102              :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     103              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     104         1702 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     105         1702 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     106              :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b
     107              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
     108              : 
     109         1702 :       CALL timeset(routineN, handle)
     110              : 
     111         1702 :       iunit = cp_logger_get_default_unit_nr()
     112              : 
     113              :       CALL get_qs_env(qs_env, &
     114              :                       qs_kind_set=qs_kind_set, &
     115              :                       atomic_kind_set=atomic_kind_set, &
     116              :                       particle_set=particle_set, &
     117              :                       cell=cell, &
     118              :                       atprop=atprop, &
     119         1702 :                       dft_control=dft_control)
     120         1702 :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
     121              : 
     122         1702 :       xtb_control => dft_control%qs_control%xtb_control
     123              : 
     124         1702 :       totalcharge = dft_control%charge
     125              : 
     126         1702 :       IF (atprop%energy) THEN
     127            0 :          CALL atprop_array_init(atprop%atecoul, natom)
     128              :       END IF
     129              : 
     130              :       ! gamma[a,b]
     131        10212 :       ALLOCATE (gab(nkind, nkind), gam(nkind))
     132        17090 :       gab = 0.0_dp
     133         5954 :       gam = 0.0_dp
     134         5954 :       DO ikind = 1, nkind
     135         4252 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     136         4252 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
     137         4252 :          IF (.NOT. defined) CYCLE
     138         4252 :          CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
     139         4252 :          gam(ikind) = gama
     140        21342 :          DO jkind = 1, nkind
     141        11136 :             CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     142        11136 :             CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
     143        11136 :             IF (.NOT. defined) CYCLE
     144        11136 :             CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
     145              :             !
     146        26524 :             gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
     147              :             !
     148              :          END DO
     149              :       END DO
     150              : 
     151              :       ! Chi[a,a]
     152         1702 :       enshift_type = xtb_control%enshift_type
     153         1702 :       IF (enshift_type == 0) THEN
     154            0 :          enshift_type = 2
     155            0 :          IF (ALL(cell%perd == 0)) enshift_type = 1
     156              :       END IF
     157         1702 :       sgamma = 8.0_dp ! see D4 for periodic systems paper
     158         1702 :       esg = 1.0_dp + EXP(sgamma)
     159         5106 :       ALLOCATE (chia(natom))
     160         1702 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
     161         9676 :       DO iatom = 1, natom
     162         7974 :          ikind = kind_of(iatom)
     163         7974 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     164         7974 :          CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
     165              :          !
     166         7974 :          IF (enshift_type == 1) THEN
     167         7974 :             scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
     168            0 :          ELSE IF (enshift_type == 2) THEN
     169            0 :             cn = cnumbers(iatom)/esg
     170            0 :             scn = LOG(esg/(esg - cnumbers(iatom)))
     171              :          ELSE
     172            0 :             CPABORT("Unknown enshift_type")
     173              :          END IF
     174        17650 :          chia(iatom) = xi - kappa*scn
     175              :          !
     176              :       END DO
     177              : 
     178         1702 :       ef_energy = 0.0_dp
     179         1702 :       IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
     180              :           dft_control%apply_efield_field) THEN
     181          984 :          ALLOCATE (efr(natom))
     182         1640 :          efr(1:natom) = 0.0_dp
     183          328 :          CALL eeq_efield_pot(qs_env, efr)
     184         3014 :          chia(1:natom) = chia(1:natom) + efr(1:natom)
     185              :       END IF
     186              : 
     187         1702 :       do_ewald = xtb_control%do_ewald
     188              : 
     189         1702 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
     190         1702 :       IF (do_ewald) THEN
     191              :          CALL get_qs_env(qs_env=qs_env, &
     192          756 :                          ewald_env=ewald_env, ewald_pw=ewald_pw)
     193              :          CALL eeq_solver(charges, lambda, eeq_energy, &
     194              :                          particle_set, kind_of, cell, chia, gam, gab, &
     195              :                          para_env, blacs_env, dft_control, eeq_sparam, &
     196              :                          totalcharge=totalcharge, ewald=do_ewald, &
     197          756 :                          ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
     198              :       ELSE
     199              :          CALL eeq_solver(charges, lambda, eeq_energy, &
     200              :                          particle_set, kind_of, cell, chia, gam, gab, &
     201              :                          para_env, blacs_env, dft_control, eeq_sparam, &
     202          946 :                          totalcharge=totalcharge, iounit=iunit)
     203              :       END IF
     204              : 
     205         1702 :       IF (dft_control%apply_period_efield .OR. dft_control%apply_efield .OR. &
     206              :           dft_control%apply_efield_field) THEN
     207          328 :          CALL eeq_efield_energy(qs_env, charges, ef_energy)
     208         1640 :          eeq_energy = eeq_energy - SUM(charges*efr)
     209          328 :          DEALLOCATE (efr)
     210              :       END IF
     211              : 
     212         1702 :       DEALLOCATE (gab, gam, chia)
     213              : 
     214         1702 :       CALL timestop(handle)
     215              : 
     216         3404 :    END SUBROUTINE xtb_eeq_calculation
     217              : 
     218              : ! **************************************************************************************************
     219              : !> \brief ...
     220              : !> \param qs_env ...
     221              : !> \param charges ...
     222              : !> \param dcharges ...
     223              : !> \param qlagrange ...
     224              : !> \param cnumbers ...
     225              : !> \param dcnum ...
     226              : !> \param eeq_sparam ...
     227              : ! **************************************************************************************************
     228           42 :    SUBROUTINE xtb_eeq_forces(qs_env, charges, dcharges, qlagrange, cnumbers, dcnum, eeq_sparam)
     229              : 
     230              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     231              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: charges, dcharges
     232              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: qlagrange
     233              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: cnumbers
     234              :       TYPE(dcnum_type), DIMENSION(:), INTENT(IN)         :: dcnum
     235              :       TYPE(eeq_solver_type), INTENT(IN)                  :: eeq_sparam
     236              : 
     237              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xtb_eeq_forces'
     238              : 
     239              :       INTEGER                                            :: atom_a, atom_b, atom_c, enshift_type, &
     240              :                                                             handle, i, ia, iatom, ikind, iunit, &
     241              :                                                             jatom, jkind, katom, kkind, natom, &
     242              :                                                             nkind
     243           42 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     244              :       LOGICAL                                            :: defined, do_ewald, use_virial
     245              :       REAL(KIND=dp)                                      :: ala, alb, alpha, cn, ctot, dr, dr2, drk, &
     246              :                                                             elag, esg, fe, gam2, gama, grc, kappa, &
     247              :                                                             qlam, qq, qq1, qq2, rcut, scn, sgamma, &
     248              :                                                             totalcharge, xi
     249              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: gam
     250           42 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: epforce, gab
     251              :       REAL(KIND=dp), DIMENSION(3)                        :: fdik, ri, rij, rik, rj
     252              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pvir
     253           42 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: chrgx, dchia, qlag
     254           42 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     255              :       TYPE(atprop_type), POINTER                         :: atprop
     256              :       TYPE(cell_type), POINTER                           :: cell
     257              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     258              :       TYPE(dft_control_type), POINTER                    :: dft_control
     259              :       TYPE(distribution_1d_type), POINTER                :: local_particles
     260              :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     261              :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     262              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     263              :       TYPE(neighbor_list_iterator_p_type), &
     264           42 :          DIMENSION(:), POINTER                           :: nl_iterator
     265              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     266           42 :          POINTER                                         :: sab_tbe
     267           42 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     268           42 :       TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
     269           42 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     270              :       TYPE(virial_type), POINTER                         :: virial
     271              :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b
     272              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
     273              : 
     274           42 :       CALL timeset(routineN, handle)
     275              : 
     276           42 :       iunit = cp_logger_get_default_unit_nr()
     277              : 
     278              :       CALL get_qs_env(qs_env, &
     279              :                       qs_kind_set=qs_kind_set, &
     280              :                       atomic_kind_set=atomic_kind_set, &
     281              :                       particle_set=particle_set, &
     282              :                       atprop=atprop, &
     283              :                       force=force, &
     284              :                       virial=virial, &
     285              :                       cell=cell, &
     286           42 :                       dft_control=dft_control)
     287           42 :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
     288           42 :       use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
     289              : 
     290           42 :       xtb_control => dft_control%qs_control%xtb_control
     291              : 
     292           42 :       totalcharge = dft_control%charge
     293              : 
     294              :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     295           42 :                                atom_of_kind=atom_of_kind, kind_of=kind_of)
     296              : 
     297              :       ! gamma[a,b]
     298          252 :       ALLOCATE (gab(nkind, nkind), gam(nkind))
     299          498 :       gab = 0.0_dp
     300          160 :       DO ikind = 1, nkind
     301          118 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     302          118 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
     303          118 :          IF (.NOT. defined) CYCLE
     304          118 :          CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
     305          118 :          gam(ikind) = gama
     306          616 :          DO jkind = 1, nkind
     307          338 :             CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     308          338 :             CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
     309          338 :             IF (.NOT. defined) CYCLE
     310          338 :             CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
     311              :             !
     312          794 :             gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
     313              :             !
     314              :          END DO
     315              :       END DO
     316              : 
     317          126 :       ALLOCATE (qlag(natom))
     318              : 
     319           42 :       do_ewald = xtb_control%do_ewald
     320              : 
     321           42 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
     322           42 :       IF (do_ewald) THEN
     323              :          CALL get_qs_env(qs_env=qs_env, &
     324           28 :                          ewald_env=ewald_env, ewald_pw=ewald_pw)
     325              :          CALL eeq_solver(qlag, qlam, elag, &
     326              :                          particle_set, kind_of, cell, -dcharges, gam, gab, &
     327              :                          para_env, blacs_env, dft_control, eeq_sparam, &
     328          172 :                          ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
     329              :       ELSE
     330              :          CALL eeq_solver(qlag, qlam, elag, &
     331              :                          particle_set, kind_of, cell, -dcharges, gam, gab, &
     332           70 :                          para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
     333              :       END IF
     334              : 
     335           42 :       enshift_type = xtb_control%enshift_type
     336           42 :       IF (enshift_type == 0) THEN
     337            0 :          enshift_type = 2
     338            0 :          IF (ALL(cell%perd == 0)) enshift_type = 1
     339              :       END IF
     340           42 :       sgamma = 8.0_dp ! see D4 for periodic systems paper
     341           42 :       esg = 1.0_dp + EXP(sgamma)
     342          168 :       ALLOCATE (chrgx(natom), dchia(natom))
     343          242 :       DO iatom = 1, natom
     344          200 :          ikind = kind_of(iatom)
     345          200 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     346          200 :          CALL get_xtb_atom_param(xtb_atom_a, xi=xi, kappa0=kappa)
     347              :          !
     348          200 :          ctot = 0.5_dp*(charges(iatom) - qlag(iatom))
     349          242 :          IF (enshift_type == 1) THEN
     350          200 :             scn = SQRT(cnumbers(iatom)) + 1.0e-14_dp
     351          200 :             dchia(iatom) = -ctot*kappa/scn
     352            0 :          ELSE IF (enshift_type == 2) THEN
     353            0 :             cn = cnumbers(iatom)
     354            0 :             scn = 1.0_dp/(esg - cn)
     355            0 :             dchia(iatom) = -ctot*kappa*scn
     356              :          ELSE
     357            0 :             CPABORT("Unknown enshift_type")
     358              :          END IF
     359              :       END DO
     360              : 
     361              :       ! Efield
     362           42 :       IF (dft_control%apply_period_efield) THEN
     363            8 :          CALL eeq_efield_force_periodic(qs_env, charges, qlag)
     364           34 :       ELSE IF (dft_control%apply_efield) THEN
     365            8 :          CALL eeq_efield_force_loc(qs_env, charges, qlag)
     366           26 :       ELSE IF (dft_control%apply_efield_field) THEN
     367            0 :          CPABORT("apply field")
     368              :       END IF
     369              : 
     370              :       ! Forces from q*X
     371              :       CALL get_qs_env(qs_env=qs_env, &
     372           42 :                       local_particles=local_particles)
     373          160 :       DO ikind = 1, nkind
     374          260 :          DO ia = 1, local_particles%n_el(ikind)
     375          100 :             iatom = local_particles%list(ikind)%array(ia)
     376          100 :             atom_a = atom_of_kind(iatom)
     377          688 :             DO i = 1, dcnum(iatom)%neighbors
     378          470 :                katom = dcnum(iatom)%nlist(i)
     379          470 :                kkind = kind_of(katom)
     380          470 :                atom_c = atom_of_kind(katom)
     381         1880 :                rik = dcnum(iatom)%rik(:, i)
     382         1880 :                drk = SQRT(SUM(rik(:)**2))
     383          570 :                IF (drk > 1.e-3_dp) THEN
     384         1880 :                   fdik(:) = dchia(iatom)*dcnum(iatom)%dvals(i)*rik(:)/drk
     385         1880 :                   force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
     386         1880 :                   force(kkind)%rho_elec(:, atom_c) = force(kkind)%rho_elec(:, atom_c) + fdik(:)
     387          470 :                   IF (use_virial) THEN
     388          374 :                      CALL virial_pair_force(virial%pv_virial, -1._dp, fdik, rik)
     389              :                   END IF
     390              :                END IF
     391              :             END DO
     392              :          END DO
     393              :       END DO
     394              : 
     395              :       ! Forces from (0.5*q+l)*dA/dR*q
     396           42 :       IF (do_ewald) THEN
     397           28 :          CALL get_qs_env(qs_env, sab_tbe=sab_tbe)
     398           28 :          CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
     399           28 :          rcut = 1.0_dp*rcut
     400           28 :          CALL neighbor_list_iterator_create(nl_iterator, sab_tbe)
     401        19968 :          DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
     402              :             CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, &
     403        19940 :                                    iatom=iatom, jatom=jatom, r=rij)
     404        19940 :             atom_a = atom_of_kind(iatom)
     405        19940 :             atom_b = atom_of_kind(jatom)
     406              :             !
     407        79760 :             dr2 = SUM(rij**2)
     408        19940 :             dr = SQRT(dr2)
     409        19940 :             IF (dr > rcut .OR. dr < 1.E-6_dp) CYCLE
     410         2309 :             fe = 1.0_dp
     411         2309 :             IF (iatom == jatom) fe = 0.5_dp
     412         2309 :             qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
     413         2309 :             gama = gab(ikind, jkind)
     414         2309 :             gam2 = gama*gama
     415              :             grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2 &
     416         2309 :                   - 2._dp*alpha*EXP(-alpha**2*dr2)*oorootpi/dr + erf(alpha*dr)/dr2
     417         2309 :             qq1 = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
     418         2309 :             qq2 = (0.5_dp*charges(jatom) - qlag(jatom))*charges(iatom)
     419         9236 :             fdik(:) = -qq1*grc*rij(:)/dr
     420         9236 :             force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
     421         9236 :             force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
     422         2309 :             IF (use_virial) THEN
     423         1554 :                CALL virial_pair_force(virial%pv_virial, fe, fdik, rij)
     424              :             END IF
     425         9236 :             fdik(:) = qq2*grc*rij(:)/dr
     426         9236 :             force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) - fdik(:)
     427         9236 :             force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) + fdik(:)
     428         2337 :             IF (use_virial) THEN
     429         1554 :                CALL virial_pair_force(virial%pv_virial, -fe, fdik, rij)
     430              :             END IF
     431              :          END DO
     432           28 :          CALL neighbor_list_iterator_release(nl_iterator)
     433              :       ELSE
     434           56 :          DO ikind = 1, nkind
     435           84 :             DO ia = 1, local_particles%n_el(ikind)
     436           28 :                iatom = local_particles%list(ikind)%array(ia)
     437           28 :                atom_a = atom_of_kind(iatom)
     438          112 :                ri(1:3) = particle_set(iatom)%r(1:3)
     439          182 :                DO jatom = 1, natom
     440          112 :                   IF (iatom == jatom) CYCLE
     441           84 :                   jkind = kind_of(jatom)
     442           84 :                   atom_b = atom_of_kind(jatom)
     443           84 :                   qq = (0.5_dp*charges(iatom) - qlag(iatom))*charges(jatom)
     444          336 :                   rj(1:3) = particle_set(jatom)%r(1:3)
     445          336 :                   rij(1:3) = ri(1:3) - rj(1:3)
     446          336 :                   rij = pbc(rij, cell)
     447          336 :                   dr2 = SUM(rij**2)
     448           84 :                   dr = SQRT(dr2)
     449           84 :                   gama = gab(ikind, jkind)
     450           84 :                   gam2 = gama*gama
     451           84 :                   grc = 2._dp*gama*EXP(-gam2*dr2)*oorootpi/dr - erf(gama*dr)/dr2
     452          336 :                   fdik(:) = qq*grc*rij(:)/dr
     453          336 :                   force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + fdik(:)
     454          364 :                   force(jkind)%rho_elec(:, atom_b) = force(jkind)%rho_elec(:, atom_b) - fdik(:)
     455              :                END DO
     456              :             END DO
     457              :          END DO
     458              :       END IF
     459              : 
     460              :       ! Forces from Ewald potential: (q+l)*A*q
     461           42 :       IF (do_ewald) THEN
     462           84 :          ALLOCATE (epforce(3, natom))
     463          604 :          epforce = 0.0_dp
     464          316 :          dchia = -charges + qlag
     465          172 :          chrgx = charges
     466              :          CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
     467           28 :                           particle_set, dchia, epforce)
     468          172 :          dchia = charges
     469          316 :          chrgx = qlag
     470              :          CALL spme_forces(ewald_env, ewald_pw, cell, particle_set, chrgx, &
     471           28 :                           particle_set, dchia, epforce)
     472          172 :          DO iatom = 1, natom
     473          144 :             ikind = kind_of(iatom)
     474          144 :             i = atom_of_kind(iatom)
     475          604 :             force(ikind)%rho_elec(:, i) = force(ikind)%rho_elec(:, i) + epforce(:, iatom)
     476              :          END DO
     477           28 :          DEALLOCATE (epforce)
     478              : 
     479              :          ! virial
     480           28 :          IF (use_virial) THEN
     481          154 :             chrgx = charges - qlag
     482           10 :             CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
     483          130 :             virial%pv_virial = virial%pv_virial + pvir
     484          154 :             chrgx = qlag
     485           10 :             CALL spme_virial(ewald_env, ewald_pw, particle_set, cell, chrgx, pvir)
     486          130 :             virial%pv_virial = virial%pv_virial - pvir
     487              :          END IF
     488              :       END IF
     489              : 
     490              :       ! return Lagrange multipliers
     491          242 :       qlagrange(1:natom) = qlag(1:natom)
     492              : 
     493           42 :       DEALLOCATE (gab, chrgx, dchia, qlag)
     494              : 
     495           42 :       CALL timestop(handle)
     496              : 
     497           84 :    END SUBROUTINE xtb_eeq_forces
     498              : 
     499              : ! **************************************************************************************************
     500              : !> \brief ...
     501              : !> \param qs_env ...
     502              : !> \param dcharges ...
     503              : !> \param qlagrange ...
     504              : !> \param eeq_sparam ...
     505              : ! **************************************************************************************************
     506           56 :    SUBROUTINE xtb_eeq_lagrange(qs_env, dcharges, qlagrange, eeq_sparam)
     507              : 
     508              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     509              :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: dcharges
     510              :       REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: qlagrange
     511              :       TYPE(eeq_solver_type), INTENT(IN)                  :: eeq_sparam
     512              : 
     513              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'xtb_eeq_lagrange'
     514              : 
     515              :       INTEGER                                            :: handle, ikind, iunit, jkind, natom, nkind
     516           56 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atom_of_kind, kind_of
     517              :       LOGICAL                                            :: defined, do_ewald
     518              :       REAL(KIND=dp)                                      :: ala, alb, elag, gama, qlam, totalcharge
     519              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: gam
     520              :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gab
     521           56 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: qlag
     522           56 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     523              :       TYPE(cell_type), POINTER                           :: cell
     524              :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
     525              :       TYPE(dft_control_type), POINTER                    :: dft_control
     526              :       TYPE(ewald_environment_type), POINTER              :: ewald_env
     527              :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     528              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     529           56 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     530           56 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     531              :       TYPE(xtb_atom_type), POINTER                       :: xtb_atom_a, xtb_atom_b
     532              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
     533              : 
     534           56 :       CALL timeset(routineN, handle)
     535              : 
     536           56 :       iunit = cp_logger_get_default_unit_nr()
     537              : 
     538              :       CALL get_qs_env(qs_env, &
     539              :                       qs_kind_set=qs_kind_set, &
     540              :                       atomic_kind_set=atomic_kind_set, &
     541              :                       particle_set=particle_set, &
     542              :                       cell=cell, &
     543           56 :                       dft_control=dft_control)
     544           56 :       CALL get_qs_env(qs_env, nkind=nkind, natom=natom)
     545              : 
     546           56 :       xtb_control => dft_control%qs_control%xtb_control
     547              : 
     548           56 :       totalcharge = dft_control%charge
     549              : 
     550              :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
     551           56 :                                atom_of_kind=atom_of_kind, kind_of=kind_of)
     552              : 
     553              :       ! gamma[a,b]
     554          336 :       ALLOCATE (gab(nkind, nkind), gam(nkind))
     555          728 :       gab = 0.0_dp
     556          224 :       DO ikind = 1, nkind
     557          168 :          CALL get_qs_kind(qs_kind_set(ikind), xtb_parameter=xtb_atom_a)
     558          168 :          CALL get_xtb_atom_param(xtb_atom_a, defined=defined)
     559          168 :          IF (.NOT. defined) CYCLE
     560          168 :          CALL get_xtb_atom_param(xtb_atom_a, alpg=ala, eta=gama)
     561          168 :          gam(ikind) = gama
     562          896 :          DO jkind = 1, nkind
     563          504 :             CALL get_qs_kind(qs_kind_set(jkind), xtb_parameter=xtb_atom_b)
     564          504 :             CALL get_xtb_atom_param(xtb_atom_b, defined=defined)
     565          504 :             IF (.NOT. defined) CYCLE
     566          504 :             CALL get_xtb_atom_param(xtb_atom_b, alpg=alb)
     567              :             !
     568         1176 :             gab(ikind, jkind) = SQRT(1._dp/(ala*ala + alb*alb))
     569              :             !
     570              :          END DO
     571              :       END DO
     572              : 
     573          168 :       ALLOCATE (qlag(natom))
     574              : 
     575           56 :       do_ewald = xtb_control%do_ewald
     576              : 
     577           56 :       CALL get_qs_env(qs_env, para_env=para_env, blacs_env=blacs_env)
     578           56 :       IF (do_ewald) THEN
     579              :          CALL get_qs_env(qs_env=qs_env, &
     580           28 :                          ewald_env=ewald_env, ewald_pw=ewald_pw)
     581              :          CALL eeq_solver(qlag, qlam, elag, &
     582              :                          particle_set, kind_of, cell, -dcharges, gam, gab, &
     583              :                          para_env, blacs_env, dft_control, eeq_sparam, &
     584          140 :                          ewald=do_ewald, ewald_env=ewald_env, ewald_pw=ewald_pw, iounit=iunit)
     585              :       ELSE
     586              :          CALL eeq_solver(qlag, qlam, elag, &
     587              :                          particle_set, kind_of, cell, -dcharges, gam, gab, &
     588          140 :                          para_env, blacs_env, dft_control, eeq_sparam, iounit=iunit)
     589              :       END IF
     590              : 
     591              :       ! return Lagrange multipliers
     592          280 :       qlagrange(1:natom) = qlag(1:natom)
     593              : 
     594           56 :       DEALLOCATE (gab, qlag)
     595              : 
     596           56 :       CALL timestop(handle)
     597              : 
     598          112 :    END SUBROUTINE xtb_eeq_lagrange
     599              : 
     600              : END MODULE xtb_eeq
        

Generated by: LCOV version 2.0-1