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