(git:374b731)
Loading...
Searching...
No Matches
topology_coordinate_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Collection of subroutine needed for topology related things
10!> \par History
11!> jgh (23-05-2004) Last atom of molecule information added
12! **************************************************************************************************
13MODULE topology_coordinate_util
14 USE atomic_kind_types, ONLY: atomic_kind_type,&
15 get_atomic_kind,&
16 set_atomic_kind
17 USE cp_log_handling, ONLY: cp_get_default_logger,&
18 cp_logger_type
19 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
20 cp_print_key_unit_nr
21 USE exclusion_types, ONLY: exclusion_type
22 USE external_potential_types, ONLY: allocate_potential,&
23 fist_potential_type,&
24 get_potential,&
25 set_potential
26 USE input_constants, ONLY: do_fist,&
27 do_skip_12,&
28 do_skip_13,&
29 do_skip_14
30 USE input_section_types, ONLY: section_vals_get,&
31 section_vals_get_subs_vals,&
32 section_vals_type,&
33 section_vals_val_get
34 USE kinds, ONLY: default_string_length,&
35 dp
36 USE memory_utilities, ONLY: reallocate
37 USE molecule_kind_types, ONLY: atom_type,&
38 get_molecule_kind,&
39 molecule_kind_type,&
40 set_molecule_kind
41 USE molecule_types, ONLY: get_molecule,&
42 molecule_type
43 USE particle_types, ONLY: allocate_particle_set,&
44 particle_type
45 USE physcon, ONLY: massunit
46 USE qmmm_types_low, ONLY: qmmm_env_mm_type
47 USE string_table, ONLY: id2str,&
48 s2s,&
49 str2id
50 USE topology_types, ONLY: atom_info_type,&
51 connectivity_info_type,&
52 topology_parameters_type
53 USE topology_util, ONLY: array1_list_type,&
54 reorder_structure
55#include "./base/base_uses.f90"
56
57 IMPLICIT NONE
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_coordinate_util'
60
61 PRIVATE
62 PUBLIC :: topology_coordinate_pack
63
64CONTAINS
65
66! **************************************************************************************************
67!> \brief Take info readin from different file format and stuff it into
68!> compatible data structure in cp2k
69!> \param particle_set ...
70!> \param atomic_kind_set ...
71!> \param molecule_kind_set ...
72!> \param molecule_set ...
73!> \param topology ...
74!> \param qmmm ...
75!> \param qmmm_env ...
76!> \param subsys_section ...
77!> \param force_env_section ...
78!> \param exclusions ...
79!> \param ignore_outside_box ...
80!> \par History
81!> Teodoro Laino - modified in order to optimize the list of molecules
82!> to build the exclusion lists
83! **************************************************************************************************
84 SUBROUTINE topology_coordinate_pack(particle_set, atomic_kind_set, &
85 molecule_kind_set, molecule_set, topology, qmmm, qmmm_env, &
86 subsys_section, force_env_section, exclusions, ignore_outside_box)
87 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
88 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
89 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
90 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
91 TYPE(topology_parameters_type), INTENT(INOUT) :: topology
92 LOGICAL, INTENT(IN), OPTIONAL :: qmmm
93 TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env
94 TYPE(section_vals_type), POINTER :: subsys_section, force_env_section
95 TYPE(exclusion_type), DIMENSION(:), OPTIONAL, &
96 POINTER :: exclusions
97 LOGICAL, INTENT(IN), OPTIONAL :: ignore_outside_box
98
99 CHARACTER(len=*), PARAMETER :: routinen = 'topology_coordinate_pack'
100
101 CHARACTER(LEN=default_string_length) :: atmname
102 INTEGER :: atom_i, atom_j, counter, dim0, dim1, &
103 dim2, dim3, first, handle, handle2, i, &
104 iatom, ikind, iw, j, k, last, &
105 method_name_id, n, natom
106 INTEGER, DIMENSION(:), POINTER :: iatomlist, id_element, id_work, kind_of, &
107 list, list2, molecule_list, &
108 natom_of_kind, wlist
109 INTEGER, DIMENSION(:, :), POINTER :: pairs
110 LOGICAL :: autogen, check, disable_exclusion_lists, do_center, explicit, found, &
111 my_ignore_outside_box, my_qmmm, present_12_excl_ei_list, present_12_excl_vdw_list
112 REAL(kind=dp) :: bounds(2, 3), cdims(3), dims(3), qeff, &
113 vec(3)
114 REAL(kind=dp), DIMENSION(:), POINTER :: charge, cpoint, mass
115 TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bend_list, ex_bond_list, &
116 ex_bond_list_ei, ex_bond_list_vdw, &
117 ex_onfo_list
118 TYPE(atom_info_type), POINTER :: atom_info
119 TYPE(atom_type), DIMENSION(:), POINTER :: atom_list
120 TYPE(atomic_kind_type), POINTER :: atomic_kind
121 TYPE(connectivity_info_type), POINTER :: conn_info
122 TYPE(cp_logger_type), POINTER :: logger
123 TYPE(fist_potential_type), POINTER :: fist_potential
124 TYPE(molecule_kind_type), POINTER :: molecule_kind
125 TYPE(molecule_type), POINTER :: molecule
126 TYPE(section_vals_type), POINTER :: exclude_section, topology_section
127
128 NULLIFY (logger)
129 logger => cp_get_default_logger()
130 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
131 extension=".subsysLog")
132 topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
133 CALL timeset(routinen, handle)
134
135 my_qmmm = .false.
136 IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
137 atom_info => topology%atom_info
138 conn_info => topology%conn_info
139 !-----------------------------------------------------------------------------
140 !-----------------------------------------------------------------------------
141 ! 1. Determine topology%[natom_type,atom_names] and save mass(natom_type)
142 ! and element(natom_type)
143 !-----------------------------------------------------------------------------
144 CALL timeset(routinen//'_1', handle2)
145 counter = 0
146 NULLIFY (id_work, mass, id_element, charge)
147 ALLOCATE (id_work(topology%natoms))
148 ALLOCATE (mass(topology%natoms))
149 ALLOCATE (id_element(topology%natoms))
150 ALLOCATE (charge(topology%natoms))
151 id_work = str2id(s2s(""))
152 IF (iw > 0) WRITE (iw, *) "molecule_kind_set ::", SIZE(molecule_kind_set)
153 DO i = 1, SIZE(molecule_kind_set)
154 j = molecule_kind_set(i)%molecule_list(1)
155 molecule => molecule_set(j)
156 molecule_kind => molecule_set(j)%molecule_kind
157 IF (iw > 0) WRITE (iw, *) "molecule number ::", j, " has molecule kind number ::", i
158 CALL get_molecule_kind(molecule_kind=molecule_kind, &
159 natom=natom, atom_list=atom_list)
160 CALL get_molecule(molecule=molecule, &
161 first_atom=first, last_atom=last)
162 IF (iw > 0) WRITE (iw, *) "boundaries of molecules (first, last) ::", first, last
163 DO j = 1, natom
164 IF (.NOT. any(id_work(1:counter) .EQ. atom_list(j)%id_name)) THEN
165 counter = counter + 1
166 id_work(counter) = atom_list(j)%id_name
167 mass(counter) = atom_info%atm_mass(first + j - 1)
168 id_element(counter) = atom_info%id_element(first + j - 1)
169 charge(counter) = atom_info%atm_charge(first + j - 1)
170 IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
171 "NEW ATOMIC KIND", id2str(id_work(counter)), mass(counter), id2str(id_element(counter)), charge(counter)
172 ELSE
173 found = .false.
174 DO k = 1, counter
175 IF ((id_work(k) == atom_list(j)%id_name) .AND. (charge(k) == atom_info%atm_charge(first + j - 1))) THEN
176 found = .true.
177 EXIT
178 END IF
179 END DO
180 IF (.NOT. found) THEN
181 counter = counter + 1
182 id_work(counter) = atom_list(j)%id_name
183 mass(counter) = atom_info%atm_mass(first + j - 1)
184 id_element(counter) = atom_info%id_element(first + j - 1)
185 charge(counter) = atom_info%atm_charge(first + j - 1)
186 IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
187 "NEW ATOMIC KIND", id2str(id_work(counter)), mass(counter), id2str(id_element(counter)), charge(counter)
188 END IF
189 END IF
190 END DO
191 END DO
192 topology%natom_type = counter
193 ALLOCATE (atom_info%id_atom_names(topology%natom_type))
194 DO k = 1, counter
195 atom_info%id_atom_names(k) = id_work(k)
196 END DO
197 DEALLOCATE (id_work)
198 CALL reallocate(mass, 1, counter)
199 CALL reallocate(id_element, 1, counter)
200 CALL reallocate(charge, 1, counter)
201 IF (iw > 0) &
202 WRITE (iw, '(5X,A,I3)') "Total Number of Atomic Kinds = ", topology%natom_type
203 CALL timestop(handle2)
204
205 !-----------------------------------------------------------------------------
206 !-----------------------------------------------------------------------------
207 ! 2. Allocate the data structure for the atomic kind information
208 !-----------------------------------------------------------------------------
209 CALL timeset(routinen//'_2', handle2)
210 NULLIFY (atomic_kind_set)
211 ALLOCATE (atomic_kind_set(topology%natom_type))
212 CALL timestop(handle2)
213
214 !-----------------------------------------------------------------------------
215 !-----------------------------------------------------------------------------
216 ! 3. Allocate the data structure for the atomic information
217 !-----------------------------------------------------------------------------
218 CALL timeset(routinen//'_3', handle2)
219 NULLIFY (particle_set)
220 CALL allocate_particle_set(particle_set, topology%natoms)
221 CALL timestop(handle2)
222
223 !-----------------------------------------------------------------------------
224 !-----------------------------------------------------------------------------
225 ! 4. Set the atomic_kind_set(ikind)%[name,kind_number,mass]
226 !-----------------------------------------------------------------------------
227 CALL timeset(routinen//'_4', handle2)
228 DO i = 1, topology%natom_type
229 atomic_kind => atomic_kind_set(i)
230 mass(i) = mass(i)*massunit
231 CALL set_atomic_kind(atomic_kind=atomic_kind, &
232 kind_number=i, &
233 name=id2str(atom_info%id_atom_names(i)), &
234 element_symbol=id2str(id_element(i)), &
235 mass=mass(i))
236 IF (iw > 0) THEN
237 WRITE (iw, '(A,I5,A,I5,4A)') "Atomic Kind n.:", i, " out of:", topology%natom_type, &
238 " name: ", trim(id2str(atom_info%id_atom_names(i))), " element: ", &
239 trim(id2str(id_element(i)))
240 END IF
241 END DO
242 DEALLOCATE (mass)
243 DEALLOCATE (id_element)
244 CALL timestop(handle2)
245
246 !-----------------------------------------------------------------------------
247 !-----------------------------------------------------------------------------
248 ! 5. Determine number of atom of each kind (ie natom_of_kind and kind_of)
249 !-----------------------------------------------------------------------------
250 CALL timeset(routinen//'_5', handle2)
251 ALLOCATE (kind_of(topology%natoms))
252 ALLOCATE (natom_of_kind(topology%natom_type))
253 kind_of(:) = 0
254 natom_of_kind(:) = 0
255 DO i = 1, topology%natom_type
256 DO j = 1, topology%natoms
257 IF ((atom_info%id_atom_names(i) == atom_info%id_atmname(j)) .AND. (charge(i) == atom_info%atm_charge(j))) THEN
258 natom_of_kind(i) = natom_of_kind(i) + 1
259 IF (kind_of(j) == 0) kind_of(j) = i
260 END IF
261 END DO
262 END DO
263 IF (any(kind_of == 0)) THEN
264 DO i = 1, topology%natoms
265 IF (kind_of(i) == 0) THEN
266 WRITE (*, *) i, kind_of(i)
267 WRITE (*, *) "Two molecules have been defined as identical molecules but atoms mismatch charges!!"
268 END IF
269 END DO
270 cpabort("")
271 END IF
272 CALL timestop(handle2)
273
274 !-----------------------------------------------------------------------------
275 !-----------------------------------------------------------------------------
276 ! 6. Set the atom_kind_set(ikind)%[natom,atom_list]
277 !-----------------------------------------------------------------------------
278 CALL timeset(routinen//'_6', handle2)
279 DO i = 1, topology%natom_type
280 atomic_kind => atomic_kind_set(i)
281 NULLIFY (iatomlist)
282 ALLOCATE (iatomlist(natom_of_kind(i)))
283 counter = 0
284 DO j = 1, topology%natoms
285 IF (kind_of(j) == i) THEN
286 counter = counter + 1
287 iatomlist(counter) = j
288 END IF
289 END DO
290 IF (iw > 0) THEN
291 WRITE (iw, '(A,I6,A)') " Atomic kind ", i, " contains particles"
292 DO j = 1, SIZE(iatomlist)
293 IF (mod(j, 5) .EQ. 0) THEN ! split long lines
294 WRITE (iw, '(I12)') iatomlist(j)
295 ELSE
296 WRITE (iw, '(I12)', advance="NO") iatomlist(j)
297 END IF
298 END DO
299 WRITE (iw, *)
300 END IF
301 CALL set_atomic_kind(atomic_kind=atomic_kind, &
302 natom=natom_of_kind(i), &
303 atom_list=iatomlist)
304 DEALLOCATE (iatomlist)
305 END DO
306 DEALLOCATE (natom_of_kind)
307 CALL timestop(handle2)
308
309 !-----------------------------------------------------------------------------
310 !-----------------------------------------------------------------------------
311 ! 7. Possibly center the coordinates and fill in coordinates in particle_set
312 !-----------------------------------------------------------------------------
313 CALL section_vals_val_get(subsys_section, &
314 "TOPOLOGY%CENTER_COORDINATES%_SECTION_PARAMETERS_", l_val=do_center)
315 CALL timeset(routinen//'_7a', handle2)
316 bounds(1, 1) = minval(atom_info%r(1, :))
317 bounds(2, 1) = maxval(atom_info%r(1, :))
318
319 bounds(1, 2) = minval(atom_info%r(2, :))
320 bounds(2, 2) = maxval(atom_info%r(2, :))
321
322 bounds(1, 3) = minval(atom_info%r(3, :))
323 bounds(2, 3) = maxval(atom_info%r(3, :))
324
325 dims = bounds(2, :) - bounds(1, :)
326 cdims(1) = topology%cell%hmat(1, 1)
327 cdims(2) = topology%cell%hmat(2, 2)
328 cdims(3) = topology%cell%hmat(3, 3)
329 IF (iw > 0) THEN
330 WRITE (iw, '(A,3F12.6)') "System sizes: ", dims, "Cell sizes (diagonal): ", cdims
331 END IF
332 check = .true.
333 DO i = 1, 3
334 IF (topology%cell%perd(i) == 0) THEN
335 check = check .AND. (dims(i) < cdims(i))
336 END IF
337 END DO
338 my_ignore_outside_box = .false.
339 IF (PRESENT(ignore_outside_box)) my_ignore_outside_box = ignore_outside_box
340 IF (.NOT. my_ignore_outside_box .AND. .NOT. check) &
341 CALL cp_abort(__location__, &
342 "A non-periodic calculation has been requested but the system size "// &
343 "exceeds the cell size in at least one of the non-periodic directions!")
344 IF (do_center) THEN
345 CALL section_vals_val_get(subsys_section, &
346 "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", explicit=explicit)
347 IF (explicit) THEN
348 CALL section_vals_val_get(subsys_section, &
349 "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", r_vals=cpoint)
350 vec = cpoint
351 ELSE
352 vec = cdims/2.0_dp
353 END IF
354 dims = (bounds(2, :) + bounds(1, :))/2.0_dp - vec
355 ELSE
356 dims = 0.0_dp
357 END IF
358 CALL timestop(handle2)
359 CALL timeset(routinen//'_7b', handle2)
360 DO i = 1, topology%natoms
361 ikind = kind_of(i)
362 IF (iw > 0) THEN
363 WRITE (iw, *) "atom number :: ", i, "kind number ::", ikind
364 END IF
365 particle_set(i)%atomic_kind => atomic_kind_set(ikind)
366 particle_set(i)%r(:) = atom_info%r(:, i) - dims
367 particle_set(i)%atom_index = i
368 END DO
369 CALL timestop(handle2)
370 DEALLOCATE (kind_of)
371
372 !-----------------------------------------------------------------------------
373 !-----------------------------------------------------------------------------
374 ! 8. Fill in the exclusions%list_exclude_vdw
375 ! 9. Fill in the exclusions%list_exclude_ei
376 ! 10. Fill in the exclusions%list_onfo
377 !-----------------------------------------------------------------------------
378 CALL timeset(routinen//'_89', handle2)
379 CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
380 CALL section_vals_val_get(subsys_section, "TOPOLOGY%DISABLE_EXCLUSION_LISTS", &
381 l_val=disable_exclusion_lists)
382 IF ((method_name_id == do_fist) .AND. (.NOT. disable_exclusion_lists)) THEN
383 cpassert(PRESENT(exclusions))
384 natom = topology%natoms
385 ! allocate exclusions. Most likely they would only be needed for the local_particles
386 ALLOCATE (exclusions(natom))
387 DO i = 1, natom
388 NULLIFY (exclusions(i)%list_exclude_vdw)
389 NULLIFY (exclusions(i)%list_exclude_ei)
390 NULLIFY (exclusions(i)%list_onfo)
391 END DO
392 ! Reorder bonds
393 ALLOCATE (ex_bond_list(natom))
394 DO i = 1, natom
395 ALLOCATE (ex_bond_list(i)%array1(0))
396 END DO
397 n = 0
398 IF (ASSOCIATED(conn_info%bond_a)) THEN
399 n = SIZE(conn_info%bond_a)
400 CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, n)
401 END IF
402
403 ! Check if a list of 1-2 exclusion bonds is defined.. if not use all bonds
404 NULLIFY (ex_bond_list_vdw, ex_bond_list_ei)
405 ! VdW
406 exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_VDW_LIST")
407 CALL section_vals_get(exclude_section, explicit=explicit)
408 present_12_excl_vdw_list = .false.
409 IF (explicit) present_12_excl_vdw_list = .true.
410 IF (present_12_excl_vdw_list) THEN
411 ALLOCATE (ex_bond_list_vdw(natom))
412 DO i = 1, natom
413 ALLOCATE (ex_bond_list_vdw(i)%array1(0))
414 END DO
415 CALL setup_exclusion_list(exclude_section, "BOND", ex_bond_list, ex_bond_list_vdw, &
416 particle_set)
417 ELSE
418 ex_bond_list_vdw => ex_bond_list
419 END IF
420 ! EI
421 exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_EI_LIST")
422 CALL section_vals_get(exclude_section, explicit=explicit)
423 present_12_excl_ei_list = .false.
424 IF (explicit) present_12_excl_ei_list = .true.
425 IF (present_12_excl_ei_list) THEN
426 ALLOCATE (ex_bond_list_ei(natom))
427 DO i = 1, natom
428 ALLOCATE (ex_bond_list_ei(i)%array1(0))
429 END DO
430 CALL setup_exclusion_list(exclude_section, "BOND", ex_bond_list, ex_bond_list_ei, &
431 particle_set)
432 ELSE
433 ex_bond_list_ei => ex_bond_list
434 END IF
435
436 CALL section_vals_val_get(topology_section, "AUTOGEN_EXCLUDE_LISTS", &
437 l_val=autogen)
438 ! Reorder bends
439 ALLOCATE (ex_bend_list(natom))
440 DO i = 1, natom
441 ALLOCATE (ex_bend_list(i)%array1(0))
442 END DO
443 IF (autogen) THEN
444 ! Construct autogenerated 1-3 pairs, i.e. all possible 1-3 pairs instead
445 ! of only the bends that are present in the topology.
446 ALLOCATE (pairs(0, 2))
447 n = 0
448 DO iatom = 1, natom
449 DO i = 1, SIZE(ex_bond_list(iatom)%array1)
450 ! a neighboring atom of iatom:
451 atom_i = ex_bond_list(iatom)%array1(i)
452 DO j = 1, i - 1
453 ! another neighboring atom of iatom
454 atom_j = ex_bond_list(iatom)%array1(j)
455 ! It is only a true bend if there is no shorter path.
456 ! No need to check if i and j correspond to the same atom.
457 ! Check if i and j are not involved in a bond:
458 check = .false.
459 DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
460 IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
461 check = .true.
462 EXIT
463 END IF
464 END DO
465 IF (check) cycle
466 ! Add the genuine 1-3 pair
467 n = n + 1
468 IF (SIZE(pairs, dim=1) <= n) THEN
469 CALL reallocate(pairs, 1, n + 5, 1, 2)
470 END IF
471 pairs(n, 1) = atom_i
472 pairs(n, 2) = atom_j
473 END DO
474 END DO
475 END DO
476 CALL reorder_structure(ex_bend_list, pairs(:, 1), pairs(:, 2), n)
477 DEALLOCATE (pairs)
478 ELSE
479 IF (ASSOCIATED(conn_info%theta_a)) THEN
480 n = SIZE(conn_info%theta_a)
481 CALL reorder_structure(ex_bend_list, conn_info%theta_a, conn_info%theta_c, n)
482 END IF
483 END IF
484
485 ! Reorder onfo
486 ALLOCATE (ex_onfo_list(natom))
487 DO i = 1, natom
488 ALLOCATE (ex_onfo_list(i)%array1(0))
489 END DO
490 IF (autogen) THEN
491 ! Construct autogenerated 1-4 pairs, i.e. all possible 1-4 pairs instead
492 ! of only the onfo's that are present in the topology.
493 ALLOCATE (pairs(0, 2))
494 n = 0
495 DO iatom = 1, natom
496 DO i = 1, SIZE(ex_bond_list(iatom)%array1)
497 ! a neighboring atom of iatom:
498 atom_i = ex_bond_list(iatom)%array1(i)
499 DO j = 1, SIZE(ex_bend_list(iatom)%array1)
500 ! a next neighboring atom of iatom:
501 atom_j = ex_bend_list(iatom)%array1(j)
502 ! It is only a true onfo if there is no shorter path.
503 ! check if i and j are not the same atom
504 IF (atom_i == atom_j) cycle
505 ! check if i and j are not involved in a bond
506 check = .false.
507 DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
508 IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
509 check = .true.
510 EXIT
511 END IF
512 END DO
513 IF (check) cycle
514 ! check if i and j are not involved in a bend
515 check = .false.
516 DO counter = 1, SIZE(ex_bend_list(atom_i)%array1)
517 IF (ex_bend_list(atom_i)%array1(counter) == atom_j) THEN
518 check = .true.
519 EXIT
520 END IF
521 END DO
522 IF (check) cycle
523 ! Add the true onfo.
524 n = n + 1
525 IF (SIZE(pairs, dim=1) <= n) THEN
526 CALL reallocate(pairs, 1, n + 5, 1, 2)
527 END IF
528 pairs(n, 1) = atom_i
529 pairs(n, 2) = atom_j
530 END DO
531 END DO
532 END DO
533 CALL reorder_structure(ex_onfo_list, pairs(:, 1), pairs(:, 2), n)
534 DEALLOCATE (pairs)
535 ELSE
536 IF (ASSOCIATED(conn_info%onfo_a)) THEN
537 n = SIZE(conn_info%onfo_a)
538 CALL reorder_structure(ex_onfo_list, conn_info%onfo_a, conn_info%onfo_b, n)
539 END IF
540 END IF
541
542 ! Build the exclusion (and onfo) list per atom.
543 DO iatom = 1, SIZE(particle_set)
544 ! Setup exclusion list for VDW: always exclude itself
545 dim0 = 1
546 ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
547 dim1 = 0
548 IF (topology%exclude_vdw == do_skip_12 .OR. &
549 topology%exclude_vdw == do_skip_13 .OR. &
550 topology%exclude_vdw == do_skip_14) dim1 = SIZE(ex_bond_list_vdw(iatom)%array1)
551 dim1 = dim1 + dim0
552 dim2 = 0
553 IF (topology%exclude_vdw == do_skip_13 .OR. &
554 topology%exclude_vdw == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
555 dim2 = dim1 + dim2
556 dim3 = 0
557 IF (topology%exclude_vdw == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
558 dim3 = dim2 + dim3
559 IF (dim3 /= 0) THEN
560 NULLIFY (list, wlist)
561 ALLOCATE (wlist(dim3))
562 wlist(dim0:dim0) = iatom
563 IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_vdw(iatom)%array1
564 IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
565 IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
566 ! Get a unique list
567 DO i = 1, SIZE(wlist) - 1
568 IF (wlist(i) == 0) cycle
569 DO j = i + 1, SIZE(wlist)
570 IF (wlist(j) == wlist(i)) wlist(j) = 0
571 END DO
572 END DO
573 dim3 = SIZE(wlist) - count(wlist == 0)
574 ALLOCATE (list(dim3))
575 j = 0
576 DO i = 1, SIZE(wlist)
577 IF (wlist(i) == 0) cycle
578 j = j + 1
579 list(j) = wlist(i)
580 END DO
581 DEALLOCATE (wlist)
582 ! Unique list completed
583 NULLIFY (list2)
584 IF ((topology%exclude_vdw == topology%exclude_ei) .AND. &
585 (.NOT. present_12_excl_ei_list) .AND. (.NOT. present_12_excl_vdw_list)) THEN
586 list2 => list
587 ELSE
588 ! Setup exclusion list for EI : always exclude itself
589 dim0 = 1
590 ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
591 dim1 = 0
592 IF (topology%exclude_ei == do_skip_12 .OR. &
593 topology%exclude_ei == do_skip_13 .OR. &
594 topology%exclude_ei == do_skip_14) dim1 = SIZE(ex_bond_list_ei(iatom)%array1)
595 dim1 = dim1 + dim0
596 dim2 = 0
597 IF (topology%exclude_ei == do_skip_13 .OR. &
598 topology%exclude_ei == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
599 dim2 = dim1 + dim2
600 dim3 = 0
601 IF (topology%exclude_ei == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
602 dim3 = dim2 + dim3
603
604 IF (dim3 /= 0) THEN
605 ALLOCATE (wlist(dim3))
606 wlist(dim0:dim0) = iatom
607 IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_ei(iatom)%array1
608 IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
609 IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
610 ! Get a unique list
611 DO i = 1, SIZE(wlist) - 1
612 IF (wlist(i) == 0) cycle
613 DO j = i + 1, SIZE(wlist)
614 IF (wlist(j) == wlist(i)) wlist(j) = 0
615 END DO
616 END DO
617 dim3 = SIZE(wlist) - count(wlist == 0)
618 ALLOCATE (list2(dim3))
619 j = 0
620 DO i = 1, SIZE(wlist)
621 IF (wlist(i) == 0) cycle
622 j = j + 1
623 list2(j) = wlist(i)
624 END DO
625 DEALLOCATE (wlist)
626 ! Unique list completed
627 END IF
628 END IF
629 END IF
630 exclusions(iatom)%list_exclude_vdw => list
631 exclusions(iatom)%list_exclude_ei => list2
632 ! Keep a list of onfo atoms for proper selection of specialized 1-4
633 ! potentials instead of conventional nonbonding potentials.
634 ALLOCATE (exclusions(iatom)%list_onfo(SIZE(ex_onfo_list(iatom)%array1)))
635 ! copy of data, not copy of pointer
636 exclusions(iatom)%list_onfo = ex_onfo_list(iatom)%array1
637 IF (iw > 0) THEN
638 IF (ASSOCIATED(list)) &
639 WRITE (iw, *) "exclusion list_vdw :: ", &
640 "atom num :", iatom, "exclusion list ::", &
641 list
642 IF (topology%exclude_vdw /= topology%exclude_ei) THEN
643 IF (ASSOCIATED(list2)) &
644 WRITE (iw, *) "exclusion list_ei :: ", &
645 "atom num :", iatom, "exclusion list ::", &
646 list2
647 END IF
648 IF (ASSOCIATED(exclusions(iatom)%list_onfo)) &
649 WRITE (iw, *) "onfo list :: ", &
650 "atom num :", iatom, "onfo list ::", &
651 exclusions(iatom)%list_onfo
652 END IF
653 END DO
654 ! deallocate onfo
655 DO i = 1, natom
656 DEALLOCATE (ex_onfo_list(i)%array1)
657 END DO
658 DEALLOCATE (ex_onfo_list)
659 ! deallocate bends
660 DO i = 1, natom
661 DEALLOCATE (ex_bend_list(i)%array1)
662 END DO
663 DEALLOCATE (ex_bend_list)
664 ! deallocate bonds
665 IF (present_12_excl_ei_list) THEN
666 DO i = 1, natom
667 DEALLOCATE (ex_bond_list_ei(i)%array1)
668 END DO
669 DEALLOCATE (ex_bond_list_ei)
670 ELSE
671 NULLIFY (ex_bond_list_ei)
672 END IF
673 IF (present_12_excl_vdw_list) THEN
674 DO i = 1, natom
675 DEALLOCATE (ex_bond_list_vdw(i)%array1)
676 END DO
677 DEALLOCATE (ex_bond_list_vdw)
678 ELSE
679 NULLIFY (ex_bond_list_vdw)
680 END IF
681 DO i = 1, natom
682 DEALLOCATE (ex_bond_list(i)%array1)
683 END DO
684 DEALLOCATE (ex_bond_list)
685 END IF
686 CALL timestop(handle2)
687 !-----------------------------------------------------------------------------
688 !-----------------------------------------------------------------------------
689 ! 11. Set the atomic_kind_set()%fist_potential%[qeff] (PART 1)
690 !-----------------------------------------------------------------------------
691 CALL timeset(routinen//'_10', handle2)
692 CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
693 IF (method_name_id == do_fist) THEN
694 DO i = 1, SIZE(atomic_kind_set)
695 atomic_kind => atomic_kind_set(i)
696 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
697 qeff = charge(i)
698 NULLIFY (fist_potential)
699 CALL allocate_potential(fist_potential)
700 CALL set_potential(potential=fist_potential, qeff=qeff)
701 CALL set_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
702 END DO
703 END IF
704 DEALLOCATE (charge)
705 CALL timestop(handle2)
706
707 !-----------------------------------------------------------------------------
708 !-----------------------------------------------------------------------------
709 ! 12. Set the atom_list for molecule_kind in molecule_kind_set (PART 2)
710 !-----------------------------------------------------------------------------
711 CALL timeset(routinen//'_11', handle2)
712 DO i = 1, SIZE(molecule_kind_set)
713 molecule_kind => molecule_kind_set(i)
714 CALL get_molecule_kind(molecule_kind=molecule_kind, &
715 natom=natom, molecule_list=molecule_list, &
716 atom_list=atom_list)
717 molecule => molecule_set(molecule_list(1))
718 CALL get_molecule(molecule=molecule, &
719 first_atom=first, last_atom=last)
720 DO j = 1, natom
721 DO k = 1, SIZE(atomic_kind_set)
722 atomic_kind => atomic_kind_set(k)
723 CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
724 IF (method_name_id == do_fist) THEN
725 CALL get_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
726 CALL get_potential(potential=fist_potential, qeff=qeff)
727 IF ((id2str(atom_list(j)%id_name) == atmname) .AND. (qeff == atom_info%atm_charge(first + j - 1))) THEN
728 atom_list(j)%atomic_kind => atomic_kind_set(k)
729 EXIT
730 END IF
731 ELSE
732 IF (id2str(atom_list(j)%id_name) == atmname) THEN
733 atom_list(j)%atomic_kind => atomic_kind_set(k)
734 EXIT
735 END IF
736 END IF
737 END DO
738 END DO
739 CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list)
740 END DO
741 CALL timestop(handle2)
742
743 CALL timestop(handle)
744 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
745 "PRINT%TOPOLOGY_INFO/UTIL_INFO")
746 END SUBROUTINE topology_coordinate_pack
747
748! **************************************************************************************************
749!> \brief Builds the exclusion list for VDW and EI if an explicit list of terms
750!> is provided by the user. Otherwise all possibilities are excluded
751!> \param exclude_section ...
752!> \param keyword ...
753!> \param ex_bond_list ...
754!> \param ex_bond_list_w ...
755!> \param particle_set ...
756!> \par History
757!> Teodoro Laino [tlaino] - 12.2009
758! **************************************************************************************************
759 SUBROUTINE setup_exclusion_list(exclude_section, keyword, ex_bond_list, &
760 ex_bond_list_w, particle_set)
761 TYPE(section_vals_type), POINTER :: exclude_section
762 CHARACTER(LEN=*), INTENT(IN) :: keyword
763 TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bond_list, ex_bond_list_w
764 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
765
766 CHARACTER(LEN=default_string_length) :: flag1, flag2
767 CHARACTER(LEN=default_string_length), &
768 DIMENSION(:), POINTER :: names
769 INTEGER :: i, ind, j, k, l, m, n_rep
770
771 cpassert(ASSOCIATED(ex_bond_list))
772 cpassert(ASSOCIATED(ex_bond_list_w))
773 SELECT CASE (keyword)
774 CASE ("BOND")
775 CALL section_vals_val_get(exclude_section, keyword, n_rep_val=n_rep)
776 DO j = 1, SIZE(ex_bond_list)
777 cpassert(ASSOCIATED(ex_bond_list(j)%array1))
778 cpassert(ASSOCIATED(ex_bond_list_w(j)%array1))
779
780 flag1 = particle_set(j)%atomic_kind%name
781 m = SIZE(ex_bond_list(j)%array1)
782 CALL reallocate(ex_bond_list_w(j)%array1, 1, m)
783
784 l = 0
785 DO k = 1, m
786 ind = ex_bond_list(j)%array1(k)
787 flag2 = particle_set(ind)%atomic_kind%name
788 DO i = 1, n_rep
789 CALL section_vals_val_get(exclude_section, keyword, i_rep_val=i, &
790 c_vals=names)
791 IF (((trim(names(1)) == trim(flag1)) .AND. (trim(names(2)) == trim(flag2))) .OR. &
792 ((trim(names(1)) == trim(flag2)) .AND. (trim(names(2)) == trim(flag1)))) THEN
793 l = l + 1
794 ex_bond_list_w(j)%array1(l) = ind
795 END IF
796 END DO
797 END DO
798 CALL reallocate(ex_bond_list_w(j)%array1, 1, l)
799 END DO
800 CASE DEFAULT
801 cpabort("")
802 END SELECT
803
804 END SUBROUTINE setup_exclusion_list
805
806END MODULE topology_coordinate_util
external_potential_types::allocate_potential
Definition external_potential_types.F:260
external_potential_types::get_potential
Definition external_potential_types.F:276
external_potential_types::set_potential
Definition external_potential_types.F:296
memory_utilities::reallocate
Definition memory_utilities.F:27
topology_util::reorder_structure
Definition topology_util.F:73
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::set_atomic_kind
subroutine, public set_atomic_kind(atomic_kind, element_symbol, name, mass, kind_number, natom, atom_list, fist_potential, shell, shell_active, damping)
Set the components of an atomic kind data set.
Definition atomic_kind_types.F:327
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1013
exclusion_types
an exclusion type
Definition exclusion_types.F:10
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::do_skip_13
integer, parameter, public do_skip_13
Definition input_constants.F:412
input_constants::do_fist
integer, parameter, public do_fist
Definition input_constants.F:48
input_constants::do_skip_12
integer, parameter, public do_skip_12
Definition input_constants.F:412
input_constants::do_skip_14
integer, parameter, public do_skip_14
Definition input_constants.F:412
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:697
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1040
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:510
molecule_kind_types::set_molecule_kind
subroutine, public set_molecule_kind(molecule_kind, name, mass, charge, kind_number, molecule_list, atom_list, nbond, bond_list, nbend, bend_list, nub, ub_list, nimpr, impr_list, nopbend, opbend_list, ntorsion, torsion_list, fixd_list, ncolv, colv_list, ng3x3, g3x3_list, ng4x6, nfixd, g4x6_list, nvsite, vsite_list, ng3x3_restraint, ng4x6_restraint, nfixd_restraint, nshell, shell_list, nvsite_restraint, bond_kind_set, bend_kind_set, ub_kind_set, torsion_kind_set, impr_kind_set, opbend_kind_set, nelectron, nsgf, molname_generated)
Set the components of a molecule kind.
Definition molecule_kind_types.F:779
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:277
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
particle_types::allocate_particle_set
subroutine, public allocate_particle_set(particle_set, nparticle)
Allocate a particle set.
Definition particle_types.F:69
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::massunit
real(kind=dp), parameter, public massunit
Definition physcon.F:141
qmmm_types_low
Definition qmmm_types_low.F:13
string_table
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition string_table.F:22
string_table::s2s
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition string_table.F:141
string_table::str2id
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
Definition string_table.F:72
string_table::id2str
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition string_table.F:115
topology_coordinate_util
Collection of subroutine needed for topology related things.
Definition topology_coordinate_util.F:13
topology_coordinate_util::topology_coordinate_pack
subroutine, public topology_coordinate_pack(particle_set, atomic_kind_set, molecule_kind_set, molecule_set, topology, qmmm, qmmm_env, subsys_section, force_env_section, exclusions, ignore_outside_box)
Take info readin from different file format and stuff it into compatible data structure in cp2k.
Definition topology_coordinate_util.F:87
topology_types
Definition topology_types.F:18
topology_util
Collection of subroutine needed for topology related things.
Definition topology_util.F:13
topology
Control for reading in different topologies and coordinates.
Definition topology.F:13
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
exclusion_types::exclusion_type
A type used to store lists of exclusions and onfos.
Definition exclusion_types.F:21
external_potential_types::fist_potential_type
Definition external_potential_types.F:85
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
molecule_kind_types::atom_type
Definition molecule_kind_types.F:56
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:161
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
qmmm_types_low::qmmm_env_mm_type
...
Definition qmmm_types_low.F:238
topology_types::atom_info_type
Definition topology_types.F:35
topology_types::connectivity_info_type
Definition topology_types.F:53
topology_types::topology_parameters_type
Definition topology_types.F:142
topology_util::array1_list_type
Definition topology_util.F:50