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 A collection of methods to treat the QM/MM links
      10             : !> \par History
      11             : !>      12.2004 created [tlaino]
      12             : !> \author Teodoro Laino
      13             : ! **************************************************************************************************
      14             : MODULE qmmm_links_methods
      15             : 
      16             :    USE cp_log_handling,                 ONLY: cp_to_string
      17             :    USE kinds,                           ONLY: dp
      18             :    USE particle_types,                  ONLY: particle_type
      19             :    USE qmmm_types_low,                  ONLY: add_set_type,&
      20             :                                               qmmm_env_qm_type,&
      21             :                                               qmmm_imomm_link_type,&
      22             :                                               qmmm_links_type
      23             : #include "./base/base_uses.f90"
      24             : 
      25             :    IMPLICIT NONE
      26             :    PRIVATE
      27             : 
      28             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      29             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_links_methods'
      30             :    PUBLIC ::  qmmm_link_Imomm_coord, &
      31             :              qmmm_link_Imomm_forces, &
      32             :              qmmm_added_chrg_coord, &
      33             :              qmmm_added_chrg_forces
      34             : 
      35             : CONTAINS
      36             : 
      37             : ! **************************************************************************************************
      38             : !> \brief correct the position for qm/mm IMOMM link type
      39             : !> \param qmmm_links ...
      40             : !> \param particles ...
      41             : !> \param qm_atom_index ...
      42             : !> \par History
      43             : !>      12.2004 created [tlaino]
      44             : !> \author Teodoro Laino
      45             : ! **************************************************************************************************
      46         338 :    SUBROUTINE qmmm_link_Imomm_coord(qmmm_links, particles, qm_atom_index)
      47             :       TYPE(qmmm_links_type), POINTER                     :: qmmm_links
      48             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
      49             :       INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
      50             : 
      51             :       INTEGER                                            :: ilink, ip, ip_mm, ip_qm, mm_index, &
      52             :                                                             n_imomm, qm_index
      53             :       REAL(KIND=dp)                                      :: alpha
      54             :       TYPE(qmmm_imomm_link_type), POINTER                :: my_link
      55             : 
      56         338 :       n_imomm = SIZE(qmmm_links%imomm)
      57         338 :       CPASSERT(n_imomm /= 0)
      58        1270 :       DO ilink = 1, n_imomm
      59         932 :          my_link => qmmm_links%imomm(ilink)%link
      60         932 :          qm_index = my_link%qm_index
      61         932 :          mm_index = my_link%mm_index
      62         932 :          alpha = 1.0_dp/my_link%alpha
      63       25616 :          DO ip = 1, SIZE(qm_atom_index)
      64       25616 :             IF (qm_atom_index(ip) == qm_index) EXIT
      65             :          END DO
      66         932 :          IF (ip == SIZE(qm_atom_index) + 1) &
      67             :             CALL cp_abort(__LOCATION__, &
      68             :                           "QM atom index ("//cp_to_string(qm_index)//") specified in the LINK section nr.("// &
      69           0 :                           cp_to_string(ilink)//") is not defined as a QM atom! Please inspect your QM_KIND sections. ")
      70         932 :          ip_qm = ip
      71       44324 :          DO ip = 1, SIZE(qm_atom_index)
      72       44324 :             IF (qm_atom_index(ip) == mm_index) EXIT
      73             :          END DO
      74         932 :          IF (ip == SIZE(qm_atom_index) + 1) &
      75             :             CALL cp_abort(__LOCATION__, &
      76             :                           "Error in setting up the MM atom index ("//cp_to_string(mm_index)// &
      77           0 :                           ") specified in the LINK section nr.("//cp_to_string(ilink)//"). Please report this bug! ")
      78         932 :          ip_mm = ip
      79        4066 :          particles(ip_mm)%r = alpha*particles(ip_mm)%r + (1.0_dp - alpha)*particles(ip_qm)%r
      80             :       END DO
      81             : 
      82         338 :    END SUBROUTINE qmmm_link_Imomm_coord
      83             : 
      84             : ! **************************************************************************************************
      85             : !> \brief correct the forces for qm/mm IMOMM link type
      86             : !> \param qmmm_links ...
      87             : !> \param particles_qm ...
      88             : !> \param qm_atom_index ...
      89             : !> \par History
      90             : !>      12.2004 created [tlaino]
      91             : !> \author Teodoro Laino
      92             : ! **************************************************************************************************
      93         274 :    SUBROUTINE qmmm_link_Imomm_forces(qmmm_links, particles_qm, qm_atom_index)
      94             :       TYPE(qmmm_links_type), POINTER                     :: qmmm_links
      95             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles_qm
      96             :       INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index
      97             : 
      98             :       INTEGER                                            :: ilink, ip, ip_mm, ip_qm, mm_index, &
      99             :                                                             n_imomm, qm_index
     100             :       REAL(KIND=dp)                                      :: alpha
     101             :       TYPE(qmmm_imomm_link_type), POINTER                :: my_link
     102             : 
     103         274 :       n_imomm = SIZE(qmmm_links%imomm)
     104         274 :       CPASSERT(n_imomm /= 0)
     105        1008 :       DO ilink = 1, n_imomm
     106         734 :          my_link => qmmm_links%imomm(ilink)%link
     107         734 :          qm_index = my_link%qm_index
     108         734 :          mm_index = my_link%mm_index
     109         734 :          alpha = 1.0_dp/my_link%alpha
     110       19558 :          DO ip = 1, SIZE(qm_atom_index)
     111       19558 :             IF (qm_atom_index(ip) == qm_index) EXIT
     112             :          END DO
     113         734 :          IF (ip == SIZE(qm_atom_index) + 1) &
     114             :             CALL cp_abort(__LOCATION__, &
     115             :                           "QM atom index ("//cp_to_string(qm_index)//") specified in the LINK section nr.("// &
     116           0 :                           cp_to_string(ilink)//") is not defined as a QM atom! Please inspect your QM_KIND sections. ")
     117         734 :          ip_qm = ip
     118       34278 :          DO ip = 1, SIZE(qm_atom_index)
     119       34278 :             IF (qm_atom_index(ip) == mm_index) EXIT
     120             :          END DO
     121         734 :          IF (ip == SIZE(qm_atom_index) + 1) &
     122             :             CALL cp_abort(__LOCATION__, &
     123             :                           "Error in setting up the MM atom index ("//cp_to_string(mm_index)// &
     124           0 :                           ") specified in the LINK section nr.("//cp_to_string(ilink)//"). Please report this bug! ")
     125         734 :          ip_mm = ip
     126        2936 :          particles_qm(ip_qm)%f = particles_qm(ip_qm)%f + particles_qm(ip_mm)%f*(1.0_dp - alpha)
     127        3210 :          particles_qm(ip_mm)%f = particles_qm(ip_mm)%f*alpha
     128             :       END DO
     129             : 
     130         274 :    END SUBROUTINE qmmm_link_Imomm_forces
     131             : 
     132             : ! **************************************************************************************************
     133             : !> \brief correct the position for added charges in qm/mm link scheme
     134             : !> \param qmmm_env ...
     135             : !> \param particles ...
     136             : !> \par History
     137             : !>      01.2005 created [tlaino]
     138             : !> \author Teodoro Laino
     139             : ! **************************************************************************************************
     140          32 :    SUBROUTINE qmmm_added_chrg_coord(qmmm_env, particles)
     141             :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     142             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
     143             : 
     144             :       INTEGER                                            :: I, Index1, Index2
     145             :       REAL(KIND=dp)                                      :: alpha
     146             :       TYPE(add_set_type), POINTER                        :: added_charges
     147             : 
     148          32 :       added_charges => qmmm_env%added_charges
     149             : 
     150         144 :       DO i = 1, added_charges%num_mm_atoms
     151         112 :          Index1 = added_charges%add_env(i)%Index1
     152         112 :          Index2 = added_charges%add_env(i)%Index2
     153         112 :          alpha = added_charges%add_env(i)%alpha
     154         928 :          added_charges%added_particles(i)%r = alpha*particles(Index1)%r + (1.0_dp - alpha)*particles(Index2)%r
     155             :       END DO
     156             : 
     157          32 :    END SUBROUTINE qmmm_added_chrg_coord
     158             : 
     159             : ! **************************************************************************************************
     160             : !> \brief correct the forces due to the  added charges in qm/mm link scheme
     161             : !> \param qmmm_env ...
     162             : !> \param particles ...
     163             : !> \par History
     164             : !>      01.2005 created [tlaino]
     165             : !> \author Teodoro Laino
     166             : ! **************************************************************************************************
     167          32 :    SUBROUTINE qmmm_added_chrg_forces(qmmm_env, particles)
     168             :       TYPE(qmmm_env_qm_type), POINTER                    :: qmmm_env
     169             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particles
     170             : 
     171             :       INTEGER                                            :: I, Index1, Index2
     172             :       REAL(KIND=dp)                                      :: alpha
     173             :       TYPE(add_set_type), POINTER                        :: added_charges
     174             : 
     175          32 :       added_charges => qmmm_env%added_charges
     176             : 
     177         144 :       DO i = 1, added_charges%num_mm_atoms
     178         112 :          Index1 = added_charges%add_env(i)%Index1
     179         112 :          Index2 = added_charges%add_env(i)%Index2
     180         112 :          alpha = added_charges%add_env(i)%alpha
     181         896 :          particles(Index1)%f = particles(Index1)%f + alpha*added_charges%added_particles(i)%f
     182         928 :          particles(Index2)%f = particles(Index2)%f + (1.0_dp - alpha)*added_charges%added_particles(i)%f
     183             :       END DO
     184             : 
     185          32 :    END SUBROUTINE qmmm_added_chrg_forces
     186             : 
     187             : END MODULE qmmm_links_methods