1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> Subroutine input_torsions changed (DG) 05-Dec-2000 9!> Output formats changed (DG) 05-Dec-2000 10!> JGH (26-01-2002) : force field parameters stored in tables, not in 11!> matrices. Input changed to have parameters labeled by the position 12!> and not atom pairs (triples etc) 13!> Teo (11.2005) : Moved all information on force field pair_potential to 14!> a much lighter memory structure 15!> Teo 09.2006 : Split all routines force_field I/O in a separate file 16!> 10.2008 Teodoro Laino [tlaino] : splitted from force_fields_input in order 17!> to collect in a unique module only the functions to read external FF 18!> \author CJM 19! ************************************************************************************************** 20MODULE force_fields_ext 21 USE cp_log_handling, ONLY: cp_get_default_logger,& 22 cp_logger_type 23 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 24 cp_print_key_unit_nr 25 USE cp_para_types, ONLY: cp_para_env_type 26 USE cp_parser_methods, ONLY: parser_get_next_line,& 27 parser_get_object,& 28 parser_search_string,& 29 parser_test_next_token 30 USE cp_parser_types, ONLY: cp_parser_type,& 31 parser_create,& 32 parser_release 33 USE cp_units, ONLY: cp_unit_to_cp2k 34 USE force_field_kind_types, ONLY: do_ff_g96 35 USE force_field_types, ONLY: amber_info_type,& 36 charmm_info_type,& 37 force_field_type,& 38 gromos_info_type 39 USE input_section_types, ONLY: section_vals_type 40 USE kinds, ONLY: default_string_length,& 41 dp 42 USE memory_utilities, ONLY: reallocate 43 USE particle_types, ONLY: particle_type 44 USE string_utilities, ONLY: uppercase 45 USE topology_amber, ONLY: rdparm_amber_8 46#include "./base/base_uses.f90" 47 48 IMPLICIT NONE 49 50 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'force_fields_ext' 51 52 PRIVATE 53 PUBLIC :: read_force_field_charmm, & 54 read_force_field_amber, & 55 read_force_field_gromos 56 57CONTAINS 58 59! ************************************************************************************************** 60!> \brief Reads the GROMOS force_field 61!> \param ff_type ... 62!> \param para_env ... 63!> \param mm_section ... 64!> \author ikuo 65! ************************************************************************************************** 66 SUBROUTINE read_force_field_gromos(ff_type, para_env, mm_section) 67 68 TYPE(force_field_type), INTENT(INOUT) :: ff_type 69 TYPE(cp_para_env_type), POINTER :: para_env 70 TYPE(section_vals_type), POINTER :: mm_section 71 72 CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_gromos', & 73 routineP = moduleN//':'//routineN 74 REAL(KIND=dp), PARAMETER :: ekt = 2.5_dp 75 76 CHARACTER(LEN=default_string_length) :: label 77 CHARACTER(LEN=default_string_length), & 78 DIMENSION(21) :: avail_section 79 CHARACTER(LEN=default_string_length), POINTER :: namearray(:) 80 INTEGER :: handle, iatom, icon, itemp, itype, iw, & 81 jatom, ncon, ntype, offset 82 LOGICAL :: found 83 REAL(KIND=dp) :: cosphi0, cost2, csq, sdet 84 TYPE(cp_logger_type), POINTER :: logger 85 TYPE(cp_parser_type), POINTER :: parser 86 TYPE(gromos_info_type), POINTER :: gro_info 87 88 CALL timeset(routineN, handle) 89 NULLIFY (logger, parser) 90 logger => cp_get_default_logger() 91 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", & 92 extension=".mmLog") 93 94 avail_section(1) = "TITLE" 95 avail_section(2) = "TOPPHYSCON" 96 avail_section(3) = "TOPVERSION" 97 avail_section(4) = "ATOMTYPENAME" 98 avail_section(5) = "RESNAME" 99 avail_section(6) = "SOLUTEATOM" 100 avail_section(7) = "BONDTYPE" 101 avail_section(8) = "BONDH" 102 avail_section(9) = "BOND" 103 avail_section(10) = "BONDANGLETYPE" 104 avail_section(11) = "BONDANGLEH" 105 avail_section(12) = "BONDANGLE" 106 avail_section(13) = "IMPDIHEDRALTYPE" 107 avail_section(14) = "IMPDIHEDRALH" 108 avail_section(15) = "IMPDIHEDRAL" 109 avail_section(16) = "DIHEDRALTYPE" 110 avail_section(17) = "DIHEDRALH" 111 avail_section(18) = "DIHEDRAL" 112 avail_section(19) = "LJPARAMETERS" 113 avail_section(20) = "SOLVENTATOM" 114 avail_section(21) = "SOLVENTCONSTR" 115 116 gro_info => ff_type%gro_info 117 gro_info%ff_gromos_type = ff_type%ff_type 118 NULLIFY (namearray) 119 ! ATOMTYPENAME SECTION 120 IF (iw > 0) WRITE (iw, '(T2,A)') 'GTOP_INFO| Parsing the ATOMTYPENAME section' 121 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 122 label = TRIM(avail_section(4)) 123 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 124 IF (found) THEN 125 CALL parser_get_next_line(parser, 1) 126 CALL parser_get_object(parser, ntype) 127 CALL reallocate(namearray, 1, ntype) 128 DO itype = 1, ntype 129 CALL parser_get_next_line(parser, 1) 130 CALL parser_get_object(parser, namearray(itype), lower_to_upper=.TRUE.) 131 IF (iw > 0) WRITE (iw, *) "GTOP_INFO| ", TRIM(namearray(itype)) 132 END DO 133 END IF 134 CALL parser_release(parser) 135 136 ! SOLVENTCONSTR SECTION 137 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the SOLVENTATOM section' 138 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 139 140 label = TRIM(avail_section(21)) 141 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 142 IF (found) THEN 143 CALL parser_get_next_line(parser, 1) 144 CALL parser_get_object(parser, ncon) 145 CALL reallocate(gro_info%solvent_k, 1, ncon) 146 CALL reallocate(gro_info%solvent_r0, 1, ncon) 147 DO icon = 1, ncon 148 CALL parser_get_next_line(parser, 1) 149 CALL parser_get_object(parser, itemp) 150 CALL parser_get_object(parser, itemp) 151 CALL parser_get_object(parser, gro_info%solvent_r0(icon)) 152 gro_info%solvent_k(icon) = 0.0_dp 153 END DO 154 END IF 155 CALL parser_release(parser) 156 157 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 158 ! BONDTYPE SECTION 159 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the BONDTYPE section' 160 label = TRIM(avail_section(7)) 161 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 162 IF (found) THEN 163 CALL parser_get_next_line(parser, 1) 164 CALL parser_get_object(parser, ntype) 165 CALL reallocate(gro_info%bond_k, 1, ntype) 166 CALL reallocate(gro_info%bond_r0, 1, ntype) 167 DO itype = 1, ntype 168 CALL parser_get_next_line(parser, 1) 169 CALL parser_get_object(parser, gro_info%bond_k(itype)) 170 CALL parser_get_object(parser, gro_info%bond_r0(itype)) 171 IF (ff_type%ff_type == do_ff_g96) THEN 172 gro_info%bond_k(itype) = cp_unit_to_cp2k(gro_info%bond_k(itype), "kjmol*nm^-4") 173 ELSE ! Assume its G87 174 gro_info%bond_k(itype) = (2.0_dp)*gro_info%bond_k(itype)*gro_info%bond_r0(itype)**2 175 gro_info%bond_k(itype) = cp_unit_to_cp2k(gro_info%bond_k(itype), "kjmol*nm^-2") 176 END IF 177 gro_info%bond_r0(itype) = cp_unit_to_cp2k(gro_info%bond_r0(itype), "nm") 178 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT BONDTYPE INFO HERE!!!!" 179 END DO 180 END IF 181 182 ! BONDANGLETYPE SECTION 183 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the BONDANGLETYPE section' 184 label = TRIM(avail_section(10)) 185 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 186 IF (found) THEN 187 CALL parser_get_next_line(parser, 1) 188 CALL parser_get_object(parser, ntype) 189 CALL reallocate(gro_info%bend_k, 1, ntype) 190 CALL reallocate(gro_info%bend_theta0, 1, ntype) 191 DO itype = 1, ntype 192 CALL parser_get_next_line(parser, 1) 193 CALL parser_get_object(parser, gro_info%bend_k(itype)) 194 CALL parser_get_object(parser, gro_info%bend_theta0(itype)) 195 gro_info%bend_theta0(itype) = cp_unit_to_cp2k(gro_info%bend_theta0(itype), "deg") 196 IF (ff_type%ff_type == do_ff_g96) THEN 197 gro_info%bend_theta0(itype) = COS(gro_info%bend_theta0(itype)) 198 ELSE ! Assume its G87 199 cost2 = COS(gro_info%bend_theta0(itype))*COS(gro_info%bend_theta0(itype)) 200 sdet = cost2*cost2 - (2.0_dp*cost2 - 1.0_dp)*(1.0_dp - ekt/gro_info%bend_k(itype)) 201 csq = (cost2 - SQRT(sdet))/(2.0_dp*cost2 - 1.0_dp) 202 gro_info%bend_k(itype) = ekt/ACOS(csq)**2 203 END IF 204 gro_info%bend_k(itype) = cp_unit_to_cp2k(gro_info%bend_k(itype), "kjmol") 205 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT BONDANGLETYPE INFO HERE!!!!" 206 END DO 207 END IF 208 209 ! IMPDIHEDRALTYPE SECTION 210 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the IMPDIHEDRALTYPE section' 211 label = TRIM(avail_section(13)) 212 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 213 IF (found) THEN 214 CALL parser_get_next_line(parser, 1) 215 CALL parser_get_object(parser, ntype) 216 CALL reallocate(gro_info%impr_k, 1, ntype) 217 CALL reallocate(gro_info%impr_phi0, 1, ntype) 218 DO itype = 1, ntype 219 CALL parser_get_next_line(parser, 1) 220 CALL parser_get_object(parser, gro_info%impr_k(itype)) 221 CALL parser_get_object(parser, gro_info%impr_phi0(itype)) 222 gro_info%impr_phi0(itype) = cp_unit_to_cp2k(gro_info%impr_phi0(itype), "deg") 223 gro_info%impr_k(itype) = cp_unit_to_cp2k(gro_info%impr_k(itype), "kjmol*deg^-2") 224 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT IMPDIHEDRALTYPE INFO HERE!!!!" 225 END DO 226 END IF 227 228 ! DIHEDRALTYPE SECTION 229 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the DIHEDRALTYPE section' 230 label = TRIM(avail_section(16)) 231 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 232 IF (found) THEN 233 CALL parser_get_next_line(parser, 1) 234 CALL parser_get_object(parser, ntype) 235 CALL reallocate(gro_info%torsion_k, 1, ntype) 236 CALL reallocate(gro_info%torsion_m, 1, ntype) 237 CALL reallocate(gro_info%torsion_phi0, 1, ntype) 238 DO itype = 1, ntype 239 CALL parser_get_next_line(parser, 1) 240 CALL parser_get_object(parser, gro_info%torsion_k(itype)) 241 CALL parser_get_object(parser, cosphi0) 242 CALL parser_get_object(parser, gro_info%torsion_m(itype)) 243 gro_info%torsion_phi0(itype) = ACOS(cosphi0) 244 gro_info%torsion_k(itype) = cp_unit_to_cp2k(gro_info%torsion_k(itype), "kjmol") 245 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT DIHEDRALTYPE INFO HERE!!!!" 246 END DO 247 END IF 248 249 CALL parser_release(parser) 250 251 ! LJPARAMETERS SECTION 252 IF (iw > 0) WRITE (iw, '(T2,A)') 'GROMOS_FF| Parsing the LJPARAMETERS section' 253 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 254 label = TRIM(avail_section(19)) 255 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 256 IF (found) THEN 257 CALL parser_get_next_line(parser, 1) 258 CALL parser_get_object(parser, ntype) 259 offset = 0 260 IF (ASSOCIATED(gro_info%nonbond_a)) offset = SIZE(gro_info%nonbond_a) 261 ntype = SIZE(namearray) 262 CALL reallocate(gro_info%nonbond_a, 1, ntype) 263 CALL reallocate(gro_info%nonbond_a_14, 1, ntype) 264 CALL reallocate(gro_info%nonbond_c6, 1, ntype, 1, ntype) 265 CALL reallocate(gro_info%nonbond_c12, 1, ntype, 1, ntype) 266 CALL reallocate(gro_info%nonbond_c6_14, 1, ntype, 1, ntype) 267 CALL reallocate(gro_info%nonbond_c12_14, 1, ntype, 1, ntype) 268 269 gro_info%nonbond_c12 = 0._dp 270 gro_info%nonbond_c6 = 0._dp 271 gro_info%nonbond_c12_14 = 0._dp 272 gro_info%nonbond_c6_14 = 0._dp 273 274 DO itype = 1, ntype*(ntype + 1)/2 275 CALL parser_get_next_line(parser, 1) 276 CALL parser_get_object(parser, iatom) 277 CALL parser_get_object(parser, jatom) 278 IF (iatom == jatom) THEN 279 gro_info%nonbond_a(iatom) = namearray(iatom) 280 gro_info%nonbond_a_14(iatom) = namearray(iatom) 281 END IF 282 CALL parser_get_object(parser, gro_info%nonbond_c12(iatom, jatom)) 283 CALL parser_get_object(parser, gro_info%nonbond_c6(iatom, jatom)) 284 CALL parser_get_object(parser, gro_info%nonbond_c12_14(iatom, jatom)) 285 CALL parser_get_object(parser, gro_info%nonbond_c6_14(iatom, jatom)) 286 gro_info%nonbond_c6(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c6(iatom, jatom), "kjmol*nm^6") 287 gro_info%nonbond_c12(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c12(iatom, jatom), "kjmol*nm^12") 288 gro_info%nonbond_c6_14(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c6_14(iatom, jatom), "kjmol*nm^6") 289 gro_info%nonbond_c12_14(iatom, jatom) = cp_unit_to_cp2k(gro_info%nonbond_c12_14(iatom, jatom), "kjmol*nm^12") 290 291 gro_info%nonbond_c6_14(jatom, iatom) = gro_info%nonbond_c6_14(iatom, jatom) 292 gro_info%nonbond_c12_14(jatom, iatom) = gro_info%nonbond_c12_14(iatom, jatom) 293 gro_info%nonbond_c6(jatom, iatom) = gro_info%nonbond_c6(iatom, jatom) 294 gro_info%nonbond_c12(jatom, iatom) = gro_info%nonbond_c12(iatom, jatom) 295 IF (iw > 0) WRITE (iw, *) "GROMOS_FF| PUT LJPARAMETERS INFO HERE!!!!" 296 END DO 297 END IF 298 CALL parser_release(parser) 299 300 CALL cp_print_key_finished_output(iw, logger, mm_section, & 301 "PRINT%FF_INFO") 302 CALL timestop(handle) 303 304 DEALLOCATE (namearray) 305 END SUBROUTINE read_force_field_gromos 306 307! ************************************************************************************************** 308!> \brief Reads the charmm force_field 309!> \param ff_type ... 310!> \param para_env ... 311!> \param mm_section ... 312!> \author ikuo 313! ************************************************************************************************** 314 SUBROUTINE read_force_field_charmm(ff_type, para_env, mm_section) 315 316 TYPE(force_field_type), INTENT(INOUT) :: ff_type 317 TYPE(cp_para_env_type), POINTER :: para_env 318 TYPE(section_vals_type), POINTER :: mm_section 319 320 CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_charmm', & 321 routineP = moduleN//':'//routineN 322 323 CHARACTER(LEN=default_string_length) :: label, string, string2, string3, string4 324 CHARACTER(LEN=default_string_length), DIMENSION(1) :: bond_section 325 CHARACTER(LEN=default_string_length), & 326 DIMENSION(19) :: avail_section 327 CHARACTER(LEN=default_string_length), DIMENSION(2) :: angl_section, impr_section, & 328 nbon_section, thet_section 329 INTEGER :: dummy, handle, ilab, iw, nbend, nbond, & 330 nimpr, nnonbond, nonfo, ntorsion, nub 331 LOGICAL :: found 332 TYPE(charmm_info_type), POINTER :: chm_info 333 TYPE(cp_logger_type), POINTER :: logger 334 TYPE(cp_parser_type), POINTER :: parser 335 336 CALL timeset(routineN, handle) 337 NULLIFY (logger, parser) 338 logger => cp_get_default_logger() 339 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", & 340 extension=".mmLog") 341 342 avail_section(1) = "BOND"; bond_section(1) = avail_section(1) 343 avail_section(11) = "BONDS" 344 avail_section(2) = "ANGL"; angl_section(1) = avail_section(2) 345 avail_section(3) = "THETA"; angl_section(2) = avail_section(3) 346 avail_section(12) = "THETAS" 347 avail_section(13) = "ANGLE" 348 avail_section(14) = "ANGLES" 349 avail_section(4) = "DIHEDRAL"; thet_section(1) = avail_section(4) 350 avail_section(15) = "DIHEDRALS" 351 avail_section(5) = "PHI"; thet_section(2) = avail_section(5) 352 avail_section(6) = "IMPROPER"; impr_section(1) = avail_section(6) 353 avail_section(7) = "IMPH"; impr_section(2) = avail_section(7) 354 avail_section(16) = "IMPHI" 355 avail_section(8) = "NONBONDED"; nbon_section(1) = avail_section(8) 356 avail_section(9) = "NBOND"; nbon_section(2) = avail_section(9) 357 avail_section(10) = "HBOND" 358 avail_section(17) = "NBFIX" 359 avail_section(18) = "CMAP" 360 avail_section(19) = "END" 361 362 chm_info => ff_type%chm_info 363 !----------------------------------------------------------------------------- 364 !----------------------------------------------------------------------------- 365 ! 1. Read in all the Bonds info from the param file here 366 ! Vbond = Kb(b-b0)^2 367 ! UNITS for Kb: [(kcal/mol)/(A^2)] to [Eh/(AU^2)] 368 !----------------------------------------------------------------------------- 369 nbond = 0 370 DO ilab = 1, SIZE(bond_section) 371 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 372 label = TRIM(bond_section(ilab)) 373 DO 374 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 375 IF (found) THEN 376 CALL parser_get_object(parser, string) 377 IF (INDEX(string, TRIM(label)) /= 1) CYCLE 378 CALL charmm_get_next_line(parser, 1) 379 DO 380 CALL parser_get_object(parser, string) 381 CALL uppercase(string) 382 IF (ANY(string == avail_section)) EXIT 383 CALL parser_get_object(parser, string2) 384 CALL uppercase(string2) 385 nbond = nbond + 1 386 CALL reallocate(chm_info%bond_a, 1, nbond) 387 CALL reallocate(chm_info%bond_b, 1, nbond) 388 CALL reallocate(chm_info%bond_k, 1, nbond) 389 CALL reallocate(chm_info%bond_r0, 1, nbond) 390 chm_info%bond_a(nbond) = string 391 chm_info%bond_b(nbond) = string2 392 CALL parser_get_object(parser, chm_info%bond_k(nbond)) 393 CALL parser_get_object(parser, chm_info%bond_r0(nbond)) 394 IF (iw > 0) WRITE (iw, *) " CHM BOND ", nbond, & 395 TRIM(chm_info%bond_a(nbond)), " ", & 396 TRIM(chm_info%bond_b(nbond)), " ", & 397 chm_info%bond_k(nbond), & 398 chm_info%bond_r0(nbond) 399 ! Do some units conversion into internal atomic units 400 chm_info%bond_r0(nbond) = cp_unit_to_cp2k(chm_info%bond_r0(nbond), "angstrom") 401 chm_info%bond_k(nbond) = cp_unit_to_cp2k(chm_info%bond_k(nbond), "kcalmol*angstrom^-2") 402 CALL charmm_get_next_line(parser, 1) 403 END DO 404 ELSE 405 EXIT 406 END IF 407 END DO 408 CALL parser_release(parser) 409 END DO 410 !----------------------------------------------------------------------------- 411 !----------------------------------------------------------------------------- 412 ! 2. Read in all the Bends and UB info from the param file here 413 ! Vangle = Ktheta(theta-theta0)^2 414 ! UNITS for Ktheta: [(kcal/mol)/(rad^2)] to [Eh/(rad^2)] 415 ! FACTOR of "2" rolled into Ktheta 416 ! Vub = Kub(S-S0)^2 417 ! UNITS for Kub: [(kcal/mol)/(A^2)] to [Eh/(AU^2)] 418 !----------------------------------------------------------------------------- 419 nbend = 0 420 nub = 0 421 DO ilab = 1, SIZE(angl_section) 422 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 423 label = TRIM(angl_section(ilab)) 424 DO 425 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 426 IF (found) THEN 427 CALL parser_get_object(parser, string) 428 IF (INDEX(string, TRIM(label)) /= 1) CYCLE 429 CALL charmm_get_next_line(parser, 1) 430 DO 431 CALL parser_get_object(parser, string) 432 CALL uppercase(string) 433 IF (ANY(string == avail_section)) EXIT 434 CALL parser_get_object(parser, string2) 435 CALL parser_get_object(parser, string3) 436 CALL uppercase(string2) 437 CALL uppercase(string3) 438 nbend = nbend + 1 439 CALL reallocate(chm_info%bend_a, 1, nbend) 440 CALL reallocate(chm_info%bend_b, 1, nbend) 441 CALL reallocate(chm_info%bend_c, 1, nbend) 442 CALL reallocate(chm_info%bend_k, 1, nbend) 443 CALL reallocate(chm_info%bend_theta0, 1, nbend) 444 chm_info%bend_a(nbend) = string 445 chm_info%bend_b(nbend) = string2 446 chm_info%bend_c(nbend) = string3 447 CALL parser_get_object(parser, chm_info%bend_k(nbend)) 448 CALL parser_get_object(parser, chm_info%bend_theta0(nbend)) 449 IF (iw > 0) WRITE (iw, *) " CHM BEND ", nbend, & 450 TRIM(chm_info%bend_a(nbend)), " ", & 451 TRIM(chm_info%bend_b(nbend)), " ", & 452 TRIM(chm_info%bend_c(nbend)), " ", & 453 chm_info%bend_k(nbend), & 454 chm_info%bend_theta0(nbend) 455 ! Do some units conversion into internal atomic units 456 chm_info%bend_theta0(nbend) = cp_unit_to_cp2k(chm_info%bend_theta0(nbend), "deg") 457 chm_info%bend_k(nbend) = cp_unit_to_cp2k(chm_info%bend_k(nbend), "kcalmol*rad^-2") 458 IF (parser_test_next_token(parser) == "FLT") THEN 459 nub = nub + 1 460 CALL reallocate(chm_info%ub_a, 1, nub) 461 CALL reallocate(chm_info%ub_b, 1, nub) 462 CALL reallocate(chm_info%ub_c, 1, nub) 463 CALL reallocate(chm_info%ub_k, 1, nub) 464 CALL reallocate(chm_info%ub_r0, 1, nub) 465 chm_info%ub_a(nub) = string 466 chm_info%ub_b(nub) = string2 467 chm_info%ub_c(nub) = string3 468 CALL parser_get_object(parser, chm_info%ub_k(nub)) 469 CALL parser_get_object(parser, chm_info%ub_r0(nub)) 470 IF (iw > 0) WRITE (iw, *) " CHM UB ", nub, & 471 TRIM(chm_info%ub_a(nub)), " ", & 472 TRIM(chm_info%ub_b(nub)), " ", & 473 TRIM(chm_info%ub_c(nub)), " ", & 474 chm_info%ub_k(nub), & 475 chm_info%ub_r0(nub) 476 ! Do some units conversion into internal atomic units 477 chm_info%ub_r0(nub) = cp_unit_to_cp2k(chm_info%ub_r0(nub), "angstrom") 478 chm_info%ub_k(nub) = cp_unit_to_cp2k(chm_info%ub_k(nub), "kcalmol*angstrom^-2") 479 END IF 480 CALL charmm_get_next_line(parser, 1) 481 END DO 482 ELSE 483 EXIT 484 END IF 485 END DO 486 CALL parser_release(parser) 487 END DO 488 !----------------------------------------------------------------------------- 489 !----------------------------------------------------------------------------- 490 ! 3. Read in all the Dihedrals info from the param file here 491 ! Vtorsion = Kphi(1+COS(n(phi)-delta)) 492 ! UNITS for Kphi: [(kcal/mol)] to [Eh] 493 !----------------------------------------------------------------------------- 494 ntorsion = 0 495 DO ilab = 1, SIZE(thet_section) 496 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 497 label = TRIM(thet_section(ilab)) 498 DO 499 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 500 IF (found) THEN 501 CALL parser_get_object(parser, string) 502 IF (INDEX(string, TRIM(label)) /= 1) CYCLE 503 CALL charmm_get_next_line(parser, 1) 504 DO 505 CALL parser_get_object(parser, string) 506 CALL uppercase(string) 507 IF (ANY(string == avail_section)) EXIT 508 CALL parser_get_object(parser, string2) 509 CALL parser_get_object(parser, string3) 510 CALL parser_get_object(parser, string4) 511 CALL uppercase(string2) 512 CALL uppercase(string3) 513 CALL uppercase(string4) 514 ntorsion = ntorsion + 1 515 CALL reallocate(chm_info%torsion_a, 1, ntorsion) 516 CALL reallocate(chm_info%torsion_b, 1, ntorsion) 517 CALL reallocate(chm_info%torsion_c, 1, ntorsion) 518 CALL reallocate(chm_info%torsion_d, 1, ntorsion) 519 CALL reallocate(chm_info%torsion_k, 1, ntorsion) 520 CALL reallocate(chm_info%torsion_m, 1, ntorsion) 521 CALL reallocate(chm_info%torsion_phi0, 1, ntorsion) 522 chm_info%torsion_a(ntorsion) = string 523 chm_info%torsion_b(ntorsion) = string2 524 chm_info%torsion_c(ntorsion) = string3 525 chm_info%torsion_d(ntorsion) = string4 526 CALL parser_get_object(parser, chm_info%torsion_k(ntorsion)) 527 CALL parser_get_object(parser, chm_info%torsion_m(ntorsion)) 528 CALL parser_get_object(parser, chm_info%torsion_phi0(ntorsion)) 529 IF (iw > 0) WRITE (iw, *) " CHM TORSION ", ntorsion, & 530 TRIM(chm_info%torsion_a(ntorsion)), " ", & 531 TRIM(chm_info%torsion_b(ntorsion)), " ", & 532 TRIM(chm_info%torsion_c(ntorsion)), " ", & 533 TRIM(chm_info%torsion_d(ntorsion)), " ", & 534 chm_info%torsion_k(ntorsion), & 535 chm_info%torsion_m(ntorsion), & 536 chm_info%torsion_phi0(ntorsion) 537 ! Do some units conversion into internal atomic units 538 chm_info%torsion_phi0(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_phi0(ntorsion), & 539 "deg") 540 chm_info%torsion_k(ntorsion) = cp_unit_to_cp2k(chm_info%torsion_k(ntorsion), "kcalmol") 541 CALL charmm_get_next_line(parser, 1) 542 END DO 543 ELSE 544 EXIT 545 END IF 546 END DO 547 CALL parser_release(parser) 548 END DO 549 !----------------------------------------------------------------------------- 550 !----------------------------------------------------------------------------- 551 ! 4. Read in all the Improper info from the param file here 552 ! Vimpr = Kpsi(psi-psi0)^2 553 ! UNITS for Kpsi: [(kcal/mol)/(rad^2)] to [Eh/(rad^2)] 554 !----------------------------------------------------------------------------- 555 nimpr = 0 556 DO ilab = 1, SIZE(impr_section) 557 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 558 label = TRIM(impr_section(ilab)) 559 DO 560 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 561 IF (found) THEN 562 CALL parser_get_object(parser, string) 563 IF (INDEX(string, TRIM(label)) /= 1) CYCLE 564 CALL charmm_get_next_line(parser, 1) 565 DO 566 CALL parser_get_object(parser, string) 567 CALL uppercase(string) 568 IF (ANY(string == avail_section)) EXIT 569 CALL parser_get_object(parser, string2) 570 CALL parser_get_object(parser, string3) 571 CALL parser_get_object(parser, string4) 572 CALL uppercase(string2) 573 CALL uppercase(string3) 574 CALL uppercase(string4) 575 nimpr = nimpr + 1 576 CALL reallocate(chm_info%impr_a, 1, nimpr) 577 CALL reallocate(chm_info%impr_b, 1, nimpr) 578 CALL reallocate(chm_info%impr_c, 1, nimpr) 579 CALL reallocate(chm_info%impr_d, 1, nimpr) 580 CALL reallocate(chm_info%impr_k, 1, nimpr) 581 CALL reallocate(chm_info%impr_phi0, 1, nimpr) 582 chm_info%impr_a(nimpr) = string 583 chm_info%impr_b(nimpr) = string2 584 chm_info%impr_c(nimpr) = string3 585 chm_info%impr_d(nimpr) = string4 586 CALL parser_get_object(parser, chm_info%impr_k(nimpr)) 587 CALL parser_get_object(parser, dummy) 588 CALL parser_get_object(parser, chm_info%impr_phi0(nimpr)) 589 IF (iw > 0) WRITE (iw, *) " CHM IMPROPERS ", nimpr, & 590 TRIM(chm_info%impr_a(nimpr)), " ", & 591 TRIM(chm_info%impr_b(nimpr)), " ", & 592 TRIM(chm_info%impr_c(nimpr)), " ", & 593 TRIM(chm_info%impr_d(nimpr)), " ", & 594 chm_info%impr_k(nimpr), & 595 chm_info%impr_phi0(nimpr) 596 ! Do some units conversion into internal atomic units 597 chm_info%impr_phi0(nimpr) = cp_unit_to_cp2k(chm_info%impr_phi0(nimpr), "deg") 598 chm_info%impr_k(nimpr) = cp_unit_to_cp2k(chm_info%impr_k(nimpr), "kcalmol") 599 CALL charmm_get_next_line(parser, 1) 600 END DO 601 ELSE 602 EXIT 603 END IF 604 END DO 605 CALL parser_release(parser) 606 END DO 607 !----------------------------------------------------------------------------- 608 !----------------------------------------------------------------------------- 609 ! 5. Read in all the Nonbonded info from the param file here 610 !----------------------------------------------------------------------------- 611 nnonbond = 0 612 nonfo = 0 613 DO ilab = 1, SIZE(nbon_section) 614 CALL parser_create(parser, ff_type%ff_file_name, para_env=para_env) 615 label = TRIM(nbon_section(ilab)) 616 DO 617 CALL parser_search_string(parser, label, .TRUE., found, begin_line=.TRUE.) 618 IF (found) THEN 619 CALL parser_get_object(parser, string) 620 IF (INDEX(string, TRIM(label)) /= 1) CYCLE 621 CALL charmm_get_next_line(parser, 1) 622 DO 623 CALL parser_get_object(parser, string) 624 CALL uppercase(string) 625 IF (ANY(string == avail_section)) EXIT 626 nnonbond = nnonbond + 1 627 CALL reallocate(chm_info%nonbond_a, 1, nnonbond) 628 CALL reallocate(chm_info%nonbond_eps, 1, nnonbond) 629 CALL reallocate(chm_info%nonbond_rmin2, 1, nnonbond) 630 chm_info%nonbond_a(nnonbond) = string 631 CALL parser_get_object(parser, chm_info%nonbond_eps(nnonbond)) 632 CALL parser_get_object(parser, chm_info%nonbond_eps(nnonbond)) 633 CALL parser_get_object(parser, chm_info%nonbond_rmin2(nnonbond)) 634 IF (iw > 0) WRITE (iw, *) " CHM NONBOND ", nnonbond, & 635 TRIM(chm_info%nonbond_a(nnonbond)), " ", & 636 chm_info%nonbond_eps(nnonbond), & 637 chm_info%nonbond_rmin2(nnonbond) 638 chm_info%nonbond_rmin2(nnonbond) = cp_unit_to_cp2k(chm_info%nonbond_rmin2(nnonbond), & 639 "angstrom") 640 chm_info%nonbond_eps(nnonbond) = cp_unit_to_cp2k(chm_info%nonbond_eps(nnonbond), & 641 "kcalmol") 642 IF (parser_test_next_token(parser) == "FLT") THEN 643 nonfo = nonfo + 1 644 CALL reallocate(chm_info%nonbond_a_14, 1, nonfo) 645 CALL reallocate(chm_info%nonbond_eps_14, 1, nonfo) 646 CALL reallocate(chm_info%nonbond_rmin2_14, 1, nonfo) 647 chm_info%nonbond_a_14(nonfo) = chm_info%nonbond_a(nnonbond) 648 CALL parser_get_object(parser, chm_info%nonbond_eps_14(nonfo)) 649 CALL parser_get_object(parser, chm_info%nonbond_eps_14(nonfo)) 650 CALL parser_get_object(parser, chm_info%nonbond_rmin2_14(nonfo)) 651 IF (iw > 0) WRITE (iw, *) " CHM ONFO ", nonfo, & 652 TRIM(chm_info%nonbond_a_14(nonfo)), " ", & 653 chm_info%nonbond_eps_14(nonfo), & 654 chm_info%nonbond_rmin2_14(nonfo) 655 chm_info%nonbond_rmin2_14(nonfo) = cp_unit_to_cp2k(chm_info%nonbond_rmin2_14(nonfo), & 656 "angstrom") 657 chm_info%nonbond_eps_14(nonfo) = cp_unit_to_cp2k(chm_info%nonbond_eps_14(nonfo), & 658 "kcalmol") 659 END IF 660 CALL charmm_get_next_line(parser, 1) 661 END DO 662 ELSE 663 EXIT 664 END IF 665 END DO 666 CALL parser_release(parser) 667 END DO 668 CALL cp_print_key_finished_output(iw, logger, mm_section, & 669 "PRINT%FF_INFO") 670 CALL timestop(handle) 671 672 END SUBROUTINE read_force_field_charmm 673 674! ************************************************************************************************** 675!> \brief Reads the AMBER force_field 676!> \param ff_type ... 677!> \param para_env ... 678!> \param mm_section ... 679!> \param particle_set ... 680!> \author Teodoro Laino [tlaino, teodoro.laino-AT-gmail.com] - 11.2008 681! ************************************************************************************************** 682 SUBROUTINE read_force_field_amber(ff_type, para_env, mm_section, particle_set) 683 684 TYPE(force_field_type), INTENT(INOUT) :: ff_type 685 TYPE(cp_para_env_type), POINTER :: para_env 686 TYPE(section_vals_type), POINTER :: mm_section 687 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 688 689 CHARACTER(len=*), PARAMETER :: routineN = 'read_force_field_amber', & 690 routineP = moduleN//':'//routineN 691 692 INTEGER :: handle, i, iw 693 TYPE(amber_info_type), POINTER :: amb_info 694 TYPE(cp_logger_type), POINTER :: logger 695 696 CALL timeset(routineN, handle) 697 NULLIFY (logger) 698 logger => cp_get_default_logger() 699 iw = cp_print_key_unit_nr(logger, mm_section, "PRINT%FF_INFO", & 700 extension=".mmLog") 701 702 amb_info => ff_type%amb_info 703 704 ! Read the Amber topology file to gather the forcefield information 705 CALL rdparm_amber_8(ff_type%ff_file_name, iw, para_env, do_connectivity=.FALSE., & 706 do_forcefield=.TRUE., particle_set=particle_set, amb_info=amb_info) 707 708 !----------------------------------------------------------------------------- 709 ! 1. Converts all the Bonds info from the param file here 710 ! Vbond = Kb(b-b0)^2 711 ! UNITS for Kb: [(kcal/mol)/(A^2)] to [Eh/(AU^2)] 712 !----------------------------------------------------------------------------- 713 DO i = 1, SIZE(amb_info%bond_a) 714 IF (iw > 0) WRITE (iw, *) " AMB BOND ", i, & 715 TRIM(amb_info%bond_a(i)), " ", & 716 TRIM(amb_info%bond_b(i)), " ", & 717 amb_info%bond_k(i), & 718 amb_info%bond_r0(i) 719 720 ! Do some units conversion into internal atomic units 721 amb_info%bond_r0(i) = cp_unit_to_cp2k(amb_info%bond_r0(i), "angstrom") 722 amb_info%bond_k(i) = cp_unit_to_cp2k(amb_info%bond_k(i), "kcalmol*angstrom^-2") 723 END DO 724 725 !----------------------------------------------------------------------------- 726 ! 2. Converts all the Bends info from the param file here 727 ! Vangle = Ktheta(theta-theta0)^2 728 ! UNITS for Ktheta: [(kcal/mol)/(rad^2)] to [Eh/(rad^2)] 729 ! FACTOR of "2" rolled into Ktheta 730 ! Vub = Kub(S-S0)^2 731 ! UNITS for Kub: [(kcal/mol)/(A^2)] to [Eh/(AU^2)] 732 !----------------------------------------------------------------------------- 733 DO i = 1, SIZE(amb_info%bend_a) 734 IF (iw > 0) WRITE (iw, *) " AMB BEND ", i, & 735 TRIM(amb_info%bend_a(i)), " ", & 736 TRIM(amb_info%bend_b(i)), " ", & 737 TRIM(amb_info%bend_c(i)), " ", & 738 amb_info%bend_k(i), & 739 amb_info%bend_theta0(i) 740 741 ! Do some units conversion into internal atomic units 742 amb_info%bend_theta0(i) = cp_unit_to_cp2k(amb_info%bend_theta0(i), "rad") 743 amb_info%bend_k(i) = cp_unit_to_cp2k(amb_info%bend_k(i), "kcalmol*rad^-2") 744 END DO 745 746 !----------------------------------------------------------------------------- 747 ! 3. Converts all the Dihedrals info from the param file here 748 ! Vtorsion = Kphi(1+COS(n(phi)-delta)) 749 ! UNITS for Kphi: [(kcal/mol)] to [Eh] 750 !----------------------------------------------------------------------------- 751 DO i = 1, SIZE(amb_info%torsion_a) 752 IF (iw > 0) WRITE (iw, *) " AMB TORSION ", i, & 753 TRIM(amb_info%torsion_a(i)), " ", & 754 TRIM(amb_info%torsion_b(i)), " ", & 755 TRIM(amb_info%torsion_c(i)), " ", & 756 TRIM(amb_info%torsion_d(i)), " ", & 757 amb_info%torsion_k(i), & 758 amb_info%torsion_m(i), & 759 amb_info%torsion_phi0(i) 760 761 ! Do some units conversion into internal atomic units 762 amb_info%torsion_phi0(i) = cp_unit_to_cp2k(amb_info%torsion_phi0(i), "rad") 763 amb_info%torsion_k(i) = cp_unit_to_cp2k(amb_info%torsion_k(i), "kcalmol") 764 END DO 765 766 !----------------------------------------------------------------------------- 767 ! 4. Converts all the Nonbonded info from the param file here 768 !----------------------------------------------------------------------------- 769 DO i = 1, SIZE(amb_info%nonbond_eps) 770 IF (iw > 0) WRITE (iw, *) " AMB NONBOND ", i, & 771 TRIM(amb_info%nonbond_a(i)), " ", & 772 amb_info%nonbond_eps(i), & 773 amb_info%nonbond_rmin2(i) 774 775 ! Do some units conversion into internal atomic units 776 amb_info%nonbond_rmin2(i) = cp_unit_to_cp2k(amb_info%nonbond_rmin2(i), "angstrom") 777 amb_info%nonbond_eps(i) = cp_unit_to_cp2k(amb_info%nonbond_eps(i), "kcalmol") 778 END DO 779 CALL cp_print_key_finished_output(iw, logger, mm_section, "PRINT%FF_INFO") 780 CALL timestop(handle) 781 END SUBROUTINE read_force_field_amber 782 783! ************************************************************************************************** 784!> \brief This function is simply a wrap to the parser_get_next_line.. 785!> Comments: This routine would not be necessary if the continuation 786!> char for CHARMM would not be the "-".. How can you choose this 787!> character in a file of numbers as a continuation char???? 788!> This sounds simply crazy.... 789!> \param parser ... 790!> \param nline ... 791!> \author Teodoro Laino - Zurich University - 06.2007 792! ************************************************************************************************** 793 SUBROUTINE charmm_get_next_line(parser, nline) 794 TYPE(cp_parser_type), POINTER :: parser 795 INTEGER, INTENT(IN) :: nline 796 797 CHARACTER(len=*), PARAMETER :: routineN = 'charmm_get_next_line', & 798 routineP = moduleN//':'//routineN 799 CHARACTER(LEN=1), PARAMETER :: continuation_char = "-" 800 801 INTEGER :: i, len_line 802 803 DO i = 1, nline 804 len_line = LEN_TRIM(parser%input_line) 805 DO WHILE (parser%input_line(len_line:len_line) == continuation_char) 806 CALL parser_get_next_line(parser, 1) 807 len_line = LEN_TRIM(parser%input_line) 808 END DO 809 CALL parser_get_next_line(parser, 1) 810 END DO 811 812 END SUBROUTINE charmm_get_next_line 813 814END MODULE force_fields_ext 815