1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Functionality to read in PSF topologies and convert it into local 8!> data structures 9!> \author ikuo 10!> tlaino 10.2006 11! ************************************************************************************************** 12MODULE topology_psf 13 USE cp_log_handling, ONLY: cp_get_default_logger,& 14 cp_logger_get_default_io_unit,& 15 cp_logger_type,& 16 cp_to_string 17 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 18 cp_print_key_generate_filename,& 19 cp_print_key_unit_nr 20 USE cp_para_types, ONLY: cp_para_env_type 21 USE cp_parser_methods, ONLY: parser_get_next_line,& 22 parser_get_object,& 23 parser_search_string,& 24 parser_test_next_token 25 USE cp_parser_types, ONLY: cp_parser_type,& 26 parser_create,& 27 parser_release 28 USE force_fields_input, ONLY: read_chrg_section 29 USE input_constants, ONLY: do_conn_psf,& 30 do_conn_psf_u 31 USE input_section_types, ONLY: section_vals_get,& 32 section_vals_get_subs_vals,& 33 section_vals_type,& 34 section_vals_val_get 35 USE kinds, ONLY: default_string_length,& 36 dp 37 USE memory_utilities, ONLY: reallocate 38 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 39 USE string_table, ONLY: id2str,& 40 s2s,& 41 str2id 42 USE string_utilities, ONLY: uppercase 43 USE topology_types, ONLY: atom_info_type,& 44 connectivity_info_type,& 45 topology_parameters_type 46 USE topology_util, ONLY: array1_list_type,& 47 reorder_structure,& 48 tag_molecule 49 USE util, ONLY: sort 50#include "./base/base_uses.f90" 51 52 IMPLICIT NONE 53 54 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_psf' 55 56 PRIVATE 57 PUBLIC :: read_topology_psf, & 58 write_topology_psf, & 59 psf_post_process, & 60 idm_psf 61 62CONTAINS 63 64! ************************************************************************************************** 65!> \brief Read PSF topology file 66!> Teodoro Laino - Introduced CHARMM31 EXT PSF standard format 67!> \param filename ... 68!> \param topology ... 69!> \param para_env ... 70!> \param subsys_section ... 71!> \param psf_type ... 72!> \par History 73!> 04-2007 Teodoro Laino - Zurich University [tlaino] 74!> This routine should contain only information read from the PSF format 75!> and all post_process should be performef in the psf_post_process 76! ************************************************************************************************** 77 SUBROUTINE read_topology_psf(filename, topology, para_env, subsys_section, psf_type) 78 CHARACTER(LEN=*), INTENT(IN) :: filename 79 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 80 TYPE(cp_para_env_type), POINTER :: para_env 81 TYPE(section_vals_type), POINTER :: subsys_section 82 INTEGER, INTENT(IN) :: psf_type 83 84 CHARACTER(len=*), PARAMETER :: routineN = 'read_topology_psf', & 85 routineP = moduleN//':'//routineN 86 87 CHARACTER(LEN=2*default_string_length) :: psf_format 88 CHARACTER(LEN=3) :: c_int 89 CHARACTER(LEN=default_string_length) :: dummy_field, field, label, strtmp1, & 90 strtmp2, strtmp3 91 INTEGER :: handle, i, iatom, ibond, idum, index_now, iphi, itheta, iw, natom, natom_prev, & 92 nbond, nbond_prev, nphi, nphi_prev, ntheta, ntheta_prev, output_unit 93 LOGICAL :: found 94 TYPE(atom_info_type), POINTER :: atom_info 95 TYPE(connectivity_info_type), POINTER :: conn_info 96 TYPE(cp_logger_type), POINTER :: logger 97 TYPE(cp_parser_type), POINTER :: parser 98 99 NULLIFY (parser, logger) 100 logger => cp_get_default_logger() 101 output_unit = cp_logger_get_default_io_unit(logger) 102 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", & 103 extension=".subsysLog") 104 CALL timeset(routineN, handle) 105 CALL parser_create(parser, filename, para_env=para_env) 106 107 atom_info => topology%atom_info 108 conn_info => topology%conn_info 109 natom_prev = 0 110 IF (ASSOCIATED(atom_info%id_molname)) natom_prev = SIZE(atom_info%id_molname) 111 c_int = 'I8' 112 label = 'PSF' 113 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 114 IF (.NOT. found) THEN 115 IF (output_unit > 0) THEN 116 WRITE (output_unit, '(A)') "ERROR| Missing PSF specification line" 117 END IF 118 CPABORT("") 119 END IF 120 DO WHILE (parser_test_next_token(parser) /= "EOL") 121 CALL parser_get_object(parser, field) 122 SELECT CASE (field(1:3)) 123 CASE ("PSF") 124 IF (psf_type == do_conn_psf) THEN 125 ! X-PLOR PSF format "similar" to the plain CHARMM PSF format 126 psf_format = '(I8,1X,A4,I5,1X,A4,1X,A4,1X,A4,1X,2G14.6,I8)' 127 ENDIF 128 CASE ("EXT") 129 IF (psf_type == do_conn_psf) THEN 130 ! EXTEnded CHARMM31 format 131 psf_format = '(I10,T12,A7,T21,I8,T30,A7,T39,A6,T47,A6,T53,F10.6,T69,F8.3,T88,I1)' 132 c_int = 'I10' 133 ELSE 134 CPABORT("PSF_INFO| "//field(1:3)//" :: not available for UPSF format!") 135 ENDIF 136 CASE DEFAULT 137 CPABORT("PSF_INFO| "//field(1:3)//" :: Unimplemented keyword in CP2K PSF/UPSF format!") 138 END SELECT 139 END DO 140 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NATOM section' 141 ! 142 ! ATOM section 143 ! 144 label = '!NATOM' 145 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 146 IF (.NOT. found) THEN 147 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NATOM section ' 148 natom = 0 149 ELSE 150 CALL parser_get_object(parser, natom) 151 IF (natom_prev + natom > topology%natoms) & 152 CALL cp_abort(__LOCATION__, & 153 "Number of atoms in connectivity control is larger than the "// & 154 "number of atoms in coordinate control. check coordinates and "// & 155 "connectivity. ") 156 IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NATOM = ', natom 157 !malloc the memory that we need 158 CALL reallocate(atom_info%id_molname, 1, natom_prev + natom) 159 CALL reallocate(atom_info%resid, 1, natom_prev + natom) 160 CALL reallocate(atom_info%id_resname, 1, natom_prev + natom) 161 CALL reallocate(atom_info%id_atmname, 1, natom_prev + natom) 162 CALL reallocate(atom_info%atm_charge, 1, natom_prev + natom) 163 CALL reallocate(atom_info%atm_mass, 1, natom_prev + natom) 164 !Read in the atom info 165 IF (psf_type == do_conn_psf_u) THEN 166 DO iatom = 1, natom 167 index_now = iatom + natom_prev 168 CALL parser_get_next_line(parser, 1) 169 READ (parser%input_line, FMT=*, ERR=9) i, & 170 strtmp1, & 171 atom_info%resid(index_now), & 172 strtmp2, & 173 dummy_field, & 174 strtmp3, & 175 atom_info%atm_charge(index_now), & 176 atom_info%atm_mass(index_now) 177 atom_info%id_molname(index_now) = str2id(s2s(strtmp1)) 178 atom_info%id_resname(index_now) = str2id(s2s(strtmp2)) 179 atom_info%id_atmname(index_now) = str2id(s2s(strtmp3)) 180 END DO 181 ELSE 182 DO iatom = 1, natom 183 index_now = iatom + natom_prev 184 CALL parser_get_next_line(parser, 1) 185 READ (parser%input_line, FMT=psf_format) & 186 i, & 187 strtmp1, & 188 atom_info%resid(index_now), & 189 strtmp2, & 190 dummy_field, & 191 strtmp3, & 192 atom_info%atm_charge(index_now), & 193 atom_info%atm_mass(index_now), & 194 idum 195 atom_info%id_molname(index_now) = str2id(s2s(strtmp1)) 196 atom_info%id_resname(index_now) = str2id(s2s(strtmp2)) 197 atom_info%id_atmname(index_now) = str2id(s2s(ADJUSTL(strtmp3))) 198 END DO 199 END IF 200 END IF 201 202 ! 203 ! BOND section 204 ! 205 nbond_prev = 0 206 IF (ASSOCIATED(conn_info%bond_a)) nbond_prev = SIZE(conn_info%bond_a) 207 208 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NBOND section' 209 IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated BOND: ', nbond_prev 210 label = '!NBOND' 211 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 212 IF (.NOT. found) THEN 213 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NBOND section ' 214 nbond = 0 215 ELSE 216 CALL parser_get_object(parser, nbond) 217 IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NBOND = ', nbond 218 !malloc the memory that we need 219 CALL reallocate(conn_info%bond_a, 1, nbond_prev + nbond) 220 CALL reallocate(conn_info%bond_b, 1, nbond_prev + nbond) 221 !Read in the bond info 222 IF (psf_type == do_conn_psf_u) THEN 223 DO ibond = 1, nbond, 4 224 CALL parser_get_next_line(parser, 1) 225 index_now = nbond_prev + ibond - 1 226 READ (parser%input_line, FMT=*, ERR=9) (conn_info%bond_a(index_now + i), & 227 conn_info%bond_b(index_now + i), & 228 i=1, MIN(4, (nbond - ibond + 1))) 229 END DO 230 ELSE 231 DO ibond = 1, nbond, 4 232 CALL parser_get_next_line(parser, 1) 233 index_now = nbond_prev + ibond - 1 234 READ (parser%input_line, FMT='(8'//TRIM(c_int)//')') & 235 (conn_info%bond_a(index_now + i), & 236 conn_info%bond_b(index_now + i), & 237 i=1, MIN(4, (nbond - ibond + 1))) 238 END DO 239 END IF 240 IF (ANY(conn_info%bond_a(nbond_prev + 1:) <= 0) .OR. & 241 ANY(conn_info%bond_a(nbond_prev + 1:) > natom) .OR. & 242 ANY(conn_info%bond_b(nbond_prev + 1:) <= 0) .OR. & 243 ANY(conn_info%bond_b(nbond_prev + 1:) > natom)) THEN 244 CPABORT("topology_read, invalid bond") 245 END IF 246 conn_info%bond_a(nbond_prev + 1:) = conn_info%bond_a(nbond_prev + 1:) + natom_prev 247 conn_info%bond_b(nbond_prev + 1:) = conn_info%bond_b(nbond_prev + 1:) + natom_prev 248 END IF 249 ! 250 ! THETA section 251 ! 252 ntheta_prev = 0 253 IF (ASSOCIATED(conn_info%theta_a)) ntheta_prev = SIZE(conn_info%theta_a) 254 255 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NTHETA section' 256 IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated THETA: ', ntheta_prev 257 label = '!NTHETA' 258 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 259 IF (.NOT. found) THEN 260 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NTHETA section ' 261 ntheta = 0 262 ELSE 263 CALL parser_get_object(parser, ntheta) 264 IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NTHETA = ', ntheta 265 !malloc the memory that we need 266 CALL reallocate(conn_info%theta_a, 1, ntheta_prev + ntheta) 267 CALL reallocate(conn_info%theta_b, 1, ntheta_prev + ntheta) 268 CALL reallocate(conn_info%theta_c, 1, ntheta_prev + ntheta) 269 !Read in the bend info 270 IF (psf_type == do_conn_psf_u) THEN 271 DO itheta = 1, ntheta, 3 272 CALL parser_get_next_line(parser, 1) 273 index_now = ntheta_prev + itheta - 1 274 READ (parser%input_line, FMT=*, ERR=9) (conn_info%theta_a(index_now + i), & 275 conn_info%theta_b(index_now + i), & 276 conn_info%theta_c(index_now + i), & 277 i=1, MIN(3, (ntheta - itheta + 1))) 278 END DO 279 ELSE 280 DO itheta = 1, ntheta, 3 281 CALL parser_get_next_line(parser, 1) 282 index_now = ntheta_prev + itheta - 1 283 READ (parser%input_line, FMT='(9'//TRIM(c_int)//')') & 284 (conn_info%theta_a(index_now + i), & 285 conn_info%theta_b(index_now + i), & 286 conn_info%theta_c(index_now + i), & 287 i=1, MIN(3, (ntheta - itheta + 1))) 288 END DO 289 END IF 290 conn_info%theta_a(ntheta_prev + 1:) = conn_info%theta_a(ntheta_prev + 1:) + natom_prev 291 conn_info%theta_b(ntheta_prev + 1:) = conn_info%theta_b(ntheta_prev + 1:) + natom_prev 292 conn_info%theta_c(ntheta_prev + 1:) = conn_info%theta_c(ntheta_prev + 1:) + natom_prev 293 END IF 294 ! 295 ! PHI section 296 ! 297 nphi_prev = 0 298 IF (ASSOCIATED(conn_info%phi_a)) nphi_prev = SIZE(conn_info%phi_a) 299 300 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NPHI section' 301 IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated PHI: ', nphi_prev 302 label = '!NPHI' 303 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 304 IF (.NOT. found) THEN 305 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NPHI section ' 306 nphi = 0 307 ELSE 308 CALL parser_get_object(parser, nphi) 309 IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NPHI = ', nphi 310 !malloc the memory that we need 311 CALL reallocate(conn_info%phi_a, 1, nphi_prev + nphi) 312 CALL reallocate(conn_info%phi_b, 1, nphi_prev + nphi) 313 CALL reallocate(conn_info%phi_c, 1, nphi_prev + nphi) 314 CALL reallocate(conn_info%phi_d, 1, nphi_prev + nphi) 315 !Read in the torsion info 316 IF (psf_type == do_conn_psf_u) THEN 317 DO iphi = 1, nphi, 2 318 CALL parser_get_next_line(parser, 1) 319 index_now = nphi_prev + iphi - 1 320 READ (parser%input_line, FMT=*, ERR=9) (conn_info%phi_a(index_now + i), & 321 conn_info%phi_b(index_now + i), & 322 conn_info%phi_c(index_now + i), & 323 conn_info%phi_d(index_now + i), & 324 i=1, MIN(2, (nphi - iphi + 1))) 325 END DO 326 ELSE 327 DO iphi = 1, nphi, 2 328 CALL parser_get_next_line(parser, 1) 329 index_now = nphi_prev + iphi - 1 330 READ (parser%input_line, FMT='(8'//TRIM(c_int)//')') & 331 (conn_info%phi_a(index_now + i), & 332 conn_info%phi_b(index_now + i), & 333 conn_info%phi_c(index_now + i), & 334 conn_info%phi_d(index_now + i), & 335 i=1, MIN(2, (nphi - iphi + 1))) 336 END DO 337 END IF 338 conn_info%phi_a(nphi_prev + 1:) = conn_info%phi_a(nphi_prev + 1:) + natom_prev 339 conn_info%phi_b(nphi_prev + 1:) = conn_info%phi_b(nphi_prev + 1:) + natom_prev 340 conn_info%phi_c(nphi_prev + 1:) = conn_info%phi_c(nphi_prev + 1:) + natom_prev 341 conn_info%phi_d(nphi_prev + 1:) = conn_info%phi_d(nphi_prev + 1:) + natom_prev 342 END IF 343 ! 344 ! IMPHI section 345 ! 346 nphi_prev = 0 347 IF (ASSOCIATED(conn_info%impr_a)) nphi_prev = SIZE(conn_info%impr_a) 348 349 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| Parsing the NIMPHI section' 350 IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Previous number of allocated IMPHI: ', nphi_prev 351 label = '!NIMPHI' 352 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 353 IF (.NOT. found) THEN 354 IF (iw > 0) WRITE (iw, '(T2,A)') 'PSF_INFO| No NIMPHI section ' 355 nphi = 0 356 ELSE 357 CALL parser_get_object(parser, nphi) 358 IF (iw > 0) WRITE (iw, '(T2,A,'//TRIM(c_int)//')') 'PSF_INFO| NIMPR = ', nphi 359 !malloc the memory that we need 360 CALL reallocate(conn_info%impr_a, 1, nphi_prev + nphi) 361 CALL reallocate(conn_info%impr_b, 1, nphi_prev + nphi) 362 CALL reallocate(conn_info%impr_c, 1, nphi_prev + nphi) 363 CALL reallocate(conn_info%impr_d, 1, nphi_prev + nphi) 364 !Read in the improper torsion info 365 IF (psf_type == do_conn_psf_u) THEN 366 DO iphi = 1, nphi, 2 367 CALL parser_get_next_line(parser, 1) 368 index_now = nphi_prev + iphi - 1 369 READ (parser%input_line, FMT=*, ERR=9) (conn_info%impr_a(index_now + i), & 370 conn_info%impr_b(index_now + i), & 371 conn_info%impr_c(index_now + i), & 372 conn_info%impr_d(index_now + i), & 373 i=1, MIN(2, (nphi - iphi + 1))) 374 END DO 375 ELSE 376 DO iphi = 1, nphi, 2 377 CALL parser_get_next_line(parser, 1) 378 index_now = nphi_prev + iphi - 1 379 READ (parser%input_line, FMT='(8'//TRIM(c_int)//')') & 380 (conn_info%impr_a(index_now + i), & 381 conn_info%impr_b(index_now + i), & 382 conn_info%impr_c(index_now + i), & 383 conn_info%impr_d(index_now + i), & 384 i=1, MIN(2, (nphi - iphi + 1))) 385 END DO 386 END IF 387 conn_info%impr_a(nphi_prev + 1:) = conn_info%impr_a(nphi_prev + 1:) + natom_prev 388 conn_info%impr_b(nphi_prev + 1:) = conn_info%impr_b(nphi_prev + 1:) + natom_prev 389 conn_info%impr_c(nphi_prev + 1:) = conn_info%impr_c(nphi_prev + 1:) + natom_prev 390 conn_info%impr_d(nphi_prev + 1:) = conn_info%impr_d(nphi_prev + 1:) + natom_prev 391 END IF 392 393 CALL parser_release(parser) 394 CALL timestop(handle) 395 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 396 "PRINT%TOPOLOGY_INFO/PSF_INFO") 397 RETURN 3989 CONTINUE 399 ! Print error and exit 400 IF (output_unit > 0) THEN 401 WRITE (output_unit, '(T2,A)') & 402 "PSF_INFO| Error while reading PSF using the unformatted PSF reading option!", & 403 "PSF_INFO| Try using PSF instead of UPSF." 404 END IF 405 406 CPABORT("Error while reading PSF data!") 407 408 END SUBROUTINE read_topology_psf 409 410! ************************************************************************************************** 411!> \brief Post processing of PSF informations 412!> \param topology ... 413!> \param subsys_section ... 414! ************************************************************************************************** 415 SUBROUTINE psf_post_process(topology, subsys_section) 416 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 417 TYPE(section_vals_type), POINTER :: subsys_section 418 419 CHARACTER(len=*), PARAMETER :: routineN = 'psf_post_process', & 420 routineP = moduleN//':'//routineN 421 422 INTEGER :: handle, i, iatom, ibond, ionfo, iw, & 423 jatom, N, natom, nbond, nonfo, nphi, & 424 ntheta 425 TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bend_list, ex_bond_list 426 TYPE(atom_info_type), POINTER :: atom_info 427 TYPE(connectivity_info_type), POINTER :: conn_info 428 TYPE(cp_logger_type), POINTER :: logger 429 430 NULLIFY (logger) 431 logger => cp_get_default_logger() 432 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", & 433 extension=".subsysLog") 434 CALL timeset(routineN, handle) 435 atom_info => topology%atom_info 436 conn_info => topology%conn_info 437 ! 438 ! PARA_RES structure 439 ! 440 natom = 0 441 nbond = 0 442 i = 0 443 IF (ASSOCIATED(atom_info%id_molname)) natom = SIZE(atom_info%id_molname) 444 IF (ASSOCIATED(conn_info%bond_a)) nbond = SIZE(conn_info%bond_a) 445 IF (ASSOCIATED(conn_info%c_bond_a)) i = SIZE(conn_info%c_bond_a) 446 DO ibond = 1, nbond 447 iatom = conn_info%bond_a(ibond) 448 jatom = conn_info%bond_b(ibond) 449 IF (topology%para_res) THEN 450 IF ((atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) .OR. & 451 (atom_info%resid(iatom) /= atom_info%resid(jatom)) .OR. & 452 (atom_info%id_resname(iatom) /= atom_info%id_resname(jatom))) THEN 453 IF (iw > 0) WRITE (iw, '(T2,A,2I6)') "PSF_INFO| PARA_RES, bond between molecules atom ", & 454 iatom, jatom 455 i = i + 1 456 CALL reallocate(conn_info%c_bond_a, 1, i) 457 CALL reallocate(conn_info%c_bond_b, 1, i) 458 conn_info%c_bond_a(i) = iatom 459 conn_info%c_bond_b(i) = jatom 460 END IF 461 ELSE 462 IF (atom_info%id_molname(iatom) /= atom_info%id_molname(jatom)) THEN 463 CPABORT("") 464 END IF 465 END IF 466 END DO 467 ! 468 ! UB structure 469 ! 470 ntheta = 0 471 IF (ASSOCIATED(conn_info%theta_a)) ntheta = SIZE(conn_info%theta_a) 472 CALL reallocate(conn_info%ub_a, 1, ntheta) 473 CALL reallocate(conn_info%ub_b, 1, ntheta) 474 CALL reallocate(conn_info%ub_c, 1, ntheta) 475 conn_info%ub_a(:) = conn_info%theta_a(:) 476 conn_info%ub_b(:) = conn_info%theta_b(:) 477 conn_info%ub_c(:) = conn_info%theta_c(:) 478 ! 479 ! ONFO structure 480 ! 481 nphi = 0 482 nonfo = 0 483 IF (ASSOCIATED(conn_info%phi_a)) nphi = SIZE(conn_info%phi_a) 484 CALL reallocate(conn_info%onfo_a, 1, nphi) 485 CALL reallocate(conn_info%onfo_b, 1, nphi) 486 conn_info%onfo_a(1:) = conn_info%phi_a(1:) 487 conn_info%onfo_b(1:) = conn_info%phi_d(1:) 488 ! Reorder bonds 489 ALLOCATE (ex_bond_list(natom)) 490 DO I = 1, natom 491 ALLOCATE (ex_bond_list(I)%array1(0)) 492 ENDDO 493 N = 0 494 IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a) 495 CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, N) 496 ! Reorder bends 497 ALLOCATE (ex_bend_list(natom)) 498 DO I = 1, natom 499 ALLOCATE (ex_bend_list(I)%array1(0)) 500 ENDDO 501 N = 0 502 IF (ASSOCIATED(conn_info%theta_a)) N = SIZE(conn_info%theta_a) 503 CALL reorder_structure(ex_bend_list, conn_info%theta_a, conn_info%theta_c, N) 504 DO ionfo = 1, nphi 505 ! Check if the torsion is not shared between angles or bonds 506 IF (ANY(ex_bond_list(conn_info%onfo_a(ionfo))%array1 == conn_info%onfo_b(ionfo)) .OR. & 507 ANY(ex_bend_list(conn_info%onfo_a(ionfo))%array1 == conn_info%onfo_b(ionfo))) CYCLE 508 nonfo = nonfo + 1 509 conn_info%onfo_a(nonfo) = conn_info%onfo_a(ionfo) 510 conn_info%onfo_b(nonfo) = conn_info%onfo_b(ionfo) 511 END DO 512 ! deallocate bends 513 DO I = 1, natom 514 DEALLOCATE (ex_bend_list(I)%array1) 515 ENDDO 516 DEALLOCATE (ex_bend_list) 517 ! deallocate bonds 518 DO I = 1, natom 519 DEALLOCATE (ex_bond_list(I)%array1) 520 ENDDO 521 DEALLOCATE (ex_bond_list) 522 ! Get unique onfo 523 ALLOCATE (ex_bond_list(natom)) 524 DO I = 1, natom 525 ALLOCATE (ex_bond_list(I)%array1(0)) 526 ENDDO 527 N = 0 528 IF (ASSOCIATED(conn_info%onfo_a)) N = nonfo 529 CALL reorder_structure(ex_bond_list, conn_info%onfo_a, conn_info%onfo_b, N) 530 nonfo = 0 531 DO I = 1, natom 532 DO ionfo = 1, SIZE(ex_bond_list(I)%array1) 533 IF (COUNT(ex_bond_list(I)%array1 == ex_bond_list(I)%array1(ionfo)) /= 1) THEN 534 ex_bond_list(I)%array1(ionfo) = 0 535 ELSE 536 IF (ex_bond_list(I)%array1(ionfo) <= I) CYCLE 537 nonfo = nonfo + 1 538 conn_info%onfo_a(nonfo) = I 539 conn_info%onfo_b(nonfo) = ex_bond_list(I)%array1(ionfo) 540 END IF 541 END DO 542 END DO 543 DO I = 1, natom 544 DEALLOCATE (ex_bond_list(I)%array1) 545 ENDDO 546 DEALLOCATE (ex_bond_list) 547 CALL reallocate(conn_info%onfo_a, 1, nonfo) 548 CALL reallocate(conn_info%onfo_b, 1, nonfo) 549 550 CALL timestop(handle) 551 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 552 "PRINT%TOPOLOGY_INFO/PSF_INFO") 553 END SUBROUTINE psf_post_process 554 555! ************************************************************************************************** 556!> \brief Input driven modification (IDM) of PSF defined structures 557!> \param topology ... 558!> \param section ... 559!> \param subsys_section ... 560!> \author Teodoro Laino - Zurich University 04.2007 561! ************************************************************************************************** 562 SUBROUTINE idm_psf(topology, section, subsys_section) 563 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 564 TYPE(section_vals_type), POINTER :: section, subsys_section 565 566 CHARACTER(len=*), PARAMETER :: routineN = 'idm_psf', routineP = moduleN//':'//routineN 567 568 INTEGER :: handle, i, iend, iend1, istart, istart1, & 569 item, iw, j, mol_id, n_rep, natom, & 570 nbond, nimpr, noe, nphi, ntheta 571 INTEGER, DIMENSION(:), POINTER :: tag_mols, tmp, wrk 572 LOGICAL :: explicit 573 TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bond_list 574 TYPE(atom_info_type), POINTER :: atom_info 575 TYPE(connectivity_info_type), POINTER :: conn_info 576 TYPE(cp_logger_type), POINTER :: logger 577 TYPE(section_vals_type), POINTER :: subsection 578 579 NULLIFY (logger) 580 logger => cp_get_default_logger() 581 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", & 582 extension=".subsysLog") 583 CALL timeset(routineN, handle) 584 CALL section_vals_get(section, explicit=explicit) 585 IF (explicit) THEN 586 atom_info => topology%atom_info 587 conn_info => topology%conn_info 588 natom = 0 589 IF (ASSOCIATED(atom_info%id_molname)) natom = SIZE(atom_info%id_molname) 590 nbond = 0 591 IF (ASSOCIATED(conn_info%bond_a)) nbond = SIZE(conn_info%bond_a) 592 ntheta = 0 593 IF (ASSOCIATED(conn_info%theta_a)) ntheta = SIZE(conn_info%theta_a) 594 nphi = 0 595 IF (ASSOCIATED(conn_info%phi_a)) nphi = SIZE(conn_info%phi_a) 596 nimpr = 0 597 IF (ASSOCIATED(conn_info%impr_a)) nimpr = SIZE(conn_info%impr_a) 598 ! Any new defined bond 599 subsection => section_vals_get_subs_vals(section, "BONDS") 600 CALL section_vals_get(subsection, explicit=explicit) 601 IF (explicit) THEN 602 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep) 603 CALL reallocate(conn_info%bond_a, 1, n_rep + nbond) 604 CALL reallocate(conn_info%bond_b, 1, n_rep + nbond) 605 DO i = 1, n_rep 606 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp) 607 conn_info%bond_a(nbond + i) = tmp(1) 608 conn_info%bond_b(nbond + i) = tmp(2) 609 END DO 610 ! And now modify the molecule name if two molecules have been bridged 611 ALLOCATE (ex_bond_list(natom)) 612 ALLOCATE (tag_mols(natom)) 613 ALLOCATE (wrk(natom)) 614 DO j = 1, natom 615 ALLOCATE (ex_bond_list(j)%array1(0)) 616 ENDDO 617 CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, nbond + n_rep) 618 ! Loop over atoms to possiblyt change molecule name 619 tag_mols = -1 620 mol_id = 1 621 DO i = 1, natom 622 IF (tag_mols(i) /= -1) CYCLE 623 CALL tag_molecule(tag_mols, ex_bond_list, i, mol_id) 624 mol_id = mol_id + 1 625 END DO 626 mol_id = mol_id - 1 627 IF (iw > 0) WRITE (iw, '(T2,A,I8)') 'PSF_INFO| Number of molecules detected after merging: ', mol_id 628 ! Now simply check about the contiguousness of molecule definition 629 CALL sort(tag_mols, natom, wrk) 630 item = tag_mols(1) 631 istart = 1 632 DO i = 2, natom 633 IF (tag_mols(i) == item) CYCLE 634 iend = i - 1 635 noe = iend - istart + 1 636 istart1 = MINVAL(wrk(istart:iend)) 637 iend1 = MAXVAL(wrk(istart:iend)) 638 CPASSERT(iend1 - istart1 + 1 == noe) 639 atom_info%id_molname(istart1:iend1) = str2id(s2s("MOL"//cp_to_string(item))) 640 item = tag_mols(i) 641 istart = i 642 END DO 643 iend = i - 1 644 noe = iend - istart + 1 645 istart1 = MINVAL(wrk(istart:iend)) 646 iend1 = MAXVAL(wrk(istart:iend)) 647 CPASSERT(iend1 - istart1 + 1 == noe) 648 atom_info%id_molname(istart1:iend1) = str2id(s2s("MOL"//cp_to_string(item))) 649 ! Deallocate bonds 650 DO i = 1, natom 651 DEALLOCATE (ex_bond_list(i)%array1) 652 ENDDO 653 DEALLOCATE (ex_bond_list) 654 DEALLOCATE (tag_mols) 655 DEALLOCATE (wrk) 656 END IF 657 ! Any new defined angle 658 subsection => section_vals_get_subs_vals(section, "ANGLES") 659 CALL section_vals_get(subsection, explicit=explicit) 660 IF (explicit) THEN 661 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep) 662 CALL reallocate(conn_info%theta_a, 1, n_rep + ntheta) 663 CALL reallocate(conn_info%theta_b, 1, n_rep + ntheta) 664 CALL reallocate(conn_info%theta_c, 1, n_rep + ntheta) 665 DO i = 1, n_rep 666 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp) 667 conn_info%theta_a(ntheta + i) = tmp(1) 668 conn_info%theta_b(ntheta + i) = tmp(2) 669 conn_info%theta_c(ntheta + i) = tmp(3) 670 END DO 671 END IF 672 ! Any new defined torsion 673 subsection => section_vals_get_subs_vals(section, "TORSIONS") 674 CALL section_vals_get(subsection, explicit=explicit) 675 IF (explicit) THEN 676 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep) 677 CALL reallocate(conn_info%phi_a, 1, n_rep + nphi) 678 CALL reallocate(conn_info%phi_b, 1, n_rep + nphi) 679 CALL reallocate(conn_info%phi_c, 1, n_rep + nphi) 680 CALL reallocate(conn_info%phi_d, 1, n_rep + nphi) 681 DO i = 1, n_rep 682 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp) 683 conn_info%phi_a(nphi + i) = tmp(1) 684 conn_info%phi_b(nphi + i) = tmp(2) 685 conn_info%phi_c(nphi + i) = tmp(3) 686 conn_info%phi_d(nphi + i) = tmp(4) 687 END DO 688 END IF 689 ! Any new defined improper 690 subsection => section_vals_get_subs_vals(section, "IMPROPERS") 691 CALL section_vals_get(subsection, explicit=explicit) 692 IF (explicit) THEN 693 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", n_rep_val=n_rep) 694 CALL reallocate(conn_info%impr_a, 1, n_rep + nimpr) 695 CALL reallocate(conn_info%impr_b, 1, n_rep + nimpr) 696 CALL reallocate(conn_info%impr_c, 1, n_rep + nimpr) 697 CALL reallocate(conn_info%impr_d, 1, n_rep + nimpr) 698 DO i = 1, n_rep 699 CALL section_vals_val_get(subsection, "_DEFAULT_KEYWORD_", i_rep_val=i, i_vals=tmp) 700 conn_info%impr_a(nimpr + i) = tmp(1) 701 conn_info%impr_b(nimpr + i) = tmp(2) 702 conn_info%impr_c(nimpr + i) = tmp(3) 703 conn_info%impr_d(nimpr + i) = tmp(4) 704 END DO 705 END IF 706 END IF 707 CALL timestop(handle) 708 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 709 "PRINT%TOPOLOGY_INFO/PSF_INFO") 710 711 END SUBROUTINE idm_psf 712 713! ************************************************************************************************** 714!> \brief Teodoro Laino - 01.2006 715!> Write PSF topology file in the CHARMM31 EXT standard format 716!> \param file_unit ... 717!> \param topology ... 718!> \param subsys_section ... 719!> \param force_env_section ... 720! ************************************************************************************************** 721 SUBROUTINE write_topology_psf(file_unit, topology, subsys_section, force_env_section) 722 INTEGER, INTENT(IN) :: file_unit 723 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 724 TYPE(section_vals_type), POINTER :: subsys_section, force_env_section 725 726 CHARACTER(len=*), PARAMETER :: routineN = 'write_topology_psf', & 727 routineP = moduleN//':'//routineN 728 729 CHARACTER(LEN=2*default_string_length) :: psf_format 730 CHARACTER(LEN=default_string_length) :: c_int, my_tag1, my_tag2, my_tag3, record 731 CHARACTER(LEN=default_string_length), & 732 DIMENSION(:), POINTER :: charge_atm 733 INTEGER :: handle, i, iw, j, my_index, nchg 734 LOGICAL :: explicit, ldum 735 REAL(KIND=dp), DIMENSION(:), POINTER :: charge_inp, charges 736 TYPE(atom_info_type), POINTER :: atom_info 737 TYPE(connectivity_info_type), POINTER :: conn_info 738 TYPE(cp_logger_type), POINTER :: logger 739 TYPE(section_vals_type), POINTER :: print_key, tmp_section 740 741 NULLIFY (logger) 742 logger => cp_get_default_logger() 743 print_key => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%DUMP_PSF") 744 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/PSF_INFO", & 745 extension=".subsysLog") 746 CALL timeset(routineN, handle) 747 748 atom_info => topology%atom_info 749 conn_info => topology%conn_info 750 751 ! Check for charges.. (need to dump them in the PSF..) 752 ALLOCATE (charges(topology%natoms)) 753 charges = atom_info%atm_charge 754 ! Collect charges from Input file.. 755 NULLIFY (tmp_section) 756 tmp_section => section_vals_get_subs_vals(force_env_section, "MM%FORCEFIELD%CHARGE") 757 CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=nchg) 758 IF (explicit) THEN 759 ALLOCATE (charge_atm(nchg)) 760 ALLOCATE (charge_inp(nchg)) 761 CALL read_chrg_section(charge_atm, charge_inp, section=tmp_section, start=0) 762 DO j = 1, topology%natoms 763 record = id2str(atom_info%id_atmname(j)) 764 ldum = qmmm_ff_precond_only_qm(record) 765 CALL uppercase(record) 766 DO i = 1, nchg 767 IF (record == charge_atm(i)) THEN 768 charges(j) = charge_inp(i) 769 EXIT 770 END IF 771 END DO 772 END DO 773 DEALLOCATE (charge_atm) 774 DEALLOCATE (charge_inp) 775 END IF 776 ! fixup for topology output 777 DO j = 1, topology%natoms 778 IF (charges(j) .EQ. -HUGE(0.0_dp)) charges(j) = -99.0_dp 779 ENDDO 780 record = cp_print_key_generate_filename(logger, print_key, & 781 extension=".psf", my_local=.FALSE.) 782 ! build the EXT format 783 c_int = "I10" 784 psf_format = '(I10,T12,A,T21,I0,T30,A,T39,A,T47,A,T53,F10.6,T69,F8.3,T88,I1)' 785 IF (iw > 0) WRITE (iw, '(T2,A)') & 786 "PSF_WRITE| Writing out PSF file with CHARMM31 EXTErnal format: ", TRIM(record) 787 788 WRITE (file_unit, FMT='(A)') "PSF EXT" 789 WRITE (file_unit, FMT='(A)') "" 790 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 1, " !NTITLE" 791 WRITE (file_unit, FMT='(A)') " CP2K generated DUMP of connectivity" 792 WRITE (file_unit, FMT='(A)') "" 793 794 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') topology%natoms, " !NATOM" 795 my_index = 1 796 i = 1 797 my_tag1 = id2str(atom_info%id_molname(i)) 798 my_tag2 = id2str(atom_info%id_resname(i)) 799 my_tag3 = id2str(atom_info%id_atmname(i)) 800 ldum = qmmm_ff_precond_only_qm(my_tag1) 801 ldum = qmmm_ff_precond_only_qm(my_tag2) 802 ldum = qmmm_ff_precond_only_qm(my_tag3) 803 WRITE (file_unit, FMT=psf_format) & 804 i, & 805 TRIM(my_tag1), & 806 my_index, & 807 TRIM(my_tag2), & 808 TRIM(my_tag3), & 809 TRIM(my_tag3), & 810 charges(i), & 811 atom_info%atm_mass(i), & 812 0 813 DO i = 2, topology%natoms 814 IF ((atom_info%map_mol_num(i) /= atom_info%map_mol_num(i - 1)) .OR. & 815 (atom_info%map_mol_res(i) /= atom_info%map_mol_res(i - 1))) my_index = my_index + 1 816 my_tag1 = id2str(atom_info%id_molname(i)) 817 my_tag2 = id2str(atom_info%id_resname(i)) 818 my_tag3 = id2str(atom_info%id_atmname(i)) 819 ldum = qmmm_ff_precond_only_qm(my_tag1) 820 ldum = qmmm_ff_precond_only_qm(my_tag2) 821 ldum = qmmm_ff_precond_only_qm(my_tag3) 822 WRITE (file_unit, FMT=psf_format) & 823 i, & 824 TRIM(my_tag1), & 825 my_index, & 826 TRIM(my_tag2), & 827 TRIM(my_tag3), & 828 TRIM(my_tag3), & 829 charges(i), & 830 atom_info%atm_mass(i), & 831 0 832 END DO 833 WRITE (file_unit, FMT='(/)') 834 DEALLOCATE (charges) 835 836 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%bond_a), " !NBOND" 837 DO i = 1, SIZE(conn_info%bond_a), 4 838 j = 0 839 DO WHILE ((j < 4) .AND. ((i + j) <= SIZE(conn_info%bond_a))) 840 WRITE (file_unit, FMT='(2('//TRIM(c_int)//'))', ADVANCE="NO") & 841 conn_info%bond_a(i + j), conn_info%bond_b(i + j) 842 j = j + 1 843 END DO 844 WRITE (file_unit, FMT='(/)', ADVANCE="NO") 845 END DO 846 WRITE (file_unit, FMT='(/)') 847 848 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%theta_a), " !NTHETA" 849 DO i = 1, SIZE(conn_info%theta_a), 3 850 j = 0 851 DO WHILE ((j < 3) .AND. ((i + j) <= SIZE(conn_info%theta_a))) 852 WRITE (file_unit, FMT='(3('//TRIM(c_int)//'))', ADVANCE="NO") & 853 conn_info%theta_a(i + j), conn_info%theta_b(i + j), & 854 conn_info%theta_c(i + j) 855 j = j + 1 856 END DO 857 WRITE (file_unit, FMT='(/)', ADVANCE="NO") 858 END DO 859 WRITE (file_unit, FMT='(/)') 860 861 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%phi_a), " !NPHI" 862 DO i = 1, SIZE(conn_info%phi_a), 2 863 j = 0 864 DO WHILE ((j < 2) .AND. ((i + j) <= SIZE(conn_info%phi_a))) 865 WRITE (file_unit, FMT='(4('//TRIM(c_int)//'))', ADVANCE="NO") & 866 conn_info%phi_a(i + j), conn_info%phi_b(i + j), & 867 conn_info%phi_c(i + j), conn_info%phi_d(i + j) 868 j = j + 1 869 END DO 870 WRITE (file_unit, FMT='(/)', ADVANCE="NO") 871 END DO 872 WRITE (file_unit, FMT='(/)') 873 874 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') SIZE(conn_info%impr_a), " !NIMPHI" 875 DO i = 1, SIZE(conn_info%impr_a), 2 876 j = 0 877 DO WHILE ((j < 2) .AND. ((i + j) <= SIZE(conn_info%impr_a))) 878 WRITE (file_unit, FMT='(4('//TRIM(c_int)//'))', ADVANCE="NO") & 879 conn_info%impr_a(i + j), conn_info%impr_b(i + j), & 880 conn_info%impr_c(i + j), conn_info%impr_d(i + j) 881 j = j + 1 882 END DO 883 WRITE (file_unit, FMT='(/)', ADVANCE="NO") 884 END DO 885 WRITE (file_unit, FMT='(/)') 886 887 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 0, " !NDON" 888 WRITE (file_unit, FMT='(/)') 889 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 0, " !NACC" 890 WRITE (file_unit, FMT='(/)') 891 WRITE (file_unit, FMT='('//TRIM(c_int)//',A)') 0, " !NNB" 892 WRITE (file_unit, FMT='(/)') 893 894 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 895 "PRINT%TOPOLOGY_INFO/PSF_INFO") 896 CALL timestop(handle) 897 898 END SUBROUTINE write_topology_psf 899 900END MODULE topology_psf 901 902