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