LCOV - code coverage report
Current view: top level - src - qmmm_se_forces.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.9 % 91 90
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 2 2

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Calculation of the derivative of the QMMM Hamiltonian integral
      10              : !>      matrix <a|\sum_i q_i|b> for semi-empirical methods
      11              : !> \author Teodoro Laino - 04.2007 [tlaino]
      12              : ! **************************************************************************************************
      13              : MODULE qmmm_se_forces
      14              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15              :                                               get_atomic_kind
      16              :    USE cell_types,                      ONLY: cell_type,&
      17              :                                               pbc
      18              :    USE cp_control_types,                ONLY: dft_control_type
      19              :    USE cp_dbcsr_api,                    ONLY: dbcsr_get_block_p,&
      20              :                                               dbcsr_p_type
      21              :    USE input_constants,                 ONLY: &
      22              :         do_method_am1, do_method_mndo, do_method_mndod, do_method_pchg, do_method_pdg, &
      23              :         do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1
      24              :    USE kinds,                           ONLY: dp
      25              :    USE message_passing,                 ONLY: mp_para_env_type
      26              :    USE multipole_types,                 ONLY: do_multipole_none
      27              :    USE particle_types,                  ONLY: particle_type
      28              :    USE qmmm_types_low,                  ONLY: qmmm_env_qm_type,&
      29              :                                               qmmm_pot_p_type,&
      30              :                                               qmmm_pot_type
      31              :    USE qmmm_util,                       ONLY: spherical_cutoff_factor
      32              :    USE qs_environment_types,            ONLY: get_qs_env,&
      33              :                                               qs_environment_type
      34              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      35              :                                               qs_kind_type
      36              :    USE qs_ks_qmmm_types,                ONLY: qs_ks_qmmm_env_type
      37              :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      38              :                                               qs_rho_type
      39              :    USE semi_empirical_int_arrays,       ONLY: se_orbital_pointer
      40              :    USE semi_empirical_integrals,        ONLY: dcorecore,&
      41              :                                               drotnuc
      42              :    USE semi_empirical_types,            ONLY: get_se_param,&
      43              :                                               se_int_control_type,&
      44              :                                               se_taper_type,&
      45              :                                               semi_empirical_create,&
      46              :                                               semi_empirical_release,&
      47              :                                               semi_empirical_type,&
      48              :                                               setup_se_int_control_type
      49              :    USE semi_empirical_utils,            ONLY: get_se_type,&
      50              :                                               se_param_set_default
      51              : #include "./base/base_uses.f90"
      52              : 
      53              :    IMPLICIT NONE
      54              : 
      55              :    PRIVATE
      56              : 
      57              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_se_forces'
      58              :    LOGICAL, PARAMETER, PRIVATE          :: debug_this_module = .FALSE.
      59              :    PUBLIC :: deriv_se_qmmm_matrix
      60              : 
      61              : CONTAINS
      62              : 
      63              : ! **************************************************************************************************
      64              : !> \brief Constructs the derivative w.r.t. 1-el semi-empirical hamiltonian
      65              : !>      QMMM terms
      66              : !> \param qs_env ...
      67              : !> \param qmmm_env ...
      68              : !> \param particles_mm ...
      69              : !> \param mm_cell ...
      70              : !> \param para_env ...
      71              : !> \param calc_force ...
      72              : !> \param Forces ...
      73              : !> \param Forces_added_charges ...
      74              : !> \author Teodoro Laino 04.2007 [created]
      75              : ! **************************************************************************************************
      76          936 :    SUBROUTINE deriv_se_qmmm_matrix(qs_env, qmmm_env, particles_mm, mm_cell, para_env, &
      77              :                                    calc_force, Forces, Forces_added_charges)
      78              : 
      79              :       TYPE(qs_environment_type), POINTER                 :: qs_env
      80              :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
      81              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm
      82              :       TYPE(cell_type), POINTER                           :: mm_cell
      83              :       TYPE(mp_para_env_type), POINTER                    :: para_env
      84              :       LOGICAL, INTENT(in), OPTIONAL                      :: calc_force
      85              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: Forces, Forces_added_charges
      86              : 
      87              :       CHARACTER(len=*), PARAMETER :: routineN = 'deriv_se_qmmm_matrix'
      88              : 
      89              :       INTEGER                                            :: handle, i, iatom, ikind, iqm, ispin, &
      90              :                                                             itype, natom, natorb_a, nkind, &
      91              :                                                             number_qm_atoms
      92          936 :       INTEGER, DIMENSION(:), POINTER                     :: list
      93              :       LOGICAL                                            :: anag, defined, found
      94              :       REAL(KIND=dp)                                      :: delta, enuclear
      95          936 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: Forces_QM, p_block_a
      96          936 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      97          936 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p
      98              :       TYPE(dft_control_type), POINTER                    :: dft_control
      99          936 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles_qm
     100          936 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     101              :       TYPE(qs_ks_qmmm_env_type), POINTER                 :: ks_qmmm_env_loc
     102              :       TYPE(qs_rho_type), POINTER                         :: rho
     103              :       TYPE(se_int_control_type)                          :: se_int_control
     104              :       TYPE(se_taper_type), POINTER                       :: se_taper
     105              :       TYPE(semi_empirical_type), POINTER                 :: se_kind_a, se_kind_mm
     106              : 
     107          936 :       CALL timeset(routineN, handle)
     108          936 :       IF (calc_force) THEN
     109          788 :          NULLIFY (rho, atomic_kind_set, qs_kind_set, se_taper)
     110          788 :          NULLIFY (se_kind_a, se_kind_mm, particles_qm)
     111              :          CALL get_qs_env(qs_env=qs_env, &
     112              :                          rho=rho, &
     113              :                          se_taper=se_taper, &
     114              :                          atomic_kind_set=atomic_kind_set, &
     115              :                          qs_kind_set=qs_kind_set, &
     116              :                          ks_qmmm_env=ks_qmmm_env_loc, &
     117              :                          dft_control=dft_control, &
     118              :                          particle_set=particles_qm, &
     119          788 :                          natom=number_qm_atoms)
     120          788 :          SELECT CASE (dft_control%qs_control%method_id)
     121              :          CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, &
     122              :                do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
     123              :             ! Go on with the calculation..
     124              :          CASE DEFAULT
     125              :             ! Otherwise stop..
     126          788 :             CPABORT("Method not available")
     127              :          END SELECT
     128          788 :          anag = dft_control%qs_control%se_control%analytical_gradients
     129          788 :          delta = dft_control%qs_control%se_control%delta
     130              :          ! Setup SE integral control type
     131              :          CALL setup_se_int_control_type( &
     132              :             se_int_control, shortrange=.FALSE., do_ewald_r3=.FALSE., &
     133              :             do_ewald_gks=.FALSE., integral_screening=dft_control%qs_control%se_control%integral_screening, &
     134          788 :             max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
     135              : 
     136              :          ! Create a fake semi-empirical type to handle the classical atom
     137         2364 :          ALLOCATE (Forces_QM(3, number_qm_atoms))
     138          788 :          CALL semi_empirical_create(se_kind_mm)
     139          788 :          CALL se_param_set_default(se_kind_mm, 0, do_method_pchg)
     140          788 :          itype = get_se_type(se_kind_mm%typ)
     141          788 :          nkind = SIZE(atomic_kind_set)
     142          788 :          enuclear = 0.0_dp
     143        17300 :          Forces_QM = 0.0_dp
     144          788 :          CALL qs_rho_get(rho, rho_ao=matrix_p)
     145              : 
     146         1576 :          DO ispin = 1, dft_control%nspins
     147              :             iqm = 0
     148         3204 :             Kinds: DO ikind = 1, nkind
     149         1628 :                CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=list)
     150         1628 :                CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
     151              :                CALL get_se_param(se_kind_a, &
     152              :                                  defined=defined, &
     153         1628 :                                  natorb=natorb_a)
     154         1628 :                IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     155         9800 :                Atoms: DO i = 1, SIZE(list)
     156         4128 :                   iqm = iqm + 1
     157         4128 :                   iatom = list(i)
     158              :                   ! Give back block
     159         4128 :                   NULLIFY (p_block_a)
     160              :                   CALL dbcsr_get_block_p(matrix=matrix_p(ispin)%matrix, &
     161         4128 :                                          row=iatom, col=iatom, BLOCK=p_block_a, found=found)
     162              : 
     163         5756 :                   IF (ASSOCIATED(p_block_a)) THEN
     164              :                      ! Expand derivative of geometrical factors
     165              :                      CALL deriv_se_qmmm_matrix_low(p_block_a, &
     166              :                                                    se_kind_a, &
     167              :                                                    se_kind_mm, &
     168              :                                                    qmmm_env%Potentials, &
     169              :                                                    particles_mm, &
     170              :                                                    qmmm_env%mm_atom_chrg, &
     171              :                                                    qmmm_env%mm_atom_index, &
     172              :                                                    mm_cell, &
     173              :                                                    iatom, &
     174              :                                                    itype, &
     175              :                                                    Forces, &
     176              :                                                    Forces_QM(:, iqm), &
     177              :                                                    se_taper, &
     178              :                                                    se_int_control, &
     179              :                                                    anag, &
     180              :                                                    delta, &
     181              :                                                    qmmm_env%spherical_cutoff, &
     182         2064 :                                                    particles_qm)
     183              :                      ! Possibly added charges
     184         2064 :                      IF (qmmm_env%move_mm_charges .OR. qmmm_env%add_mm_charges) THEN
     185              :                         CALL deriv_se_qmmm_matrix_low(p_block_a, &
     186              :                                                       se_kind_a, &
     187              :                                                       se_kind_mm, &
     188              :                                                       qmmm_env%added_charges%potentials, &
     189              :                                                       qmmm_env%added_charges%added_particles, &
     190              :                                                       qmmm_env%added_charges%mm_atom_chrg, &
     191              :                                                       qmmm_env%added_charges%mm_atom_index, &
     192              :                                                       mm_cell, &
     193              :                                                       iatom, &
     194              :                                                       itype, &
     195              :                                                       Forces_added_charges, &
     196              :                                                       Forces_QM(:, iqm), &
     197              :                                                       se_taper, &
     198              :                                                       se_int_control, &
     199              :                                                       anag, &
     200              :                                                       delta, &
     201              :                                                       qmmm_env%spherical_cutoff, &
     202            0 :                                                       particles_qm)
     203              :                      END IF
     204              :                   END IF
     205              :                END DO Atoms
     206              :             END DO Kinds
     207              :          END DO
     208          788 :          CPASSERT(iqm == number_qm_atoms)
     209              :          ! Transfer QM gradients to the QM particles..
     210        33812 :          CALL para_env%sum(Forces_QM)
     211          788 :          iqm = 0
     212         2416 :          DO ikind = 1, nkind
     213         1628 :             CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=list)
     214         1628 :             CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
     215              :             CALL get_se_param(se_kind_a, &
     216              :                               defined=defined, &
     217         1628 :                               natorb=natorb_a)
     218         1628 :             IF (.NOT. defined .OR. natorb_a < 1) CYCLE
     219         8172 :             DO i = 1, SIZE(list)
     220         4128 :                iqm = iqm + 1
     221         4128 :                iatom = qmmm_env%qm_atom_index(list(i))
     222        30524 :                particles_mm(iatom)%f(:) = particles_mm(iatom)%f(:) + Forces_QM(:, iqm)
     223              :             END DO
     224              :          END DO
     225              :          ! MM forces will be handled directly from the QMMM module in the same way
     226              :          ! as for GPW/GAPW methods
     227          788 :          DEALLOCATE (Forces_QM)
     228          788 :          CALL semi_empirical_release(se_kind_mm)
     229              : 
     230              :       END IF
     231          936 :       CALL timestop(handle)
     232          936 :    END SUBROUTINE deriv_se_qmmm_matrix
     233              : 
     234              : ! **************************************************************************************************
     235              : !> \brief Low Level : Computes derivatives of the 1-el semi-empirical QMMM
     236              : !>                  hamiltonian block w.r.t. MM and QM coordinates
     237              : !> \param p_block_a ...
     238              : !> \param se_kind_a ...
     239              : !> \param se_kind_mm ...
     240              : !> \param potentials ...
     241              : !> \param particles_mm ...
     242              : !> \param mm_charges ...
     243              : !> \param mm_atom_index ...
     244              : !> \param mm_cell ...
     245              : !> \param IndQM ...
     246              : !> \param itype ...
     247              : !> \param forces ...
     248              : !> \param forces_qm ...
     249              : !> \param se_taper ...
     250              : !> \param se_int_control ...
     251              : !> \param anag ...
     252              : !> \param delta ...
     253              : !> \param qmmm_spherical_cutoff ...
     254              : !> \param particles_qm ...
     255              : !> \author Teodoro Laino 04.2007 [created]
     256              : ! **************************************************************************************************
     257         4128 :    SUBROUTINE deriv_se_qmmm_matrix_low(p_block_a, se_kind_a, se_kind_mm, &
     258              :                                        potentials, particles_mm, mm_charges, mm_atom_index, &
     259         2064 :                                        mm_cell, IndQM, itype, forces, forces_qm, se_taper, &
     260              :                                        se_int_control, anag, delta, qmmm_spherical_cutoff, particles_qm)
     261              : 
     262              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: p_block_a
     263              :       TYPE(semi_empirical_type), POINTER                 :: se_kind_a, se_kind_mm
     264              :       TYPE(qmmm_pot_p_type), DIMENSION(:), POINTER       :: potentials
     265              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles_mm
     266              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: mm_charges
     267              :       INTEGER, DIMENSION(:), POINTER                     :: mm_atom_index
     268              :       TYPE(cell_type), POINTER                           :: mm_cell
     269              :       INTEGER, INTENT(IN)                                :: IndQM, itype
     270              :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: forces
     271              :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: forces_qm
     272              :       TYPE(se_taper_type), POINTER                       :: se_taper
     273              :       TYPE(se_int_control_type), INTENT(IN)              :: se_int_control
     274              :       LOGICAL, INTENT(IN)                                :: anag
     275              :       REAL(KIND=dp), INTENT(IN)                          :: delta, qmmm_spherical_cutoff(2)
     276              :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles_qm
     277              : 
     278              :       CHARACTER(len=*), PARAMETER :: routineN = 'deriv_se_qmmm_matrix_low'
     279              : 
     280              :       INTEGER                                            :: handle, i1, i1L, i2, Imm, Imp, IndMM, &
     281              :                                                             Ipot, j1, j1L
     282              :       REAL(KIND=dp)                                      :: rt1, rt2, rt3, sph_chrg_factor
     283              :       REAL(KIND=dp), DIMENSION(3)                        :: denuc, force_ab, r_pbc, rij
     284              :       REAL(KIND=dp), DIMENSION(3, 45)                    :: de1b
     285              :       TYPE(qmmm_pot_type), POINTER                       :: Pot
     286              : 
     287         2064 :       CALL timeset(routineN, handle)
     288              :       ! Loop Over MM atoms - parallelization over MM atoms...
     289              :       ! Loop over Pot stores atoms with the same charge
     290         5802 :       MainLoopPot: DO Ipot = 1, SIZE(Potentials)
     291         3738 :          Pot => Potentials(Ipot)%Pot
     292              :          ! Loop over atoms belonging to this type
     293      8436363 :          LoopMM: DO Imp = 1, SIZE(Pot%mm_atom_index)
     294      8430561 :             Imm = Pot%mm_atom_index(Imp)
     295      8430561 :             IndMM = mm_atom_index(Imm)
     296     33722244 :             r_pbc = pbc(particles_mm(IndMM)%r - particles_qm(IndQM)%r, mm_cell)
     297      8430561 :             rt1 = r_pbc(1)
     298      8430561 :             rt2 = r_pbc(2)
     299      8430561 :             rt3 = r_pbc(3)
     300     33722244 :             rij = (/rt1, rt2, rt3/)
     301      8430561 :             se_kind_mm%zeff = mm_charges(Imm)
     302              :             ! Computes the screening factor for the spherical cutoff
     303      8430561 :             IF (qmmm_spherical_cutoff(1) > 0.0_dp) THEN
     304       910272 :                CALL spherical_cutoff_factor(qmmm_spherical_cutoff, rij, sph_chrg_factor)
     305       910272 :                se_kind_mm%zeff = se_kind_mm%zeff*sph_chrg_factor
     306              :             END IF
     307      8430561 :             IF (ABS(se_kind_mm%zeff) <= EPSILON(0.0_dp)) CYCLE
     308              :             ! Integrals derivatives involving QM - MM atoms
     309              :             CALL drotnuc(se_kind_a, se_kind_mm, rij, itype=itype, de1b=de1b, &
     310              :                          se_int_control=se_int_control, anag=anag, delta=delta, &
     311      7671782 :                          se_taper=se_taper)
     312              :             CALL dcorecore(se_kind_a, se_kind_mm, rij, itype=itype, denuc=denuc, &
     313              :                            se_int_control=se_int_control, anag=anag, delta=delta, &
     314      7671782 :                            se_taper=se_taper)
     315              :             ! Nucler - Nuclear term
     316     30687128 :             force_ab(1:3) = -denuc(1:3)
     317              :             ! Force contribution from the QMMM Hamiltonian
     318      7671782 :             i2 = 0
     319     25586785 :             DO i1L = 1, se_kind_a%natorb
     320     17915003 :                i1 = se_orbital_pointer(i1L)
     321     38401445 :                DO j1L = 1, i1L - 1
     322     20486442 :                   j1 = se_orbital_pointer(j1L)
     323     20486442 :                   i2 = i2 + 1
     324     99860771 :                   force_ab = force_ab - 2.0_dp*de1b(:, i2)*p_block_a(i1, j1)
     325              :                END DO
     326     17915003 :                j1 = se_orbital_pointer(j1L)
     327     17915003 :                i2 = i2 + 1
     328     79331794 :                force_ab = force_ab - de1b(:, i2)*p_block_a(i1, j1)
     329              :             END DO
     330              :             ! The array of QM forces are really the forces
     331     30687128 :             forces_qm(:) = forces_qm(:) - force_ab
     332              :             ! The one of MM atoms are instead gradients
     333     30690866 :             forces(:, Imm) = forces(:, Imm) - force_ab
     334              :          END DO LoopMM
     335              :       END DO MainLoopPot
     336         2064 :       CALL timestop(handle)
     337         2064 :    END SUBROUTINE deriv_se_qmmm_matrix_low
     338              : 
     339              : END MODULE qmmm_se_forces
        

Generated by: LCOV version 2.0-1