Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 11500 : 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, err
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 11500 : INTEGER, DIMENSION(:), POINTER :: iatomlist, id_element, id_work, kind_of, &
107 11500 : list, list2, molecule_list, &
108 11500 : natom_of_kind, wlist
109 11500 : 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 11500 : REAL(KIND=dp), DIMENSION(:), POINTER :: charge, cpoint, mass
115 11500 : TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bend_list, ex_bond_list, &
116 11500 : ex_bond_list_ei, ex_bond_list_vdw, &
117 11500 : ex_onfo_list
118 : TYPE(atom_info_type), POINTER :: atom_info
119 11500 : 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 11500 : NULLIFY (logger)
129 23000 : logger => cp_get_default_logger()
130 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
131 11500 : extension=".subsysLog")
132 11500 : topology_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY")
133 11500 : CALL timeset(routineN, handle)
134 :
135 11500 : my_qmmm = .FALSE.
136 11500 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env)) my_qmmm = qmmm
137 11500 : atom_info => topology%atom_info
138 11500 : 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 11500 : CALL timeset(routineN//'_1', handle2)
145 11500 : counter = 0
146 11500 : NULLIFY (id_work, mass, id_element, charge)
147 34500 : ALLOCATE (id_work(topology%natoms))
148 34500 : ALLOCATE (mass(topology%natoms))
149 23000 : ALLOCATE (id_element(topology%natoms))
150 23000 : ALLOCATE (charge(topology%natoms))
151 777463 : id_work = str2id(s2s(""))
152 11500 : IF (iw > 0) WRITE (iw, *) "molecule_kind_set ::", SIZE(molecule_kind_set)
153 158459 : DO i = 1, SIZE(molecule_kind_set)
154 146959 : j = molecule_kind_set(i)%molecule_list(1)
155 146959 : molecule => molecule_set(j)
156 146959 : molecule_kind => molecule_set(j)%molecule_kind
157 146959 : IF (iw > 0) WRITE (iw, *) "molecule number ::", j, " has molecule kind number ::", i
158 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
159 146959 : natom=natom, atom_list=atom_list)
160 : CALL get_molecule(molecule=molecule, &
161 146959 : first_atom=first, last_atom=last)
162 146959 : IF (iw > 0) WRITE (iw, *) "boundaries of molecules (first, last) ::", first, last
163 566132 : DO j = 1, natom
164 18762951 : IF (.NOT. ANY(id_work(1:counter) == atom_list(j)%id_name)) THEN
165 27695 : counter = counter + 1
166 27695 : id_work(counter) = atom_list(j)%id_name
167 27695 : mass(counter) = atom_info%atm_mass(first + j - 1)
168 27695 : id_element(counter) = atom_info%id_element(first + j - 1)
169 27695 : charge(counter) = atom_info%atm_charge(first + j - 1)
170 27695 : IF (iw > 0) WRITE (iw, '(7X,A,1X,A5,F10.5,5X,A2,5X,F10.5)') &
171 84 : "NEW ATOMIC KIND", id2str(id_work(counter)), mass(counter), id2str(id_element(counter)), charge(counter)
172 : ELSE
173 18439423 : found = .FALSE.
174 18439423 : DO k = 1, counter
175 18439423 : 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 233019 : 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 11500 : topology%natom_type = counter
193 34500 : ALLOCATE (atom_info%id_atom_names(topology%natom_type))
194 39719 : DO k = 1, counter
195 39719 : atom_info%id_atom_names(k) = id_work(k)
196 : END DO
197 11500 : DEALLOCATE (id_work)
198 11500 : CALL reallocate(mass, 1, counter)
199 11500 : CALL reallocate(id_element, 1, counter)
200 11500 : CALL reallocate(charge, 1, counter)
201 11500 : IF (iw > 0) THEN
202 27 : WRITE (iw, '(5X,A,I3)') "Total Number of Atomic Kinds = ", topology%natom_type
203 : END IF
204 11500 : CALL timestop(handle2)
205 :
206 : !-----------------------------------------------------------------------------
207 : !-----------------------------------------------------------------------------
208 : ! 2. Allocate the data structure for the atomic kind information
209 : !-----------------------------------------------------------------------------
210 11500 : CALL timeset(routineN//'_2', handle2)
211 11500 : NULLIFY (atomic_kind_set)
212 62719 : ALLOCATE (atomic_kind_set(topology%natom_type))
213 11500 : CALL timestop(handle2)
214 :
215 : !-----------------------------------------------------------------------------
216 : !-----------------------------------------------------------------------------
217 : ! 3. Allocate the data structure for the atomic information
218 : !-----------------------------------------------------------------------------
219 11500 : CALL timeset(routineN//'_3', handle2)
220 11500 : NULLIFY (particle_set)
221 11500 : CALL allocate_particle_set(particle_set, topology%natoms)
222 11500 : CALL timestop(handle2)
223 :
224 : !-----------------------------------------------------------------------------
225 : !-----------------------------------------------------------------------------
226 : ! 4. Set the atomic_kind_set(ikind)%[name,kind_number,mass]
227 : !-----------------------------------------------------------------------------
228 11500 : CALL timeset(routineN//'_4', handle2)
229 39719 : DO i = 1, topology%natom_type
230 28219 : atomic_kind => atomic_kind_set(i)
231 28219 : mass(i) = mass(i)*massunit
232 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
233 : kind_number=i, &
234 : name=id2str(atom_info%id_atom_names(i)), &
235 : element_symbol=id2str(id_element(i)), &
236 28219 : mass=mass(i))
237 39719 : IF (iw > 0) THEN
238 84 : WRITE (iw, '(A,I5,A,I5,4A)') "Atomic Kind n.:", i, " out of:", topology%natom_type, &
239 84 : " name: ", TRIM(id2str(atom_info%id_atom_names(i))), " element: ", &
240 168 : TRIM(id2str(id_element(i)))
241 : END IF
242 : END DO
243 11500 : DEALLOCATE (mass)
244 11500 : DEALLOCATE (id_element)
245 11500 : CALL timestop(handle2)
246 :
247 : !-----------------------------------------------------------------------------
248 : !-----------------------------------------------------------------------------
249 : ! 5. Determine number of atom of each kind (ie natom_of_kind and kind_of)
250 : !-----------------------------------------------------------------------------
251 11500 : CALL timeset(routineN//'_5', handle2)
252 34500 : ALLOCATE (kind_of(topology%natoms))
253 34500 : ALLOCATE (natom_of_kind(topology%natom_type))
254 777463 : kind_of(:) = 0
255 39719 : natom_of_kind(:) = 0
256 39719 : DO i = 1, topology%natom_type
257 37108489 : DO j = 1, topology%natoms
258 37096989 : IF ((atom_info%id_atom_names(i) == atom_info%id_atmname(j)) .AND. (charge(i) == atom_info%atm_charge(j))) THEN
259 765963 : natom_of_kind(i) = natom_of_kind(i) + 1
260 765963 : IF (kind_of(j) == 0) kind_of(j) = i
261 : END IF
262 : END DO
263 : END DO
264 777463 : IF (ANY(kind_of == 0)) THEN
265 0 : DO i = 1, topology%natoms
266 0 : IF (kind_of(i) == 0) THEN
267 0 : WRITE (err, "(1X,I5)") i
268 : END IF
269 : END DO
270 : CALL cp_abort(__LOCATION__, &
271 : "Two molecules have been defined as identical molecules but atoms "// &
272 0 : "mismatch charges! Check these atoms:"//TRIM(err))
273 : END IF
274 11500 : CALL timestop(handle2)
275 :
276 : !-----------------------------------------------------------------------------
277 : !-----------------------------------------------------------------------------
278 : ! 6. Set the atom_kind_set(ikind)%[natom,atom_list]
279 : !-----------------------------------------------------------------------------
280 11500 : CALL timeset(routineN//'_6', handle2)
281 39719 : DO i = 1, topology%natom_type
282 28219 : atomic_kind => atomic_kind_set(i)
283 : NULLIFY (iatomlist)
284 84657 : ALLOCATE (iatomlist(natom_of_kind(i)))
285 28219 : counter = 0
286 37096989 : DO j = 1, topology%natoms
287 37096989 : IF (kind_of(j) == i) THEN
288 765963 : counter = counter + 1
289 765963 : iatomlist(counter) = j
290 : END IF
291 : END DO
292 28219 : IF (iw > 0) THEN
293 84 : WRITE (iw, '(A,I6,A)') " Atomic kind ", i, " contains particles"
294 1616 : DO J = 1, SIZE(iatomlist)
295 1616 : IF (MOD(J, 5) == 0) THEN ! split long lines
296 271 : WRITE (iw, '(I12)') iatomlist(J)
297 : ELSE
298 1261 : WRITE (iw, '(I12)', ADVANCE="NO") iatomlist(J)
299 : END IF
300 : END DO
301 84 : WRITE (iw, *)
302 : END IF
303 : CALL set_atomic_kind(atomic_kind=atomic_kind, &
304 : natom=natom_of_kind(i), &
305 28219 : atom_list=iatomlist)
306 39719 : DEALLOCATE (iatomlist)
307 : END DO
308 11500 : DEALLOCATE (natom_of_kind)
309 11500 : CALL timestop(handle2)
310 :
311 : !-----------------------------------------------------------------------------
312 : !-----------------------------------------------------------------------------
313 : ! 7. Possibly center the coordinates and fill in coordinates in particle_set
314 : !-----------------------------------------------------------------------------
315 : CALL section_vals_val_get(subsys_section, &
316 11500 : "TOPOLOGY%CENTER_COORDINATES%_SECTION_PARAMETERS_", l_val=do_center)
317 11500 : CALL timeset(routineN//'_7a', handle2)
318 777463 : bounds(1, 1) = MINVAL(atom_info%r(1, :))
319 777463 : bounds(2, 1) = MAXVAL(atom_info%r(1, :))
320 :
321 777463 : bounds(1, 2) = MINVAL(atom_info%r(2, :))
322 777463 : bounds(2, 2) = MAXVAL(atom_info%r(2, :))
323 :
324 777463 : bounds(1, 3) = MINVAL(atom_info%r(3, :))
325 777463 : bounds(2, 3) = MAXVAL(atom_info%r(3, :))
326 :
327 46000 : dims = bounds(2, :) - bounds(1, :)
328 11500 : cdims(1) = topology%cell%hmat(1, 1)
329 11500 : cdims(2) = topology%cell%hmat(2, 2)
330 11500 : cdims(3) = topology%cell%hmat(3, 3)
331 11500 : IF (iw > 0) THEN
332 27 : WRITE (iw, '(A,3F12.6)') "System sizes: ", dims, "Cell sizes (diagonal): ", cdims
333 : END IF
334 11500 : check = .TRUE.
335 46000 : DO i = 1, 3
336 46000 : IF (topology%cell%perd(i) == 0) THEN
337 12440 : check = check .AND. (dims(i) < cdims(i))
338 : END IF
339 : END DO
340 11500 : my_ignore_outside_box = .FALSE.
341 11500 : IF (PRESENT(ignore_outside_box)) my_ignore_outside_box = ignore_outside_box
342 11500 : IF (.NOT. my_ignore_outside_box .AND. .NOT. check) THEN
343 : CALL cp_abort(__LOCATION__, &
344 : "A non-periodic calculation has been requested but the system size "// &
345 0 : "exceeds the cell size in at least one of the non-periodic directions!")
346 : END IF
347 11500 : IF (do_center) THEN
348 : CALL section_vals_val_get(subsys_section, &
349 2440 : "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", explicit=explicit)
350 2440 : IF (explicit) THEN
351 : CALL section_vals_val_get(subsys_section, &
352 0 : "TOPOLOGY%CENTER_COORDINATES%CENTER_POINT", r_vals=cpoint)
353 0 : vec = cpoint
354 : ELSE
355 9760 : vec = cdims/2.0_dp
356 : END IF
357 9760 : dims = (bounds(2, :) + bounds(1, :))/2.0_dp - vec
358 : ELSE
359 9060 : dims = 0.0_dp
360 : END IF
361 11500 : CALL timestop(handle2)
362 11500 : CALL timeset(routineN//'_7b', handle2)
363 777463 : DO i = 1, topology%natoms
364 765963 : ikind = kind_of(i)
365 765963 : IF (iw > 0) THEN
366 1532 : WRITE (iw, *) "atom number :: ", i, "kind number ::", ikind
367 : END IF
368 765963 : particle_set(i)%atomic_kind => atomic_kind_set(ikind)
369 5361741 : particle_set(i)%r(:) = atom_info%r(:, i) - dims
370 777463 : particle_set(i)%atom_index = i
371 : END DO
372 11500 : CALL timestop(handle2)
373 11500 : DEALLOCATE (kind_of)
374 :
375 : !-----------------------------------------------------------------------------
376 : !-----------------------------------------------------------------------------
377 : ! 8. Fill in the exclusions%list_exclude_vdw
378 : ! 9. Fill in the exclusions%list_exclude_ei
379 : ! 10. Fill in the exclusions%list_onfo
380 : !-----------------------------------------------------------------------------
381 11500 : CALL timeset(routineN//'_89', handle2)
382 11500 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
383 : CALL section_vals_val_get(subsys_section, "TOPOLOGY%DISABLE_EXCLUSION_LISTS", &
384 11500 : l_val=disable_exclusion_lists)
385 11500 : IF ((method_name_id == do_fist) .AND. (.NOT. disable_exclusion_lists)) THEN
386 2486 : CPASSERT(PRESENT(exclusions))
387 2486 : natom = topology%natoms
388 : ! allocate exclusions. Most likely they would only be needed for the local_particles
389 640348 : ALLOCATE (exclusions(natom))
390 635376 : DO I = 1, natom
391 632890 : NULLIFY (exclusions(i)%list_exclude_vdw)
392 632890 : NULLIFY (exclusions(i)%list_exclude_ei)
393 635376 : NULLIFY (exclusions(i)%list_onfo)
394 : END DO
395 : ! Reorder bonds
396 640348 : ALLOCATE (ex_bond_list(natom))
397 635376 : DO I = 1, natom
398 635376 : ALLOCATE (ex_bond_list(I)%array1(0))
399 : END DO
400 2486 : N = 0
401 2486 : IF (ASSOCIATED(conn_info%bond_a)) THEN
402 2486 : N = SIZE(conn_info%bond_a)
403 2486 : CALL reorder_structure(ex_bond_list, conn_info%bond_a, conn_info%bond_b, N)
404 : END IF
405 :
406 : ! Check if a list of 1-2 exclusion bonds is defined.. if not use all bonds
407 : NULLIFY (ex_bond_list_vdw, ex_bond_list_ei)
408 : ! VdW
409 2486 : exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_VDW_LIST")
410 2486 : CALL section_vals_get(exclude_section, explicit=explicit)
411 2486 : present_12_excl_vdw_list = .FALSE.
412 2486 : IF (explicit) present_12_excl_vdw_list = .TRUE.
413 : IF (present_12_excl_vdw_list) THEN
414 40 : ALLOCATE (ex_bond_list_vdw(natom))
415 32 : DO I = 1, natom
416 32 : ALLOCATE (ex_bond_list_vdw(I)%array1(0))
417 : END DO
418 : CALL setup_exclusion_list(exclude_section, "BOND", ex_bond_list, ex_bond_list_vdw, &
419 8 : particle_set)
420 : ELSE
421 2478 : ex_bond_list_vdw => ex_bond_list
422 : END IF
423 : ! EI
424 2486 : exclude_section => section_vals_get_subs_vals(topology_section, "EXCLUDE_EI_LIST")
425 2486 : CALL section_vals_get(exclude_section, explicit=explicit)
426 2486 : present_12_excl_ei_list = .FALSE.
427 2486 : IF (explicit) present_12_excl_ei_list = .TRUE.
428 : IF (present_12_excl_ei_list) THEN
429 50 : ALLOCATE (ex_bond_list_ei(natom))
430 40 : DO I = 1, natom
431 40 : ALLOCATE (ex_bond_list_ei(I)%array1(0))
432 : END DO
433 : CALL setup_exclusion_list(exclude_section, "BOND", ex_bond_list, ex_bond_list_ei, &
434 10 : particle_set)
435 : ELSE
436 2476 : ex_bond_list_ei => ex_bond_list
437 : END IF
438 :
439 : CALL section_vals_val_get(topology_section, "AUTOGEN_EXCLUDE_LISTS", &
440 2486 : l_val=autogen)
441 : ! Reorder bends
442 637862 : ALLOCATE (ex_bend_list(natom))
443 635376 : DO I = 1, natom
444 635376 : ALLOCATE (ex_bend_list(I)%array1(0))
445 : END DO
446 2486 : IF (autogen) THEN
447 : ! Construct autogenerated 1-3 pairs, i.e. all possible 1-3 pairs instead
448 : ! of only the bends that are present in the topology.
449 4 : ALLOCATE (pairs(0, 2))
450 4 : N = 0
451 26 : DO iatom = 1, natom
452 62 : DO i = 1, SIZE(ex_bond_list(iatom)%array1)
453 : ! a neighboring atom of iatom:
454 36 : atom_i = ex_bond_list(iatom)%array1(i)
455 92 : DO j = 1, i - 1
456 : ! another neighboring atom of iatom
457 34 : atom_j = ex_bond_list(iatom)%array1(j)
458 : ! It is only a true bend if there is no shorter path.
459 : ! No need to check if i and j correspond to the same atom.
460 : ! Check if i and j are not involved in a bond:
461 34 : check = .FALSE.
462 70 : DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
463 70 : IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
464 : check = .TRUE.
465 : EXIT
466 : END IF
467 : END DO
468 34 : IF (check) CYCLE
469 : ! Add the genuine 1-3 pair
470 34 : N = N + 1
471 34 : IF (SIZE(pairs, dim=1) <= N) THEN
472 8 : CALL reallocate(pairs, 1, N + 5, 1, 2)
473 : END IF
474 34 : pairs(N, 1) = atom_i
475 70 : pairs(N, 2) = atom_j
476 : END DO
477 : END DO
478 : END DO
479 4 : CALL reorder_structure(ex_bend_list, pairs(:, 1), pairs(:, 2), N)
480 4 : DEALLOCATE (pairs)
481 : ELSE
482 2482 : IF (ASSOCIATED(conn_info%theta_a)) THEN
483 2482 : N = SIZE(conn_info%theta_a)
484 2482 : CALL reorder_structure(ex_bend_list, conn_info%theta_a, conn_info%theta_c, N)
485 : END IF
486 : END IF
487 :
488 : ! Reorder onfo
489 637862 : ALLOCATE (ex_onfo_list(natom))
490 635376 : DO I = 1, natom
491 635376 : ALLOCATE (ex_onfo_list(I)%array1(0))
492 : END DO
493 2486 : IF (autogen) THEN
494 : ! Construct autogenerated 1-4 pairs, i.e. all possible 1-4 pairs instead
495 : ! of only the onfo's that are present in the topology.
496 4 : ALLOCATE (pairs(0, 2))
497 4 : N = 0
498 26 : DO iatom = 1, natom
499 62 : DO i = 1, SIZE(ex_bond_list(iatom)%array1)
500 : ! a neighboring atom of iatom:
501 36 : atom_i = ex_bond_list(iatom)%array1(i)
502 130 : DO j = 1, SIZE(ex_bend_list(iatom)%array1)
503 : ! a next neighboring atom of iatom:
504 72 : atom_j = ex_bend_list(iatom)%array1(j)
505 : ! It is only a true onfo if there is no shorter path.
506 : ! check if i and j are not the same atom
507 72 : IF (atom_i == atom_j) CYCLE
508 : ! check if i and j are not involved in a bond
509 72 : check = .FALSE.
510 230 : DO counter = 1, SIZE(ex_bond_list(atom_i)%array1)
511 230 : IF (ex_bond_list(atom_i)%array1(counter) == atom_j) THEN
512 : check = .TRUE.
513 : EXIT
514 : END IF
515 : END DO
516 72 : IF (check) CYCLE
517 : ! check if i and j are not involved in a bend
518 4 : check = .FALSE.
519 8 : DO counter = 1, SIZE(ex_bend_list(atom_i)%array1)
520 8 : IF (ex_bend_list(atom_i)%array1(counter) == atom_j) THEN
521 : check = .TRUE.
522 : EXIT
523 : END IF
524 : END DO
525 4 : IF (check) CYCLE
526 : ! Add the true onfo.
527 4 : N = N + 1
528 4 : IF (SIZE(pairs, dim=1) <= N) THEN
529 2 : CALL reallocate(pairs, 1, N + 5, 1, 2)
530 : END IF
531 4 : pairs(N, 1) = atom_i
532 108 : pairs(N, 2) = atom_j
533 : END DO
534 : END DO
535 : END DO
536 4 : CALL reorder_structure(ex_onfo_list, pairs(:, 1), pairs(:, 2), N)
537 4 : DEALLOCATE (pairs)
538 : ELSE
539 2482 : IF (ASSOCIATED(conn_info%onfo_a)) THEN
540 2476 : N = SIZE(conn_info%onfo_a)
541 2476 : CALL reorder_structure(ex_onfo_list, conn_info%onfo_a, conn_info%onfo_b, N)
542 : END IF
543 : END IF
544 :
545 : ! Build the exclusion (and onfo) list per atom.
546 635376 : DO iatom = 1, SIZE(particle_set)
547 : ! Setup exclusion list for VDW: always exclude itself
548 632890 : dim0 = 1
549 : ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
550 632890 : dim1 = 0
551 : IF (topology%exclude_vdw == do_skip_12 .OR. &
552 632890 : topology%exclude_vdw == do_skip_13 .OR. &
553 632296 : topology%exclude_vdw == do_skip_14) dim1 = SIZE(ex_bond_list_vdw(iatom)%array1)
554 632890 : dim1 = dim1 + dim0
555 632890 : dim2 = 0
556 632890 : IF (topology%exclude_vdw == do_skip_13 .OR. &
557 631690 : topology%exclude_vdw == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
558 632890 : dim2 = dim1 + dim2
559 632890 : dim3 = 0
560 632890 : IF (topology%exclude_vdw == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
561 632890 : dim3 = dim2 + dim3
562 632890 : IF (dim3 /= 0) THEN
563 632890 : NULLIFY (list, wlist)
564 1898670 : ALLOCATE (wlist(dim3))
565 1265780 : wlist(dim0:dim0) = iatom
566 1604576 : IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_vdw(iatom)%array1
567 1156480 : IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
568 634438 : IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
569 : ! Get a unique list
570 2129714 : DO i = 1, SIZE(wlist) - 1
571 1496824 : IF (wlist(i) == 0) CYCLE
572 5431388 : DO j = i + 1, SIZE(wlist)
573 4798568 : IF (wlist(j) == wlist(i)) wlist(j) = 0
574 : END DO
575 : END DO
576 2762604 : dim3 = SIZE(wlist) - COUNT(wlist == 0)
577 1898670 : ALLOCATE (list(dim3))
578 632890 : j = 0
579 2762604 : DO i = 1, SIZE(wlist)
580 2129714 : IF (wlist(i) == 0) CYCLE
581 2043830 : j = j + 1
582 2762604 : list(j) = wlist(i)
583 : END DO
584 632890 : DEALLOCATE (wlist)
585 : ! Unique list completed
586 632890 : NULLIFY (list2)
587 : IF ((topology%exclude_vdw == topology%exclude_ei) .AND. &
588 632890 : (.NOT. present_12_excl_ei_list) .AND. (.NOT. present_12_excl_vdw_list)) THEN
589 : list2 => list
590 : ELSE
591 : ! Setup exclusion list for EI : always exclude itself
592 1770 : dim0 = 1
593 : ! exclude bond-neighbors (only if do_skip_12 .OR. do_skip_13 .OR. do_skip_14)
594 1770 : dim1 = 0
595 : IF (topology%exclude_ei == do_skip_12 .OR. &
596 1770 : topology%exclude_ei == do_skip_13 .OR. &
597 1326 : topology%exclude_ei == do_skip_14) dim1 = SIZE(ex_bond_list_ei(iatom)%array1)
598 1770 : dim1 = dim1 + dim0
599 1770 : dim2 = 0
600 1770 : IF (topology%exclude_ei == do_skip_13 .OR. &
601 864 : topology%exclude_ei == do_skip_14) dim2 = SIZE(ex_bend_list(iatom)%array1)
602 1770 : dim2 = dim1 + dim2
603 1770 : dim3 = 0
604 1770 : IF (topology%exclude_ei == do_skip_14) dim3 = SIZE(ex_onfo_list(iatom)%array1)
605 1770 : dim3 = dim2 + dim3
606 :
607 1770 : IF (dim3 /= 0) THEN
608 5310 : ALLOCATE (wlist(dim3))
609 3540 : wlist(dim0:dim0) = iatom
610 4098 : IF (dim1 > dim0) wlist(dim0 + 1:dim1) = ex_bond_list_ei(iatom)%array1
611 4266 : IF (dim2 > dim1) wlist(dim1 + 1:dim2) = ex_bend_list(iatom)%array1
612 2922 : IF (dim3 > dim2) wlist(dim2 + 1:dim3) = ex_onfo_list(iatom)%array1
613 : ! Get a unique list
614 7746 : DO i = 1, SIZE(wlist) - 1
615 5976 : IF (wlist(i) == 0) CYCLE
616 28896 : DO j = i + 1, SIZE(wlist)
617 27126 : IF (wlist(j) == wlist(i)) wlist(j) = 0
618 : END DO
619 : END DO
620 9516 : dim3 = SIZE(wlist) - COUNT(wlist == 0)
621 5310 : ALLOCATE (list2(dim3))
622 1770 : j = 0
623 9516 : DO i = 1, SIZE(wlist)
624 7746 : IF (wlist(i) == 0) CYCLE
625 7746 : j = j + 1
626 9516 : list2(j) = wlist(i)
627 : END DO
628 1770 : DEALLOCATE (wlist)
629 : ! Unique list completed
630 : END IF
631 : END IF
632 : END IF
633 632890 : exclusions(iatom)%list_exclude_vdw => list
634 632890 : exclusions(iatom)%list_exclude_ei => list2
635 : ! Keep a list of onfo atoms for proper selection of specialized 1-4
636 : ! potentials instead of conventional nonbonding potentials.
637 1340895 : ALLOCATE (exclusions(iatom)%list_onfo(SIZE(ex_onfo_list(iatom)%array1)))
638 : ! copy of data, not copy of pointer
639 1253926 : exclusions(iatom)%list_onfo = ex_onfo_list(iatom)%array1
640 635376 : IF (iw > 0) THEN
641 1244 : IF (ASSOCIATED(list)) THEN
642 1244 : WRITE (iw, *) "exclusion list_vdw :: ", &
643 1244 : "atom num :", iatom, "exclusion list ::", &
644 6006 : list
645 : END IF
646 1244 : IF (topology%exclude_vdw /= topology%exclude_ei) THEN
647 9 : IF (ASSOCIATED(list2)) THEN
648 9 : WRITE (iw, *) "exclusion list_ei :: ", &
649 9 : "atom num :", iatom, "exclusion list ::", &
650 35 : list2
651 : END IF
652 : END IF
653 1244 : IF (ASSOCIATED(exclusions(iatom)%list_onfo)) THEN
654 1244 : WRITE (iw, *) "onfo list :: ", &
655 1244 : "atom num :", iatom, "onfo list ::", &
656 3322 : exclusions(iatom)%list_onfo
657 : END IF
658 : END IF
659 : END DO
660 : ! deallocate onfo
661 635376 : DO I = 1, natom
662 635376 : DEALLOCATE (ex_onfo_list(I)%array1)
663 : END DO
664 2486 : DEALLOCATE (ex_onfo_list)
665 : ! deallocate bends
666 635376 : DO I = 1, natom
667 635376 : DEALLOCATE (ex_bend_list(I)%array1)
668 : END DO
669 2486 : DEALLOCATE (ex_bend_list)
670 : ! deallocate bonds
671 2486 : IF (present_12_excl_ei_list) THEN
672 40 : DO I = 1, natom
673 40 : DEALLOCATE (ex_bond_list_ei(I)%array1)
674 : END DO
675 10 : DEALLOCATE (ex_bond_list_ei)
676 : ELSE
677 : NULLIFY (ex_bond_list_ei)
678 : END IF
679 2486 : IF (present_12_excl_vdw_list) THEN
680 32 : DO I = 1, natom
681 32 : DEALLOCATE (ex_bond_list_vdw(I)%array1)
682 : END DO
683 8 : DEALLOCATE (ex_bond_list_vdw)
684 : ELSE
685 : NULLIFY (ex_bond_list_vdw)
686 : END IF
687 635376 : DO I = 1, natom
688 635376 : DEALLOCATE (ex_bond_list(I)%array1)
689 : END DO
690 9944 : DEALLOCATE (ex_bond_list)
691 : END IF
692 11500 : CALL timestop(handle2)
693 : !-----------------------------------------------------------------------------
694 : !-----------------------------------------------------------------------------
695 : ! 11. Set the atomic_kind_set()%fist_potential%[qeff] (PART 1)
696 : !-----------------------------------------------------------------------------
697 11500 : CALL timeset(routineN//'_10', handle2)
698 11500 : CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
699 11500 : IF (method_name_id == do_fist) THEN
700 13906 : DO i = 1, SIZE(atomic_kind_set)
701 11266 : atomic_kind => atomic_kind_set(i)
702 11266 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
703 11266 : qeff = charge(i)
704 11266 : NULLIFY (fist_potential)
705 11266 : CALL allocate_potential(fist_potential)
706 11266 : CALL set_potential(potential=fist_potential, qeff=qeff)
707 13906 : CALL set_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
708 : END DO
709 : END IF
710 11500 : DEALLOCATE (charge)
711 11500 : CALL timestop(handle2)
712 :
713 : !-----------------------------------------------------------------------------
714 : !-----------------------------------------------------------------------------
715 : ! 12. Set the atom_list for molecule_kind in molecule_kind_set (PART 2)
716 : !-----------------------------------------------------------------------------
717 11500 : CALL timeset(routineN//'_11', handle2)
718 158459 : DO i = 1, SIZE(molecule_kind_set)
719 146959 : molecule_kind => molecule_kind_set(i)
720 : CALL get_molecule_kind(molecule_kind=molecule_kind, &
721 : natom=natom, molecule_list=molecule_list, &
722 146959 : atom_list=atom_list)
723 146959 : molecule => molecule_set(molecule_list(1))
724 : CALL get_molecule(molecule=molecule, &
725 146959 : first_atom=first, last_atom=last)
726 407673 : DO j = 1, natom
727 18855127 : DO k = 1, SIZE(atomic_kind_set)
728 18708168 : atomic_kind => atomic_kind_set(k)
729 18708168 : CALL get_atomic_kind(atomic_kind=atomic_kind, name=atmname)
730 18708168 : IF (method_name_id == do_fist) THEN
731 18366637 : CALL get_atomic_kind(atomic_kind=atomic_kind, fist_potential=fist_potential)
732 18366637 : CALL get_potential(potential=fist_potential, qeff=qeff)
733 18366637 : IF ((id2str(atom_list(j)%id_name) == atmname) .AND. (qeff == atom_info%atm_charge(first + j - 1))) THEN
734 183303 : atom_list(j)%atomic_kind => atomic_kind_set(k)
735 18549940 : EXIT
736 : END IF
737 : ELSE
738 341531 : IF (id2str(atom_list(j)%id_name) == atmname) THEN
739 77411 : atom_list(j)%atomic_kind => atomic_kind_set(k)
740 418942 : EXIT
741 : END IF
742 : END IF
743 : END DO
744 : END DO
745 305418 : CALL set_molecule_kind(molecule_kind=molecule_kind, atom_list=atom_list)
746 : END DO
747 11500 : CALL timestop(handle2)
748 :
749 11500 : CALL timestop(handle)
750 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
751 11500 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
752 138000 : END SUBROUTINE topology_coordinate_pack
753 :
754 : ! **************************************************************************************************
755 : !> \brief Builds the exclusion list for VDW and EI if an explicit list of terms
756 : !> is provided by the user. Otherwise all possibilities are excluded
757 : !> \param exclude_section ...
758 : !> \param keyword ...
759 : !> \param ex_bond_list ...
760 : !> \param ex_bond_list_w ...
761 : !> \param particle_set ...
762 : !> \par History
763 : !> Teodoro Laino [tlaino] - 12.2009
764 : ! **************************************************************************************************
765 18 : SUBROUTINE setup_exclusion_list(exclude_section, keyword, ex_bond_list, &
766 : ex_bond_list_w, particle_set)
767 : TYPE(section_vals_type), POINTER :: exclude_section
768 : CHARACTER(LEN=*), INTENT(IN) :: keyword
769 : TYPE(array1_list_type), DIMENSION(:), POINTER :: ex_bond_list, ex_bond_list_w
770 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
771 :
772 : CHARACTER(LEN=default_string_length) :: flag1, flag2
773 : CHARACTER(LEN=default_string_length), &
774 18 : DIMENSION(:), POINTER :: names
775 : INTEGER :: i, ind, j, k, l, m, n_rep
776 :
777 0 : CPASSERT(ASSOCIATED(ex_bond_list))
778 18 : CPASSERT(ASSOCIATED(ex_bond_list_w))
779 18 : SELECT CASE (keyword)
780 : CASE ("BOND")
781 18 : CALL section_vals_val_get(exclude_section, keyword, n_rep_val=n_rep)
782 72 : DO j = 1, SIZE(ex_bond_list)
783 54 : CPASSERT(ASSOCIATED(ex_bond_list(j)%array1))
784 54 : CPASSERT(ASSOCIATED(ex_bond_list_w(j)%array1))
785 :
786 54 : flag1 = particle_set(j)%atomic_kind%name
787 54 : m = SIZE(ex_bond_list(j)%array1)
788 54 : CALL reallocate(ex_bond_list_w(j)%array1, 1, m)
789 :
790 54 : l = 0
791 126 : DO k = 1, m
792 72 : ind = ex_bond_list(j)%array1(k)
793 72 : flag2 = particle_set(ind)%atomic_kind%name
794 150 : DO i = 1, n_rep
795 : CALL section_vals_val_get(exclude_section, keyword, i_rep_val=i, &
796 24 : c_vals=names)
797 24 : IF (((TRIM(names(1)) == TRIM(flag1)) .AND. (TRIM(names(2)) == TRIM(flag2))) .OR. &
798 72 : ((TRIM(names(1)) == TRIM(flag2)) .AND. (TRIM(names(2)) == TRIM(flag1)))) THEN
799 24 : l = l + 1
800 24 : ex_bond_list_w(j)%array1(l) = ind
801 : END IF
802 : END DO
803 : END DO
804 72 : CALL reallocate(ex_bond_list_w(j)%array1, 1, l)
805 : END DO
806 : CASE DEFAULT
807 : CALL cp_abort(__LOCATION__, &
808 : "<BOND> is supported as the <keyword> for "// &
809 : "setup_exclusion_list, found unknown option "// &
810 18 : "<"//TRIM(keyword)//">")
811 : END SELECT
812 :
813 18 : END SUBROUTINE setup_exclusion_list
814 :
815 : END MODULE topology_coordinate_util
|
