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             : !> \author teo
      10             : ! **************************************************************************************************
      11             : MODULE qmmm_topology_util
      12             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      13             :                                               cp_logger_get_default_io_unit,&
      14             :                                               cp_logger_type
      15             :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      16             :                                               cp_print_key_unit_nr
      17             :    USE input_section_types,             ONLY: section_vals_type
      18             :    USE kinds,                           ONLY: default_string_length
      19             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      20             :    USE molecule_types,                  ONLY: get_molecule,&
      21             :                                               molecule_type
      22             :    USE qmmm_types_low,                  ONLY: qmmm_env_mm_type
      23             :    USE string_table,                    ONLY: id2str,&
      24             :                                               s2s,&
      25             :                                               str2id
      26             :    USE string_utilities,                ONLY: compress
      27             :    USE topology_types,                  ONLY: topology_parameters_type
      28             : #include "./base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             :    PRIVATE
      32             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qmmm_topology_util'
      33             : 
      34             :    PUBLIC :: qmmm_coordinate_control, &
      35             :              qmmm_connectivity_control
      36             : 
      37             : CONTAINS
      38             : 
      39             : ! **************************************************************************************************
      40             : !> \brief Modifies the atom_info%id_atmname
      41             : !> \param topology ...
      42             : !> \param qmmm_env ...
      43             : !> \param subsys_section ...
      44             : !> \par History
      45             : !>      11.2004 created [tlaino]
      46             : !> \author Teodoro Laino
      47             : ! **************************************************************************************************
      48         394 :    SUBROUTINE qmmm_coordinate_control(topology, qmmm_env, subsys_section)
      49             : 
      50             :       TYPE(topology_parameters_type), INTENT(INOUT)      :: topology
      51             :       TYPE(qmmm_env_mm_type), POINTER                    :: qmmm_env
      52             :       TYPE(section_vals_type), POINTER                   :: subsys_section
      53             : 
      54             :       CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_coordinate_control'
      55             : 
      56             :       CHARACTER(LEN=default_string_length)               :: prefix_lnk
      57             :       INTEGER                                            :: handle, iatm, iw
      58             :       LOGICAL                                            :: qmmm_index_in_range
      59             :       TYPE(cp_logger_type), POINTER                      :: logger
      60             : 
      61         394 :       CALL timeset(routineN, handle)
      62         394 :       NULLIFY (logger)
      63         394 :       logger => cp_get_default_logger()
      64             :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
      65         394 :                                 extension=".subsysLog")
      66         394 :       IF (iw > 0) WRITE (iw, *) "  Entering qmmm_coordinate_control"
      67             :       !
      68             :       ! setting ilast and ifirst for QM molecule
      69             :       !
      70         394 :       CPASSERT(SIZE(qmmm_env%qm_atom_index) /= 0)
      71             :       qmmm_index_in_range = (MAXVAL(qmmm_env%qm_atom_index) <= SIZE(topology%atom_info%id_atmname)) &
      72        6178 :                             .AND. (MINVAL(qmmm_env%qm_atom_index) > 0)
      73           0 :       CPASSERT(qmmm_index_in_range)
      74        3286 :       DO iatm = 1, SIZE(qmmm_env%qm_atom_index)
      75             :          topology%atom_info%id_atmname(qmmm_env%qm_atom_index(iatm)) = &
      76        2892 :             str2id(s2s("_QM_"//TRIM(id2str(topology%atom_info%id_atmname(qmmm_env%qm_atom_index(iatm))))))
      77             :          topology%atom_info%id_resname(qmmm_env%qm_atom_index(iatm)) = &
      78        3286 :             str2id(s2s("_QM_"//TRIM(id2str(topology%atom_info%id_resname(qmmm_env%qm_atom_index(iatm))))))
      79             :       END DO
      80             :       !
      81             :       ! Modify type for MM link atoms
      82             :       !
      83         394 :       IF (ASSOCIATED(qmmm_env%mm_link_atoms)) THEN
      84         256 :          DO iatm = 1, SIZE(qmmm_env%mm_link_atoms)
      85         194 :             prefix_lnk = "_LNK000"
      86         194 :             WRITE (prefix_lnk(5:), '(I20)') iatm
      87         194 :             CALL compress(prefix_lnk, .TRUE.)
      88             :             topology%atom_info%id_atmname(qmmm_env%mm_link_atoms(iatm)) = &
      89         194 :                str2id(s2s(TRIM(prefix_lnk)//TRIM(id2str(topology%atom_info%id_atmname(qmmm_env%mm_link_atoms(iatm))))))
      90             :             topology%atom_info%id_resname(qmmm_env%mm_link_atoms(iatm)) = &
      91         256 :                str2id(s2s(TRIM(prefix_lnk)//TRIM(id2str(topology%atom_info%id_resname(qmmm_env%mm_link_atoms(iatm))))))
      92             :          END DO
      93             :       END IF
      94             :       !
      95         394 :       IF (iw > 0) WRITE (iw, *) "  Exiting  qmmm_coordinate_control"
      96             :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
      97         394 :                                         "PRINT%TOPOLOGY_INFO/UTIL_INFO")
      98         394 :       CALL timestop(handle)
      99         394 :    END SUBROUTINE qmmm_coordinate_control
     100             : 
     101             : ! **************************************************************************************************
     102             : !> \brief Set up the connectivity for QM/MM calculations
     103             : !> \param molecule_set ...
     104             : !> \param qmmm_env ...
     105             : !> \param subsys_section ...
     106             : !> \par History
     107             : !>      12.2004 created [tlaino]
     108             : !> \author Teodoro Laino
     109             : ! **************************************************************************************************
     110         394 :    SUBROUTINE qmmm_connectivity_control(molecule_set, &
     111             :                                         qmmm_env, subsys_section)
     112             : 
     113             :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     114             :       TYPE(qmmm_env_mm_type), POINTER                    :: qmmm_env
     115             :       TYPE(section_vals_type), POINTER                   :: subsys_section
     116             : 
     117             :       CHARACTER(len=*), PARAMETER :: routineN = 'qmmm_connectivity_control'
     118             : 
     119             :       INTEGER                                            :: first_atom, handle, i, imolecule, iw, &
     120             :                                                             last_atom, natom, output_unit, &
     121             :                                                             qm_mol_num
     122         394 :       INTEGER, DIMENSION(:), POINTER                     :: qm_atom_index, qm_molecule_index
     123             :       LOGICAL                                            :: detected_link
     124             :       TYPE(cp_logger_type), POINTER                      :: logger
     125             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind
     126             :       TYPE(molecule_type), POINTER                       :: molecule
     127             : 
     128         394 :       NULLIFY (qm_atom_index, qm_molecule_index, molecule, molecule_kind)
     129         394 :       detected_link = .FALSE.
     130         788 :       logger => cp_get_default_logger()
     131         394 :       output_unit = cp_logger_get_default_io_unit(logger)
     132             :       iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
     133         394 :                                 extension=".subsysLog")
     134         394 :       CALL timeset(routineN, handle)
     135         394 :       qm_mol_num = 0
     136         394 :       qm_atom_index => qmmm_env%qm_atom_index
     137       62090 :       DO imolecule = 1, SIZE(molecule_set)
     138       61696 :          IF (iw > 0) WRITE (iw, *) "Entering molecule number ::", imolecule
     139       61696 :          molecule => molecule_set(imolecule)
     140             :          CALL get_molecule(molecule, molecule_kind=molecule_kind, &
     141       61696 :                            first_atom=first_atom, last_atom=last_atom)
     142     2272628 :          IF (ANY(qm_atom_index >= first_atom .AND. qm_atom_index <= last_atom)) &
     143        1390 :             qm_mol_num = qm_mol_num + 1
     144             :       END DO
     145             :       !
     146        1182 :       ALLOCATE (qm_molecule_index(qm_mol_num))
     147         394 :       qm_mol_num = 0
     148       62090 :       DO imolecule = 1, SIZE(molecule_set)
     149       61696 :          IF (iw > 0) WRITE (iw, *) "Entering molecule number ::", imolecule
     150       61696 :          molecule => molecule_set(imolecule)
     151             :          CALL get_molecule(molecule, molecule_kind=molecule_kind, &
     152       61696 :                            first_atom=first_atom, last_atom=last_atom)
     153       61696 :          natom = last_atom - first_atom + 1
     154     2273022 :          IF (ANY(qm_atom_index >= first_atom .AND. qm_atom_index <= last_atom)) THEN
     155         996 :             qm_mol_num = qm_mol_num + 1
     156             :             !
     157             :             ! Check if all atoms of the molecule are QM or if a QM/MM link scheme
     158             :             ! need to be used...
     159             :             !
     160         996 :             detected_link = .FALSE.
     161        4682 :             DO i = first_atom, last_atom
     162       97482 :                IF (.NOT. ANY(qm_atom_index == i)) detected_link = .TRUE.
     163             :             END DO
     164         996 :             IF (detected_link) THEN
     165          78 :                IF (iw > 0) WRITE (iw, fmt='(A)', ADVANCE="NO") " QM/MM link detected..."
     166          78 :                IF (.NOT. qmmm_env%qmmm_link) THEN
     167           0 :                   IF (iw > 0) WRITE (iw, fmt='(A)') " Missing LINK section in input file!"
     168             :                   WRITE (output_unit, '(T2,"QMMM_CONNECTIVITY|",A)') &
     169           0 :                      " ERROR in the QM/MM connectivity. A QM/MM LINK was detected but", &
     170           0 :                      " no LINK section was provided in the Input file!", &
     171           0 :                      " This very probably can be identified as an error in the specified QM", &
     172           0 :                      " indexes or in a missing LINK section. Check your structure!"
     173           0 :                   CPABORT("")
     174             :                END IF
     175             :             END IF
     176         996 :             qm_molecule_index(qm_mol_num) = imolecule
     177             :          END IF
     178             :       END DO
     179         394 :       IF (ASSOCIATED(qmmm_env%qm_molecule_index)) DEALLOCATE (qmmm_env%qm_molecule_index)
     180         394 :       qmmm_env%qm_molecule_index => qm_molecule_index
     181         394 :       IF (iw > 0) WRITE (iw, *) "    QM molecule index ::", qm_molecule_index
     182             :       CALL cp_print_key_finished_output(iw, logger, subsys_section, &
     183         394 :                                         "PRINT%TOPOLOGY_INFO/UTIL_INFO")
     184         394 :       CALL timestop(handle)
     185             : 
     186         788 :    END SUBROUTINE qmmm_connectivity_control
     187             : 
     188             : END MODULE qmmm_topology_util