1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Collection of subroutine needed for topology related things 8!> \par History 9!> jgh (23-05-2004) Last atom of molecule information added 10! ************************************************************************************************** 11MODULE topology_constraint_util 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind,& 14 is_hydrogen 15 USE cell_types, ONLY: use_perd_x,& 16 use_perd_xy,& 17 use_perd_xyz,& 18 use_perd_xz,& 19 use_perd_y,& 20 use_perd_yz,& 21 use_perd_z 22 USE colvar_methods, ONLY: colvar_eval_mol_f 23 USE colvar_types, ONLY: & 24 colvar_clone, colvar_counters, colvar_create, colvar_p_reallocate, colvar_release, & 25 colvar_setup, colvar_type, dist_colvar_id, torsion_colvar_id, xyz_diag_colvar_id, & 26 xyz_outerdiag_colvar_id 27 USE colvar_utils, ONLY: post_process_colvar 28 USE cp_log_handling, ONLY: cp_get_default_logger,& 29 cp_logger_type,& 30 cp_to_string 31 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 32 cp_print_key_unit_nr 33 USE input_constants, ONLY: do_constr_atomic,& 34 do_constr_molec 35 USE input_section_types, ONLY: section_vals_get,& 36 section_vals_get_subs_vals,& 37 section_vals_type,& 38 section_vals_val_get,& 39 section_vals_val_set 40 USE kinds, ONLY: default_string_length,& 41 dp 42 USE memory_utilities, ONLY: reallocate 43 USE molecule_kind_types, ONLY: & 44 atom_type, bond_type, colvar_constraint_type, fixd_constraint_type, g3x3_constraint_type, & 45 g4x6_constraint_type, get_molecule_kind, molecule_kind_type, set_molecule_kind, & 46 setup_colvar_counters, vsite_constraint_type 47 USE molecule_types, ONLY: get_molecule,& 48 global_constraint_type,& 49 local_colvar_constraint_type,& 50 local_constraint_type,& 51 local_g3x3_constraint_type,& 52 local_g4x6_constraint_type,& 53 molecule_type,& 54 set_molecule 55 USE particle_types, ONLY: particle_type 56 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 57 USE qmmm_types_low, ONLY: qmmm_env_mm_type 58 USE topology_types, ONLY: constr_list_type,& 59 constraint_info_type,& 60 topology_parameters_type 61#include "./base/base_uses.f90" 62 63 IMPLICIT NONE 64 65 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_constraint_util' 66 67 PRIVATE 68 PUBLIC :: topology_constraint_pack 69 70CONTAINS 71 72! ************************************************************************************************** 73!> \brief Pack in all the information needed for the constraints 74!> \param molecule_kind_set ... 75!> \param molecule_set ... 76!> \param topology ... 77!> \param qmmm_env ... 78!> \param particle_set ... 79!> \param input_file ... 80!> \param subsys_section ... 81!> \param gci ... 82! ************************************************************************************************** 83 SUBROUTINE topology_constraint_pack(molecule_kind_set, molecule_set, & 84 topology, qmmm_env, particle_set, input_file, subsys_section, gci) 85 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 86 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 87 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 88 TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env 89 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 90 TYPE(section_vals_type), POINTER :: input_file, subsys_section 91 TYPE(global_constraint_type), POINTER :: gci 92 93 CHARACTER(len=*), PARAMETER :: routineN = 'topology_constraint_pack', & 94 routineP = moduleN//':'//routineN 95 96 CHARACTER(LEN=2) :: element_symbol 97 CHARACTER(LEN=default_string_length) :: molname, name 98 CHARACTER(LEN=default_string_length), & 99 DIMENSION(:), POINTER :: atom_typeh, cnds 100 INTEGER :: cind, first, first_atom, gind, handle, handle2, i, ii, itype, iw, j, k, k1loc, & 101 k2loc, kk, last, last_atom, m, n_start_colv, natom, nbond, ncolv_glob, ncolv_mol, & 102 nfixd_list_gci, nfixd_restart, nfixd_restraint, nfixed_atoms, ng3x3, ng3x3_restraint, & 103 ng4x6, ng4x6_restraint, nhdist, nmolecule, nrep, nvsite, nvsite_restraint, offset 104 INTEGER, DIMENSION(:), POINTER :: constr_x_glob, inds, molecule_list 105 LOGICAL :: exclude_mm, exclude_qm, fix_atom_mm, fix_atom_molname, fix_atom_qm, & 106 fix_atom_qmmm, fix_fixed_atom, found_molname, is_qm, ishbond, ldummy, & 107 restart_restraint_clv, restart_restraint_pos, use_clv_info 108 LOGICAL, ALLOCATABLE, DIMENSION(:) :: missed_molname 109 REAL(KIND=dp) :: rmod, rvec(3) 110 REAL(KIND=dp), DIMENSION(:), POINTER :: hdist, r 111 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list 112 TYPE(atomic_kind_type), POINTER :: atomic_kind 113 TYPE(bond_type), DIMENSION(:), POINTER :: bond_list 114 TYPE(colvar_constraint_type), DIMENSION(:), & 115 POINTER :: colv_list 116 TYPE(colvar_counters) :: ncolv 117 TYPE(constr_list_type), DIMENSION(:), POINTER :: constr_x_mol 118 TYPE(constraint_info_type), POINTER :: cons_info 119 TYPE(cp_logger_type), POINTER :: logger 120 TYPE(fixd_constraint_type), DIMENSION(:), POINTER :: fixd_list, fixd_list_gci 121 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list 122 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list 123 TYPE(local_colvar_constraint_type), DIMENSION(:), & 124 POINTER :: lcolv 125 TYPE(local_constraint_type), POINTER :: lci 126 TYPE(local_g3x3_constraint_type), DIMENSION(:), & 127 POINTER :: lg3x3 128 TYPE(local_g4x6_constraint_type), DIMENSION(:), & 129 POINTER :: lg4x6 130 TYPE(molecule_kind_type), POINTER :: molecule_kind 131 TYPE(molecule_type), POINTER :: molecule 132 TYPE(section_vals_type), POINTER :: colvar_func_info, colvar_rest, & 133 fixd_restr_rest, hbonds_section 134 TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list 135 136 NULLIFY (logger, constr_x_mol, constr_x_glob) 137 logger => cp_get_default_logger() 138 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", & 139 extension=".subsysLog") 140 CALL timeset(routineN, handle) 141 CALL timeset(routineN//"_1", handle2) 142 143 cons_info => topology%cons_info 144 hbonds_section => section_vals_get_subs_vals(input_file, & 145 "MOTION%CONSTRAINT%HBONDS") 146 fixd_restr_rest => section_vals_get_subs_vals(input_file, & 147 "MOTION%CONSTRAINT%FIX_ATOM_RESTART") 148 CALL section_vals_get(fixd_restr_rest, explicit=restart_restraint_pos) 149 colvar_rest => section_vals_get_subs_vals(input_file, & 150 "MOTION%CONSTRAINT%COLVAR_RESTART") 151 CALL section_vals_get(colvar_rest, explicit=restart_restraint_clv) 152 colvar_func_info => section_vals_get_subs_vals(subsys_section, & 153 "COLVAR%COLVAR_FUNC_INFO") 154 CALL section_vals_get(colvar_func_info, explicit=use_clv_info) 155 !----------------------------------------------------------------------------- 156 !----------------------------------------------------------------------------- 157 ! 1. NULLIFY the molecule_set(imol)%lci via set_molecule_set 158 !----------------------------------------------------------------------------- 159 DO i = 1, topology%nmol 160 molecule => molecule_set(i) 161 NULLIFY (lci) 162 ! only allocate the lci if constraints are active. Can this stuff be distributed ? 163 IF (topology%const_atom .OR. topology%const_hydr .OR. & 164 topology%const_33 .OR. topology%const_46 .OR. & 165 topology%const_colv .OR. topology%const_vsite) THEN 166 ALLOCATE (lci) 167 NULLIFY (lci%lcolv) 168 NULLIFY (lci%lg3x3) 169 NULLIFY (lci%lg4x6) 170 ENDIF 171 CALL set_molecule(molecule, lci=lci) 172 END DO 173 ALLOCATE (gci) 174 NULLIFY (gci%lcolv, & 175 gci%lg3x3, & 176 gci%lg4x6, & 177 gci%fixd_list, & 178 gci%colv_list, & 179 gci%g3x3_list, & 180 gci%g4x6_list, & 181 gci%vsite_list) 182 gci%ntot = 0 183 gci%ng3x3 = 0 184 gci%ng4x6 = 0 185 gci%nvsite = 0 186 gci%ng3x3_restraint = 0 187 gci%ng4x6_restraint = 0 188 gci%nvsite_restraint = 0 189 CALL setup_colvar_counters(gci%colv_list, gci%ncolv) 190 gci%nrestraint = gci%ng3x3_restraint + & 191 gci%ng4x6_restraint + & 192 gci%nvsite_restraint + & 193 gci%ncolv%nrestraint 194 CALL timestop(handle2) 195 CALL timeset(routineN//"_2", handle2) 196 !----------------------------------------------------------------------------- 197 !----------------------------------------------------------------------------- 198 ! 2. Add more stuff to COLVAR constraint if constraint hydrogen is on 199 !----------------------------------------------------------------------------- 200 IF (topology%const_hydr) THEN 201 topology%const_colv = .TRUE. 202 NULLIFY (atom_typeh, hdist) 203 ALLOCATE (constr_x_mol(SIZE(molecule_kind_set))) 204 DO i = 1, SIZE(molecule_kind_set) 205 ALLOCATE (constr_x_mol(i)%constr(1)) 206 constr_x_mol(i)%constr(1) = 1 207 END DO 208 CALL section_vals_val_get(hbonds_section, "MOLECULE", n_rep_val=nrep) 209 IF (nrep /= 0) THEN 210 NULLIFY (inds) 211 DO i = 1, SIZE(molecule_kind_set) 212 constr_x_mol(i)%constr(1) = 0 213 END DO 214 CALL section_vals_val_get(hbonds_section, "MOLECULE", i_vals=inds) 215 DO i = 1, SIZE(inds) 216 constr_x_mol(inds(i))%constr(1) = 1 217 END DO 218 ELSE 219 CALL section_vals_val_get(hbonds_section, "MOLNAME", n_rep_val=nrep) 220 IF (nrep /= 0) THEN 221 NULLIFY (cnds) 222 DO i = 1, SIZE(molecule_kind_set) 223 constr_x_mol(i)%constr(1) = 0 224 END DO 225 CALL section_vals_val_get(hbonds_section, "MOLNAME", c_vals=cnds) 226 DO i = 1, SIZE(cnds) 227 found_molname = .FALSE. 228 DO k = 1, SIZE(molecule_kind_set) 229 molecule_kind => molecule_kind_set(k) 230 name = molecule_kind%name 231 ldummy = qmmm_ff_precond_only_qm(id1=name) 232 IF (cnds(i) == name) THEN 233 constr_x_mol(k)%constr(1) = 1 234 found_molname = .TRUE. 235 END IF 236 END DO 237 CALL print_warning_molname(found_molname, cnds(i)) 238 END DO 239 END IF 240 END IF 241 CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", n_rep_val=nrep) 242 IF (nrep /= 0) & 243 CALL section_vals_val_get(hbonds_section, "ATOM_TYPE", c_vals=atom_typeh) 244 CALL section_vals_val_get(hbonds_section, "TARGETS", n_rep_val=nrep) 245 IF (nrep /= 0) & 246 CALL section_vals_val_get(hbonds_section, "TARGETS", r_vals=hdist) 247 IF (ASSOCIATED(hdist)) THEN 248 CPASSERT(SIZE(hdist) == SIZE(atom_typeh)) 249 END IF 250 CALL section_vals_val_get(hbonds_section, "exclude_qm", l_val=exclude_qm) 251 CALL section_vals_val_get(hbonds_section, "exclude_mm", l_val=exclude_mm) 252 nhdist = 0 253 DO i = 1, SIZE(molecule_kind_set) 254 molecule_kind => molecule_kind_set(i) 255 IF (constr_x_mol(i)%constr(1) == 0) CYCLE 256 CALL get_molecule_kind(molecule_kind=molecule_kind, & 257 bond_list=bond_list, nbond=nbond, atom_list=atom_list, & 258 molecule_list=molecule_list) 259 ! Let's tag all requested atoms involving Hydrogen 260 ! on the first molecule of this kind 261 molecule => molecule_set(molecule_list(1)) 262 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 263 natom = last_atom - first_atom + 1 264 DO k = 1, nbond 265 ishbond = .FALSE. 266 j = bond_list(k)%a 267 IF (j < 1 .OR. j > natom) CYCLE 268 atomic_kind => atom_list(j)%atomic_kind 269 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name) 270 is_qm = qmmm_ff_precond_only_qm(id1=name) 271 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE. 272 IF (is_qm .AND. exclude_qm) ishbond = .FALSE. 273 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE. 274 IF (.NOT. ishbond) THEN 275 j = bond_list(k)%b 276 IF (j < 1 .OR. j > natom) CYCLE 277 atomic_kind => atom_list(j)%atomic_kind 278 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name) 279 is_qm = qmmm_ff_precond_only_qm(id1=name) 280 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE. 281 IF (is_qm .AND. exclude_qm) ishbond = .FALSE. 282 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE. 283 END IF 284 IF (ishbond) THEN 285 nhdist = nhdist + 1 286 END IF 287 END DO 288 END DO 289 n_start_colv = cons_info%nconst_colv 290 cons_info%nconst_colv = nhdist + n_start_colv 291 CALL reallocate(cons_info%const_colv_mol, 1, cons_info%nconst_colv) 292 CALL reallocate(cons_info%const_colv_molname, 1, cons_info%nconst_colv) 293 CALL reallocate(cons_info%const_colv_target, 1, cons_info%nconst_colv) 294 CALL reallocate(cons_info%const_colv_target_growth, 1, cons_info%nconst_colv) 295 CALL colvar_p_reallocate(cons_info%colvar_set, 1, cons_info%nconst_colv) 296 ! Fill in Restraints info 297 CALL reallocate(cons_info%colv_intermolecular, 1, cons_info%nconst_colv) 298 CALL reallocate(cons_info%colv_restraint, 1, cons_info%nconst_colv) 299 CALL reallocate(cons_info%colv_k0, 1, cons_info%nconst_colv) 300 CALL reallocate(cons_info%colv_exclude_qm, 1, cons_info%nconst_colv) 301 CALL reallocate(cons_info%colv_exclude_mm, 1, cons_info%nconst_colv) 302 ! Bonds involving hydrogens are by their nature only intramolecular 303 cons_info%colv_intermolecular(n_start_colv + 1:cons_info%nconst_colv) = .FALSE. 304 cons_info%colv_exclude_qm(n_start_colv + 1:cons_info%nconst_colv) = .FALSE. 305 cons_info%colv_exclude_mm(n_start_colv + 1:cons_info%nconst_colv) = .FALSE. 306 cons_info%colv_restraint(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_restraint 307 cons_info%colv_k0(n_start_colv + 1:cons_info%nconst_colv) = cons_info%hbonds_k0 308 ! 309 nhdist = 0 310 DO i = 1, SIZE(molecule_kind_set) 311 IF (constr_x_mol(i)%constr(1) == 0) CYCLE 312 molecule_kind => molecule_kind_set(i) 313 CALL get_molecule_kind(molecule_kind=molecule_kind, & 314 bond_list=bond_list, nbond=nbond, atom_list=atom_list, & 315 molecule_list=molecule_list) 316 molecule => molecule_set(molecule_list(1)) 317 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 318 natom = last_atom - first_atom + 1 319 offset = first_atom - 1 320 DO k = 1, nbond 321 ishbond = .FALSE. 322 j = bond_list(k)%a 323 IF (j < 1 .OR. j > natom) CYCLE 324 atomic_kind => atom_list(j)%atomic_kind 325 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name) 326 is_qm = qmmm_ff_precond_only_qm(id1=name) 327 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE. 328 IF (is_qm .AND. exclude_qm) ishbond = .FALSE. 329 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE. 330 IF (.NOT. ishbond) THEN 331 j = bond_list(k)%b 332 IF (j < 1 .OR. j > natom) CYCLE 333 atomic_kind => atom_list(j)%atomic_kind 334 CALL get_atomic_kind(atomic_kind=atomic_kind, name=name) 335 is_qm = qmmm_ff_precond_only_qm(id1=name) 336 IF ((name(1:1) == "H") .OR. is_hydrogen(atomic_kind)) ishbond = .TRUE. 337 IF (is_qm .AND. exclude_qm) ishbond = .FALSE. 338 IF (.NOT. (is_qm) .AND. exclude_mm) ishbond = .FALSE. 339 END IF 340 IF (ishbond) THEN 341 nhdist = nhdist + 1 342 rvec = particle_set(offset + bond_list(k)%a)%r - particle_set(offset + bond_list(k)%b)%r 343 rmod = SQRT(DOT_PRODUCT(rvec, rvec)) 344 IF (ASSOCIATED(hdist)) THEN 345 IF (SIZE(hdist) > 0) THEN 346 IF (bond_list(k)%a == j) atomic_kind => atom_list(bond_list(k)%b)%atomic_kind 347 IF (bond_list(k)%b == j) atomic_kind => atom_list(bond_list(k)%a)%atomic_kind 348 CALL get_atomic_kind(atomic_kind=atomic_kind, & 349 name=name, element_symbol=element_symbol) 350 ldummy = qmmm_ff_precond_only_qm(id1=name) 351 DO m = 1, SIZE(hdist) 352 IF (TRIM(name) == TRIM(atom_typeh(m))) EXIT 353 IF (TRIM(element_symbol) == TRIM(atom_typeh(m))) EXIT 354 END DO 355 IF (m <= SIZE(hdist)) THEN 356 rmod = hdist(m) 357 END IF 358 END IF 359 END IF 360 cons_info%const_colv_mol(nhdist + n_start_colv) = i 361 cons_info%const_colv_molname(nhdist + n_start_colv) = "UNDEF" 362 cons_info%const_colv_target(nhdist + n_start_colv) = rmod 363 cons_info%const_colv_target_growth(nhdist + n_start_colv) = 0.0_dp 364 CALL colvar_create(cons_info%colvar_set(nhdist + n_start_colv)%colvar, & 365 dist_colvar_id) 366 cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%i_at = bond_list(k)%a 367 cons_info%colvar_set(nhdist + n_start_colv)%colvar%dist_param%j_at = bond_list(k)%b 368 CALL colvar_setup(cons_info%colvar_set(nhdist + n_start_colv)%colvar) 369 END IF 370 END DO 371 END DO 372 DO j = 1, SIZE(constr_x_mol) 373 DEALLOCATE (constr_x_mol(j)%constr) 374 END DO 375 DEALLOCATE (constr_x_mol) 376 END IF 377 378 CALL timestop(handle2) 379 CALL timeset(routineN//"_3", handle2) 380 !----------------------------------------------------------------------------- 381 !----------------------------------------------------------------------------- 382 ! 3. Set the COLVAR constraint molecule_kind_set(ikind)%colv_list 383 !----------------------------------------------------------------------------- 384 IF (topology%const_colv) THEN 385 ! Post Process of COLVARS.. 386 DO ii = 1, SIZE(cons_info%colvar_set) 387 CALL post_process_colvar(cons_info%colvar_set(ii)%colvar, particle_set) 388 END DO 389 ! Real constraint/restraint part.. 390 CALL give_constraint_array(cons_info%const_colv_mol, & 391 cons_info%const_colv_molname, & 392 cons_info%colv_intermolecular, & 393 constr_x_mol, & 394 constr_x_glob, & 395 molecule_kind_set, & 396 cons_info%colv_exclude_qm, & 397 cons_info%colv_exclude_mm) 398 ! Intramolecular constraints 399 gind = 0 400 cind = 0 401 DO ii = 1, SIZE(molecule_kind_set) 402 molecule_kind => molecule_kind_set(ii) 403 CALL get_molecule_kind(molecule_kind=molecule_kind, & 404 nmolecule=nmolecule, molecule_list=molecule_list) 405 ncolv_mol = SIZE(constr_x_mol(ii)%constr) 406 ALLOCATE (colv_list(ncolv_mol)) 407 ! Starting index of the first molecule of this kind. 408 ! We need the index if no target is provided in the input file 409 ! for the collective variable.. The target will be computed on the 410 ! first molecule of the kind... 411 molecule => molecule_set(molecule_list(1)) 412 CALL get_molecule(molecule, first_atom=first_atom) 413 CALL setup_colv_list(colv_list, constr_x_mol(ii)%constr, gind, & 414 cons_info, topology, particle_set, restart_restraint_clv, & 415 colvar_rest, first_atom) 416 CALL setup_colvar_counters(colv_list, ncolv) 417 CALL set_molecule_kind(molecule_kind, colv_list=colv_list, ncolv=ncolv) 418 DO j = 1, nmolecule 419 molecule => molecule_set(molecule_list(j)) 420 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 421 ALLOCATE (lcolv(ncolv_mol)) 422 CALL setup_lcolv(lcolv, constr_x_mol(ii)%constr, first_atom, last_atom, & 423 cons_info, particle_set, colvar_func_info, use_clv_info, cind) 424 CALL set_molecule(molecule=molecule, lcolv=lcolv) 425 END DO 426 END DO 427 DO j = 1, SIZE(constr_x_mol) 428 DEALLOCATE (constr_x_mol(j)%constr) 429 END DO 430 DEALLOCATE (constr_x_mol) 431 ! Intermolecular constraints 432 ncolv_glob = 0 433 IF (ASSOCIATED(constr_x_glob)) THEN 434 ncolv_glob = SIZE(constr_x_glob) 435 ALLOCATE (colv_list(ncolv_glob)) 436 CALL setup_colv_list(colv_list, constr_x_glob, gind, cons_info, & 437 topology, particle_set, restart_restraint_clv, colvar_rest, & 438 first_atom=1) 439 CALL setup_colvar_counters(colv_list, ncolv) 440 ALLOCATE (lcolv(ncolv_glob)) 441 CALL setup_lcolv(lcolv, constr_x_glob, 1, SIZE(particle_set), cons_info, & 442 particle_set, colvar_func_info, use_clv_info, cind) 443 gci%colv_list => colv_list 444 gci%lcolv => lcolv 445 gci%ncolv = ncolv 446 ! Total number of Intermolecular constraints 447 gci%ntot = gci%ncolv%ntot + gci%ntot 448 DEALLOCATE (constr_x_glob) 449 END IF 450 END IF 451 452 CALL timestop(handle2) 453 CALL timeset(routineN//"_4", handle2) 454 !----------------------------------------------------------------------------- 455 !----------------------------------------------------------------------------- 456 ! 4. Set the group 3x3 constraint g3x3_list 457 !----------------------------------------------------------------------------- 458 IF (topology%const_33) THEN 459 CALL give_constraint_array(cons_info%const_g33_mol, & 460 cons_info%const_g33_molname, & 461 cons_info%g33_intermolecular, & 462 constr_x_mol, & 463 constr_x_glob, & 464 molecule_kind_set, & 465 cons_info%g33_exclude_qm, & 466 cons_info%g33_exclude_mm) 467 ! Intramolecular constraints 468 DO ii = 1, SIZE(molecule_kind_set) 469 molecule_kind => molecule_kind_set(ii) 470 CALL get_molecule_kind(molecule_kind=molecule_kind, & 471 nmolecule=nmolecule, & 472 molecule_list=molecule_list) 473 ng3x3 = SIZE(constr_x_mol(ii)%constr) 474 ALLOCATE (g3x3_list(ng3x3)) 475 CALL setup_g3x3_list(g3x3_list, constr_x_mol(ii)%constr, cons_info, ng3x3_restraint) 476 CALL set_molecule_kind(molecule_kind, ng3x3=ng3x3, ng3x3_restraint=ng3x3_restraint, g3x3_list=g3x3_list) 477 DO j = 1, nmolecule 478 molecule => molecule_set(molecule_list(j)) 479 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 480 ALLOCATE (lg3x3(ng3x3)) 481 CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom) 482 CALL set_molecule(molecule=molecule, lg3x3=lg3x3) 483 END DO 484 END DO 485 DO j = 1, SIZE(constr_x_mol) 486 DEALLOCATE (constr_x_mol(j)%constr) 487 END DO 488 DEALLOCATE (constr_x_mol) 489 ! Intermolecular constraints 490 IF (ASSOCIATED(constr_x_glob)) THEN 491 ng3x3 = SIZE(constr_x_glob) 492 ALLOCATE (g3x3_list(ng3x3)) 493 CALL setup_g3x3_list(g3x3_list, constr_x_glob, cons_info, ng3x3_restraint) 494 ALLOCATE (lg3x3(ng3x3)) 495 CALL setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom) 496 gci%g3x3_list => g3x3_list 497 gci%lg3x3 => lg3x3 498 gci%ng3x3 = ng3x3 499 gci%ng3x3_restraint = ng3x3_restraint 500 ! Total number of Intermolecular constraints 501 gci%ntot = 3*gci%ng3x3 + gci%ntot 502 DEALLOCATE (constr_x_glob) 503 END IF 504 END IF 505 506 CALL timestop(handle2) 507 CALL timeset(routineN//"_5", handle2) 508 !----------------------------------------------------------------------------- 509 !----------------------------------------------------------------------------- 510 ! 5. Set the group 4x6 constraint g4x6_list 511 !----------------------------------------------------------------------------- 512 IF (topology%const_46) THEN 513 CALL give_constraint_array(cons_info%const_g46_mol, & 514 cons_info%const_g46_molname, & 515 cons_info%g46_intermolecular, & 516 constr_x_mol, & 517 constr_x_glob, & 518 molecule_kind_set, & 519 cons_info%g46_exclude_qm, & 520 cons_info%g46_exclude_mm) 521 ! Intramolecular constraints 522 DO ii = 1, SIZE(molecule_kind_set) 523 molecule_kind => molecule_kind_set(ii) 524 CALL get_molecule_kind(molecule_kind=molecule_kind, & 525 nmolecule=nmolecule, molecule_list=molecule_list) 526 ng4x6 = SIZE(constr_x_mol(ii)%constr) 527 ALLOCATE (g4x6_list(ng4x6)) 528 CALL setup_g4x6_list(g4x6_list, constr_x_mol(ii)%constr, cons_info, ng4x6_restraint) 529 CALL set_molecule_kind(molecule_kind, ng4x6=ng4x6, ng4x6_restraint=ng4x6_restraint, g4x6_list=g4x6_list) 530 DO j = 1, nmolecule 531 molecule => molecule_set(molecule_list(j)) 532 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 533 ALLOCATE (lg4x6(ng4x6)) 534 CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom) 535 CALL set_molecule(molecule=molecule, lg4x6=lg4x6) 536 END DO 537 END DO 538 DO j = 1, SIZE(constr_x_mol) 539 DEALLOCATE (constr_x_mol(j)%constr) 540 END DO 541 DEALLOCATE (constr_x_mol) 542 ! Intermolecular constraints 543 IF (ASSOCIATED(constr_x_glob)) THEN 544 ng4x6 = SIZE(constr_x_glob) 545 ALLOCATE (g4x6_list(ng4x6)) 546 CALL setup_g4x6_list(g4x6_list, constr_x_glob, cons_info, ng4x6_restraint) 547 ALLOCATE (lg4x6(ng4x6)) 548 CALL setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom) 549 gci%g4x6_list => g4x6_list 550 gci%lg4x6 => lg4x6 551 gci%ng4x6 = ng4x6 552 gci%ng4x6_restraint = ng4x6_restraint 553 ! Total number of Intermolecular constraints 554 gci%ntot = 6*gci%ng4x6 + gci%ntot 555 DEALLOCATE (constr_x_glob) 556 END IF 557 END IF 558 559 CALL timestop(handle2) 560 CALL timeset(routineN//"_6", handle2) 561 !----------------------------------------------------------------------------- 562 !----------------------------------------------------------------------------- 563 ! 6. Set the group vsite constraint vsite_list 564 !----------------------------------------------------------------------------- 565 IF (topology%const_vsite) THEN 566 CALL give_constraint_array(cons_info%const_vsite_mol, & 567 cons_info%const_vsite_molname, & 568 cons_info%vsite_intermolecular, & 569 constr_x_mol, & 570 constr_x_glob, & 571 molecule_kind_set, & 572 cons_info%vsite_exclude_qm, & 573 cons_info%vsite_exclude_mm) 574 ! Intramolecular constraints 575 DO ii = 1, SIZE(molecule_kind_set) 576 molecule_kind => molecule_kind_set(ii) 577 CALL get_molecule_kind(molecule_kind=molecule_kind, & 578 nmolecule=nmolecule, molecule_list=molecule_list) 579 nvsite = SIZE(constr_x_mol(ii)%constr) 580 ALLOCATE (vsite_list(nvsite)) 581 CALL setup_vsite_list(vsite_list, constr_x_mol(ii)%constr, cons_info, nvsite_restraint) 582 CALL set_molecule_kind(molecule_kind, nvsite=nvsite, nvsite_restraint=nvsite_restraint, & 583 vsite_list=vsite_list) 584 END DO 585 DO j = 1, SIZE(constr_x_mol) 586 DEALLOCATE (constr_x_mol(j)%constr) 587 END DO 588 DEALLOCATE (constr_x_mol) 589 ! Intermolecular constraints 590 IF (ASSOCIATED(constr_x_glob)) THEN 591 nvsite = SIZE(constr_x_glob) 592 ALLOCATE (vsite_list(nvsite)) 593 CALL setup_vsite_list(vsite_list, constr_x_glob, cons_info, nvsite_restraint) 594 gci%vsite_list => vsite_list 595 gci%nvsite = nvsite 596 gci%nvsite_restraint = nvsite_restraint 597 ! Total number of Intermolecular constraints 598 gci%ntot = gci%nvsite + gci%ntot 599 DEALLOCATE (constr_x_glob) 600 END IF 601 END IF 602 CALL timestop(handle2) 603 CALL timeset(routineN//"_7", handle2) 604 !----------------------------------------------------------------------------- 605 !----------------------------------------------------------------------------- 606 ! 7. Set the group fixed_atom constraint fixd_list 607 !----------------------------------------------------------------------------- 608 IF (topology%const_atom) THEN 609 ALLOCATE (fixd_list_gci(SIZE(particle_set))) 610 nfixd_list_gci = 0 611 ALLOCATE (missed_molname(SIZE(cons_info%fixed_molnames, 1))) 612 missed_molname = .TRUE. 613 nfixd_restart = 0 614 DO i = 1, SIZE(molecule_kind_set) 615 molecule_kind => molecule_kind_set(i) 616 CALL get_molecule_kind(molecule_kind=molecule_kind, & 617 nmolecule=nmolecule, molecule_list=molecule_list, name=molname) 618 is_qm = qmmm_ff_precond_only_qm(id1=molname) 619 WHERE (molname .EQ. cons_info%fixed_molnames) 620 missed_molname = .FALSE. 621 END WHERE 622 ! Try to figure out how many atoms of the list belong to this molecule_kind 623 nfixed_atoms = 0 624 DO j = 1, nmolecule 625 molecule => molecule_set(molecule_list(j)) 626 CALL get_molecule(molecule, first_atom=first, last_atom=last) 627 fix_atom_molname = .FALSE. 628 IF (ASSOCIATED(cons_info%fixed_molnames)) THEN 629 DO k = 1, SIZE(cons_info%fixed_molnames) 630 IF (cons_info%fixed_molnames(k) .EQ. molname) THEN 631 fix_atom_molname = .TRUE. 632 IF (is_qm .AND. cons_info%fixed_exclude_qm(k)) fix_atom_molname = .FALSE. 633 IF ((.NOT. is_qm) .AND. cons_info%fixed_exclude_mm(k)) fix_atom_molname = .FALSE. 634 END IF 635 END DO 636 ENDIF 637 DO k = first, last 638 fix_atom_qmmm = .FALSE. 639 IF (PRESENT(qmmm_env)) THEN 640 SELECT CASE (cons_info%freeze_qm) 641 CASE (do_constr_atomic) 642 IF (ANY(qmmm_env%qm_atom_index == k)) fix_atom_qmmm = .TRUE. 643 CASE (do_constr_molec) 644 IF (ANY(qmmm_env%qm_molecule_index == molecule_list(j))) fix_atom_qmmm = .TRUE. 645 END SELECT 646 SELECT CASE (cons_info%freeze_mm) 647 CASE (do_constr_atomic) 648 IF (ALL(qmmm_env%qm_atom_index /= k)) fix_atom_qmmm = .TRUE. 649 CASE (do_constr_molec) 650 IF (ALL(qmmm_env%qm_molecule_index /= molecule_list(j))) fix_atom_qmmm = .TRUE. 651 END SELECT 652 END IF 653 IF (ANY(cons_info%fixed_atoms == k) .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN 654 nfixed_atoms = nfixed_atoms + 1 655 END IF 656 END DO 657 END DO 658 ALLOCATE (fixd_list(nfixed_atoms)) 659 kk = 0 660 nfixd_restraint = 0 661 IF (nfixed_atoms /= 0) THEN 662 DO j = 1, nmolecule 663 molecule => molecule_set(molecule_list(j)) 664 CALL get_molecule(molecule, first_atom=first, last_atom=last) 665 fix_atom_molname = .FALSE. 666 IF (ASSOCIATED(cons_info%fixed_molnames)) THEN 667 DO k1loc = 1, SIZE(cons_info%fixed_molnames) 668 IF (cons_info%fixed_molnames(k1loc) .EQ. molname) THEN 669 fix_atom_molname = .TRUE. 670 itype = cons_info%fixed_mol_type(k1loc) 671 EXIT 672 END IF 673 END DO 674 ENDIF 675 DO k = first, last 676 ! FIXED LIST ATOMS 677 fix_fixed_atom = .FALSE. 678 DO k2loc = 1, SIZE(cons_info%fixed_atoms) 679 IF (cons_info%fixed_atoms(k2loc) == k) THEN 680 fix_fixed_atom = .TRUE. 681 itype = cons_info%fixed_type(k2loc) 682 EXIT 683 END IF 684 END DO 685 ! QMMM FIXED ATOMS (QM OR MM) 686 fix_atom_qmmm = .FALSE. 687 fix_atom_mm = .FALSE. 688 fix_atom_qm = .FALSE. 689 IF (PRESENT(qmmm_env)) THEN 690 SELECT CASE (cons_info%freeze_qm) 691 CASE (do_constr_atomic) 692 IF (ANY(qmmm_env%qm_atom_index == k)) THEN 693 fix_atom_qmmm = .TRUE. 694 fix_atom_qm = .TRUE. 695 itype = cons_info%freeze_qm_type 696 END IF 697 CASE (do_constr_molec) 698 IF (ANY(qmmm_env%qm_molecule_index == molecule_list(j))) THEN 699 fix_atom_qmmm = .TRUE. 700 fix_atom_qm = .TRUE. 701 itype = cons_info%freeze_qm_type 702 END IF 703 END SELECT 704 SELECT CASE (cons_info%freeze_mm) 705 CASE (do_constr_atomic) 706 IF (ALL(qmmm_env%qm_atom_index /= k)) THEN 707 fix_atom_qmmm = .TRUE. 708 fix_atom_mm = .TRUE. 709 itype = cons_info%freeze_mm_type 710 END IF 711 CASE (do_constr_molec) 712 IF (ALL(qmmm_env%qm_molecule_index /= molecule_list(j))) THEN 713 fix_atom_qmmm = .TRUE. 714 fix_atom_mm = .TRUE. 715 itype = cons_info%freeze_mm_type 716 END IF 717 END SELECT 718 ! We should never reach this point but let's check it anyway 719 IF (fix_atom_qm .AND. fix_atom_mm) THEN 720 CALL cp_abort(__LOCATION__, & 721 "Atom number: "//cp_to_string(k)// & 722 " has been defined both QM and MM. General Error!") 723 END IF 724 END IF 725 ! Check that the fixed atom constraint/restraint is unique 726 IF ((fix_fixed_atom .AND. fix_atom_qmmm) .OR. (fix_fixed_atom .AND. fix_atom_molname) & 727 .OR. (fix_atom_qmmm .AND. fix_atom_molname)) THEN 728 CALL cp_abort(__LOCATION__, & 729 "Atom number: "//cp_to_string(k)// & 730 " has been constrained/restrained to be fixed in more than one"// & 731 " input section. Check and correct your input file!") 732 END IF 733 ! Let's store the atom index 734 IF (fix_fixed_atom .OR. fix_atom_qmmm .OR. fix_atom_molname) THEN 735 kk = kk + 1 736 fixd_list(kk)%fixd = k 737 fixd_list(kk)%coord = particle_set(k)%r 738 fixd_list(kk)%itype = itype 739 ! Possibly Restraint 740 IF (fix_fixed_atom) THEN 741 fixd_list(kk)%restraint%active = cons_info%fixed_restraint(k2loc) 742 fixd_list(kk)%restraint%k0 = cons_info%fixed_k0(k2loc) 743 ELSEIF (fix_atom_qm) THEN 744 fixd_list(kk)%restraint%active = cons_info%fixed_qm_restraint 745 fixd_list(kk)%restraint%k0 = cons_info%fixed_qm_k0 746 ELSEIF (fix_atom_mm) THEN 747 fixd_list(kk)%restraint%active = cons_info%fixed_mm_restraint 748 fixd_list(kk)%restraint%k0 = cons_info%fixed_mm_k0 749 ELSEIF (fix_atom_molname) THEN 750 fixd_list(kk)%restraint%active = cons_info%fixed_mol_restraint(k1loc) 751 fixd_list(kk)%restraint%k0 = cons_info%fixed_mol_k0(k1loc) 752 ELSE 753 ! Should never reach this point 754 CPABORT("") 755 END IF 756 IF (fixd_list(kk)%restraint%active) THEN 757 nfixd_restraint = nfixd_restraint + 1 758 nfixd_restart = nfixd_restart + 1 759 ! Check that we use the components that we really want.. 760 SELECT CASE (itype) 761 CASE (use_perd_x) 762 fixd_list(kk)%coord(2) = HUGE(0.0_dp) 763 fixd_list(kk)%coord(3) = HUGE(0.0_dp) 764 CASE (use_perd_y) 765 fixd_list(kk)%coord(1) = HUGE(0.0_dp) 766 fixd_list(kk)%coord(3) = HUGE(0.0_dp) 767 CASE (use_perd_z) 768 fixd_list(kk)%coord(1) = HUGE(0.0_dp) 769 fixd_list(kk)%coord(2) = HUGE(0.0_dp) 770 CASE (use_perd_xy) 771 fixd_list(kk)%coord(3) = HUGE(0.0_dp) 772 CASE (use_perd_xz) 773 fixd_list(kk)%coord(2) = HUGE(0.0_dp) 774 CASE (use_perd_yz) 775 fixd_list(kk)%coord(1) = HUGE(0.0_dp) 776 END SELECT 777 IF (restart_restraint_pos) THEN 778 ! Read coord0 value for restraint 779 CALL section_vals_val_get(fixd_restr_rest, "_DEFAULT_KEYWORD_", & 780 i_rep_val=nfixd_restart, r_vals=r) 781 SELECT CASE (itype) 782 CASE (use_perd_x) 783 CPASSERT(SIZE(r) == 1) 784 fixd_list(kk)%coord(1) = r(1) 785 CASE (use_perd_y) 786 CPASSERT(SIZE(r) == 1) 787 fixd_list(kk)%coord(2) = r(1) 788 CASE (use_perd_z) 789 CPASSERT(SIZE(r) == 1) 790 fixd_list(kk)%coord(3) = r(1) 791 CASE (use_perd_xy) 792 CPASSERT(SIZE(r) == 2) 793 fixd_list(kk)%coord(1) = r(1) 794 fixd_list(kk)%coord(2) = r(2) 795 CASE (use_perd_xz) 796 CPASSERT(SIZE(r) == 2) 797 fixd_list(kk)%coord(1) = r(1) 798 fixd_list(kk)%coord(3) = r(2) 799 CASE (use_perd_yz) 800 CPASSERT(SIZE(r) == 2) 801 fixd_list(kk)%coord(2) = r(1) 802 fixd_list(kk)%coord(3) = r(2) 803 CASE (use_perd_xyz) 804 CPASSERT(SIZE(r) == 3) 805 fixd_list(kk)%coord(1:3) = r(1:3) 806 END SELECT 807 ELSE 808 ! Write coord0 value for restraint 809 SELECT CASE (itype) 810 CASE (use_perd_x) 811 ALLOCATE (r(1)) 812 r(1) = fixd_list(kk)%coord(1) 813 CASE (use_perd_y) 814 ALLOCATE (r(1)) 815 r(1) = fixd_list(kk)%coord(2) 816 CASE (use_perd_z) 817 ALLOCATE (r(1)) 818 r(1) = fixd_list(kk)%coord(3) 819 CASE (use_perd_xy) 820 ALLOCATE (r(2)) 821 r(1) = fixd_list(kk)%coord(1) 822 r(2) = fixd_list(kk)%coord(2) 823 CASE (use_perd_xz) 824 ALLOCATE (r(2)) 825 r(1) = fixd_list(kk)%coord(1) 826 r(2) = fixd_list(kk)%coord(3) 827 CASE (use_perd_yz) 828 ALLOCATE (r(2)) 829 r(1) = fixd_list(kk)%coord(1) 830 r(2) = fixd_list(kk)%coord(3) 831 CASE (use_perd_xyz) 832 ALLOCATE (r(3)) 833 r(1:3) = fixd_list(kk)%coord(1:3) 834 END SELECT 835 CALL section_vals_val_set(fixd_restr_rest, "_DEFAULT_KEYWORD_", & 836 i_rep_val=nfixd_restart, r_vals_ptr=r) 837 END IF 838 END IF 839 END IF 840 END DO 841 END DO 842 END IF 843 IF (iw > 0) THEN 844 WRITE (iw, *) "MOLECULE KIND:", i, " NR. FIXED ATOMS:", SIZE(fixd_list(:)%fixd), " LIST::", fixd_list(:)%fixd 845 END IF 846 CALL set_molecule_kind(molecule_kind, nfixd=nfixed_atoms, nfixd_restraint=nfixd_restraint, & 847 fixd_list=fixd_list) 848 fixd_list_gci(nfixd_list_gci + 1:nfixd_list_gci + nfixed_atoms) = fixd_list 849 nfixd_list_gci = nfixd_list_gci + nfixed_atoms 850 END DO 851 IF (iw > 0) THEN 852 WRITE (iw, *) "TOTAL NUMBER OF FIXED ATOMS:", nfixd_list_gci 853 END IF 854 CPASSERT(COUNT(missed_molname) == 0) 855 DEALLOCATE (missed_molname) 856 ! Intermolecular constraints 857 IF (gci%ntot /= 0) THEN 858 ALLOCATE (fixd_list(nfixd_list_gci)) 859 fixd_list(1:nfixd_list_gci) = fixd_list_gci(1:nfixd_list_gci) 860 gci%fixd_list => fixd_list 861 END IF 862 DEALLOCATE (fixd_list_gci) 863 END IF 864 ! Final setup of the number of possible restraints 865 gci%nrestraint = gci%ng3x3_restraint + & 866 gci%ng4x6_restraint + & 867 gci%nvsite_restraint + & 868 gci%ncolv%nrestraint 869 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 870 "PRINT%TOPOLOGY_INFO/UTIL_INFO") 871 CALL timestop(handle2) 872 CALL timestop(handle) 873 END SUBROUTINE topology_constraint_pack 874 875! ************************************************************************************************** 876!> \brief Setup the colv_list for the packing of constraints 877!> \param colv_list ... 878!> \param ilist ... 879!> \param gind ... 880!> \param cons_info ... 881!> \param topology ... 882!> \param particle_set ... 883!> \param restart_restraint_clv ... 884!> \param colvar_rest ... 885!> \param first_atom ... 886!> \par History 887!> Updated 2007 for intermolecular constraints 888!> \author Teodoro Laino [2007] 889! ************************************************************************************************** 890 SUBROUTINE setup_colv_list(colv_list, ilist, gind, cons_info, topology, & 891 particle_set, restart_restraint_clv, colvar_rest, first_atom) 892 893 TYPE(colvar_constraint_type), DIMENSION(:), & 894 POINTER :: colv_list 895 INTEGER, DIMENSION(:), POINTER :: ilist 896 INTEGER, INTENT(INOUT) :: gind 897 TYPE(constraint_info_type), POINTER :: cons_info 898 TYPE(topology_parameters_type), INTENT(INOUT) :: topology 899 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 900 POINTER :: particle_set 901 LOGICAL, INTENT(IN) :: restart_restraint_clv 902 TYPE(section_vals_type), POINTER :: colvar_rest 903 INTEGER, INTENT(IN) :: first_atom 904 905 CHARACTER(len=*), PARAMETER :: routineN = 'setup_colv_list', & 906 routineP = moduleN//':'//routineN 907 908 INTEGER :: j, kdim, kk, ncolv_mol 909 REAL(KIND=dp) :: rmod 910 TYPE(colvar_type), POINTER :: local_colvar 911 912 ncolv_mol = 0 913 DO kk = 1, SIZE(ilist) 914 j = ilist(kk) 915 ncolv_mol = ncolv_mol + 1 916 kdim = SIZE(cons_info%colvar_set(j)%colvar%i_atom) 917 ALLOCATE (colv_list(ncolv_mol)%i_atoms(kdim)) 918 colv_list(ncolv_mol)%inp_seq_num = j 919 colv_list(ncolv_mol)%type_id = cons_info%colvar_set(j)%colvar%type_id 920 colv_list(ncolv_mol)%i_atoms = cons_info%colvar_set(j)%colvar%i_atom 921 colv_list(ncolv_mol)%use_points = cons_info%colvar_set(j)%colvar%use_points 922 ! Restraint 923 colv_list(ncolv_mol)%restraint%active = cons_info%colv_restraint(j) 924 colv_list(ncolv_mol)%restraint%k0 = cons_info%colv_k0(j) 925 IF (cons_info%const_colv_target(j) == -HUGE(0.0_dp)) THEN 926 ! Let's compute the value.. 927 NULLIFY (local_colvar) 928 CALL colvar_clone(local_colvar, cons_info%colvar_set(j)%colvar, & 929 i_atom_offset=first_atom - 1) 930 CALL colvar_eval_mol_f(local_colvar, topology%cell, particle_set) 931 colv_list(ncolv_mol)%expected_value = local_colvar%ss 932 CALL colvar_release(local_colvar) 933 ELSE 934 colv_list(ncolv_mol)%expected_value = cons_info%const_colv_target(j) 935 END IF 936 colv_list(ncolv_mol)%expected_value_growth_speed = cons_info%const_colv_target_growth(j) 937 ! In case of Restraint let's check for possible restart values 938 IF (colv_list(ncolv_mol)%restraint%active .AND. & 939 (colv_list(ncolv_mol)%expected_value_growth_speed == 0.0_dp)) THEN 940 gind = gind + 1 941 IF (restart_restraint_clv) THEN 942 CALL section_vals_val_get(colvar_rest, "_DEFAULT_KEYWORD_", & 943 i_rep_val=gind, r_val=rmod) 944 colv_list(ncolv_mol)%expected_value = rmod 945 ELSE 946 rmod = colv_list(ncolv_mol)%expected_value 947 CALL section_vals_val_set(colvar_rest, "_DEFAULT_KEYWORD_", & 948 i_rep_val=gind, r_val=rmod) 949 END IF 950 END IF 951 ! Only if torsion let's take into account the singularity in the definition 952 ! of the dihedral 953 IF (cons_info%colvar_set(j)%colvar%type_id == torsion_colvar_id) THEN 954 cons_info%colvar_set(j)%colvar%torsion_param%o0 = colv_list(ncolv_mol)%expected_value 955 END IF 956 END DO 957 END SUBROUTINE setup_colv_list 958 959! ************************************************************************************************** 960!> \brief Setup the g3x3_list for the packing of constraints 961!> \param g3x3_list ... 962!> \param ilist ... 963!> \param cons_info ... 964!> \param ng3x3_restraint ... 965!> \par History 966!> Updated 2007 for intermolecular constraints 967!> \author Teodoro Laino [2007] 968! ************************************************************************************************** 969 SUBROUTINE setup_g3x3_list(g3x3_list, ilist, cons_info, ng3x3_restraint) 970 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list 971 INTEGER, DIMENSION(:), POINTER :: ilist 972 TYPE(constraint_info_type), POINTER :: cons_info 973 INTEGER, INTENT(OUT) :: ng3x3_restraint 974 975 CHARACTER(len=*), PARAMETER :: routineN = 'setup_g3x3_list', & 976 routineP = moduleN//':'//routineN 977 978 INTEGER :: j, ng3x3 979 980 ng3x3_restraint = 0 981 DO ng3x3 = 1, SIZE(ilist) 982 j = ilist(ng3x3) 983 g3x3_list(ng3x3)%a = cons_info%const_g33_a(j) 984 g3x3_list(ng3x3)%b = cons_info%const_g33_b(j) 985 g3x3_list(ng3x3)%c = cons_info%const_g33_c(j) 986 g3x3_list(ng3x3)%dab = cons_info%const_g33_dab(j) 987 g3x3_list(ng3x3)%dac = cons_info%const_g33_dac(j) 988 g3x3_list(ng3x3)%dbc = cons_info%const_g33_dbc(j) 989 ! Restraint 990 g3x3_list(ng3x3)%restraint%active = cons_info%g33_restraint(j) 991 g3x3_list(ng3x3)%restraint%k0 = cons_info%g33_k0(j) 992 IF (g3x3_list(ng3x3)%restraint%active) ng3x3_restraint = ng3x3_restraint + 1 993 END DO 994 995 END SUBROUTINE setup_g3x3_list 996 997! ************************************************************************************************** 998!> \brief Setup the g4x6_list for the packing of constraints 999!> \param g4x6_list ... 1000!> \param ilist ... 1001!> \param cons_info ... 1002!> \param ng4x6_restraint ... 1003!> \par History 1004!> Updated 2007 for intermolecular constraints 1005!> \author Teodoro Laino [2007] 1006! ************************************************************************************************** 1007 SUBROUTINE setup_g4x6_list(g4x6_list, ilist, cons_info, ng4x6_restraint) 1008 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list 1009 INTEGER, DIMENSION(:), POINTER :: ilist 1010 TYPE(constraint_info_type), POINTER :: cons_info 1011 INTEGER, INTENT(OUT) :: ng4x6_restraint 1012 1013 CHARACTER(len=*), PARAMETER :: routineN = 'setup_g4x6_list', & 1014 routineP = moduleN//':'//routineN 1015 1016 INTEGER :: j, ng4x6 1017 1018 ng4x6 = 0 1019 ng4x6_restraint = 0 1020 DO ng4x6 = 1, SIZE(ilist) 1021 j = ilist(ng4x6) 1022 g4x6_list(ng4x6)%a = cons_info%const_g46_a(j) 1023 g4x6_list(ng4x6)%b = cons_info%const_g46_b(j) 1024 g4x6_list(ng4x6)%c = cons_info%const_g46_c(j) 1025 g4x6_list(ng4x6)%d = cons_info%const_g46_d(j) 1026 g4x6_list(ng4x6)%dab = cons_info%const_g46_dab(j) 1027 g4x6_list(ng4x6)%dac = cons_info%const_g46_dac(j) 1028 g4x6_list(ng4x6)%dbc = cons_info%const_g46_dbc(j) 1029 g4x6_list(ng4x6)%dad = cons_info%const_g46_dad(j) 1030 g4x6_list(ng4x6)%dbd = cons_info%const_g46_dbd(j) 1031 g4x6_list(ng4x6)%dcd = cons_info%const_g46_dcd(j) 1032 ! Restraint 1033 g4x6_list(ng4x6)%restraint%active = cons_info%g46_restraint(j) 1034 g4x6_list(ng4x6)%restraint%k0 = cons_info%g46_k0(j) 1035 IF (g4x6_list(ng4x6)%restraint%active) ng4x6_restraint = ng4x6_restraint + 1 1036 END DO 1037 1038 END SUBROUTINE setup_g4x6_list 1039 1040! ************************************************************************************************** 1041!> \brief Setup the vsite_list for the packing of constraints 1042!> \param vsite_list ... 1043!> \param ilist ... 1044!> \param cons_info ... 1045!> \param nvsite_restraint ... 1046!> \par History 1047!> \author Marcel Baer [2008] 1048! ************************************************************************************************** 1049 SUBROUTINE setup_vsite_list(vsite_list, ilist, cons_info, nvsite_restraint) 1050 TYPE(vsite_constraint_type), DIMENSION(:), POINTER :: vsite_list 1051 INTEGER, DIMENSION(:), POINTER :: ilist 1052 TYPE(constraint_info_type), POINTER :: cons_info 1053 INTEGER, INTENT(OUT) :: nvsite_restraint 1054 1055 CHARACTER(len=*), PARAMETER :: routineN = 'setup_vsite_list', & 1056 routineP = moduleN//':'//routineN 1057 1058 INTEGER :: j, nvsite 1059 1060 nvsite = 0 1061 nvsite_restraint = 0 1062 DO nvsite = 1, SIZE(ilist) 1063 j = ilist(nvsite) 1064 vsite_list(nvsite)%a = cons_info%const_vsite_a(j) 1065 vsite_list(nvsite)%b = cons_info%const_vsite_b(j) 1066 vsite_list(nvsite)%c = cons_info%const_vsite_c(j) 1067 vsite_list(nvsite)%d = cons_info%const_vsite_d(j) 1068 vsite_list(nvsite)%wbc = cons_info%const_vsite_wbc(j) 1069 vsite_list(nvsite)%wdc = cons_info%const_vsite_wdc(j) 1070 ! Restraint 1071 vsite_list(nvsite)%restraint%active = cons_info%vsite_restraint(j) 1072 vsite_list(nvsite)%restraint%k0 = cons_info%vsite_k0(j) 1073 IF (vsite_list(nvsite)%restraint%active) nvsite_restraint = nvsite_restraint + 1 1074 END DO 1075 1076 END SUBROUTINE setup_vsite_list 1077! ************************************************************************************************** 1078!> \brief Setup the lcolv for the packing of constraints 1079!> \param lcolv ... 1080!> \param ilist ... 1081!> \param first_atom ... 1082!> \param last_atom ... 1083!> \param cons_info ... 1084!> \param particle_set ... 1085!> \param colvar_func_info ... 1086!> \param use_clv_info ... 1087!> \param cind ... 1088!> \par History 1089!> Updated 2007 for intermolecular constraints 1090!> \author Teodoro Laino [2007] 1091! ************************************************************************************************** 1092 SUBROUTINE setup_lcolv(lcolv, ilist, first_atom, last_atom, cons_info, & 1093 particle_set, colvar_func_info, use_clv_info, & 1094 cind) 1095 TYPE(local_colvar_constraint_type), DIMENSION(:), & 1096 POINTER :: lcolv 1097 INTEGER, DIMENSION(:), POINTER :: ilist 1098 INTEGER, INTENT(IN) :: first_atom, last_atom 1099 TYPE(constraint_info_type), POINTER :: cons_info 1100 TYPE(particle_type), DIMENSION(:), OPTIONAL, & 1101 POINTER :: particle_set 1102 TYPE(section_vals_type), POINTER :: colvar_func_info 1103 LOGICAL, INTENT(IN) :: use_clv_info 1104 INTEGER, INTENT(INOUT) :: cind 1105 1106 CHARACTER(len=*), PARAMETER :: routineN = 'setup_lcolv', routineP = moduleN//':'//routineN 1107 1108 INTEGER :: ind, k, kk 1109 REAL(KIND=dp), DIMENSION(:), POINTER :: r_vals 1110 1111 DO kk = 1, SIZE(ilist) 1112 k = ilist(kk) 1113 lcolv(kk)%init = .FALSE. 1114 lcolv(kk)%lambda = 0.0_dp 1115 lcolv(kk)%sigma = 0.0_dp 1116 1117 ! Set Up colvar variable 1118 NULLIFY (lcolv(kk)%colvar, lcolv(kk)%colvar_old) 1119 ! Colvar 1120 CALL colvar_clone(lcolv(kk)%colvar, cons_info%colvar_set(k)%colvar, & 1121 i_atom_offset=first_atom - 1) 1122 1123 ! Some COLVARS may need additional information for evaluating the 1124 ! functional form: this is the case for COLVARS which depend on the 1125 ! initial position of the atoms: This information is stored in a proper 1126 ! container in the COLVAR_RESTART section.. 1127 IF ((lcolv(kk)%colvar%type_id == xyz_diag_colvar_id) .OR. & 1128 (lcolv(kk)%colvar%type_id == xyz_outerdiag_colvar_id)) THEN 1129 cind = cind + 1 1130 IF (use_clv_info) THEN 1131 CALL section_vals_val_get(colvar_func_info, "_DEFAULT_KEYWORD_", & 1132 i_rep_val=cind, r_vals=r_vals) 1133 SELECT CASE (lcolv(kk)%colvar%type_id) 1134 CASE (xyz_diag_colvar_id) 1135 CPASSERT(SIZE(r_vals) == 3) 1136 lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals 1137 CASE (xyz_outerdiag_colvar_id) 1138 CPASSERT(SIZE(r_vals) == 6) 1139 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3) 1140 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 2) = r_vals(4:6) 1141 END SELECT 1142 ELSE 1143 SELECT CASE (lcolv(kk)%colvar%type_id) 1144 CASE (xyz_diag_colvar_id) 1145 ALLOCATE (r_vals(3)) 1146 ind = first_atom - 1 + lcolv(kk)%colvar%xyz_diag_param%i_atom 1147 r_vals = particle_set(ind)%r 1148 lcolv(kk)%colvar%xyz_diag_param%r0 = r_vals 1149 CASE (xyz_outerdiag_colvar_id) 1150 ALLOCATE (r_vals(6)) 1151 ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(1) 1152 r_vals(1:3) = particle_set(ind)%r 1153 ind = first_atom - 1 + lcolv(kk)%colvar%xyz_outerdiag_param%i_atoms(2) 1154 r_vals(4:6) = particle_set(ind)%r 1155 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 1) = r_vals(1:3) 1156 lcolv(kk)%colvar%xyz_outerdiag_param%r0(:, 2) = r_vals(4:6) 1157 END SELECT 1158 CALL section_vals_val_set(colvar_func_info, "_DEFAULT_KEYWORD_", & 1159 i_rep_val=cind, r_vals_ptr=r_vals) 1160 END IF 1161 END IF 1162 1163 ! Setup Colvar_old 1164 CALL colvar_clone(lcolv(kk)%colvar_old, lcolv(kk)%colvar) 1165 1166 ! Check for consistency in the constraint definition 1167 IF (ANY(lcolv(kk)%colvar%i_atom > last_atom) .OR. & 1168 ANY(lcolv(kk)%colvar%i_atom < first_atom)) THEN 1169 WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!" 1170 WRITE (*, '(T5,"|",T8,A)') "A constraint has been defined for a molecule type", & 1171 " but the atoms specified in the constraint and the atoms defined for", & 1172 " the molecule DO NOT match!", & 1173 "This could be very probable due to a wrong connectivity, or an error", & 1174 " in the constraint specification in the input file.", & 1175 " Please check it carefully!" 1176 CPABORT("") 1177 END IF 1178 END DO 1179 END SUBROUTINE setup_lcolv 1180 1181! ************************************************************************************************** 1182!> \brief Setup the lg3x3 for the packing of constraints 1183!> \param lg3x3 ... 1184!> \param g3x3_list ... 1185!> \param first_atom ... 1186!> \param last_atom ... 1187!> \par History 1188!> Updated 2007 for intermolecular constraints 1189!> \author Teodoro Laino [2007] 1190! ************************************************************************************************** 1191 SUBROUTINE setup_lg3x3(lg3x3, g3x3_list, first_atom, last_atom) 1192 TYPE(local_g3x3_constraint_type), DIMENSION(:), & 1193 POINTER :: lg3x3 1194 TYPE(g3x3_constraint_type), DIMENSION(:), POINTER :: g3x3_list 1195 INTEGER, INTENT(IN) :: first_atom, last_atom 1196 1197 CHARACTER(len=*), PARAMETER :: routineN = 'setup_lg3x3', routineP = moduleN//':'//routineN 1198 1199 INTEGER :: kk 1200 1201 DO kk = 1, SIZE(lg3x3) 1202 lg3x3(kk)%init = .FALSE. 1203 lg3x3(kk)%scale = 0.0_dp 1204 lg3x3(kk)%scale_old = 0.0_dp 1205 lg3x3(kk)%fa = 0.0_dp 1206 lg3x3(kk)%fb = 0.0_dp 1207 lg3x3(kk)%fc = 0.0_dp 1208 lg3x3(kk)%ra_old = 0.0_dp 1209 lg3x3(kk)%rb_old = 0.0_dp 1210 lg3x3(kk)%rc_old = 0.0_dp 1211 lg3x3(kk)%va = 0.0_dp 1212 lg3x3(kk)%vb = 0.0_dp 1213 lg3x3(kk)%vc = 0.0_dp 1214 lg3x3(kk)%lambda = 0.0_dp 1215 IF ((g3x3_list(kk)%a + first_atom - 1 < first_atom) .OR. & 1216 (g3x3_list(kk)%b + first_atom - 1 < first_atom) .OR. & 1217 (g3x3_list(kk)%c + first_atom - 1 < first_atom) .OR. & 1218 (g3x3_list(kk)%a + first_atom - 1 > last_atom) .OR. & 1219 (g3x3_list(kk)%b + first_atom - 1 > last_atom) .OR. & 1220 (g3x3_list(kk)%c + first_atom - 1 > last_atom)) THEN 1221 WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!" 1222 WRITE (*, '(T5,"|",T8,A)') "A constraint has been defined for a molecule type", & 1223 " but the atoms specified in the constraint and the atoms defined for", & 1224 " the molecule DO NOT match!", & 1225 "This could be very probable due to a wrong connectivity, or an error", & 1226 " in the constraint specification in the input file.", & 1227 " Please check it carefully!" 1228 CPABORT("") 1229 END IF 1230 END DO 1231 1232 END SUBROUTINE setup_lg3x3 1233 1234! ************************************************************************************************** 1235!> \brief Setup the lg4x6 for the packing of constraints 1236!> \param lg4x6 ... 1237!> \param g4x6_list ... 1238!> \param first_atom ... 1239!> \param last_atom ... 1240!> \par History 1241!> Updated 2007 for intermolecular constraints 1242!> \author Teodoro Laino [2007] 1243! ************************************************************************************************** 1244 SUBROUTINE setup_lg4x6(lg4x6, g4x6_list, first_atom, last_atom) 1245 TYPE(local_g4x6_constraint_type), DIMENSION(:), & 1246 POINTER :: lg4x6 1247 TYPE(g4x6_constraint_type), DIMENSION(:), POINTER :: g4x6_list 1248 INTEGER, INTENT(IN) :: first_atom, last_atom 1249 1250 CHARACTER(len=*), PARAMETER :: routineN = 'setup_lg4x6', routineP = moduleN//':'//routineN 1251 1252 INTEGER :: kk 1253 1254 DO kk = 1, SIZE(lg4x6) 1255 lg4x6(kk)%init = .FALSE. 1256 lg4x6(kk)%scale = 0.0_dp 1257 lg4x6(kk)%scale_old = 0.0_dp 1258 lg4x6(kk)%fa = 0.0_dp 1259 lg4x6(kk)%fb = 0.0_dp 1260 lg4x6(kk)%fc = 0.0_dp 1261 lg4x6(kk)%fd = 0.0_dp 1262 lg4x6(kk)%fe = 0.0_dp 1263 lg4x6(kk)%ff = 0.0_dp 1264 lg4x6(kk)%ra_old = 0.0_dp 1265 lg4x6(kk)%rb_old = 0.0_dp 1266 lg4x6(kk)%rc_old = 0.0_dp 1267 lg4x6(kk)%rd_old = 0.0_dp 1268 lg4x6(kk)%re_old = 0.0_dp 1269 lg4x6(kk)%rf_old = 0.0_dp 1270 lg4x6(kk)%va = 0.0_dp 1271 lg4x6(kk)%vb = 0.0_dp 1272 lg4x6(kk)%vc = 0.0_dp 1273 lg4x6(kk)%vd = 0.0_dp 1274 lg4x6(kk)%ve = 0.0_dp 1275 lg4x6(kk)%vf = 0.0_dp 1276 lg4x6(kk)%lambda = 0.0_dp 1277 IF ((g4x6_list(kk)%a + first_atom - 1 < first_atom) .OR. & 1278 (g4x6_list(kk)%b + first_atom - 1 < first_atom) .OR. & 1279 (g4x6_list(kk)%c + first_atom - 1 < first_atom) .OR. & 1280 (g4x6_list(kk)%d + first_atom - 1 < first_atom) .OR. & 1281 (g4x6_list(kk)%a + first_atom - 1 > last_atom) .OR. & 1282 (g4x6_list(kk)%b + first_atom - 1 > last_atom) .OR. & 1283 (g4x6_list(kk)%c + first_atom - 1 > last_atom) .OR. & 1284 (g4x6_list(kk)%d + first_atom - 1 > last_atom)) THEN 1285 WRITE (*, '(T5,"|",T8,A)') "Error in constraints setup!" 1286 WRITE (*, '(T5,"|",T8,A)') "A constrained has been defined for a molecule type", & 1287 " but the atoms specified in the constraint and the atoms defined for", & 1288 " the molecule DO NOT match!", & 1289 "This could be very probable due to a wrong connectivity, or an error", & 1290 " in the constraint specification in the input file.", & 1291 " Please check it carefully!" 1292 CPABORT("") 1293 END IF 1294 END DO 1295 1296 END SUBROUTINE setup_lg4x6 1297 1298! ************************************************************************************************** 1299!> \brief Gives back a list of molecule to which apply the constraint 1300!> \param const_mol ... 1301!> \param const_molname ... 1302!> \param const_intermolecular ... 1303!> \param constr_x_mol ... 1304!> \param constr_x_glob ... 1305!> \param molecule_kind_set ... 1306!> \param exclude_qm ... 1307!> \param exclude_mm ... 1308!> \par History 1309!> Updated 2007 for intermolecular constraints 1310!> \author Teodoro Laino [2006] 1311! ************************************************************************************************** 1312 SUBROUTINE give_constraint_array(const_mol, const_molname, const_intermolecular, & 1313 constr_x_mol, constr_x_glob, molecule_kind_set, exclude_qm, exclude_mm) 1314 1315 INTEGER, DIMENSION(:), POINTER :: const_mol 1316 CHARACTER(LEN=default_string_length), & 1317 DIMENSION(:), POINTER :: const_molname 1318 LOGICAL, DIMENSION(:), POINTER :: const_intermolecular 1319 TYPE(constr_list_type), DIMENSION(:), POINTER :: constr_x_mol 1320 INTEGER, DIMENSION(:), POINTER :: constr_x_glob 1321 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 1322 LOGICAL, DIMENSION(:), POINTER :: exclude_qm, exclude_mm 1323 1324 CHARACTER(len=*), PARAMETER :: routineN = 'give_constraint_array', & 1325 routineP = moduleN//':'//routineN 1326 1327 CHARACTER(LEN=default_string_length) :: myname, name 1328 INTEGER :: handle, i, iglob, isize, k 1329 LOGICAL :: found_molname, is_qm 1330 TYPE(molecule_kind_type), POINTER :: molecule_kind 1331 1332 CALL timeset(routineN, handle) 1333 NULLIFY (molecule_kind) 1334 ALLOCATE (constr_x_mol(SIZE(molecule_kind_set))) 1335 DO i = 1, SIZE(constr_x_mol) 1336 NULLIFY (constr_x_mol(i)%constr) 1337 ALLOCATE (constr_x_mol(i)%constr(0)) 1338 END DO 1339 CPASSERT(SIZE(const_mol) == SIZE(const_molname)) 1340 iglob = 0 1341 DO i = 1, SIZE(const_mol) 1342 IF (const_intermolecular(i)) THEN 1343 ! Intermolecular constraint 1344 iglob = iglob + 1 1345 CALL reallocate(constr_x_glob, 1, iglob) 1346 constr_x_glob(iglob) = i 1347 ELSE 1348 ! Intramolecular constraint 1349 IF (const_mol(i) /= 0) THEN 1350 k = const_mol(i) 1351 IF (k > SIZE(molecule_kind_set)) & 1352 CALL cp_abort(__LOCATION__, & 1353 "A constraint has been specified providing the molecule index. But the "// & 1354 " molecule index ("//cp_to_string(k)//") is out of range of the possible"// & 1355 " molecule kinds ("//cp_to_string(SIZE(molecule_kind_set))//").") 1356 isize = SIZE(constr_x_mol(k)%constr) 1357 CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1) 1358 constr_x_mol(k)%constr(isize + 1) = i 1359 ELSE 1360 myname = const_molname(i) 1361 found_molname = .FALSE. 1362 DO k = 1, SIZE(molecule_kind_set) 1363 molecule_kind => molecule_kind_set(k) 1364 name = molecule_kind%name 1365 is_qm = qmmm_ff_precond_only_qm(id1=name) 1366 IF (is_qm .AND. exclude_qm(i)) CYCLE 1367 IF (.NOT. is_qm .AND. exclude_mm(i)) CYCLE 1368 IF (name == myname) THEN 1369 isize = SIZE(constr_x_mol(k)%constr) 1370 CALL reallocate(constr_x_mol(k)%constr, 1, isize + 1) 1371 constr_x_mol(k)%constr(isize + 1) = i 1372 found_molname = .TRUE. 1373 END IF 1374 END DO 1375 CALL print_warning_molname(found_molname, myname) 1376 END IF 1377 END IF 1378 END DO 1379 CALL timestop(handle) 1380 END SUBROUTINE give_constraint_array 1381 1382! ************************************************************************************************** 1383!> \brief Prints a warning message if undefined molnames are used to define constraints 1384!> \param found ... 1385!> \param name ... 1386!> \author Teodoro Laino [2007] - Zurich University 1387! ************************************************************************************************** 1388 SUBROUTINE print_warning_molname(found, name) 1389 LOGICAL, INTENT(IN) :: found 1390 CHARACTER(LEN=*), INTENT(IN) :: name 1391 1392 CHARACTER(len=*), PARAMETER :: routineN = 'print_warning_molname', & 1393 routineP = moduleN//':'//routineN 1394 1395 IF (.NOT. found) & 1396 CALL cp_warn(__LOCATION__, & 1397 " MOLNAME ("//TRIM(name)//") was defined for constraints, but this molecule name "// & 1398 "is not defined. Please check carefully your PDB, PSF (has priority over PDB) or "// & 1399 "input driven CP2K coordinates. In case you may not find the reason for this warning "// & 1400 "it may be a good idea to print all molecule information (including kind name) activating "// & 1401 "the print_key MOLECULES specific of the SUBSYS%PRINT section. ") 1402 1403 END SUBROUTINE print_warning_molname 1404 1405END MODULE topology_constraint_util 1406