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