(git:0de0cc2)
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 ! **************************************************************************************************
14  USE atomic_kind_types, ONLY: atomic_kind_type,&
18  cp_logger_type
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,&
32  section_vals_type,&
34  USE kinds, ONLY: default_string_length,&
35  dp
36  USE memory_utilities, ONLY: reallocate
37  USE molecule_kind_types, ONLY: atom_type,&
39  molecule_kind_type,&
41  USE molecule_types, ONLY: get_molecule,&
42  molecule_type
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,&
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 
64 CONTAINS
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 
806 END MODULE topology_coordinate_util
Define the atomic kind types and their sub types.
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.
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.
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
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)
...
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,...
an exclusion type
Definition of the atomic potential types.
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_skip_13
integer, parameter, public do_fist
integer, parameter, public do_skip_12
integer, parameter, public do_skip_14
objects that represent the structure of input sections and the data contained in an input section
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
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
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
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
Utility routines for the memory handling.
Define the molecule kind structure types and the corresponding functionality.
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.
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.
Define the data structure for the molecule information.
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.
Define the data structure for the particle information.
subroutine, public allocate_particle_set(particle_set, nparticle)
Allocate a particle set.
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public massunit
Definition: physcon.F:141
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition: string_table.F:22
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition: string_table.F:141
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
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition: string_table.F:115
Collection of subroutine needed for topology related things.
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.
Collection of subroutine needed for topology related things.
Definition: topology_util.F:13
Control for reading in different topologies and coordinates.
Definition: topology.F:13