Line data Source code
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 : ! **************************************************************************************************
13 : MODULE 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 :
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 9512 : 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 9512 : INTEGER, DIMENSION(:), POINTER :: iatomlist, id_element, id_work, kind_of, &
107 9512 : list, list2, molecule_list, &
108 9512 : natom_of_kind, wlist
109 9512 : 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 9512 : REAL(KIND=dp), DIMENSION(:), POINTER :: charge, cpoint, mass
115 9512 : TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bend_list, ex_bond_list, &
116 9512 : ex_bond_list_ei, ex_bond_list_vdw, &
117 9512 : ex_onfo_list
118 : TYPE(atom_info_type), POINTER :: atom_info
119 9512 : 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 9512 : NULLIFY (logger)
129 19024 : logger => cp_get_default_logger()
130 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
131 9512 : extension=".subsysLog")
132 9512 : topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
133 9512 : CALL timeset(routineN, handle)
134 :
135 9512 : my_qmmm = .FALSE.
136 9512 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
137 9512 : atom_info => topology%atom_info
138 9512 : 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 9512 : CALL timeset(routineN//'_1', handle2)
145 9512 : counter = 0
146 9512 : NULLIFY (id_work, mass, id_element, charge)
147 28536 : ALLOCATE (id_work(topology%natoms))
148 28536 : ALLOCATE (mass(topology%natoms))
149 19024 : ALLOCATE (id_element(topology%natoms))
150 19024 : ALLOCATE (charge(topology%natoms))
151 760063 : id_work = str2id(s2s(""))
152 9512 : IF (iw > 0) WRITE (iw, *) "molecule_kind_set ::", SIZE(molecule_kind_set)
153 141933 : DO i = 1, SIZE(molecule_kind_set)
154 132421 : j = molecule_kind_set(i)%molecule_list(1)
155 132421 : molecule => molecule_set(j)
156 132421 : molecule_kind => molecule_set(j)%molecule_kind
157 132421 : IF (iw > 0) WRITE (iw, *) "molecule number ::", j, " has molecule kind number ::", i
158 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
159 132421 : natom=natom, atom_list=atom_list)
160 : CALL get_molecule(molecule=molecule, &
161 132421 : first_atom=first, last_atom=last)
162 132421 : IF (iw > 0) WRITE (iw, *) "boundaries of molecules (first, last) ::", first, last
163 531608 : DO j = 1, natom
164 18747125 : IF (.NOT. ANY(id_work(1:counter) .EQ. atom_list(j)%id_name)) THEN
165 24071 : counter = counter + 1
166 24071 : id_work(counter) = atom_list(j)%id_name
167 24071 : mass(counter) = atom_info%atm_mass(first + j - 1)
168 24071 : id_element(counter) = atom_info%id_element(first + j - 1)
169 24071 : charge(counter) = atom_info%atm_charge(first + j - 1)
170 24071 : IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
171 83 : "NEW ATOMIC KIND", id2str(id_work(counter)), mass(counter), id2str(id_element(counter)), charge(counter)
172 : ELSE
173 18443683 : found = .FALSE.
174 18443683 : DO k = 1, counter
175 18443683 : 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 233183 : IF (.NOT. found) THEN
181 524 : counter = counter + 1
182 524 : id_work(counter) = atom_list(j)%id_name
183 524 : mass(counter) = atom_info%atm_mass(first + j - 1)
184 524 : id_element(counter) = atom_info%id_element(first + j - 1)
185 524 : charge(counter) = atom_info%atm_charge(first + j - 1)
186 524 : IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
187 0 : "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 9512 : topology%natom_type = counter
193 28536 : ALLOCATE (atom_info%id_atom_names(topology%natom_type))
194 34107 : DO k = 1, counter
195 34107 : atom_info%id_atom_names(k) = id_work(k)
196 : END DO
197 9512 : DEALLOCATE (id_work)
198 9512 : CALL reallocate(mass, 1, counter)
199 9512 : CALL reallocate(id_element, 1, counter)
200 9512 : CALL reallocate(charge, 1, counter)
201 9512 : IF (iw > 0) &
202 26 : WRITE (iw, '(5X,A,I3)') "Total Number of Atomic Kinds = ", topology%natom_type
203 9512 : CALL timestop(handle2)
204 :
205 : !-----------------------------------------------------------------------------
206 : !-----------------------------------------------------------------------------
207 : ! 2. Allocate the data structure for the atomic kind information
208 : !-----------------------------------------------------------------------------
209 9512 : CALL timeset(routineN//'_2', handle2)
210 9512 : NULLIFY (atomic_kind_set)
211 53131 : ALLOCATE (atomic_kind_set(topology%natom_type))
212 9512 : CALL timestop(handle2)
213 :
214 : !-----------------------------------------------------------------------------
215 : !-----------------------------------------------------------------------------
216 : ! 3. Allocate the data structure for the atomic information
217 : !-----------------------------------------------------------------------------
218 9512 : CALL timeset(routineN//'_3', handle2)
219 9512 : NULLIFY (particle_set)
220 9512 : CALL allocate_particle_set(particle_set, topology%natoms)
221 9512 : CALL timestop(handle2)
222 :
223 : !-----------------------------------------------------------------------------
224 : !-----------------------------------------------------------------------------
225 : ! 4. Set the atomic_kind_set(ikind)%[name,kind_number,mass]
226 : !-----------------------------------------------------------------------------
227 9512 : CALL timeset(routineN//'_4', handle2)
228 34107 : DO i = 1, topology%natom_type
229 24595 : atomic_kind => atomic_kind_set(i)
230 24595 : 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 24595 : mass=mass(i))
236 34107 : IF (iw > 0) THEN
237 83 : WRITE (iw, '(A,I5,A,I5,4A)') "Atomic Kind n.:", i, " out of:", topology%natom_type, &
238 83 : " name: ", TRIM(id2str(atom_info%id_atom_names(i))), " element: ", &
239 166 : TRIM(id2str(id_element(i)))
240 : END IF
241 : END DO
242 9512 : DEALLOCATE (mass)
243 9512 : DEALLOCATE (id_element)
244 9512 : CALL timestop(handle2)
245 :
246 : !-----------------------------------------------------------------------------
247 : !-----------------------------------------------------------------------------
248 : ! 5. Determine number of atom of each kind (ie natom_of_kind and kind_of)
249 : !-----------------------------------------------------------------------------
250 9512 : CALL timeset(routineN//'_5', handle2)
251 28536 : ALLOCATE (kind_of(topology%natoms))
252 28536 : ALLOCATE (natom_of_kind(topology%natom_type))
253 760063 : kind_of(:) = 0
254 34107 : natom_of_kind(:) = 0
255 34107 : DO i = 1, topology%natom_type
256 37078597 : DO j = 1, topology%natoms
257 37069085 : IF ((atom_info%id_atom_names(i) == atom_info%id_atmname(j)) .AND. (charge(i) == atom_info%atm_charge(j))) THEN
258 750551 : natom_of_kind(i) = natom_of_kind(i) + 1
259 750551 : IF (kind_of(j) == 0) kind_of(j) = i
260 : END IF
261 : END DO
262 : END DO
263 760063 : IF (ANY(kind_of == 0)) THEN
264 0 : DO i = 1, topology%natoms
265 0 : IF (kind_of(i) == 0) THEN
266 0 : WRITE (*, *) i, kind_of(i)
267 0 : WRITE (*, *) "Two molecules have been defined as identical molecules but atoms mismatch charges!"
268 : END IF
269 : END DO
270 0 : CPABORT("")
271 : END IF
272 9512 : CALL timestop(handle2)
273 :
274 : !-----------------------------------------------------------------------------
275 : !-----------------------------------------------------------------------------
276 : ! 6. Set the atom_kind_set(ikind)%[natom,atom_list]
277 : !-----------------------------------------------------------------------------
278 9512 : CALL timeset(routineN//'_6', handle2)
279 34107 : DO i = 1, topology%natom_type
280 24595 : atomic_kind => atomic_kind_set(i)
281 : NULLIFY (iatomlist)
282 73785 : ALLOCATE (iatomlist(natom_of_kind(i)))
283 24595 : counter = 0
284 37069085 : DO j = 1, topology%natoms
285 37069085 : IF (kind_of(j) == i) THEN
286 750551 : counter = counter + 1
287 750551 : iatomlist(counter) = j
288 : END IF
289 : END DO
290 24595 : IF (iw > 0) THEN
291 83 : WRITE (iw, '(A,I6,A)') " Atomic kind ", i, " contains particles"
292 1614 : DO J = 1, SIZE(iatomlist)
293 1614 : IF (MOD(J, 5) .EQ. 0) THEN ! split long lines
294 271 : WRITE (iw, '(I12)') iatomlist(J)
295 : ELSE
296 1260 : WRITE (iw, '(I12)', ADVANCE="NO") iatomlist(J)
297 : END IF
298 : END DO
299 83 : WRITE (iw, *)
300 : END IF
301 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
302 : natom=natom_of_kind(i), &
303 24595 : atom_list=iatomlist)
304 34107 : DEALLOCATE (iatomlist)
305 : END DO
306 9512 : DEALLOCATE (natom_of_kind)
307 9512 : 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 9512 : "TOPOLOGY%CENTER_COORDINATES%_SECTION_PARAMETERS_", l_val=do_center)
315 9512 : CALL timeset(routineN//'_7a', handle2)
316 769575 : bounds(1, 1) = MINVAL(atom_info%r(1, :))
317 769575 : bounds(2, 1) = MAXVAL(atom_info%r(1, :))
318 :
319 769575 : bounds(1, 2) = MINVAL(atom_info%r(2, :))
320 769575 : bounds(2, 2) = MAXVAL(atom_info%r(2, :))
321 :
322 769575 : bounds(1, 3) = MINVAL(atom_info%r(3, :))
323 769575 : bounds(2, 3) = MAXVAL(atom_info%r(3, :))
324 :
325 38048 : dims = bounds(2, :) - bounds(1, :)
326 9512 : cdims(1) = topology%cell%hmat(1, 1)
327 9512 : cdims(2) = topology%cell%hmat(2, 2)
328 9512 : cdims(3) = topology%cell%hmat(3, 3)
329 9512 : IF (iw > 0) THEN
330 26 : WRITE (iw, '(A,3F12.6)') "System sizes: ", dims, "Cell sizes (diagonal): ", cdims
331 : END IF
332 9512 : check = .TRUE.
333 38048 : DO i = 1, 3
334 38048 : IF (topology%cell%perd(i) == 0) THEN
335 8886 : check = check .AND. (dims(i) < cdims(i))
336 : END IF
337 : END DO
338 9512 : my_ignore_outside_box = .FALSE.
339 9512 : IF (PRESENT(ignore_outside_box)) my_ignore_outside_box = ignore_outside_box
340 9512 : 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 0 : "exceeds the cell size in at least one of the non-periodic directions!")
344 9512 : IF (do_center) THEN
345 : CALL section_vals_val_get(subsys_section, &
346 1520 : "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", explicit=explicit)
347 1520 : IF (explicit) THEN
348 : CALL section_vals_val_get(subsys_section, &
349 0 : "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", r_vals=cpoint)
350 0 : vec = cpoint
351 : ELSE
352 6080 : vec = cdims/2.0_dp
353 : END IF
354 6080 : dims = (bounds(2, :) + bounds(1, :))/2.0_dp - vec
355 : ELSE
356 7992 : dims = 0.0_dp
357 : END IF
358 9512 : CALL timestop(handle2)
359 9512 : CALL timeset(routineN//'_7b', handle2)
360 760063 : DO i = 1, topology%natoms
361 750551 : ikind = kind_of(i)
362 750551 : IF (iw > 0) THEN
363 1531 : WRITE (iw, *) "atom number :: ", i, "kind number ::", ikind
364 : END IF
365 750551 : particle_set(i)%atomic_kind => atomic_kind_set(ikind)
366 5253857 : particle_set(i)%r(:) = atom_info%r(:, i) - dims
367 760063 : particle_set(i)%atom_index = i
368 : END DO
369 9512 : CALL timestop(handle2)
370 9512 : 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 9512 : CALL timeset(routineN//'_89', handle2)
379 9512 : 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 9512 : l_val=disable_exclusion_lists)
382 9512 : IF ((method_name_id == do_fist) .AND. (.NOT. disable_exclusion_lists)) THEN
383 2494 : CPASSERT(PRESENT(exclusions))
384 2494 : natom = topology%natoms
385 : ! allocate exclusions. Most likely they would only be needed for the local_particles
386 639848 : ALLOCATE (exclusions(natom))
387 634860 : DO I = 1, natom
388 632366 : NULLIFY (exclusions(i)%list_exclude_vdw)
389 632366 : NULLIFY (exclusions(i)%list_exclude_ei)
390 634860 : NULLIFY (exclusions(i)%list_onfo)
391 : END DO
392 : ! Reorder bonds
393 639848 : ALLOCATE (ex_bond_list(natom))
394 634860 : DO I = 1, natom
395 634860 : ALLOCATE (ex_bond_list(I)%array1(0))
396 : END DO
397 2494 : N = 0
398 2494 : IF (ASSOCIATED(conn_info%bond_a)) THEN
399 2494 : N = SIZE(conn_info%bond_a)
400 2494 : 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 2494 : NULLIFY (ex_bond_list_vdw, ex_bond_list_ei)
405 : ! VdW
406 2494 : exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_VDW_LIST")
407 2494 : CALL section_vals_get(exclude_section, explicit=explicit)
408 2494 : present_12_excl_vdw_list = .FALSE.
409 2494 : IF (explicit) present_12_excl_vdw_list = .TRUE.
410 : IF (present_12_excl_vdw_list) THEN
411 40 : ALLOCATE (ex_bond_list_vdw(natom))
412 32 : DO I = 1, natom
413 32 : 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 8 : particle_set)
417 : ELSE
418 2486 : ex_bond_list_vdw => ex_bond_list
419 : END IF
420 : ! EI
421 2494 : exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_EI_LIST")
422 2494 : CALL section_vals_get(exclude_section, explicit=explicit)
423 2494 : present_12_excl_ei_list = .FALSE.
424 2494 : IF (explicit) present_12_excl_ei_list = .TRUE.
425 : IF (present_12_excl_ei_list) THEN
426 50 : ALLOCATE (ex_bond_list_ei(natom))
427 40 : DO I = 1, natom
428 40 : 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 10 : particle_set)
432 : ELSE
433 2484 : ex_bond_list_ei => ex_bond_list
434 : END IF
435 :
436 : CALL section_vals_val_get(topology_section, "AUTOGEN_EXCLUDE_LISTS", &
437 2494 : l_val=autogen)
438 : ! Reorder bends
439 637354 : ALLOCATE (ex_bend_list(natom))
440 634860 : DO I = 1, natom
441 634860 : ALLOCATE (ex_bend_list(I)%array1(0))
442 : END DO
443 2494 : 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 4 : ALLOCATE (pairs(0, 2))
447 4 : N = 0
448 26 : DO iatom = 1, natom
449 62 : DO i = 1, SIZE(ex_bond_list(iatom)%array1)
450 : ! a neighboring atom of iatom:
451 36 : atom_i = ex_bond_list(iatom)%array1(i)
452 92 : DO j = 1, i - 1
453 : ! another neighboring atom of iatom
454 34 : 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 34 : check = .FALSE.
459 70 : DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
460 70 : IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
461 : check = .TRUE.
462 : EXIT
463 : END IF
464 : END DO
465 34 : IF (check) CYCLE
466 : ! Add the genuine 1-3 pair
467 34 : N = N + 1
468 34 : IF (SIZE(pairs, dim=1) <= N) THEN
469 8 : CALL reallocate(pairs, 1, N + 5, 1, 2)
470 : END IF
471 34 : pairs(N, 1) = atom_i
472 70 : pairs(N, 2) = atom_j
473 : END DO
474 : END DO
475 : END DO
476 4 : CALL reorder_structure(ex_bend_list, pairs(:, 1), pairs(:, 2), N)
477 4 : DEALLOCATE (pairs)
478 : ELSE
479 2490 : IF (ASSOCIATED(conn_info%theta_a)) THEN
480 2490 : N = SIZE(conn_info%theta_a)
481 2490 : 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 637354 : ALLOCATE (ex_onfo_list(natom))
487 634860 : DO I = 1, natom
488 634860 : ALLOCATE (ex_onfo_list(I)%array1(0))
489 : END DO
490 2494 : 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 4 : ALLOCATE (pairs(0, 2))
494 4 : N = 0
495 26 : DO iatom = 1, natom
496 62 : DO i = 1, SIZE(ex_bond_list(iatom)%array1)
497 : ! a neighboring atom of iatom:
498 36 : atom_i = ex_bond_list(iatom)%array1(i)
499 130 : DO j = 1, SIZE(ex_bend_list(iatom)%array1)
500 : ! a next neighboring atom of iatom:
501 72 : 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 72 : IF (atom_i == atom_j) CYCLE
505 : ! check if i and j are not involved in a bond
506 72 : check = .FALSE.
507 230 : DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
508 230 : IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
509 : check = .TRUE.
510 : EXIT
511 : END IF
512 : END DO
513 72 : IF (check) CYCLE
514 : ! check if i and j are not involved in a bend
515 4 : check = .FALSE.
516 8 : DO counter = 1, SIZE(ex_bend_list(atom_i)%array1)
517 8 : IF (ex_bend_list(atom_i)%array1(counter) == atom_j) THEN
518 : check = .TRUE.
519 : EXIT
520 : END IF
521 : END DO
522 4 : IF (check) CYCLE
523 : ! Add the true onfo.
524 4 : N = N + 1
525 4 : IF (SIZE(pairs, dim=1) <= N) THEN
526 2 : CALL reallocate(pairs, 1, N + 5, 1, 2)
527 : END IF
528 4 : pairs(N, 1) = atom_i
529 108 : pairs(N, 2) = atom_j
530 : END DO
531 : END DO
532 : END DO
533 4 : CALL reorder_structure(ex_onfo_list, pairs(:, 1), pairs(:, 2), N)
534 4 : DEALLOCATE (pairs)
535 : ELSE
536 2490 : IF (ASSOCIATED(conn_info%onfo_a)) THEN
537 2484 : N = SIZE(conn_info%onfo_a)
538 2484 : 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 634860 : DO iatom = 1, SIZE(particle_set)
544 : ! Setup exclusion list for VDW: always exclude itself
545 632366 : dim0 = 1
546 : ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
547 632366 : dim1 = 0
548 : IF (topology%exclude_vdw == do_skip_12 .OR. &
549 632366 : topology%exclude_vdw == do_skip_13 .OR. &
550 631772 : topology%exclude_vdw == do_skip_14) dim1 = SIZE(ex_bond_list_vdw(iatom)%array1)
551 632366 : dim1 = dim1 + dim0
552 632366 : dim2 = 0
553 632366 : IF (topology%exclude_vdw == do_skip_13 .OR. &
554 631166 : topology%exclude_vdw == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
555 632366 : dim2 = dim1 + dim2
556 632366 : dim3 = 0
557 632366 : IF (topology%exclude_vdw == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
558 632366 : dim3 = dim2 + dim3
559 632366 : IF (dim3 /= 0) THEN
560 632366 : NULLIFY (list, wlist)
561 1897098 : ALLOCATE (wlist(dim3))
562 1264732 : wlist(dim0:dim0) = iatom
563 1602504 : IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_vdw(iatom)%array1
564 1152872 : IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
565 633914 : IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
566 : ! Get a unique list
567 2124558 : DO i = 1, SIZE(wlist) - 1
568 1492192 : IF (wlist(i) == 0) CYCLE
569 5405448 : DO j = i + 1, SIZE(wlist)
570 4773152 : IF (wlist(j) == wlist(i)) wlist(j) = 0
571 : END DO
572 : END DO
573 2756924 : dim3 = SIZE(wlist) - COUNT(wlist == 0)
574 1897098 : ALLOCATE (list(dim3))
575 632366 : j = 0
576 2756924 : DO i = 1, SIZE(wlist)
577 2124558 : IF (wlist(i) == 0) CYCLE
578 2038674 : j = j + 1
579 2756924 : list(j) = wlist(i)
580 : END DO
581 632366 : DEALLOCATE (wlist)
582 : ! Unique list completed
583 632366 : NULLIFY (list2)
584 : IF ((topology%exclude_vdw == topology%exclude_ei) .AND. &
585 632366 : (.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 1770 : dim0 = 1
590 : ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
591 1770 : dim1 = 0
592 : IF (topology%exclude_ei == do_skip_12 .OR. &
593 1770 : topology%exclude_ei == do_skip_13 .OR. &
594 1326 : topology%exclude_ei == do_skip_14) dim1 = SIZE(ex_bond_list_ei(iatom)%array1)
595 1770 : dim1 = dim1 + dim0
596 1770 : dim2 = 0
597 1770 : IF (topology%exclude_ei == do_skip_13 .OR. &
598 864 : topology%exclude_ei == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
599 1770 : dim2 = dim1 + dim2
600 1770 : dim3 = 0
601 1770 : IF (topology%exclude_ei == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
602 1770 : dim3 = dim2 + dim3
603 :
604 1770 : IF (dim3 /= 0) THEN
605 5310 : ALLOCATE (wlist(dim3))
606 3540 : wlist(dim0:dim0) = iatom
607 4098 : IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_ei(iatom)%array1
608 4266 : IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
609 2922 : IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
610 : ! Get a unique list
611 7746 : DO i = 1, SIZE(wlist) - 1
612 5976 : IF (wlist(i) == 0) CYCLE
613 28896 : DO j = i + 1, SIZE(wlist)
614 27126 : IF (wlist(j) == wlist(i)) wlist(j) = 0
615 : END DO
616 : END DO
617 9516 : dim3 = SIZE(wlist) - COUNT(wlist == 0)
618 5310 : ALLOCATE (list2(dim3))
619 1770 : j = 0
620 9516 : DO i = 1, SIZE(wlist)
621 7746 : IF (wlist(i) == 0) CYCLE
622 7746 : j = j + 1
623 9516 : list2(j) = wlist(i)
624 : END DO
625 1770 : DEALLOCATE (wlist)
626 : ! Unique list completed
627 : END IF
628 : END IF
629 : END IF
630 632366 : exclusions(iatom)%list_exclude_vdw => list
631 632366 : 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 1339271 : ALLOCATE (exclusions(iatom)%list_onfo(SIZE(ex_onfo_list(iatom)%array1)))
635 : ! copy of data, not copy of pointer
636 1244186 : exclusions(iatom)%list_onfo = ex_onfo_list(iatom)%array1
637 634860 : IF (iw > 0) THEN
638 1243 : IF (ASSOCIATED(list)) &
639 1243 : WRITE (iw, *) "exclusion list_vdw :: ", &
640 1243 : "atom num :", iatom, "exclusion list ::", &
641 6003 : list
642 1243 : IF (topology%exclude_vdw /= topology%exclude_ei) THEN
643 9 : IF (ASSOCIATED(list2)) &
644 9 : WRITE (iw, *) "exclusion list_ei :: ", &
645 9 : "atom num :", iatom, "exclusion list ::", &
646 35 : list2
647 : END IF
648 1243 : IF (ASSOCIATED(exclusions(iatom)%list_onfo)) &
649 1243 : WRITE (iw, *) "onfo list :: ", &
650 1243 : "atom num :", iatom, "onfo list ::", &
651 3320 : exclusions(iatom)%list_onfo
652 : END IF
653 : END DO
654 : ! deallocate onfo
655 634860 : DO I = 1, natom
656 634860 : DEALLOCATE (ex_onfo_list(I)%array1)
657 : END DO
658 2494 : DEALLOCATE (ex_onfo_list)
659 : ! deallocate bends
660 634860 : DO I = 1, natom
661 634860 : DEALLOCATE (ex_bend_list(I)%array1)
662 : END DO
663 2494 : DEALLOCATE (ex_bend_list)
664 : ! deallocate bonds
665 2494 : IF (present_12_excl_ei_list) THEN
666 40 : DO I = 1, natom
667 40 : DEALLOCATE (ex_bond_list_ei(I)%array1)
668 : END DO
669 10 : DEALLOCATE (ex_bond_list_ei)
670 : ELSE
671 2484 : NULLIFY (ex_bond_list_ei)
672 : END IF
673 2494 : IF (present_12_excl_vdw_list) THEN
674 32 : DO I = 1, natom
675 32 : DEALLOCATE (ex_bond_list_vdw(I)%array1)
676 : END DO
677 8 : DEALLOCATE (ex_bond_list_vdw)
678 : ELSE
679 2486 : NULLIFY (ex_bond_list_vdw)
680 : END IF
681 634860 : DO I = 1, natom
682 634860 : DEALLOCATE (ex_bond_list(I)%array1)
683 : END DO
684 9976 : DEALLOCATE (ex_bond_list)
685 : END IF
686 9512 : CALL timestop(handle2)
687 : !-----------------------------------------------------------------------------
688 : !-----------------------------------------------------------------------------
689 : ! 11. Set the atomic_kind_set()%fist_potential%[qeff] (PART 1)
690 : !-----------------------------------------------------------------------------
691 9512 : CALL timeset(routineN//'_10', handle2)
692 9512 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
693 9512 : IF (method_name_id == do_fist) THEN
694 13912 : DO i = 1, SIZE(atomic_kind_set)
695 11270 : atomic_kind => atomic_kind_set(i)
696 11270 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
697 11270 : qeff = charge(i)
698 11270 : NULLIFY (fist_potential)
699 11270 : CALL allocate_potential(fist_potential)
700 11270 : CALL set_potential(potential=fist_potential, qeff=qeff)
701 13912 : CALL set_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
702 : END DO
703 : END IF
704 9512 : DEALLOCATE (charge)
705 9512 : CALL timestop(handle2)
706 :
707 : !-----------------------------------------------------------------------------
708 : !-----------------------------------------------------------------------------
709 : ! 12. Set the atom_list for molecule_kind in molecule_kind_set (PART 2)
710 : !-----------------------------------------------------------------------------
711 9512 : CALL timeset(routineN//'_11', handle2)
712 141933 : DO i = 1, SIZE(molecule_kind_set)
713 132421 : molecule_kind => molecule_kind_set(i)
714 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
715 : natom=natom, molecule_list=molecule_list, &
716 132421 : atom_list=atom_list)
717 132421 : molecule => molecule_set(molecule_list(1))
718 : CALL get_molecule(molecule=molecule, &
719 132421 : first_atom=first, last_atom=last)
720 389675 : DO j = 1, natom
721 18839301 : DO k = 1, SIZE(atomic_kind_set)
722 18706880 : atomic_kind => atomic_kind_set(k)
723 18706880 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
724 18706880 : IF (method_name_id == do_fist) THEN
725 18385419 : CALL get_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
726 18385419 : CALL get_potential(potential=fist_potential, qeff=qeff)
727 18385419 : IF ((id2str(atom_list(j)%id_name) == atmname) .AND. (qeff == atom_info%atm_charge(first + j - 1))) THEN
728 194785 : atom_list(j)%atomic_kind => atomic_kind_set(k)
729 18580204 : EXIT
730 : END IF
731 : ELSE
732 321461 : IF (id2str(atom_list(j)%id_name) == atmname) THEN
733 62469 : atom_list(j)%atomic_kind => atomic_kind_set(k)
734 383930 : EXIT
735 : END IF
736 : END IF
737 : END DO
738 : END DO
739 274354 : CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list)
740 : END DO
741 9512 : CALL timestop(handle2)
742 :
743 9512 : CALL timestop(handle)
744 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
745 9512 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
746 114144 : 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 18 : 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 18 : DIMENSION(:), POINTER :: names
769 : INTEGER :: i, ind, j, k, l, m, n_rep
770 :
771 0 : CPASSERT(ASSOCIATED(ex_bond_list))
772 18 : CPASSERT(ASSOCIATED(ex_bond_list_w))
773 18 : SELECT CASE (keyword)
774 : CASE ("BOND")
775 18 : CALL section_vals_val_get(exclude_section, keyword, n_rep_val=n_rep)
776 72 : DO j = 1, SIZE(ex_bond_list)
777 54 : CPASSERT(ASSOCIATED(ex_bond_list(j)%array1))
778 54 : CPASSERT(ASSOCIATED(ex_bond_list_w(j)%array1))
779 :
780 54 : flag1 = particle_set(j)%atomic_kind%name
781 54 : m = SIZE(ex_bond_list(j)%array1)
782 54 : CALL reallocate(ex_bond_list_w(j)%array1, 1, m)
783 :
784 54 : l = 0
785 126 : DO k = 1, m
786 72 : ind = ex_bond_list(j)%array1(k)
787 72 : flag2 = particle_set(ind)%atomic_kind%name
788 150 : DO i = 1, n_rep
789 : CALL section_vals_val_get(exclude_section, keyword, i_rep_val=i, &
790 24 : c_vals=names)
791 24 : IF (((TRIM(names(1)) == TRIM(flag1)) .AND. (TRIM(names(2)) == TRIM(flag2))) .OR. &
792 72 : ((TRIM(names(1)) == TRIM(flag2)) .AND. (TRIM(names(2)) == TRIM(flag1)))) THEN
793 24 : l = l + 1
794 24 : ex_bond_list_w(j)%array1(l) = ind
795 : END IF
796 : END DO
797 : END DO
798 72 : CALL reallocate(ex_bond_list_w(j)%array1, 1, l)
799 : END DO
800 : CASE DEFAULT
801 18 : CPABORT("")
802 : END SELECT
803 :
804 18 : END SUBROUTINE setup_exclusion_list
805 :
806 : END MODULE topology_coordinate_util
|
