1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Reads the input sections "topology" 8!> \par History 9!> JGH (26-01-2002) Added read_topology_section 10!> \author JGH 11! ************************************************************************************************** 12MODULE topology_input 13 USE colvar_types, ONLY: colvar_clone,& 14 colvar_p_type 15 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit,& 16 cp_to_string 17 USE input_constants, ONLY: do_conn_generate,& 18 do_conn_mol_set,& 19 do_conn_off,& 20 do_conn_user,& 21 do_constr_none,& 22 do_coord_off 23 USE input_section_types, ONLY: section_vals_get,& 24 section_vals_get_subs_vals,& 25 section_vals_type,& 26 section_vals_val_get,& 27 section_vals_val_unset 28 USE kinds, ONLY: default_string_length,& 29 dp 30 USE memory_utilities, ONLY: reallocate 31 USE topology_types, ONLY: constraint_info_type,& 32 topology_parameters_type 33#include "./base/base_uses.f90" 34 35 IMPLICIT NONE 36 37 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_input' 38 39 PRIVATE 40 PUBLIC :: read_topology_section, read_constraints_section 41 42CONTAINS 43 44! ************************************************************************************************** 45!> \brief reads the input section topology 46!> \param topology ... 47!> \param topology_section ... 48!> \par History 49!> none 50!> \author JGH (26-01-2002) 51! ************************************************************************************************** 52 SUBROUTINE read_topology_section(topology, topology_section) 53 TYPE(topology_parameters_type) :: topology 54 TYPE(section_vals_type), POINTER :: topology_section 55 56 CHARACTER(len=*), PARAMETER :: routineN = 'read_topology_section', & 57 routineP = moduleN//':'//routineN 58 59 INTEGER :: handle, ival 60 61 CALL timeset(routineN, handle) 62 CALL section_vals_val_get(topology_section, "CHARGE_OCCUP", l_val=topology%charge_occup) 63 CALL section_vals_val_get(topology_section, "CHARGE_BETA", l_val=topology%charge_beta) 64 CALL section_vals_val_get(topology_section, "CHARGE_EXTENDED", l_val=topology%charge_extended) 65 ival = COUNT((/topology%charge_occup, topology%charge_beta, topology%charge_extended/)) 66 IF (ival > 1) & 67 CPABORT("Only one between <CHARGE_OCCUP,CHARGE_BETA,CHARGE_EXTENDED> can be defined! ") 68 CALL section_vals_val_get(topology_section, "PARA_RES", l_val=topology%para_res) 69 CALL section_vals_val_get(topology_section, "GENERATE%REORDER", l_val=topology%reorder_atom) 70 CALL section_vals_val_get(topology_section, "GENERATE%CREATE_MOLECULES", l_val=topology%create_molecules) 71 CALL section_vals_val_get(topology_section, "MOL_CHECK", l_val=topology%molecules_check) 72 CALL section_vals_val_get(topology_section, "USE_G96_VELOCITY", l_val=topology%use_g96_velocity) 73 CALL section_vals_val_get(topology_section, "COORD_FILE_FORMAT", i_val=topology%coord_type) 74 SELECT CASE (topology%coord_type) 75 CASE (do_coord_off) 76 ! Do Nothing 77 CASE DEFAULT 78 topology%coordinate = .TRUE. 79 CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", c_val=topology%coord_file_name) 80 END SELECT 81 CALL section_vals_val_get(topology_section, "CONN_FILE_FORMAT", i_val=topology%conn_type) 82 SELECT CASE (topology%conn_type) 83 CASE (do_conn_off, do_conn_generate, do_conn_mol_set, do_conn_user) 84 ! Do Nothing 85 CASE DEFAULT 86 CALL section_vals_val_get(topology_section, "CONN_FILE_NAME", c_val=topology%conn_file_name) 87 END SELECT 88 CALL section_vals_val_get(topology_section, "EXCLUDE_VDW", i_val=topology%exclude_vdw) 89 CALL section_vals_val_get(topology_section, "EXCLUDE_EI", i_val=topology%exclude_ei) 90 CALL section_vals_val_get(topology_section, "GENERATE%BONDPARM", i_val=topology%bondparm_type) 91 CALL section_vals_val_get(topology_section, "GENERATE%BONDPARM_FACTOR", r_val=topology%bondparm_factor) 92 CALL timestop(handle) 93 END SUBROUTINE read_topology_section 94 95! ************************************************************************************************** 96!> \brief Read all the distance parameters. Put them in the 97!> constraint_distance array. 98!> \param topology ... 99!> \param colvar_p ... 100!> \param constraint_section ... 101!> \par History 102!> JGH (26-01-2002) Distance parameters are now stored in tables. The position 103!> within the table is used as handle for the topology 104!> teo Read the CONSTRAINT section within the new input style 105!> \author teo 106! ************************************************************************************************** 107 SUBROUTINE read_constraints_section(topology, colvar_p, constraint_section) 108 109 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 110 TYPE(colvar_p_type), DIMENSION(:), POINTER :: colvar_p 111 TYPE(section_vals_type), POINTER :: constraint_section 112 113 CHARACTER(len=*), PARAMETER :: routineN = 'read_constraints_section', & 114 routineP = moduleN//':'//routineN 115 116 CHARACTER(LEN=default_string_length), & 117 DIMENSION(:), POINTER :: tmpstringlist 118 INTEGER :: icolvar, ig, isize, isize_old, itype, & 119 jg, msize, msize_old, n_rep, ncons, & 120 nrep 121 INTEGER, DIMENSION(:), POINTER :: ilist, tmplist 122 LOGICAL :: explicit 123 REAL(KIND=dp), DIMENSION(:), POINTER :: rlist 124 TYPE(constraint_info_type), POINTER :: cons_info 125 TYPE(section_vals_type), POINTER :: collective_section, fix_atom_section, & 126 g3x3_section, g4x6_section, & 127 hbonds_section, vsite_section 128 129 cons_info => topology%cons_info 130 IF (ASSOCIATED(constraint_section)) THEN 131 hbonds_section => section_vals_get_subs_vals(constraint_section, "HBONDS") 132 g3x3_section => section_vals_get_subs_vals(constraint_section, "G3X3") 133 g4x6_section => section_vals_get_subs_vals(constraint_section, "G4X6") 134 vsite_section => section_vals_get_subs_vals(constraint_section, "VIRTUAL_SITE") 135 fix_atom_section => section_vals_get_subs_vals(constraint_section, "FIXED_ATOMS") 136 collective_section => section_vals_get_subs_vals(constraint_section, "COLLECTIVE") 137 ! HBONDS 138 CALL section_vals_get(hbonds_section, explicit=topology%const_hydr) 139 CALL check_restraint(hbonds_section, & 140 is_restraint=cons_info%hbonds_restraint, & 141 k0=cons_info%hbonds_k0, & 142 label="HBONDS") 143 ! G3X3 144 CALL section_vals_get(g3x3_section, explicit=explicit, n_repetition=ncons) 145 IF (explicit) THEN 146 topology%const_33 = .TRUE. 147 cons_info%nconst_g33 = ncons 148 ! 149 ALLOCATE (cons_info%const_g33_mol(ncons)) 150 ALLOCATE (cons_info%const_g33_molname(ncons)) 151 ALLOCATE (cons_info%const_g33_a(ncons)) 152 ALLOCATE (cons_info%const_g33_b(ncons)) 153 ALLOCATE (cons_info%const_g33_c(ncons)) 154 ALLOCATE (cons_info%const_g33_dab(ncons)) 155 ALLOCATE (cons_info%const_g33_dac(ncons)) 156 ALLOCATE (cons_info%const_g33_dbc(ncons)) 157 ALLOCATE (cons_info%g33_intermolecular(ncons)) 158 ALLOCATE (cons_info%g33_restraint(ncons)) 159 ALLOCATE (cons_info%g33_k0(ncons)) 160 ALLOCATE (cons_info%g33_exclude_qm(ncons)) 161 ALLOCATE (cons_info%g33_exclude_mm(ncons)) 162 DO ig = 1, ncons 163 CALL check_restraint(g3x3_section, & 164 is_restraint=cons_info%g33_restraint(ig), & 165 k0=cons_info%g33_k0(ig), & 166 i_rep_section=ig, & 167 label="G3X3") 168 cons_info%const_g33_mol(ig) = 0 169 cons_info%const_g33_molname(ig) = "UNDEF" 170 ! Exclude QM or MM 171 CALL section_vals_val_get(g3x3_section, "EXCLUDE_QM", i_rep_section=ig, & 172 l_val=cons_info%g33_exclude_qm(ig)) 173 CALL section_vals_val_get(g3x3_section, "EXCLUDE_MM", i_rep_section=ig, & 174 l_val=cons_info%g33_exclude_mm(ig)) 175 ! Intramolecular restraint 176 CALL section_vals_val_get(g3x3_section, "INTERMOLECULAR", i_rep_section=ig, & 177 l_val=cons_info%g33_intermolecular(ig)) 178 ! If it is intramolecular let's unset (in case user did it) 179 ! the molecule and molname field 180 IF (cons_info%g33_intermolecular(ig)) THEN 181 CALL section_vals_val_unset(g3x3_section, "MOLECULE", i_rep_section=ig) 182 CALL section_vals_val_unset(g3x3_section, "MOLNAME", i_rep_section=ig) 183 END IF 184 ! Let's tag to which molecule we want to apply constraints 185 CALL section_vals_val_get(g3x3_section, "MOLECULE", i_rep_section=ig, & 186 n_rep_val=nrep) 187 IF (nrep /= 0) THEN 188 CALL section_vals_val_get(g3x3_section, "MOLECULE", i_rep_section=ig, & 189 i_val=cons_info%const_g33_mol(ig)) 190 END IF 191 CALL section_vals_val_get(g3x3_section, "MOLNAME", i_rep_section=ig, & 192 n_rep_val=nrep) 193 IF (nrep /= 0) THEN 194 CALL section_vals_val_get(g3x3_section, "MOLNAME", i_rep_section=ig, & 195 c_val=cons_info%const_g33_molname(ig)) 196 END IF 197 IF ((cons_info%const_g33_mol(ig) /= 0) .AND. (cons_info%const_g33_molname(ig) /= "UNDEF")) THEN 198 CPABORT("") 199 END IF 200 IF ((cons_info%const_g33_mol(ig) == 0) .AND. (cons_info%const_g33_molname(ig) == "UNDEF") .AND. & 201 (.NOT. cons_info%g33_intermolecular(ig))) THEN 202 CPABORT("") 203 END IF 204 CALL section_vals_val_get(g3x3_section, "ATOMS", i_rep_section=ig, & 205 i_vals=ilist) 206 CALL section_vals_val_get(g3x3_section, "DISTANCES", i_rep_section=ig, & 207 r_vals=rlist) 208 cons_info%const_g33_a(ig) = ilist(1) 209 cons_info%const_g33_b(ig) = ilist(2) 210 cons_info%const_g33_c(ig) = ilist(3) 211 212 cons_info%const_g33_dab(ig) = rlist(1) 213 cons_info%const_g33_dac(ig) = rlist(2) 214 cons_info%const_g33_dbc(ig) = rlist(3) 215 END DO 216 END IF 217 ! G4X6 218 CALL section_vals_get(g4x6_section, explicit=explicit, n_repetition=ncons) 219 IF (explicit) THEN 220 topology%const_46 = .TRUE. 221 cons_info%nconst_g46 = ncons 222 ! 223 ALLOCATE (cons_info%const_g46_mol(ncons)) 224 ALLOCATE (cons_info%const_g46_molname(ncons)) 225 ALLOCATE (cons_info%const_g46_a(ncons)) 226 ALLOCATE (cons_info%const_g46_b(ncons)) 227 ALLOCATE (cons_info%const_g46_c(ncons)) 228 ALLOCATE (cons_info%const_g46_d(ncons)) 229 ALLOCATE (cons_info%const_g46_dab(ncons)) 230 ALLOCATE (cons_info%const_g46_dac(ncons)) 231 ALLOCATE (cons_info%const_g46_dbc(ncons)) 232 ALLOCATE (cons_info%const_g46_dad(ncons)) 233 ALLOCATE (cons_info%const_g46_dbd(ncons)) 234 ALLOCATE (cons_info%const_g46_dcd(ncons)) 235 ALLOCATE (cons_info%g46_intermolecular(ncons)) 236 ALLOCATE (cons_info%g46_restraint(ncons)) 237 ALLOCATE (cons_info%g46_k0(ncons)) 238 ALLOCATE (cons_info%g46_exclude_qm(ncons)) 239 ALLOCATE (cons_info%g46_exclude_mm(ncons)) 240 DO ig = 1, ncons 241 CALL check_restraint(g4x6_section, & 242 is_restraint=cons_info%g46_restraint(ig), & 243 k0=cons_info%g46_k0(ig), & 244 i_rep_section=ig, & 245 label="G4X6") 246 cons_info%const_g46_mol(ig) = 0 247 cons_info%const_g46_molname(ig) = "UNDEF" 248 ! Exclude QM or MM 249 CALL section_vals_val_get(g4x6_section, "EXCLUDE_QM", i_rep_section=ig, & 250 l_val=cons_info%g46_exclude_qm(ig)) 251 CALL section_vals_val_get(g4x6_section, "EXCLUDE_MM", i_rep_section=ig, & 252 l_val=cons_info%g46_exclude_mm(ig)) 253 ! Intramolecular restraint 254 CALL section_vals_val_get(g4x6_section, "INTERMOLECULAR", i_rep_section=ig, & 255 l_val=cons_info%g46_intermolecular(ig)) 256 ! If it is intramolecular let's unset (in case user did it) 257 ! the molecule and molname field 258 IF (cons_info%g46_intermolecular(ig)) THEN 259 CALL section_vals_val_unset(g4x6_section, "MOLECULE", i_rep_section=ig) 260 CALL section_vals_val_unset(g4x6_section, "MOLNAME", i_rep_section=ig) 261 END IF 262 ! Let's tag to which molecule we want to apply constraints 263 CALL section_vals_val_get(g4x6_section, "MOLECULE", i_rep_section=ig, & 264 n_rep_val=nrep) 265 IF (nrep /= 0) THEN 266 CALL section_vals_val_get(g4x6_section, "MOLECULE", i_rep_section=ig, & 267 i_val=cons_info%const_g46_mol(ig)) 268 END IF 269 CALL section_vals_val_get(g4x6_section, "MOLNAME", i_rep_section=ig, & 270 n_rep_val=nrep) 271 IF (nrep /= 0) THEN 272 CALL section_vals_val_get(g4x6_section, "MOLNAME", i_rep_section=ig, & 273 c_val=cons_info%const_g46_molname(ig)) 274 END IF 275 IF ((cons_info%const_g46_mol(ig) /= 0) .AND. (cons_info%const_g46_molname(ig) /= "UNDEF")) THEN 276 CPABORT("") 277 END IF 278 IF ((cons_info%const_g46_mol(ig) == 0) .AND. (cons_info%const_g46_molname(ig) == "UNDEF") .AND. & 279 (.NOT. cons_info%g46_intermolecular(ig))) THEN 280 CPABORT("") 281 END IF 282 CALL section_vals_val_get(g4x6_section, "ATOMS", i_rep_section=ig, & 283 i_vals=ilist) 284 CALL section_vals_val_get(g4x6_section, "DISTANCES", i_rep_section=ig, & 285 r_vals=rlist) 286 cons_info%const_g46_a(ig) = ilist(1) 287 cons_info%const_g46_b(ig) = ilist(2) 288 cons_info%const_g46_c(ig) = ilist(3) 289 cons_info%const_g46_d(ig) = ilist(4) 290 cons_info%const_g46_dab(ig) = rlist(1) 291 cons_info%const_g46_dac(ig) = rlist(2) 292 cons_info%const_g46_dad(ig) = rlist(3) 293 cons_info%const_g46_dbc(ig) = rlist(4) 294 cons_info%const_g46_dbd(ig) = rlist(5) 295 cons_info%const_g46_dcd(ig) = rlist(6) 296 END DO 297 END IF 298 ! virtual 299 CALL section_vals_get(vsite_section, explicit=explicit, n_repetition=ncons) 300 IF (explicit) THEN 301 topology%const_vsite = .TRUE. 302 cons_info%nconst_vsite = ncons 303 ! 304 ALLOCATE (cons_info%const_vsite_mol(ncons)) 305 ALLOCATE (cons_info%const_vsite_molname(ncons)) 306 ALLOCATE (cons_info%const_vsite_a(ncons)) 307 ALLOCATE (cons_info%const_vsite_b(ncons)) 308 ALLOCATE (cons_info%const_vsite_c(ncons)) 309 ALLOCATE (cons_info%const_vsite_d(ncons)) 310 ALLOCATE (cons_info%const_vsite_wbc(ncons)) 311 ALLOCATE (cons_info%const_vsite_wdc(ncons)) 312 ALLOCATE (cons_info%vsite_intermolecular(ncons)) 313 ALLOCATE (cons_info%vsite_restraint(ncons)) 314 ALLOCATE (cons_info%vsite_k0(ncons)) 315 ALLOCATE (cons_info%vsite_exclude_qm(ncons)) 316 ALLOCATE (cons_info%vsite_exclude_mm(ncons)) 317 DO ig = 1, ncons 318 CALL check_restraint(vsite_section, & 319 is_restraint=cons_info%vsite_restraint(ig), & 320 k0=cons_info%vsite_k0(ig), & 321 i_rep_section=ig, & 322 label="Virtual_SITE") 323 cons_info%const_vsite_mol(ig) = 0 324 cons_info%const_vsite_molname(ig) = "UNDEF" 325 ! Exclude QM or MM 326 CALL section_vals_val_get(vsite_section, "EXCLUDE_QM", i_rep_section=ig, & 327 l_val=cons_info%vsite_exclude_qm(ig)) 328 CALL section_vals_val_get(vsite_section, "EXCLUDE_MM", i_rep_section=ig, & 329 l_val=cons_info%vsite_exclude_mm(ig)) 330 ! Intramolecular restraint 331 CALL section_vals_val_get(vsite_section, "INTERMOLECULAR", i_rep_section=ig, & 332 l_val=cons_info%vsite_intermolecular(ig)) 333 ! If it is intramolecular let's unset (in case user did it) 334 ! the molecule and molname field 335 IF (cons_info%vsite_intermolecular(ig)) THEN 336 CALL section_vals_val_unset(vsite_section, "MOLECULE", i_rep_section=ig) 337 CALL section_vals_val_unset(vsite_section, "MOLNAME", i_rep_section=ig) 338 END IF 339 ! Let's tag to which molecule we want to apply constraints 340 CALL section_vals_val_get(vsite_section, "MOLECULE", i_rep_section=ig, & 341 n_rep_val=nrep) 342 IF (nrep /= 0) THEN 343 CALL section_vals_val_get(vsite_section, "MOLECULE", i_rep_section=ig, & 344 i_val=cons_info%const_vsite_mol(ig)) 345 END IF 346 CALL section_vals_val_get(vsite_section, "MOLNAME", i_rep_section=ig, & 347 n_rep_val=nrep) 348 IF (nrep /= 0) THEN 349 CALL section_vals_val_get(vsite_section, "MOLNAME", i_rep_section=ig, & 350 c_val=cons_info%const_vsite_molname(ig)) 351 END IF 352 IF ((cons_info%const_vsite_mol(ig) /= 0) .AND. (cons_info%const_vsite_molname(ig) /= "UNDEF")) THEN 353 CPABORT("") 354 END IF 355 IF ((cons_info%const_vsite_mol(ig) == 0) .AND. (cons_info%const_vsite_molname(ig) == "UNDEF") .AND. & 356 (.NOT. cons_info%vsite_intermolecular(ig))) THEN 357 CPABORT("") 358 END IF 359 CALL section_vals_val_get(vsite_section, "ATOMS", i_rep_section=ig, & 360 i_vals=ilist) 361 CALL section_vals_val_get(vsite_section, "PARAMETERS", i_rep_section=ig, & 362 r_vals=rlist) 363 cons_info%const_vsite_a(ig) = ilist(1) 364 cons_info%const_vsite_b(ig) = ilist(2) 365 cons_info%const_vsite_c(ig) = ilist(3) 366 cons_info%const_vsite_d(ig) = ilist(4) 367 cons_info%const_vsite_wbc(ig) = rlist(1) 368 cons_info%const_vsite_wdc(ig) = rlist(2) 369 END DO 370 END IF 371 ! FIXED ATOMS 372 CALL section_vals_get(fix_atom_section, explicit=explicit, n_repetition=ncons) 373 IF (explicit) THEN 374 NULLIFY (tmplist, tmpstringlist) 375 isize = 0 376 msize = 0 377 ALLOCATE (cons_info%fixed_atoms(isize)) 378 ALLOCATE (cons_info%fixed_type(isize)) 379 ALLOCATE (cons_info%fixed_restraint(isize)) 380 ALLOCATE (cons_info%fixed_k0(isize)) 381 ALLOCATE (cons_info%fixed_molnames(msize)) 382 ALLOCATE (cons_info%fixed_mol_type(isize)) 383 ALLOCATE (cons_info%fixed_mol_restraint(msize)) 384 ALLOCATE (cons_info%fixed_mol_k0(msize)) 385 ALLOCATE (cons_info%fixed_exclude_qm(ncons)) 386 ALLOCATE (cons_info%fixed_exclude_mm(ncons)) 387 DO ig = 1, ncons 388 isize_old = isize 389 msize_old = msize 390 CALL section_vals_val_get(fix_atom_section, "COMPONENTS_TO_FIX", i_rep_section=ig, & 391 i_val=itype) 392 CALL section_vals_val_get(fix_atom_section, "LIST", i_rep_section=ig, & 393 n_rep_val=n_rep) 394 DO jg = 1, n_rep 395 CALL section_vals_val_get(fix_atom_section, "LIST", i_rep_section=ig, & 396 i_rep_val=jg, i_vals=tmplist) 397 CALL reallocate(cons_info%fixed_atoms, 1, isize + SIZE(tmplist)) 398 cons_info%fixed_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist 399 CALL reallocate(cons_info%fixed_restraint, 1, isize + SIZE(tmplist)) 400 CALL reallocate(cons_info%fixed_k0, 1, isize + SIZE(tmplist)) 401 CALL reallocate(cons_info%fixed_type, 1, isize + SIZE(tmplist)) 402 cons_info%fixed_type(isize + 1:isize + SIZE(tmplist)) = itype 403 isize = SIZE(cons_info%fixed_atoms) 404 END DO 405 !Check for restraints 406 IF ((isize - isize_old) > 0) THEN 407 CALL check_restraint(fix_atom_section, & 408 is_restraint=cons_info%fixed_restraint(isize_old + 1), & 409 k0=cons_info%fixed_k0(isize_old + 1), & 410 i_rep_section=ig, & 411 label="FIXED ATOM") 412 cons_info%fixed_restraint(isize_old + 1:isize) = cons_info%fixed_restraint(isize_old + 1) 413 cons_info%fixed_k0(isize_old + 1:isize) = cons_info%fixed_k0(isize_old + 1) 414 END IF 415 CALL section_vals_val_get(fix_atom_section, "MOLNAME", i_rep_section=ig, & 416 n_rep_val=n_rep) 417 IF (n_rep /= 0) THEN 418 DO jg = 1, n_rep 419 CALL section_vals_val_get(fix_atom_section, "MOLNAME", i_rep_section=ig, & 420 i_rep_val=jg, c_vals=tmpstringlist) 421 CALL reallocate(cons_info%fixed_molnames, 1, msize + SIZE(tmpstringlist, 1)) 422 CALL reallocate(cons_info%fixed_mol_type, 1, msize + SIZE(tmpstringlist, 1)) 423 CALL reallocate(cons_info%fixed_mol_restraint, 1, msize + SIZE(tmpstringlist, 1)) 424 CALL reallocate(cons_info%fixed_mol_k0, 1, msize + SIZE(tmpstringlist, 1)) 425 cons_info%fixed_molnames(msize + 1:msize + SIZE(tmpstringlist, 1)) = tmpstringlist 426 cons_info%fixed_mol_type(msize + 1:msize + SIZE(tmpstringlist, 1)) = itype 427 msize = SIZE(cons_info%fixed_molnames) 428 END DO 429 ! Exclude QM or MM work only if defined MOLNAME 430 CALL reallocate(cons_info%fixed_exclude_qm, 1, msize) 431 CALL reallocate(cons_info%fixed_exclude_mm, 1, msize) 432 CALL section_vals_val_get(fix_atom_section, "EXCLUDE_QM", i_rep_section=ig, & 433 l_val=cons_info%fixed_exclude_qm(msize_old + 1)) 434 CALL section_vals_val_get(fix_atom_section, "EXCLUDE_MM", i_rep_section=ig, & 435 l_val=cons_info%fixed_exclude_mm(msize_old + 1)) 436 cons_info%fixed_exclude_qm(msize_old + 1:msize) = cons_info%fixed_exclude_qm(msize_old + 1) 437 cons_info%fixed_exclude_mm(msize_old + 1:msize) = cons_info%fixed_exclude_mm(msize_old + 1) 438 END IF 439 !Check for restraints 440 IF (n_rep /= 0) THEN 441 CALL check_restraint(fix_atom_section, & 442 is_restraint=cons_info%fixed_mol_restraint(msize_old + 1), & 443 k0=cons_info%fixed_mol_k0(msize_old + 1), & 444 i_rep_section=ig, & 445 label="FIXED ATOM") 446 cons_info%fixed_mol_restraint(msize_old + 1:msize) = cons_info%fixed_mol_restraint(msize_old + 1) 447 cons_info%fixed_mol_k0(msize_old + 1:msize) = cons_info%fixed_mol_k0(msize_old + 1) 448 END IF 449 CALL section_vals_val_get(fix_atom_section, "MM_SUBSYS", i_rep_section=ig, & 450 n_rep_val=nrep, explicit=explicit) 451 IF (nrep == 1 .AND. explicit) THEN 452 CPASSERT(cons_info%freeze_mm == do_constr_none) 453 CALL section_vals_val_get(fix_atom_section, "MM_SUBSYS", i_val=cons_info%freeze_mm, & 454 i_rep_section=ig) 455 cons_info%freeze_mm_type = itype 456 END IF 457 CALL section_vals_val_get(fix_atom_section, "QM_SUBSYS", i_rep_section=ig, & 458 n_rep_val=nrep, explicit=explicit) 459 IF (nrep == 1 .AND. explicit) THEN 460 CPASSERT(cons_info%freeze_qm == do_constr_none) 461 CALL section_vals_val_get(fix_atom_section, "QM_SUBSYS", i_val=cons_info%freeze_qm, & 462 i_rep_section=ig) 463 cons_info%freeze_qm_type = itype 464 END IF 465 IF (cons_info%freeze_mm /= do_constr_none) THEN 466 CALL check_restraint(fix_atom_section, & 467 is_restraint=cons_info%fixed_mm_restraint, & 468 k0=cons_info%fixed_mm_k0, & 469 i_rep_section=ig, & 470 label="FIXED ATOM") 471 END IF 472 IF (cons_info%freeze_qm /= do_constr_none) THEN 473 CALL check_restraint(fix_atom_section, & 474 is_restraint=cons_info%fixed_qm_restraint, & 475 k0=cons_info%fixed_qm_k0, & 476 i_rep_section=ig, & 477 label="FIXED ATOM") 478 END IF 479 480 END DO 481 IF ((isize /= 0) .OR. (msize /= 0) .OR. & 482 (cons_info%freeze_mm /= do_constr_none) .OR. & 483 (cons_info%freeze_qm /= do_constr_none)) THEN 484 topology%const_atom = .TRUE. 485 END IF 486 END IF 487 ! Collective Constraints 488 CALL section_vals_get(collective_section, explicit=explicit, n_repetition=ncons) 489 IF (explicit) THEN 490 topology%const_colv = .TRUE. 491 DO ig = 1, ncons 492 CALL section_vals_val_get(collective_section, "COLVAR", i_rep_section=ig, i_val=icolvar) 493 CPASSERT(icolvar <= SIZE(colvar_p)) 494 END DO 495 cons_info%nconst_colv = ncons 496 ALLOCATE (cons_info%const_colv_mol(ncons)) 497 ALLOCATE (cons_info%const_colv_molname(ncons)) 498 ALLOCATE (cons_info%const_colv_target(ncons)) 499 ALLOCATE (cons_info%const_colv_target_growth(ncons)) 500 ALLOCATE (cons_info%colvar_set(ncons)) 501 ALLOCATE (cons_info%colv_intermolecular(ncons)) 502 ALLOCATE (cons_info%colv_restraint(ncons)) 503 ALLOCATE (cons_info%colv_k0(ncons)) 504 ALLOCATE (cons_info%colv_exclude_qm(ncons)) 505 ALLOCATE (cons_info%colv_exclude_mm(ncons)) 506 DO ig = 1, ncons 507 CALL check_restraint(collective_section, & 508 is_restraint=cons_info%colv_restraint(ig), & 509 k0=cons_info%colv_k0(ig), & 510 i_rep_section=ig, & 511 label="COLLECTIVE") 512 cons_info%const_colv_mol(ig) = 0 513 cons_info%const_colv_molname(ig) = "UNDEF" 514 ! Exclude QM or MM 515 CALL section_vals_val_get(collective_section, "EXCLUDE_QM", i_rep_section=ig, & 516 l_val=cons_info%colv_exclude_qm(ig)) 517 CALL section_vals_val_get(collective_section, "EXCLUDE_MM", i_rep_section=ig, & 518 l_val=cons_info%colv_exclude_mm(ig)) 519 ! Intramolecular restraint 520 CALL section_vals_val_get(collective_section, "INTERMOLECULAR", i_rep_section=ig, & 521 l_val=cons_info%colv_intermolecular(ig)) 522 ! If it is intramolecular let's unset (in case user did it) 523 ! the molecule and molname field 524 IF (cons_info%colv_intermolecular(ig)) THEN 525 CALL section_vals_val_unset(collective_section, "MOLECULE", i_rep_section=ig) 526 CALL section_vals_val_unset(collective_section, "MOLNAME", i_rep_section=ig) 527 END IF 528 ! Let's tag to which molecule we want to apply constraints 529 CALL section_vals_val_get(collective_section, "MOLECULE", i_rep_section=ig, & 530 n_rep_val=nrep) 531 IF (nrep /= 0) THEN 532 CALL section_vals_val_get(collective_section, "MOLECULE", i_rep_section=ig, & 533 i_val=cons_info%const_colv_mol(ig)) 534 END IF 535 CALL section_vals_val_get(collective_section, "MOLNAME", i_rep_section=ig, & 536 n_rep_val=nrep) 537 IF (nrep /= 0) THEN 538 CALL section_vals_val_get(collective_section, "MOLNAME", i_rep_section=ig, & 539 c_val=cons_info%const_colv_molname(ig)) 540 END IF 541 IF (((cons_info%const_colv_mol(ig) /= 0) .AND. (cons_info%const_colv_molname(ig) /= "UNDEF"))) THEN 542 CPABORT("Both MOLNAME and MOLECULE specified for CONSTRAINT section. ") 543 END IF 544 IF ((cons_info%const_colv_mol(ig) == 0) .AND. (cons_info%const_colv_molname(ig) == "UNDEF") .AND. & 545 (.NOT. cons_info%colv_intermolecular(ig))) THEN 546 CALL cp_abort(__LOCATION__, & 547 "Constraint section error: you have to specify at least one of the "// & 548 "following keywords: MOLECULE, MOLNAME or INTERMOLECULAR! ") 549 END IF 550 NULLIFY (cons_info%colvar_set(ig)%colvar) 551 CALL section_vals_val_get(collective_section, "COLVAR", i_rep_section=ig, & 552 i_val=icolvar) 553 CALL colvar_clone(cons_info%colvar_set(ig)%colvar, & 554 colvar_p(icolvar)%colvar) 555 CALL section_vals_val_get(collective_section, "TARGET", & 556 n_rep_val=n_rep, i_rep_section=ig) 557 IF (n_rep /= 0) THEN 558 CALL section_vals_val_get(collective_section, "TARGET", & 559 r_val=cons_info%const_colv_target(ig), i_rep_section=ig) 560 ELSE 561 cons_info%const_colv_target(ig) = -HUGE(0.0_dp) 562 END IF 563 CALL section_vals_val_get(collective_section, "TARGET_GROWTH", & 564 r_val=cons_info%const_colv_target_growth(ig), i_rep_section=ig) 565 END DO 566 END IF 567 END IF 568 569 END SUBROUTINE read_constraints_section 570 571! ************************************************************************************************** 572!> \brief Reads input and decides if apply restraints instead of constraints 573!> \param cons_section ... 574!> \param is_restraint ... 575!> \param k0 ... 576!> \param i_rep_section ... 577!> \param label ... 578!> \author teo 579! ************************************************************************************************** 580 SUBROUTINE check_restraint(cons_section, is_restraint, k0, i_rep_section, label) 581 TYPE(section_vals_type), POINTER :: cons_section 582 LOGICAL, INTENT(OUT) :: is_restraint 583 REAL(KIND=dp), INTENT(OUT) :: k0 584 INTEGER, INTENT(IN), OPTIONAL :: i_rep_section 585 CHARACTER(LEN=*), INTENT(IN) :: label 586 587 CHARACTER(len=*), PARAMETER :: routineN = 'check_restraint', & 588 routineP = moduleN//':'//routineN 589 590 CHARACTER(LEN=default_string_length) :: nlabel 591 INTEGER :: output_unit 592 LOGICAL :: explicit 593 TYPE(section_vals_type), POINTER :: restraint_section 594 595 is_restraint = .FALSE. 596 output_unit = cp_logger_get_default_io_unit() 597 CALL section_vals_get(cons_section, explicit=explicit) 598 IF (explicit) THEN 599 restraint_section => section_vals_get_subs_vals(cons_section, "RESTRAINT", & 600 i_rep_section=i_rep_section) 601 CALL section_vals_get(restraint_section, explicit=is_restraint) 602 IF (is_restraint) THEN 603 CALL section_vals_val_get(restraint_section, "K", r_val=k0) 604 IF (output_unit > 0) THEN 605 nlabel = cp_to_string(i_rep_section) 606 WRITE (output_unit, FMT='(T2,"RESTRAINT|",1X,A,F9.6)') & 607 "Active restraint on "//label//" section Nr."// & 608 TRIM(nlabel)//". K [a.u.]=", k0 609 END IF 610 END IF 611 END IF 612 END SUBROUTINE check_restraint 613 614END MODULE topology_input 615 616