1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7MODULE topology_gromos 8 USE cp_log_handling, ONLY: cp_get_default_logger,& 9 cp_logger_type 10 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 11 cp_print_key_unit_nr 12 USE cp_para_types, ONLY: cp_para_env_type 13 USE cp_parser_methods, ONLY: parser_get_next_line,& 14 parser_get_object,& 15 parser_search_string 16 USE cp_parser_types, ONLY: cp_parser_type,& 17 parser_create,& 18 parser_release 19 USE cp_units, ONLY: cp_unit_to_cp2k 20 USE input_cp2k_restarts_util, ONLY: section_velocity_val_set 21 USE input_section_types, ONLY: section_get_rval,& 22 section_vals_get_subs_vals,& 23 section_vals_type 24 USE kinds, ONLY: default_string_length,& 25 dp 26 USE memory_utilities, ONLY: reallocate 27 USE string_table, ONLY: s2s,& 28 str2id 29 USE topology_types, ONLY: atom_info_type,& 30 connectivity_info_type,& 31 topology_parameters_type 32#include "./base/base_uses.f90" 33 34 IMPLICIT NONE 35 36 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_gromos' 37 38 PRIVATE 39 PUBLIC :: read_topology_gromos, read_coordinate_g96 40 41CONTAINS 42 43! ************************************************************************************************** 44!> \brief Read GROMOS topology file 45!> \param file_name ... 46!> \param topology ... 47!> \param para_env ... 48!> \param subsys_section ... 49! ************************************************************************************************** 50 SUBROUTINE read_topology_gromos(file_name, topology, para_env, subsys_section) 51 CHARACTER(LEN=*), INTENT(IN) :: file_name 52 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 53 TYPE(cp_para_env_type), POINTER :: para_env 54 TYPE(section_vals_type), POINTER :: subsys_section 55 56 CHARACTER(len=*), PARAMETER :: routineN = 'read_topology_gromos', & 57 routineP = moduleN//':'//routineN 58 59 CHARACTER :: ctemp 60 CHARACTER(LEN=default_string_length) :: label, string 61 CHARACTER(LEN=default_string_length), & 62 DIMENSION(21) :: avail_section 63 CHARACTER(LEN=default_string_length), POINTER :: namearray1(:), namearray2(:) 64 INTEGER :: begin, handle, i, iatom, ibond, icon, ii(50), index_now, iresid, isolvent, itemp, & 65 itype, iw, jatom, natom, natom_prev, nbond, nbond_prev, ncon, nsolute, nsolvent, ntype, & 66 offset, offset2, stat 67 INTEGER, POINTER :: ba(:), bb(:), na(:) 68 LOGICAL :: found 69 REAL(KIND=dp) :: ftemp 70 REAL(KIND=dp), POINTER :: ac(:), am(:) 71 TYPE(atom_info_type), POINTER :: atom_info 72 TYPE(connectivity_info_type), POINTER :: conn_info 73 TYPE(cp_logger_type), POINTER :: logger 74 TYPE(cp_parser_type), POINTER :: parser 75 76 NULLIFY (parser, logger) 77 NULLIFY (namearray1, namearray2) 78 logger => cp_get_default_logger() 79 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/GTOP_INFO", & 80 extension=".subsysLog") 81 CALL timeset(routineN, handle) 82 83 avail_section(1) = "TITLE" 84 avail_section(2) = "TOPPHYSCON" 85 avail_section(3) = "TOPVERSION" 86 avail_section(4) = "ATOMTYPENAME" 87 avail_section(5) = "RESNAME" 88 avail_section(6) = "SOLUTEATOM" 89 avail_section(7) = "BONDTYPE" 90 avail_section(8) = "BONDH" 91 avail_section(9) = "BOND" 92 avail_section(10) = "BONDANGLETYPE" 93 avail_section(11) = "BONDANGLEH" 94 avail_section(12) = "BONDANGLE" 95 avail_section(13) = "IMPDIHEDRALTYPE" 96 avail_section(14) = "IMPDIHEDRALH" 97 avail_section(15) = "IMPDIHEDRAL" 98 avail_section(16) = "DIHEDRALTYPE" 99 avail_section(17) = "DIHEDRALH" 100 avail_section(18) = "DIHEDRAL" 101 avail_section(19) = "LJPARAMETERS" 102 avail_section(20) = "SOLVENTATOM" 103 avail_section(21) = "SOLVENTCONSTR" 104 105 atom_info => topology%atom_info 106 conn_info => topology%conn_info 107 108 natom_prev = 0 109 IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname) 110 ! TITLE SECTION 111 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TITLE section' 112 CALL parser_create(parser, file_name, para_env=para_env) 113 label = TRIM(avail_section(1)) 114 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 115 IF (found) THEN 116 DO 117 CALL parser_get_next_line(parser, 1) 118 CALL parser_get_object(parser, string, string_length=80) 119 IF (string == TRIM("END")) EXIT 120 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(string) 121 END DO 122 END IF 123 CALL parser_release(parser) 124 ! TOPPHYSCON SECTION 125 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TOPPHYSCON section' 126 CALL parser_create(parser, file_name, para_env=para_env) 127 label = TRIM(avail_section(2)) 128 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 129 IF (found) THEN 130 DO 131 CALL parser_get_next_line(parser, 1) 132 CALL parser_get_object(parser, string) 133 IF (string == TRIM("END")) EXIT 134 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(string) 135 END DO 136 END IF 137 CALL parser_release(parser) 138 ! TOPVERSION SECTION 139 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the TOPVERSION section' 140 CALL parser_create(parser, file_name, para_env=para_env) 141 label = TRIM(avail_section(3)) 142 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 143 IF (found) THEN 144 DO 145 CALL parser_get_next_line(parser, 1) 146 CALL parser_get_object(parser, string) 147 IF (string == TRIM("END")) EXIT 148 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(string) 149 END DO 150 END IF 151 CALL parser_release(parser) 152 ! ATOMTYPENAME SECTION 153 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the ATOMTYPENAME section' 154 CALL parser_create(parser, file_name, para_env=para_env) 155 label = TRIM(avail_section(4)) 156 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 157 IF (found) THEN 158 CALL parser_get_next_line(parser, 1) 159 CALL parser_get_object(parser, ntype) 160 CALL reallocate(namearray1, 1, ntype) 161 DO itype = 1, ntype 162 CALL parser_get_next_line(parser, 1) 163 CALL parser_get_object(parser, namearray1(itype)) 164 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(namearray1(itype)) 165 END DO 166 END IF 167 CALL parser_release(parser) 168 ! RESNAME SECTION 169 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the RESNAME section' 170 CALL parser_create(parser, file_name, para_env=para_env) 171 label = TRIM(avail_section(5)) 172 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 173 IF (found) THEN 174 CALL parser_get_next_line(parser, 1) 175 CALL parser_get_object(parser, ntype) 176 CALL reallocate(namearray2, 1, ntype) 177 DO itype = 1, ntype 178 CALL parser_get_next_line(parser, 1) 179 CALL parser_get_object(parser, namearray2(itype)) 180 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(namearray2(itype)) 181 END DO 182 END IF 183 CALL parser_release(parser) 184 ! SOLUTEATOM SECTION 185 iresid = 1 186 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the SOLUTEATOM section' 187 CALL parser_create(parser, file_name, para_env=para_env) 188 label = TRIM(avail_section(6)) 189 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 190 IF (found) THEN 191 CALL parser_get_next_line(parser, 1) 192 CALL parser_get_object(parser, natom) 193 CALL reallocate(atom_info%id_molname, 1, natom_prev + natom) 194 CALL reallocate(atom_info%resid, 1, natom_prev + natom) 195 CALL reallocate(atom_info%id_resname, 1, natom_prev + natom) 196 CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom) 197 CALL reallocate(atom_info%id_element, 1, natom_prev + natom) 198 CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom) 199 CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom) 200 CALL parser_get_next_line(parser, 1) 201 DO iatom = 1, natom 202 index_now = iatom + natom_prev 203 CALL parser_get_object(parser, itemp) 204 CALL parser_get_object(parser, itemp) 205 atom_info%resid(index_now) = itemp 206 atom_info%id_molname(index_now) = str2id(s2s(namearray2(itemp))) 207 atom_info%id_resname(index_now) = str2id(s2s(namearray2(itemp))) 208 CALL parser_get_object(parser, string) 209 CALL parser_get_object(parser, itemp) 210 atom_info%id_atmname(index_now) = str2id(s2s(namearray1(itemp))) 211 atom_info%id_element(index_now) = str2id(s2s(namearray1(itemp))) 212 CALL parser_get_object(parser, atom_info%atm_mass(index_now)) 213 CALL parser_get_object(parser, atom_info%atm_charge(index_now)) 214 CALL parser_get_object(parser, itemp) 215 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT SOLUTEATOM INFO HERE!!!!" 216 CALL parser_get_object(parser, ntype) 217 DO i = 1, 50 218 ii(i) = -1 219 END DO 220 stat = -1 221 begin = 1 222 IF (ntype /= 0) THEN 223 DO 224 IF (begin .EQ. 1) THEN 225 READ (parser%input_line, iostat=stat, FMT=*) itemp, itemp, ctemp, itemp, ftemp, ftemp, & 226 itemp, itemp, (ii(i), i=begin, ntype) 227 ELSE IF (begin .GT. 1) THEN 228 CALL parser_get_next_line(parser, 1) 229 READ (parser%input_line, iostat=stat, FMT=*) (ii(i), i=begin, ntype) 230 END IF 231 DO i = ntype, 1, -1 232 IF (ii(i) .LT. 0) THEN 233 begin = i 234 END IF 235 END DO 236 IF (stat .EQ. 0) EXIT 237 END DO 238 END IF 239 ! 1-4 list 240 CALL parser_get_next_line(parser, 1) 241 CALL parser_get_object(parser, ntype) 242 IF (ntype /= 0) THEN 243 itemp = (itemp - 1)/6 + 1 244 offset = 0 245 IF (ASSOCIATED(conn_info%onfo_a)) offset = SIZE(conn_info%onfo_a) 246 CALL reallocate(conn_info%onfo_a, 1, offset + ntype) 247 CALL reallocate(conn_info%onfo_b, 1, offset + ntype) 248 conn_info%onfo_a(offset + 1:offset + ntype) = index_now 249 DO i = 1, 50 250 ii(i) = -1 251 END DO 252 stat = -1 253 begin = 1 254 IF (ntype /= 0) THEN 255 DO 256 IF (begin .EQ. 1) THEN 257 READ (parser%input_line, iostat=stat, FMT=*) itemp, (ii(i), i=begin, ntype) 258 ELSE IF (begin .GT. 1) THEN 259 CALL parser_get_next_line(parser, 1) 260 READ (parser%input_line, iostat=stat, FMT=*) (ii(i), i=begin, ntype) 261 END IF 262 DO i = ntype, 1, -1 263 IF (ii(i) .LT. 0) THEN 264 begin = i 265 END IF 266 END DO 267 IF (stat .EQ. 0) EXIT 268 END DO 269 DO i = 1, ntype 270 conn_info%onfo_b(offset + i) = ii(i) 271 END DO 272 END IF 273 CALL parser_get_next_line(parser, 1) 274 ELSE 275 CALL parser_get_next_line(parser, 1) 276 END IF 277 END DO 278 END IF 279 nsolute = natom 280 CALL parser_release(parser) 281 282 CALL parser_create(parser, file_name, para_env=para_env) 283 ! BONDH SECTION 284 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDH section' 285 label = TRIM(avail_section(8)) 286 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 287 IF (found) THEN 288 CALL parser_get_next_line(parser, 1) 289 CALL parser_get_object(parser, ntype) 290 offset = 0 291 IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a) 292 CALL reallocate(conn_info%bond_a, 1, offset + ntype) 293 CALL reallocate(conn_info%bond_b, 1, offset + ntype) 294 CALL reallocate(conn_info%bond_type, 1, offset + ntype) 295 DO itype = 1, ntype 296 CALL parser_get_next_line(parser, 1) 297 CALL parser_get_object(parser, conn_info%bond_a(offset + itype)) 298 CALL parser_get_object(parser, conn_info%bond_b(offset + itype)) 299 CALL parser_get_object(parser, itemp) 300 conn_info%bond_type(offset + itype) = itemp 301 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDH INFO HERE!!!!" 302 END DO 303 conn_info%bond_a(offset + 1:offset + ntype) = conn_info%bond_a(offset + 1:offset + ntype) + natom_prev 304 conn_info%bond_b(offset + 1:offset + ntype) = conn_info%bond_b(offset + 1:offset + ntype) + natom_prev 305 END IF 306 ! BOND SECTION 307 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BOND section' 308 label = TRIM(avail_section(9)) 309 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 310 IF (found) THEN 311 CALL parser_get_next_line(parser, 1) 312 CALL parser_get_object(parser, ntype) 313 offset = 0 314 IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a) 315 CALL reallocate(conn_info%bond_a, 1, offset + ntype) 316 CALL reallocate(conn_info%bond_b, 1, offset + ntype) 317 CALL reallocate(conn_info%bond_type, 1, offset + ntype) 318 DO itype = 1, ntype 319 CALL parser_get_next_line(parser, 1) 320 CALL parser_get_object(parser, conn_info%bond_a(offset + itype)) 321 CALL parser_get_object(parser, conn_info%bond_b(offset + itype)) 322 CALL parser_get_object(parser, itemp) 323 conn_info%bond_type(offset + itype) = itemp 324 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BOND INFO HERE!!!!" 325 END DO 326 conn_info%bond_a(offset + 1:offset + ntype) = conn_info%bond_a(offset + 1:offset + ntype) + natom_prev 327 conn_info%bond_b(offset + 1:offset + ntype) = conn_info%bond_b(offset + 1:offset + ntype) + natom_prev 328 END IF 329 ! BONDANGLEH SECTION 330 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDANGLEH section' 331 label = TRIM(avail_section(11)) 332 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 333 IF (found) THEN 334 CALL parser_get_next_line(parser, 1) 335 CALL parser_get_object(parser, ntype) 336 offset = 0 337 IF (ASSOCIATED(conn_info%theta_a)) offset = SIZE(conn_info%theta_a) 338 CALL reallocate(conn_info%theta_a, 1, offset + ntype) 339 CALL reallocate(conn_info%theta_b, 1, offset + ntype) 340 CALL reallocate(conn_info%theta_c, 1, offset + ntype) 341 CALL reallocate(conn_info%theta_type, 1, offset + ntype) 342 DO itype = 1, ntype 343 CALL parser_get_next_line(parser, 1) 344 CALL parser_get_object(parser, conn_info%theta_a(offset + itype)) 345 CALL parser_get_object(parser, conn_info%theta_b(offset + itype)) 346 CALL parser_get_object(parser, conn_info%theta_c(offset + itype)) 347 CALL parser_get_object(parser, itemp) 348 conn_info%theta_type(offset + itype) = itemp 349 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDANGLEH INFO HERE!!!!" 350 END DO 351 conn_info%theta_a(offset + 1:offset + ntype) = conn_info%theta_a(offset + 1:offset + ntype) + natom_prev 352 conn_info%theta_b(offset + 1:offset + ntype) = conn_info%theta_b(offset + 1:offset + ntype) + natom_prev 353 conn_info%theta_c(offset + 1:offset + ntype) = conn_info%theta_c(offset + 1:offset + ntype) + natom_prev 354 END IF 355 ! BONDANGLE SECTION 356 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the BONDANGLE section' 357 label = TRIM(avail_section(12)) 358 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 359 IF (found) THEN 360 CALL parser_get_next_line(parser, 1) 361 CALL parser_get_object(parser, ntype) 362 offset = 0 363 IF (ASSOCIATED(conn_info%theta_a)) offset = SIZE(conn_info%theta_a) 364 CALL reallocate(conn_info%theta_a, 1, offset + ntype) 365 CALL reallocate(conn_info%theta_b, 1, offset + ntype) 366 CALL reallocate(conn_info%theta_c, 1, offset + ntype) 367 CALL reallocate(conn_info%theta_type, 1, offset + ntype) 368 DO itype = 1, ntype 369 CALL parser_get_next_line(parser, 1) 370 CALL parser_get_object(parser, conn_info%theta_a(offset + itype)) 371 CALL parser_get_object(parser, conn_info%theta_b(offset + itype)) 372 CALL parser_get_object(parser, conn_info%theta_c(offset + itype)) 373 CALL parser_get_object(parser, itemp) 374 conn_info%theta_type(offset + itype) = itemp 375 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT BONDANGLE INFO HERE!!!!" 376 END DO 377 conn_info%theta_a(offset + 1:offset + ntype) = conn_info%theta_a(offset + 1:offset + ntype) + natom_prev 378 conn_info%theta_b(offset + 1:offset + ntype) = conn_info%theta_b(offset + 1:offset + ntype) + natom_prev 379 conn_info%theta_c(offset + 1:offset + ntype) = conn_info%theta_c(offset + 1:offset + ntype) + natom_prev 380 END IF 381 ! IMPDIHEDRALH SECTION 382 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the IMPDIHEDRALH section' 383 label = TRIM(avail_section(14)) 384 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 385 IF (found) THEN 386 CALL parser_get_next_line(parser, 1) 387 CALL parser_get_object(parser, ntype) 388 offset = 0 389 IF (ASSOCIATED(conn_info%impr_a)) offset = SIZE(conn_info%impr_a) 390 CALL reallocate(conn_info%impr_a, 1, offset + ntype) 391 CALL reallocate(conn_info%impr_b, 1, offset + ntype) 392 CALL reallocate(conn_info%impr_c, 1, offset + ntype) 393 CALL reallocate(conn_info%impr_d, 1, offset + ntype) 394 CALL reallocate(conn_info%impr_type, 1, offset + ntype) 395 DO itype = 1, ntype 396 CALL parser_get_next_line(parser, 1) 397 CALL parser_get_object(parser, conn_info%impr_a(offset + itype)) 398 CALL parser_get_object(parser, conn_info%impr_b(offset + itype)) 399 CALL parser_get_object(parser, conn_info%impr_c(offset + itype)) 400 CALL parser_get_object(parser, conn_info%impr_d(offset + itype)) 401 CALL parser_get_object(parser, itemp) 402 conn_info%impr_type(offset + itype) = itemp 403 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT IMPDIHEDRALH INFO HERE!!!!" 404 END DO 405 conn_info%impr_a(offset + 1:offset + ntype) = conn_info%impr_a(offset + 1:offset + ntype) + natom_prev 406 conn_info%impr_b(offset + 1:offset + ntype) = conn_info%impr_b(offset + 1:offset + ntype) + natom_prev 407 conn_info%impr_c(offset + 1:offset + ntype) = conn_info%impr_c(offset + 1:offset + ntype) + natom_prev 408 conn_info%impr_d(offset + 1:offset + ntype) = conn_info%impr_d(offset + 1:offset + ntype) + natom_prev 409 END IF 410 ! IMPDIHEDRAL SECTION 411 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the IMPDIHEDRAL section' 412 label = TRIM(avail_section(15)) 413 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 414 IF (found) THEN 415 CALL parser_get_next_line(parser, 1) 416 CALL parser_get_object(parser, ntype) 417 offset = 0 418 IF (ASSOCIATED(conn_info%impr_a)) offset = SIZE(conn_info%impr_a) 419 CALL reallocate(conn_info%impr_a, 1, offset + ntype) 420 CALL reallocate(conn_info%impr_b, 1, offset + ntype) 421 CALL reallocate(conn_info%impr_c, 1, offset + ntype) 422 CALL reallocate(conn_info%impr_d, 1, offset + ntype) 423 CALL reallocate(conn_info%impr_type, 1, offset + ntype) 424 DO itype = 1, ntype 425 CALL parser_get_next_line(parser, 1) 426 CALL parser_get_object(parser, conn_info%impr_a(offset + itype)) 427 CALL parser_get_object(parser, conn_info%impr_b(offset + itype)) 428 CALL parser_get_object(parser, conn_info%impr_c(offset + itype)) 429 CALL parser_get_object(parser, conn_info%impr_d(offset + itype)) 430 CALL parser_get_object(parser, itemp) 431 conn_info%impr_type(offset + itype) = itemp 432 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT IMPDIHEDRAL INFO HERE!!!!" 433 END DO 434 conn_info%impr_a(offset + 1:offset + ntype) = conn_info%impr_a(offset + 1:offset + ntype) + natom_prev 435 conn_info%impr_b(offset + 1:offset + ntype) = conn_info%impr_b(offset + 1:offset + ntype) + natom_prev 436 conn_info%impr_c(offset + 1:offset + ntype) = conn_info%impr_c(offset + 1:offset + ntype) + natom_prev 437 conn_info%impr_d(offset + 1:offset + ntype) = conn_info%impr_d(offset + 1:offset + ntype) + natom_prev 438 END IF 439 ! DIHEDRALH SECTION 440 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the DIHEDRALH section' 441 label = TRIM(avail_section(17)) 442 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 443 IF (found) THEN 444 CALL parser_get_next_line(parser, 1) 445 CALL parser_get_object(parser, ntype) 446 offset = 0 447 IF (ASSOCIATED(conn_info%phi_a)) offset = SIZE(conn_info%phi_a) 448 CALL reallocate(conn_info%phi_a, 1, offset + ntype) 449 CALL reallocate(conn_info%phi_b, 1, offset + ntype) 450 CALL reallocate(conn_info%phi_c, 1, offset + ntype) 451 CALL reallocate(conn_info%phi_d, 1, offset + ntype) 452 CALL reallocate(conn_info%phi_type, 1, offset + ntype) 453 DO itype = 1, ntype 454 CALL parser_get_next_line(parser, 1) 455 CALL parser_get_object(parser, conn_info%phi_a(offset + itype)) 456 CALL parser_get_object(parser, conn_info%phi_b(offset + itype)) 457 CALL parser_get_object(parser, conn_info%phi_c(offset + itype)) 458 CALL parser_get_object(parser, conn_info%phi_d(offset + itype)) 459 CALL parser_get_object(parser, itemp) 460 conn_info%phi_type(offset + itype) = itemp 461 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT DIHEDRALH INFO HERE!!!!" 462 END DO 463 conn_info%phi_a(offset + 1:offset + ntype) = conn_info%phi_a(offset + 1:offset + ntype) + natom_prev 464 conn_info%phi_b(offset + 1:offset + ntype) = conn_info%phi_b(offset + 1:offset + ntype) + natom_prev 465 conn_info%phi_c(offset + 1:offset + ntype) = conn_info%phi_c(offset + 1:offset + ntype) + natom_prev 466 conn_info%phi_d(offset + 1:offset + ntype) = conn_info%phi_d(offset + 1:offset + ntype) + natom_prev 467 END IF 468 ! DIHEDRAL SECTION 469 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the DIHEDRAL section' 470 label = TRIM(avail_section(18)) 471 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 472 IF (found) THEN 473 CALL parser_get_next_line(parser, 1) 474 CALL parser_get_object(parser, ntype) 475 offset = 0 476 IF (ASSOCIATED(conn_info%phi_a)) offset = SIZE(conn_info%phi_a) 477 CALL reallocate(conn_info%phi_a, 1, offset + ntype) 478 CALL reallocate(conn_info%phi_b, 1, offset + ntype) 479 CALL reallocate(conn_info%phi_c, 1, offset + ntype) 480 CALL reallocate(conn_info%phi_d, 1, offset + ntype) 481 CALL reallocate(conn_info%phi_type, 1, offset + ntype) 482 DO itype = 1, ntype 483 CALL parser_get_next_line(parser, 1) 484 CALL parser_get_object(parser, conn_info%phi_a(offset + itype)) 485 CALL parser_get_object(parser, conn_info%phi_b(offset + itype)) 486 CALL parser_get_object(parser, conn_info%phi_c(offset + itype)) 487 CALL parser_get_object(parser, conn_info%phi_d(offset + itype)) 488 CALL parser_get_object(parser, itemp) 489 conn_info%phi_type(offset + itype) = itemp 490 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT DIHEDRAL INFO HERE!!!!" 491 END DO 492 conn_info%phi_a(offset + 1:offset + ntype) = conn_info%phi_a(offset + 1:offset + ntype) + natom_prev 493 conn_info%phi_b(offset + 1:offset + ntype) = conn_info%phi_b(offset + 1:offset + ntype) + natom_prev 494 conn_info%phi_c(offset + 1:offset + ntype) = conn_info%phi_c(offset + 1:offset + ntype) + natom_prev 495 conn_info%phi_d(offset + 1:offset + ntype) = conn_info%phi_d(offset + 1:offset + ntype) + natom_prev 496 END IF 497 CALL parser_release(parser) 498 499 ! SOLVENTATOM and SOLVENTCONSTR SECTION 500 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the SOLVENTATOM section' 501 nsolvent = (SIZE(atom_info%r(1, :)) - nsolute)/3 502 503 NULLIFY (na, am, ac, ba, bb) 504 CALL parser_create(parser, file_name, para_env=para_env) 505 label = TRIM(avail_section(20)) 506 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 507 IF (found) THEN 508 CALL parser_get_next_line(parser, 1) 509 CALL parser_get_object(parser, natom) 510 CALL reallocate(na, 1, natom) 511 CALL reallocate(am, 1, natom) 512 CALL reallocate(ac, 1, natom) 513 DO iatom = 1, natom 514 CALL parser_get_next_line(parser, 1) 515 CALL parser_get_object(parser, itemp) 516 CALL parser_get_object(parser, string) 517 CALL parser_get_object(parser, na(iatom)) 518 CALL parser_get_object(parser, am(iatom)) 519 CALL parser_get_object(parser, ac(iatom)) 520 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| PUT SOLVENTATOM INFO HERE!!!!" 521 END DO 522 END IF 523 label = TRIM(avail_section(21)) 524 ncon = 0 525 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 526 IF (found) THEN 527 CALL parser_get_next_line(parser, 1) 528 CALL parser_get_object(parser, ncon) 529 CALL reallocate(ba, 1, ncon) 530 CALL reallocate(bb, 1, ncon) 531 DO icon = 1, ncon 532 CALL parser_get_next_line(parser, 1) 533 CALL parser_get_object(parser, ba(icon)) 534 CALL parser_get_object(parser, bb(icon)) 535 END DO 536 END IF 537 CALL parser_release(parser) 538 539 offset = 0 540 IF (ASSOCIATED(atom_info%id_molname)) offset = SIZE(atom_info%id_molname) 541 CALL reallocate(atom_info%id_molname, 1, offset + nsolvent*natom) 542 CALL reallocate(atom_info%resid, 1, offset + nsolvent*natom) 543 CALL reallocate(atom_info%id_resname, 1, offset + nsolvent*natom) 544 CALL reallocate(atom_info%id_atmname, 1, offset + nsolvent*natom) 545 CALL reallocate(atom_info%id_element, 1, offset + nsolvent*natom) 546 CALL reallocate(atom_info%atm_charge, 1, offset + nsolvent*natom) 547 CALL reallocate(atom_info%atm_mass, 1, offset + nsolvent*natom) 548 DO isolvent = 1, nsolvent 549 offset = nsolute + natom*isolvent - natom 550 DO iatom = 1, natom 551 index_now = offset + iatom 552 atom_info%id_atmname(index_now) = str2id(s2s(namearray1(na(iatom)))) 553 atom_info%id_element(index_now) = str2id(s2s(namearray1(na(iatom)))) 554 atom_info%id_molname(index_now) = str2id(s2s("SOL")) 555 atom_info%id_resname(index_now) = str2id(s2s("SOL")) 556 atom_info%resid(index_now) = isolvent 557 atom_info%atm_mass(index_now) = am(iatom) 558 atom_info%atm_charge(index_now) = ac(iatom) 559 END DO 560 END DO 561 562 offset = 0 563 IF (ASSOCIATED(conn_info%bond_a)) offset = SIZE(conn_info%bond_a) 564 offset2 = MAXVAL(conn_info%bond_type(:)) 565 CALL reallocate(conn_info%bond_a, 1, offset + ncon*nsolvent) 566 CALL reallocate(conn_info%bond_b, 1, offset + ncon*nsolvent) 567 CALL reallocate(conn_info%bond_type, 1, offset + ncon*nsolvent) 568 offset = offset - ncon 569 DO isolvent = 1, nsolvent 570 offset = offset + ncon 571 DO icon = 1, ncon 572 conn_info%bond_a(offset + icon) = nsolute + isolvent*ncon - ncon + ba(icon) 573 conn_info%bond_b(offset + icon) = nsolute + isolvent*ncon - ncon + bb(icon) 574 conn_info%bond_type(offset + icon) = offset2 + isolvent*ncon - ncon + icon 575 END DO 576 END DO 577 ! PARA_RES structure 578 i = 0 579 nbond_prev = 0 580 IF (ASSOCIATED(conn_info%c_bond_a)) i = SIZE(conn_info%c_bond_a) 581 nbond = SIZE(conn_info%bond_a) 582 DO ibond = 1 + nbond_prev, nbond + nbond_prev 583 iatom = conn_info%bond_a(ibond) 584 jatom = conn_info%bond_b(ibond) 585 IF (topology%para_res) THEN 586 IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. & 587 (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. & 588 (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN 589 IF (iw > 0) WRITE (iw, '(T2,A,2I3)') "GTOP_INFO| PARA_RES, bond between molecules atom ", & 590 iatom, jatom 591 i = i + 1 592 CALL reallocate(conn_info%c_bond_a, 1, i) 593 CALL reallocate(conn_info%c_bond_b, 1, i) 594 CALL reallocate(conn_info%c_bond_type, 1, i) 595 conn_info%c_bond_a(i) = iatom 596 conn_info%c_bond_b(i) = jatom 597 conn_info%c_bond_type(i) = conn_info%bond_type(ibond) 598 END IF 599 ELSE 600 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) THEN 601 CPABORT("") 602 END IF 603 END IF 604 END DO 605 606 DEALLOCATE (namearray1) 607 DEALLOCATE (namearray2) 608 DEALLOCATE (na) 609 DEALLOCATE (am) 610 DEALLOCATE (ac) 611 IF (ASSOCIATED(ba)) & 612 DEALLOCATE (ba) 613 IF (ASSOCIATED(bb)) & 614 DEALLOCATE (bb) 615 616 CALL timestop(handle) 617 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 618 "PRINT%TOPOLOGY_INFO/GTOP_INFO") 619 620 END SUBROUTINE read_topology_gromos 621 622! ************************************************************************************************** 623!> \brief ... 624!> \param topology ... 625!> \param para_env ... 626!> \param subsys_section ... 627! ************************************************************************************************** 628 SUBROUTINE read_coordinate_g96(topology, para_env, subsys_section) 629 TYPE(topology_parameters_type) :: topology 630 TYPE(cp_para_env_type), POINTER :: para_env 631 TYPE(section_vals_type), POINTER :: subsys_section 632 633 CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_g96', & 634 routineP = moduleN//':'//routineN 635 INTEGER, PARAMETER :: nblock = 1000 636 637 CHARACTER(LEN=default_string_length) :: label, string, strtmp, strtmp2 638 CHARACTER(LEN=default_string_length), DIMENSION(5) :: avail_section 639 INTEGER :: handle, itemp, iw, natom, newsize 640 LOGICAL :: found 641 REAL(KIND=dp) :: pfactor 642 REAL(KIND=dp), DIMENSION(:, :), POINTER :: velocity 643 TYPE(atom_info_type), POINTER :: atom_info 644 TYPE(cp_logger_type), POINTER :: logger 645 TYPE(cp_parser_type), POINTER :: parser 646 TYPE(section_vals_type), POINTER :: velocity_section 647 648 NULLIFY (parser, logger, velocity) 649 logger => cp_get_default_logger() 650 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/G96_INFO", & 651 extension=".subsysLog") 652 CALL timeset(routineN, handle) 653 654 pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR") 655 atom_info => topology%atom_info 656 IF (iw > 0) WRITE (iw, *) " Reading in G96 file ", TRIM(topology%coord_file_name) 657 avail_section(1) = "TITLE" 658 avail_section(2) = "TIMESTEP" 659 avail_section(3) = "POSITION" 660 avail_section(4) = "VELOCITY" 661 avail_section(5) = "BOX" 662 663 ! Element is assigned on the basis of the atm_name 664 topology%aa_element = .TRUE. 665 666 CALL reallocate(atom_info%id_molname, 1, nblock) 667 CALL reallocate(atom_info%id_resname, 1, nblock) 668 CALL reallocate(atom_info%resid, 1, nblock) 669 CALL reallocate(atom_info%id_atmname, 1, nblock) 670 CALL reallocate(atom_info%id_element, 1, nblock) 671 CALL reallocate(atom_info%r, 1, 3, 1, nblock) 672 CALL reallocate(atom_info%atm_mass, 1, nblock) 673 CALL reallocate(atom_info%atm_charge, 1, nblock) 674 CALL reallocate(atom_info%occup, 1, nblock) 675 CALL reallocate(atom_info%beta, 1, nblock) 676 ! TITLE SECTION 677 IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the TITLE section' 678 CALL parser_create(parser, topology%coord_file_name, para_env=para_env) 679 label = TRIM(avail_section(1)) 680 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 681 IF (found) THEN 682 DO 683 CALL parser_get_next_line(parser, 1) 684 CALL parser_get_object(parser, string, string_length=default_string_length) 685 IF (string == TRIM("END")) EXIT 686 IF (iw > 0) WRITE (iw, *) "G96_INFO| ", TRIM(string) 687 END DO 688 END IF 689 CALL parser_release(parser) 690 ! POSITION SECTION 691 IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the POSITION section' 692 CALL parser_create(parser, topology%coord_file_name, para_env=para_env) 693 label = TRIM(avail_section(3)) 694 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 695 IF (found) THEN 696 natom = 0 697 CALL parser_get_next_line(parser, 1) 698 CALL parser_get_object(parser, string, string_length=default_string_length) 699 DO 700 IF (string == TRIM("END")) EXIT 701 natom = natom + 1 702 IF (natom > SIZE(atom_info%id_molname)) THEN 703 newsize = INT(pfactor*natom) 704 CALL reallocate(atom_info%id_molname, 1, newsize) 705 CALL reallocate(atom_info%id_resname, 1, newsize) 706 CALL reallocate(atom_info%resid, 1, newsize) 707 CALL reallocate(atom_info%id_atmname, 1, newsize) 708 CALL reallocate(atom_info%id_element, 1, newsize) 709 CALL reallocate(atom_info%r, 1, 3, 1, newsize) 710 CALL reallocate(atom_info%atm_mass, 1, newsize) 711 CALL reallocate(atom_info%atm_charge, 1, newsize) 712 CALL reallocate(atom_info%occup, 1, newsize) 713 CALL reallocate(atom_info%beta, 1, newsize) 714 END IF 715 READ (string, *) & 716 atom_info%resid(natom), strtmp, strtmp2, & 717 itemp, atom_info%r(1, natom), atom_info%r(2, natom), atom_info%r(3, natom) 718 atom_info%id_resname(natom) = str2id(s2s(strtmp)) 719 atom_info%id_atmname(natom) = str2id(s2s(strtmp2)) 720 atom_info%id_molname(natom) = atom_info%id_resname(natom) 721 atom_info%id_element(natom) = atom_info%id_atmname(natom) 722 atom_info%r(1, natom) = cp_unit_to_cp2k(atom_info%r(1, natom), "nm") 723 atom_info%r(2, natom) = cp_unit_to_cp2k(atom_info%r(2, natom), "nm") 724 atom_info%r(3, natom) = cp_unit_to_cp2k(atom_info%r(3, natom), "nm") 725 IF (iw > 0) WRITE (iw, *) "G96_INFO| PUT POSITION INFO HERE!!!!" 726 CALL parser_get_next_line(parser, 1) 727 CALL parser_get_object(parser, string, string_length=default_string_length) 728 END DO 729 END IF 730 CALL parser_release(parser) 731 CALL reallocate(velocity, 1, 3, 1, natom) 732 ! VELOCITY SECTION 733 IF (topology%use_g96_velocity) THEN 734 IF (iw > 0) WRITE (iw, '(T2,A)') 'G96_INFO| Parsing the VELOCITY section' 735 CALL parser_create(parser, topology%coord_file_name, para_env=para_env) 736 label = TRIM(avail_section(4)) 737 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 738 IF (found) THEN 739 natom = 0 740 CALL parser_get_next_line(parser, 1) 741 CALL parser_get_object(parser, string, string_length=default_string_length) 742 DO 743 IF (string == TRIM("END")) EXIT 744 natom = natom + 1 745 READ (string, *) & 746 atom_info%resid(natom), strtmp, strtmp2, & 747 itemp, velocity(1, natom), velocity(2, natom), velocity(3, natom) 748 atom_info%id_resname(natom) = str2id(strtmp) 749 atom_info%id_atmname(natom) = str2id(strtmp2) 750 atom_info%id_molname(natom) = atom_info%id_resname(natom) 751 atom_info%id_element(natom) = atom_info%id_atmname(natom) 752 velocity(1, natom) = cp_unit_to_cp2k(velocity(1, natom), "nm*ps^-1") 753 velocity(2, natom) = cp_unit_to_cp2k(velocity(2, natom), "nm*ps^-1") 754 velocity(3, natom) = cp_unit_to_cp2k(velocity(3, natom), "nm*ps^-1") 755 IF (iw > 0) WRITE (iw, *) "G96_INFO| PUT VELOCITY INFO HERE!!!!" 756 CALL parser_get_next_line(parser, 1) 757 CALL parser_get_object(parser, string, string_length=default_string_length) 758 END DO 759 CALL parser_release(parser) 760 velocity_section => section_vals_get_subs_vals(subsys_section, "VELOCITY") 761 CALL section_velocity_val_set(velocity_section, velocity=velocity, & 762 conv_factor=1.0_dp) 763 END IF 764 END IF 765 DEALLOCATE (velocity) 766 767 CALL reallocate(atom_info%id_molname, 1, natom) 768 CALL reallocate(atom_info%id_resname, 1, natom) 769 CALL reallocate(atom_info%resid, 1, natom) 770 CALL reallocate(atom_info%id_atmname, 1, natom) 771 CALL reallocate(atom_info%id_element, 1, natom) 772 CALL reallocate(atom_info%r, 1, 3, 1, natom) 773 CALL reallocate(atom_info%atm_mass, 1, natom) 774 CALL reallocate(atom_info%atm_charge, 1, natom) 775 CALL reallocate(atom_info%occup, 1, natom) 776 CALL reallocate(atom_info%beta, 1, natom) 777 778 IF (.NOT. topology%para_res) atom_info%resid(:) = 1 779 780 topology%natoms = natom 781 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 782 "PRINT%TOPOLOGY_INFO/G96_INFO") 783 CALL timestop(handle) 784 785 END SUBROUTINE read_coordinate_g96 786 787END MODULE topology_gromos 788 789