Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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_util
14 : USE cp_log_handling, ONLY: cp_get_default_logger,&
15 : cp_logger_get_default_io_unit,&
16 : cp_logger_type,&
17 : cp_to_string
18 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
19 : cp_print_key_unit_nr
20 : USE graphcon, ONLY: graph_type,&
21 : hash_molecule,&
22 : reorder_graph,&
23 : vertex
24 : USE input_section_types, ONLY: section_vals_get,&
25 : section_vals_get_subs_vals,&
26 : section_vals_type,&
27 : section_vals_val_get,&
28 : section_vals_val_set
29 : USE kinds, ONLY: default_string_length,&
30 : dp
31 : USE mm_mapping_library, ONLY: amber_map,&
32 : charmm_map,&
33 : gromos_map
34 : USE periodic_table, ONLY: get_ptable_info
35 : USE qmmm_types_low, ONLY: qmmm_env_mm_type
36 : USE string_table, ONLY: id2str,&
37 : str2id
38 : USE string_utilities, ONLY: uppercase
39 : USE topology_types, ONLY: atom_info_type,&
40 : connectivity_info_type,&
41 : topology_parameters_type
42 : USE util, ONLY: sort
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_util'
48 :
49 : TYPE array1_list_type
50 : INTEGER, POINTER, DIMENSION(:) :: array1 => NULL()
51 : END TYPE array1_list_type
52 :
53 : TYPE array2_list_type
54 : INTEGER, POINTER, DIMENSION(:) :: array1 => NULL()
55 : INTEGER, POINTER, DIMENSION(:) :: array2 => NULL()
56 : END TYPE array2_list_type
57 :
58 : PRIVATE
59 : PUBLIC :: topology_set_atm_mass, &
60 : topology_reorder_atoms, &
61 : topology_molecules_check, &
62 : check_subsys_element, &
63 : reorder_structure, &
64 : find_molecule, &
65 : array1_list_type, &
66 : array2_list_type, &
67 : give_back_molecule, &
68 : reorder_list_array, &
69 : tag_molecule
70 :
71 : INTERFACE reorder_structure
72 : MODULE PROCEDURE reorder_structure1d, reorder_structure2d
73 : END INTERFACE
74 :
75 : CONTAINS
76 :
77 : ! **************************************************************************************************
78 : !> \brief ...
79 : !> \param topology ...
80 : !> \param qmmm ...
81 : !> \param qmmm_env_mm ...
82 : !> \param subsys_section ...
83 : !> \param force_env_section ...
84 : !> \par History
85 : !> Teodoro Laino 09.2006 - Rewritten with a graph matching algorithm
86 : ! **************************************************************************************************
87 314 : SUBROUTINE topology_reorder_atoms(topology, qmmm, qmmm_env_mm, subsys_section, force_env_section)
88 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
89 : LOGICAL, INTENT(in), OPTIONAL :: qmmm
90 : TYPE(qmmm_env_mm_type), OPTIONAL, POINTER :: qmmm_env_mm
91 : TYPE(section_vals_type), POINTER :: subsys_section, force_env_section
92 :
93 : CHARACTER(len=*), PARAMETER :: routineN = 'topology_reorder_atoms'
94 :
95 : CHARACTER(LEN=default_string_length) :: mol_id
96 314 : CHARACTER(LEN=default_string_length), POINTER :: molname(:), telement(:), &
97 314 : tlabel_atmname(:), tlabel_molname(:), &
98 314 : tlabel_resname(:)
99 : INTEGER :: handle, i, iatm, iindex, ikind, imol, imol_ref, iref, iund, iw, j, k, location, &
100 : max_mol_num, mm_index, n, n_rep, n_var, natom, natom_loc, nkind, nlinks, old_hash, &
101 : old_mol, output_unit, qm_index, unique_mol
102 314 : INTEGER, DIMENSION(:), POINTER :: mm_link_atoms, qm_atom_index
103 314 : INTEGER, POINTER :: atm_map1(:), atm_map2(:), atm_map3(:), map_atom_type(:), &
104 314 : map_mol_hash(:), mm_indexes_n(:), mm_indexes_v(:), mol_bnd(:, :), mol_hash(:), &
105 314 : mol_num(:), new_position(:), order(:), tmp_n(:), tmp_v(:), wrk(:)
106 : LOGICAL :: explicit, matches, my_qmmm
107 314 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: tr
108 314 : REAL(KIND=dp), POINTER :: tatm_charge(:), tatm_mass(:)
109 314 : TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list
110 : TYPE(atom_info_type), POINTER :: atom_info
111 : TYPE(connectivity_info_type), POINTER :: conn_info
112 : TYPE(cp_logger_type), POINTER :: logger
113 314 : TYPE(graph_type), DIMENSION(:), POINTER :: reference_set
114 : TYPE(section_vals_type), POINTER :: generate_section, isolated_section, &
115 : qm_kinds, qmmm_link_section, &
116 : qmmm_section
117 314 : TYPE(vertex), DIMENSION(:), POINTER :: reference, unordered
118 :
119 314 : NULLIFY (logger, generate_section, isolated_section, tmp_v, tmp_n)
120 628 : logger => cp_get_default_logger()
121 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
122 314 : extension=".subsysLog")
123 314 : CALL timeset(routineN, handle)
124 314 : output_unit = cp_logger_get_default_io_unit(logger)
125 314 : IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER | ")')
126 :
127 314 : my_qmmm = .FALSE.
128 314 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env_mm)) my_qmmm = qmmm
129 :
130 314 : atom_info => topology%atom_info
131 314 : conn_info => topology%conn_info
132 314 : natom = topology%natoms
133 :
134 314 : NULLIFY (new_position, reference_set)
135 314 : NULLIFY (tlabel_atmname, telement, mol_num, tlabel_molname, tlabel_resname)
136 314 : NULLIFY (tr, tatm_charge, tatm_mass, atm_map1, atm_map2, atm_map3)
137 : ! This routine can be called only at a very high level where these structures are still
138 : ! not even taken into account...
139 314 : CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_num))
140 314 : CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_typ))
141 314 : CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_res))
142 : !ALLOCATE all the temporary arrays needed for this routine
143 942 : ALLOCATE (new_position(natom))
144 628 : ALLOCATE (mol_num(natom))
145 942 : ALLOCATE (molname(natom))
146 628 : ALLOCATE (tlabel_atmname(natom))
147 628 : ALLOCATE (tlabel_molname(natom))
148 628 : ALLOCATE (tlabel_resname(natom))
149 942 : ALLOCATE (tr(3, natom))
150 942 : ALLOCATE (tatm_charge(natom))
151 628 : ALLOCATE (tatm_mass(natom))
152 628 : ALLOCATE (telement(natom))
153 628 : ALLOCATE (atm_map1(natom))
154 628 : ALLOCATE (atm_map2(natom))
155 :
156 : ! The only information we have at this level is the connectivity of the system.
157 : ! 0. Build a list of mapping atom types
158 628 : ALLOCATE (map_atom_type(natom))
159 : ! 1. Build a list of bonds
160 9500 : ALLOCATE (atom_bond_list(natom))
161 8872 : DO I = 1, natom
162 8558 : map_atom_type(I) = atom_info%id_atmname(i)
163 8872 : ALLOCATE (atom_bond_list(I)%array1(0))
164 : END DO
165 314 : N = 0
166 314 : IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
167 314 : CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
168 628 : ALLOCATE (atom_info%map_mol_num(natom))
169 8872 : atom_info%map_mol_num = -1
170 314 : CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
171 8872 : max_mol_num = MAXVAL(atom_info%map_mol_num)
172 : ! In atom_info%map_mol_num have already been mapped molecules
173 942 : ALLOCATE (mol_bnd(2, max_mol_num))
174 942 : ALLOCATE (mol_hash(max_mol_num))
175 628 : ALLOCATE (map_mol_hash(max_mol_num))
176 : ! 2. Sort the map_mol_num array.. atm_map1 contains now the mapped index
177 : ! of the reordered array
178 314 : CALL sort(atom_info%map_mol_num, natom, atm_map1)
179 314 : old_mol = 0
180 314 : iindex = 0
181 314 : imol = 0
182 8872 : DO i = 1, natom
183 8558 : IF (old_mol /= atom_info%map_mol_num(I)) THEN
184 1420 : old_mol = atom_info%map_mol_num(I)
185 1420 : iindex = 0
186 1420 : IF (imol > 0) THEN
187 1106 : mol_bnd(2, imol) = i - 1
188 : END IF
189 1420 : imol = imol + 1
190 1420 : mol_bnd(1, imol) = i
191 : END IF
192 8558 : iindex = iindex + 1
193 8872 : atm_map2(atm_map1(i)) = iindex
194 : END DO
195 314 : mol_bnd(2, imol) = natom
196 : ! Indexes of the two molecules to check
197 314 : iref = 1
198 314 : iund = max_mol_num/2 + 1
199 : ! Allocate reference and unordered
200 314 : NULLIFY (reference, unordered)
201 : ! This is the real matching of graphs
202 1734 : DO j = 1, max_mol_num
203 : CALL setup_graph(j, reference, map_atom_type, &
204 1420 : atom_bond_list, mol_bnd, atm_map1, atm_map2)
205 :
206 4260 : ALLOCATE (order(SIZE(reference)))
207 1420 : CALL hash_molecule(reference, order, mol_hash(j))
208 :
209 1420 : DEALLOCATE (order)
210 9978 : DO I = 1, SIZE(reference)
211 9978 : DEALLOCATE (reference(I)%bonds)
212 : END DO
213 1734 : DEALLOCATE (reference)
214 : END DO
215 : ! Reorder molecules hashes
216 314 : CALL sort(mol_hash, max_mol_num, map_mol_hash)
217 : ! Now find unique molecules and reorder atoms too (if molecules match)
218 314 : old_hash = -1
219 314 : unique_mol = 0
220 314 : natom_loc = 0
221 314 : IF (output_unit > 0) THEN
222 : WRITE (output_unit, '(T2,"REORDER | ",A)') &
223 157 : "Reordering Molecules. The Reordering of molecules will override all", &
224 157 : "information regarding molecule names and residue names.", &
225 314 : "New ones will be provided on the basis of the connectivity!"
226 : END IF
227 1734 : DO j = 1, max_mol_num
228 1734 : IF (mol_hash(j) /= old_hash) THEN
229 326 : unique_mol = unique_mol + 1
230 326 : old_hash = mol_hash(j)
231 : CALL setup_graph_set(reference_set, unique_mol, map_mol_hash(j), &
232 : map_atom_type, atom_bond_list, mol_bnd, atm_map1, atm_map2, &
233 326 : atm_map3)
234 : ! Reorder Last added reference
235 326 : mol_id = TRIM(ADJUSTL(cp_to_string(unique_mol)))
236 5636 : DO i = 1, SIZE(atm_map3)
237 5310 : natom_loc = natom_loc + 1
238 5310 : new_position(natom_loc) = atm_map3(i)
239 5310 : molname(natom_loc) = mol_id
240 5636 : mol_num(natom_loc) = unique_mol
241 : END DO
242 326 : DEALLOCATE (atm_map3)
243 : ELSE
244 : CALL setup_graph(map_mol_hash(j), unordered, map_atom_type, &
245 1094 : atom_bond_list, mol_bnd, atm_map1, atm_map2, atm_map3)
246 1150 : DO imol_ref = 1, unique_mol
247 : !
248 1150 : reference => reference_set(imol_ref)%graph
249 3450 : ALLOCATE (order(SIZE(reference)))
250 1150 : CALL reorder_graph(reference, unordered, order, matches)
251 1150 : IF (matches) EXIT
252 2300 : DEALLOCATE (order)
253 : END DO
254 1094 : IF (matches) THEN
255 : ! Reorder according the correct index
256 3282 : ALLOCATE (wrk(SIZE(order)))
257 1094 : CALL sort(order, SIZE(order), wrk)
258 4342 : DO i = 1, SIZE(order)
259 3248 : natom_loc = natom_loc + 1
260 3248 : new_position(natom_loc) = atm_map3(wrk(i))
261 3248 : molname(natom_loc) = mol_id
262 4342 : mol_num(natom_loc) = unique_mol
263 : END DO
264 : !
265 1094 : DEALLOCATE (order)
266 1094 : DEALLOCATE (wrk)
267 : ELSE
268 0 : unique_mol = unique_mol + 1
269 : CALL setup_graph_set(reference_set, unique_mol, map_mol_hash(j), &
270 : map_atom_type, atom_bond_list, mol_bnd, atm_map1, atm_map2, &
271 0 : atm_map3)
272 : ! Reorder Last added reference
273 0 : mol_id = TRIM(ADJUSTL(cp_to_string(unique_mol)))
274 0 : DO i = 1, SIZE(atm_map3)
275 0 : natom_loc = natom_loc + 1
276 0 : new_position(natom_loc) = atm_map3(i)
277 0 : molname(natom_loc) = mol_id
278 0 : mol_num(natom_loc) = unique_mol
279 : END DO
280 0 : DEALLOCATE (atm_map3)
281 : END IF
282 4342 : DO I = 1, SIZE(unordered)
283 4342 : DEALLOCATE (unordered(I)%bonds)
284 : END DO
285 1094 : DEALLOCATE (unordered)
286 1094 : DEALLOCATE (atm_map3)
287 : END IF
288 : END DO
289 314 : IF (output_unit > 0) THEN
290 157 : WRITE (output_unit, '(T2,"REORDER | ",A,I7,A)') "Number of unique molecules found:", unique_mol, "."
291 : END IF
292 314 : CPASSERT(natom_loc == natom)
293 314 : DEALLOCATE (map_atom_type)
294 314 : DEALLOCATE (atm_map1)
295 314 : DEALLOCATE (atm_map2)
296 314 : DEALLOCATE (mol_bnd)
297 314 : DEALLOCATE (mol_hash)
298 314 : DEALLOCATE (map_mol_hash)
299 : ! Deallocate working arrays
300 8872 : DO I = 1, natom
301 8872 : DEALLOCATE (atom_bond_list(I)%array1)
302 : END DO
303 314 : DEALLOCATE (atom_bond_list)
304 314 : DEALLOCATE (atom_info%map_mol_num)
305 : ! Deallocate reference_set
306 640 : DO i = 1, SIZE(reference_set)
307 5636 : DO j = 1, SIZE(reference_set(i)%graph)
308 5636 : DEALLOCATE (reference_set(i)%graph(j)%bonds)
309 : END DO
310 640 : DEALLOCATE (reference_set(i)%graph)
311 : END DO
312 314 : DEALLOCATE (reference_set)
313 : !Lets swap the atoms now
314 8872 : DO iatm = 1, natom
315 8558 : location = new_position(iatm)
316 8558 : tlabel_molname(iatm) = id2str(atom_info%id_molname(location))
317 8558 : tlabel_resname(iatm) = id2str(atom_info%id_resname(location))
318 8558 : tlabel_atmname(iatm) = id2str(atom_info%id_atmname(location))
319 8558 : telement(iatm) = id2str(atom_info%id_element(location))
320 8558 : tr(1, iatm) = atom_info%r(1, location)
321 8558 : tr(2, iatm) = atom_info%r(2, location)
322 8558 : tr(3, iatm) = atom_info%r(3, location)
323 8558 : tatm_charge(iatm) = atom_info%atm_charge(location)
324 8872 : tatm_mass(iatm) = atom_info%atm_mass(location)
325 : END DO
326 314 : IF (topology%create_molecules) THEN
327 8858 : DO iatm = 1, natom
328 8546 : tlabel_molname(iatm) = "MOL"//TRIM(molname(iatm))
329 8858 : tlabel_resname(iatm) = "R"//TRIM(molname(iatm))
330 : END DO
331 312 : topology%create_molecules = .FALSE.
332 : END IF
333 8872 : DO iatm = 1, natom
334 8558 : atom_info%id_molname(iatm) = str2id(tlabel_molname(iatm))
335 8558 : atom_info%id_resname(iatm) = str2id(tlabel_resname(iatm))
336 8558 : atom_info%id_atmname(iatm) = str2id(tlabel_atmname(iatm))
337 8558 : atom_info%id_element(iatm) = str2id(telement(iatm))
338 8558 : atom_info%resid(iatm) = mol_num(iatm)
339 8558 : atom_info%r(1, iatm) = tr(1, iatm)
340 8558 : atom_info%r(2, iatm) = tr(2, iatm)
341 8558 : atom_info%r(3, iatm) = tr(3, iatm)
342 8558 : atom_info%atm_charge(iatm) = tatm_charge(iatm)
343 8872 : atom_info%atm_mass(iatm) = tatm_mass(iatm)
344 : END DO
345 :
346 : ! Let's reorder all the list provided in the input..
347 942 : ALLOCATE (wrk(SIZE(new_position)))
348 314 : CALL sort(new_position, SIZE(new_position), wrk)
349 :
350 : ! NOTE: In the future the whole list of possible integers should be updated at this level..
351 : ! Let's update all the integer lists in the qmmm_env_mm section and in the input sections
352 : ! from where qmmm_env_qm will be read.
353 314 : IF (my_qmmm) THEN
354 : ! Update the qmmm_env_mm
355 2 : CPASSERT(SIZE(qmmm_env_mm%qm_atom_index) /= 0)
356 8 : CPASSERT(ALL(qmmm_env_mm%qm_atom_index /= 0))
357 6 : ALLOCATE (qm_atom_index(SIZE(qmmm_env_mm%qm_atom_index)))
358 8 : DO iatm = 1, SIZE(qmmm_env_mm%qm_atom_index)
359 8 : qm_atom_index(iatm) = wrk(qmmm_env_mm%qm_atom_index(iatm))
360 : END DO
361 16 : qmmm_env_mm%qm_atom_index = qm_atom_index
362 8 : CPASSERT(ALL(qmmm_env_mm%qm_atom_index /= 0))
363 2 : DEALLOCATE (qm_atom_index)
364 : ! Update the qmmm_section: MM_INDEX of the QM_KIND
365 2 : qmmm_section => section_vals_get_subs_vals(force_env_section, "QMMM")
366 2 : qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
367 2 : CALL section_vals_get(qm_kinds, n_repetition=nkind)
368 6 : DO ikind = 1, nkind
369 4 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
370 10 : DO k = 1, n_var
371 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
372 4 : i_vals=mm_indexes_v)
373 12 : ALLOCATE (mm_indexes_n(SIZE(mm_indexes_v)))
374 8 : DO j = 1, SIZE(mm_indexes_v)
375 8 : mm_indexes_n(j) = wrk(mm_indexes_v(j))
376 : END DO
377 : CALL section_vals_val_set(qm_kinds, "MM_INDEX", i_rep_section=ikind, &
378 8 : i_vals_ptr=mm_indexes_n, i_rep_val=k)
379 : END DO
380 : END DO
381 : ! Handle the link atoms
382 4 : IF (qmmm_env_mm%qmmm_link) THEN
383 : ! Update the qmmm_env_mm
384 2 : CPASSERT(SIZE(qmmm_env_mm%mm_link_atoms) > 0)
385 6 : ALLOCATE (mm_link_atoms(SIZE(qmmm_env_mm%mm_link_atoms)))
386 4 : DO iatm = 1, SIZE(qmmm_env_mm%mm_link_atoms)
387 4 : mm_link_atoms(iatm) = wrk(qmmm_env_mm%mm_link_atoms(iatm))
388 : END DO
389 8 : qmmm_env_mm%mm_link_atoms = mm_link_atoms
390 4 : CPASSERT(ALL(qmmm_env_mm%mm_link_atoms /= 0))
391 2 : DEALLOCATE (mm_link_atoms)
392 : ! Update the qmmm_section: MM_INDEX,QM_INDEX of the LINK atom list
393 2 : qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
394 2 : CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
395 2 : CPASSERT(nlinks /= 0)
396 6 : DO ikind = 1, nlinks
397 2 : CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
398 2 : CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
399 2 : mm_index = wrk(mm_index)
400 2 : qm_index = wrk(qm_index)
401 2 : CALL section_vals_val_set(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
402 6 : CALL section_vals_val_set(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
403 : END DO
404 : END IF
405 : END IF
406 : !
407 : !LIST of ISOLATED atoms
408 : !
409 314 : generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE")
410 314 : isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS")
411 314 : CALL section_vals_get(isolated_section, explicit=explicit)
412 314 : IF (explicit) THEN
413 4 : CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep)
414 10 : DO i = 1, n_rep
415 6 : CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i)
416 18 : ALLOCATE (tmp_n(SIZE(tmp_v)))
417 50 : DO j = 1, SIZE(tmp_v)
418 50 : tmp_n(j) = wrk(tmp_v(j))
419 : END DO
420 10 : CALL section_vals_val_set(isolated_section, "LIST", i_vals_ptr=tmp_n, i_rep_val=i)
421 : END DO
422 : END IF
423 314 : DEALLOCATE (wrk)
424 : !DEALLOCATE all the temporary arrays needed for this routine
425 314 : DEALLOCATE (new_position)
426 314 : DEALLOCATE (tlabel_atmname)
427 314 : DEALLOCATE (tlabel_molname)
428 314 : DEALLOCATE (tlabel_resname)
429 314 : DEALLOCATE (telement)
430 314 : DEALLOCATE (tr)
431 314 : DEALLOCATE (tatm_charge)
432 314 : DEALLOCATE (tatm_mass)
433 314 : DEALLOCATE (molname)
434 314 : DEALLOCATE (mol_num)
435 :
436 : ! DEALLOCATE the bond structures in the connectivity
437 314 : DEALLOCATE (conn_info%bond_a)
438 314 : DEALLOCATE (conn_info%bond_b)
439 314 : IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER | ")')
440 314 : CALL timestop(handle)
441 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
442 314 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
443 1256 : END SUBROUTINE topology_reorder_atoms
444 :
445 : ! **************************************************************************************************
446 : !> \brief Set up a SET of graph kind
447 : !> \param graph_set ...
448 : !> \param idim ...
449 : !> \param ind ...
450 : !> \param array2 ...
451 : !> \param atom_bond_list ...
452 : !> \param map_mol ...
453 : !> \param atm_map1 ...
454 : !> \param atm_map2 ...
455 : !> \param atm_map3 ...
456 : !> \author Teodoro Laino 10.2006
457 : ! **************************************************************************************************
458 326 : SUBROUTINE setup_graph_set(graph_set, idim, ind, array2, atom_bond_list, map_mol, &
459 : atm_map1, atm_map2, atm_map3)
460 : TYPE(graph_type), DIMENSION(:), POINTER :: graph_set
461 : INTEGER, INTENT(IN) :: idim, ind
462 : INTEGER, DIMENSION(:), POINTER :: array2
463 : TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
464 : INTEGER, DIMENSION(:, :), POINTER :: map_mol
465 : INTEGER, DIMENSION(:), POINTER :: atm_map1, atm_map2, atm_map3
466 :
467 : INTEGER :: ldim
468 326 : TYPE(graph_type), DIMENSION(:), POINTER :: tmp_graph_set
469 :
470 326 : ldim = 0
471 326 : NULLIFY (tmp_graph_set)
472 326 : IF (ASSOCIATED(graph_set)) THEN
473 12 : ldim = SIZE(graph_set)
474 12 : CPASSERT(ldim + 1 == idim)
475 : MARK_USED(idim)
476 : NULLIFY (tmp_graph_set)
477 12 : CALL allocate_graph_set(graph_set, tmp_graph_set)
478 : END IF
479 326 : CALL allocate_graph_set(tmp_graph_set, graph_set, ldim, ldim + 1)
480 : CALL setup_graph(ind, graph_set(ldim + 1)%graph, array2, &
481 326 : atom_bond_list, map_mol, atm_map1, atm_map2, atm_map3)
482 :
483 326 : END SUBROUTINE setup_graph_set
484 :
485 : ! **************************************************************************************************
486 : !> \brief Allocate a new graph_set deallocating an old one..
487 : !> \param graph_set_in ...
488 : !> \param graph_set_out ...
489 : !> \param ldim_in ...
490 : !> \param ldim_out ...
491 : !> \author Teodoro Laino 10.2006
492 : ! **************************************************************************************************
493 338 : SUBROUTINE allocate_graph_set(graph_set_in, graph_set_out, ldim_in, ldim_out)
494 : TYPE(graph_type), DIMENSION(:), POINTER :: graph_set_in, graph_set_out
495 : INTEGER, INTENT(IN), OPTIONAL :: ldim_in, ldim_out
496 :
497 : INTEGER :: b_dim, i, j, mydim_in, mydim_out, v_dim
498 :
499 338 : CPASSERT(.NOT. ASSOCIATED(graph_set_out))
500 338 : mydim_in = 0
501 338 : mydim_out = 0
502 338 : IF (ASSOCIATED(graph_set_in)) THEN
503 24 : mydim_in = SIZE(graph_set_in)
504 24 : mydim_out = SIZE(graph_set_in)
505 : END IF
506 338 : IF (PRESENT(ldim_in)) mydim_in = ldim_in
507 338 : IF (PRESENT(ldim_out)) mydim_out = ldim_out
508 1372 : ALLOCATE (graph_set_out(mydim_out))
509 696 : DO i = 1, mydim_out
510 696 : NULLIFY (graph_set_out(i)%graph)
511 : END DO
512 : ! Copy graph structure into the temporary array
513 370 : DO i = 1, mydim_in
514 32 : v_dim = SIZE(graph_set_in(i)%graph)
515 332 : ALLOCATE (graph_set_out(i)%graph(v_dim))
516 268 : DO j = 1, v_dim
517 236 : graph_set_out(i)%graph(j)%kind = graph_set_in(i)%graph(j)%kind
518 236 : b_dim = SIZE(graph_set_in(i)%graph(j)%bonds)
519 700 : ALLOCATE (graph_set_out(i)%graph(j)%bonds(b_dim))
520 1304 : graph_set_out(i)%graph(j)%bonds = graph_set_in(i)%graph(j)%bonds
521 268 : DEALLOCATE (graph_set_in(i)%graph(j)%bonds)
522 : END DO
523 370 : DEALLOCATE (graph_set_in(i)%graph)
524 : END DO
525 338 : IF (ASSOCIATED(graph_set_in)) THEN
526 24 : DEALLOCATE (graph_set_in)
527 : END IF
528 :
529 338 : END SUBROUTINE allocate_graph_set
530 :
531 : ! **************************************************************************************************
532 : !> \brief Set up a graph kind
533 : !> \param ind ...
534 : !> \param graph ...
535 : !> \param array2 ...
536 : !> \param atom_bond_list ...
537 : !> \param map_mol ...
538 : !> \param atm_map1 ...
539 : !> \param atm_map2 ...
540 : !> \param atm_map3 ...
541 : !> \author Teodoro Laino 09.2006
542 : ! **************************************************************************************************
543 2840 : SUBROUTINE setup_graph(ind, graph, array2, atom_bond_list, map_mol, &
544 : atm_map1, atm_map2, atm_map3)
545 : INTEGER, INTENT(IN) :: ind
546 : TYPE(vertex), DIMENSION(:), POINTER :: graph
547 : INTEGER, DIMENSION(:), POINTER :: array2
548 : TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
549 : INTEGER, DIMENSION(:, :), POINTER :: map_mol
550 : INTEGER, DIMENSION(:), POINTER :: atm_map1, atm_map2
551 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: atm_map3
552 :
553 : INTEGER :: i, idim, ifirst, ilast, j, nbonds, &
554 : nelement
555 :
556 2840 : IF (PRESENT(atm_map3)) THEN
557 1420 : CPASSERT(.NOT. ASSOCIATED(atm_map3))
558 : END IF
559 2840 : CPASSERT(.NOT. ASSOCIATED(graph))
560 : ! Setup reference graph
561 2840 : idim = 0
562 2840 : ifirst = map_mol(1, ind)
563 2840 : ilast = map_mol(2, ind)
564 2840 : nelement = ilast - ifirst + 1
565 25636 : ALLOCATE (graph(nelement))
566 2840 : IF (PRESENT(atm_map3)) THEN
567 4260 : ALLOCATE (atm_map3(nelement))
568 : END IF
569 19956 : DO i = ifirst, ilast
570 17116 : idim = idim + 1
571 17116 : graph(idim)%kind = array2(atm_map1(i))
572 17116 : nbonds = SIZE(atom_bond_list(atm_map1(i))%array1)
573 51260 : ALLOCATE (graph(idim)%bonds(nbonds))
574 49220 : DO j = 1, nbonds
575 49220 : graph(idim)%bonds(j) = atm_map2(atom_bond_list(atm_map1(i))%array1(j))
576 : END DO
577 19956 : IF (PRESENT(atm_map3)) THEN
578 8558 : atm_map3(idim) = atm_map1(i)
579 : END IF
580 : END DO
581 :
582 2840 : END SUBROUTINE setup_graph
583 :
584 : ! **************************************************************************************************
585 : !> \brief Order arrays of lists
586 : !> \param Ilist1 ...
587 : !> \param Ilist2 ...
588 : !> \param Ilist3 ...
589 : !> \param Ilist4 ...
590 : !> \param nsize ...
591 : !> \param ndim ...
592 : !> \author Teodoro Laino 09.2006
593 : ! **************************************************************************************************
594 261436 : RECURSIVE SUBROUTINE reorder_list_array(Ilist1, Ilist2, Ilist3, Ilist4, nsize, ndim)
595 : INTEGER, DIMENSION(:), POINTER :: Ilist1
596 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Ilist2, Ilist3, Ilist4
597 : INTEGER, INTENT(IN) :: nsize, ndim
598 :
599 : INTEGER :: i, iend, istart, ldim
600 261436 : INTEGER, DIMENSION(:), POINTER :: tmp, wrk
601 :
602 0 : CPASSERT(nsize > 0)
603 759414 : ALLOCATE (wrk(Ndim))
604 261436 : CALL sort(Ilist1, Ndim, wrk)
605 261436 : IF (nsize /= 1) THEN
606 212277 : ALLOCATE (tmp(Ndim))
607 1054990 : tmp = Ilist2(1:Ndim)
608 527495 : DO i = 1, Ndim
609 527495 : Ilist2(i) = tmp(wrk(i))
610 : END DO
611 19724 : SELECT CASE (nsize)
612 : CASE (3)
613 293172 : tmp = Ilist3(1:Ndim)
614 146586 : DO i = 1, Ndim
615 146586 : Ilist3(i) = tmp(wrk(i))
616 : END DO
617 : CASE (4)
618 96364 : tmp = Ilist3(1:Ndim)
619 48182 : DO i = 1, Ndim
620 48182 : Ilist3(i) = tmp(wrk(i))
621 : END DO
622 96364 : tmp = Ilist4(1:Ndim)
623 162149 : DO i = 1, Ndim
624 48182 : Ilist4(i) = tmp(wrk(i))
625 : END DO
626 : END SELECT
627 113967 : DEALLOCATE (tmp)
628 113967 : istart = 1
629 527495 : DO i = 1, Ndim
630 527495 : IF (Ilist1(i) /= Ilist1(istart)) THEN
631 133780 : iend = i - 1
632 133780 : ldim = iend - istart + 1
633 : CALL reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
634 133780 : ldim, istart, iend)
635 133780 : istart = i
636 : END IF
637 : END DO
638 : ! Last term to sort
639 113967 : iend = Ndim
640 113967 : ldim = iend - istart + 1
641 : CALL reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
642 113967 : ldim, istart, iend)
643 : END IF
644 261436 : DEALLOCATE (wrk)
645 261436 : END SUBROUTINE reorder_list_array
646 :
647 : ! **************************************************************************************************
648 : !> \brief Low level routine for ordering arrays of lists
649 : !> \param Ilist2 ...
650 : !> \param Ilist3 ...
651 : !> \param Ilist4 ...
652 : !> \param nsize ...
653 : !> \param ldim ...
654 : !> \param istart ...
655 : !> \param iend ...
656 : !> \author Teodoro Laino 09.2006
657 : ! **************************************************************************************************
658 247747 : RECURSIVE SUBROUTINE reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
659 : ldim, istart, iend)
660 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Ilist2, Ilist3, Ilist4
661 : INTEGER, INTENT(IN) :: nsize, ldim, istart, iend
662 :
663 247747 : INTEGER, DIMENSION(:), POINTER :: tmp_2, tmp_3, tmp_4
664 :
665 395216 : SELECT CASE (nsize)
666 : CASE (2)
667 433170 : ALLOCATE (tmp_2(ldim))
668 779298 : tmp_2(:) = Ilist2(istart:iend)
669 147469 : CALL reorder_list_array(tmp_2, nsize=nsize - 1, ndim=ldim)
670 779298 : Ilist2(istart:iend) = tmp_2(:)
671 147469 : DEALLOCATE (tmp_2)
672 : CASE (3)
673 243952 : ALLOCATE (tmp_2(ldim))
674 240698 : ALLOCATE (tmp_3(ldim))
675 418528 : tmp_2(:) = Ilist2(istart:iend)
676 497676 : tmp_3(:) = Ilist3(istart:iend)
677 82402 : CALL reorder_list_array(tmp_2, tmp_3, nsize=nsize - 1, ndim=ldim)
678 418528 : Ilist2(istart:iend) = tmp_2(:)
679 418528 : Ilist3(istart:iend) = tmp_3(:)
680 82402 : DEALLOCATE (tmp_2)
681 82402 : DEALLOCATE (tmp_3)
682 : CASE (4)
683 50462 : ALLOCATE (tmp_2(ldim))
684 47296 : ALLOCATE (tmp_3(ldim))
685 47296 : ALLOCATE (tmp_4(ldim))
686 124724 : tmp_2(:) = Ilist2(istart:iend)
687 139434 : tmp_3(:) = Ilist3(istart:iend)
688 139434 : tmp_4(:) = Ilist4(istart:iend)
689 17876 : CALL reorder_list_array(tmp_2, tmp_3, tmp_4, nsize=nsize - 1, ndim=ldim)
690 124724 : Ilist2(istart:iend) = tmp_2(:)
691 124724 : Ilist3(istart:iend) = tmp_3(:)
692 124724 : Ilist4(istart:iend) = tmp_4(:)
693 17876 : DEALLOCATE (tmp_2)
694 17876 : DEALLOCATE (tmp_3)
695 17876 : DEALLOCATE (tmp_4)
696 : END SELECT
697 :
698 247747 : END SUBROUTINE reorder_list_array_low
699 :
700 : ! **************************************************************************************************
701 : !> \brief ...
702 : !> \param icheck ...
703 : !> \param bond_list ...
704 : !> \param i ...
705 : !> \param mol_natom ...
706 : !> \param mol_map ...
707 : !> \param my_mol ...
708 : !> \author Teodoro Laino 09.2006
709 : ! **************************************************************************************************
710 209112 : RECURSIVE SUBROUTINE give_back_molecule(icheck, bond_list, i, mol_natom, mol_map, my_mol)
711 : LOGICAL, DIMENSION(:), POINTER :: icheck
712 : TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
713 : INTEGER, INTENT(IN) :: i
714 : INTEGER, INTENT(INOUT) :: mol_natom
715 : INTEGER, DIMENSION(:), POINTER :: mol_map
716 : INTEGER, INTENT(IN) :: my_mol
717 :
718 : INTEGER :: j, k
719 :
720 209112 : IF (mol_map(i) == my_mol) THEN
721 209104 : icheck(i) = .TRUE.
722 423728 : DO j = 1, SIZE(bond_list(i)%array1)
723 214624 : k = bond_list(i)%array1(j)
724 214624 : IF (icheck(k)) CYCLE
725 106312 : mol_natom = mol_natom + 1
726 423728 : CALL give_back_molecule(icheck, bond_list, k, mol_natom, mol_map, my_mol)
727 : END DO
728 : ELSE
729 : ! Do nothing means only that bonds were found between molecules
730 : ! as we defined them.. This could easily be a bond detected but not
731 : ! physically present.. so just skip this part and go on..
732 : END IF
733 209112 : END SUBROUTINE give_back_molecule
734 :
735 : ! **************************************************************************************************
736 : !> \brief gives back a mapping of molecules.. icheck needs to be initialized with -1
737 : !> \param icheck ...
738 : !> \param bond_list ...
739 : !> \param i ...
740 : !> \param my_mol ...
741 : !> \author Teodoro Laino 04.2007 - Zurich University
742 : ! **************************************************************************************************
743 6706 : RECURSIVE SUBROUTINE tag_molecule(icheck, bond_list, i, my_mol)
744 : INTEGER, DIMENSION(:), POINTER :: icheck
745 : TYPE(array1_list_type), DIMENSION(:), POINTER :: bond_list
746 : INTEGER, INTENT(IN) :: i, my_mol
747 :
748 : INTEGER :: j, k
749 :
750 6706 : icheck(i) = my_mol
751 17898 : DO j = 1, SIZE(bond_list(i)%array1)
752 11192 : k = bond_list(i)%array1(j)
753 11192 : IF (k <= i) CYCLE
754 17898 : CALL tag_molecule(icheck, bond_list, k, my_mol)
755 : END DO
756 :
757 6706 : END SUBROUTINE tag_molecule
758 :
759 : ! **************************************************************************************************
760 : !> \brief Given two lists of corresponding element a complex type is built in
761 : !> which each element of the type has a list of mapping elements
762 : !> \param work ...
763 : !> \param list1 ...
764 : !> \param list2 ...
765 : !> \param N ...
766 : !> \author Teodoro Laino 08.2006
767 : ! **************************************************************************************************
768 71993 : SUBROUTINE reorder_structure1d(work, list1, list2, N)
769 : TYPE(array1_list_type), DIMENSION(:), &
770 : INTENT(INOUT) :: work
771 : INTEGER, DIMENSION(:), INTENT(IN) :: list1, list2
772 : INTEGER, INTENT(IN) :: N
773 :
774 : INTEGER :: I, index1, index2, Nsize
775 71993 : INTEGER, DIMENSION(:), POINTER :: wrk_tmp
776 :
777 3374455 : DO I = 1, N
778 3302462 : index1 = list1(I)
779 3302462 : index2 = list2(I)
780 :
781 3302462 : wrk_tmp => work(index1)%array1
782 3302462 : Nsize = SIZE(wrk_tmp)
783 9907386 : ALLOCATE (work(index1)%array1(Nsize + 1))
784 6455025 : work(index1)%array1(1:Nsize) = wrk_tmp
785 3302462 : work(index1)%array1(Nsize + 1) = index2
786 3302462 : DEALLOCATE (wrk_tmp)
787 :
788 3302462 : wrk_tmp => work(index2)%array1
789 3302462 : Nsize = SIZE(wrk_tmp)
790 9907386 : ALLOCATE (work(index2)%array1(Nsize + 1))
791 4519961 : work(index2)%array1(1:Nsize) = wrk_tmp
792 3302462 : work(index2)%array1(Nsize + 1) = index1
793 3374455 : DEALLOCATE (wrk_tmp)
794 : END DO
795 :
796 71993 : END SUBROUTINE reorder_structure1d
797 :
798 : ! **************************************************************************************************
799 : !> \brief Given two lists of corresponding element a complex type is built in
800 : !> which each element of the type has a list of mapping elements
801 : !> \param work ...
802 : !> \param list1 ...
803 : !> \param list2 ...
804 : !> \param list3 ...
805 : !> \param N ...
806 : !> \author Teodoro Laino 09.2006
807 : ! **************************************************************************************************
808 8143 : SUBROUTINE reorder_structure2d(work, list1, list2, list3, N)
809 : TYPE(array2_list_type), DIMENSION(:), &
810 : INTENT(INOUT) :: work
811 : INTEGER, DIMENSION(:), INTENT(IN) :: list1, list2, list3
812 : INTEGER, INTENT(IN) :: N
813 :
814 : INTEGER :: I, index1, index2, index3, Nsize
815 8143 : INTEGER, DIMENSION(:), POINTER :: wrk_tmp
816 :
817 1147844 : DO I = 1, N
818 1139701 : index1 = list1(I)
819 1139701 : index2 = list2(I)
820 1139701 : index3 = list3(I)
821 :
822 1139701 : wrk_tmp => work(index1)%array1
823 1139701 : Nsize = SIZE(wrk_tmp)
824 3419103 : ALLOCATE (work(index1)%array1(Nsize + 1))
825 23874175 : work(index1)%array1(1:Nsize) = wrk_tmp
826 1139701 : work(index1)%array1(Nsize + 1) = index2
827 1139701 : DEALLOCATE (wrk_tmp)
828 :
829 1139701 : wrk_tmp => work(index2)%array1
830 1139701 : Nsize = SIZE(wrk_tmp)
831 3419103 : ALLOCATE (work(index2)%array1(Nsize + 1))
832 15755001 : work(index2)%array1(1:Nsize) = wrk_tmp
833 1139701 : work(index2)%array1(Nsize + 1) = index1
834 1139701 : DEALLOCATE (wrk_tmp)
835 :
836 1139701 : wrk_tmp => work(index1)%array2
837 1139701 : Nsize = SIZE(wrk_tmp)
838 3419103 : ALLOCATE (work(index1)%array2(Nsize + 1))
839 23874175 : work(index1)%array2(1:Nsize) = wrk_tmp
840 1139701 : work(index1)%array2(Nsize + 1) = index3
841 1139701 : DEALLOCATE (wrk_tmp)
842 :
843 1139701 : wrk_tmp => work(index2)%array2
844 1139701 : Nsize = SIZE(wrk_tmp)
845 3419103 : ALLOCATE (work(index2)%array2(Nsize + 1))
846 15755001 : work(index2)%array2(1:Nsize) = wrk_tmp
847 1139701 : work(index2)%array2(Nsize + 1) = -index3
848 1147844 : DEALLOCATE (wrk_tmp)
849 : END DO
850 :
851 8143 : END SUBROUTINE reorder_structure2d
852 :
853 : ! **************************************************************************************************
854 : !> \brief each atom will be assigned a molecule number based on bonded fragments
855 : !> The array mol_info should be initialized with -1 before calling the
856 : !> find_molecule routine
857 : !> \param atom_bond_list ...
858 : !> \param mol_info ...
859 : !> \param mol_name ...
860 : !> \author Joost 05.2006
861 : ! **************************************************************************************************
862 10588 : SUBROUTINE find_molecule(atom_bond_list, mol_info, mol_name)
863 : TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
864 : INTEGER, DIMENSION(:), POINTER :: mol_info, mol_name
865 :
866 : INTEGER :: I, my_mol_name, N, nmol
867 :
868 10588 : N = SIZE(atom_bond_list)
869 10588 : nmol = 0
870 773715 : DO I = 1, N
871 773715 : IF (mol_info(I) == -1) THEN
872 330784 : nmol = nmol + 1
873 330784 : my_mol_name = mol_name(I)
874 : CALL spread_mol(atom_bond_list, mol_info, i, nmol, my_mol_name, &
875 330784 : mol_name)
876 : END IF
877 : END DO
878 10588 : END SUBROUTINE find_molecule
879 :
880 : ! **************************************************************************************************
881 : !> \brief spreads the molnumber over the bonded list
882 : !> \param atom_bond_list ...
883 : !> \param mol_info ...
884 : !> \param iatom ...
885 : !> \param imol ...
886 : !> \param my_mol_name ...
887 : !> \param mol_name ...
888 : !> \author Joost 05.2006
889 : ! **************************************************************************************************
890 763127 : RECURSIVE SUBROUTINE spread_mol(atom_bond_list, mol_info, iatom, imol, &
891 : my_mol_name, mol_name)
892 : TYPE(array1_list_type), DIMENSION(:), INTENT(IN) :: atom_bond_list
893 : INTEGER, DIMENSION(:), POINTER :: mol_info
894 : INTEGER, INTENT(IN) :: iatom, imol, my_mol_name
895 : INTEGER, DIMENSION(:), POINTER :: mol_name
896 :
897 : INTEGER :: atom_b, i
898 :
899 763127 : mol_info(iatom) = imol
900 1795793 : DO I = 1, SIZE(atom_bond_list(iatom)%array1)
901 1032666 : atom_b = atom_bond_list(iatom)%array1(I)
902 : ! In this way we're really sure that all atoms belong to the same
903 : ! molecule. This should correct possible errors in the generation of
904 : ! the bond list..
905 1795793 : IF (mol_info(atom_b) == -1 .AND. my_mol_name == mol_name(atom_b)) THEN
906 432343 : CALL spread_mol(atom_bond_list, mol_info, atom_b, imol, my_mol_name, mol_name)
907 432343 : IF (mol_info(atom_b) /= imol) CPABORT("internal error")
908 : END IF
909 : END DO
910 763127 : END SUBROUTINE spread_mol
911 :
912 : ! **************************************************************************************************
913 : !> \brief Use info from periodic table and set atm_mass
914 : !> \param topology ...
915 : !> \param subsys_section ...
916 : ! **************************************************************************************************
917 9791 : SUBROUTINE topology_set_atm_mass(topology, subsys_section)
918 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
919 : TYPE(section_vals_type), POINTER :: subsys_section
920 :
921 : CHARACTER(len=*), PARAMETER :: routineN = 'topology_set_atm_mass'
922 :
923 : CHARACTER(LEN=2) :: upper_sym_1
924 : CHARACTER(LEN=default_string_length) :: atmname_upper
925 : CHARACTER(LEN=default_string_length), &
926 9791 : DIMENSION(:), POINTER :: keyword
927 : INTEGER :: handle, i, i_rep, iatom, ielem_found, &
928 : iw, n_rep, natom
929 : LOGICAL :: user_defined
930 9791 : REAL(KIND=dp), DIMENSION(:), POINTER :: mass
931 : TYPE(atom_info_type), POINTER :: atom_info
932 : TYPE(cp_logger_type), POINTER :: logger
933 : TYPE(section_vals_type), POINTER :: kind_section
934 :
935 9791 : NULLIFY (logger)
936 19582 : logger => cp_get_default_logger()
937 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
938 9791 : extension=".subsysLog")
939 9791 : CALL timeset(routineN, handle)
940 :
941 9791 : atom_info => topology%atom_info
942 9791 : natom = topology%natoms
943 :
944 : ! Available external info
945 9791 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
946 9791 : CALL section_vals_get(kind_section, n_repetition=n_rep)
947 25657 : ALLOCATE (keyword(n_rep))
948 25657 : ALLOCATE (mass(n_rep))
949 21998 : mass = HUGE(0.0_dp)
950 21998 : DO i_rep = 1, n_rep
951 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
952 12207 : c_val=keyword(i_rep), i_rep_section=i_rep)
953 12207 : CALL uppercase(keyword(i_rep))
954 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
955 12207 : keyword_name="MASS", n_rep_val=i)
956 12207 : IF (i > 0) CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
957 22096 : keyword_name="MASS", r_val=mass(i_rep))
958 : END DO
959 : !
960 288929 : DO iatom = 1, natom
961 : !If we reach this point then we've definitely identified the element..
962 : !Let's look if an external mass has been defined..
963 279138 : user_defined = .FALSE.
964 392136 : DO i = 1, SIZE(keyword)
965 113596 : atmname_upper = id2str(atom_info%id_atmname(iatom))
966 113596 : CALL uppercase(atmname_upper)
967 392136 : IF (TRIM(atmname_upper) == TRIM(keyword(i)) .AND. mass(i) /= HUGE(0.0_dp)) THEN
968 598 : atom_info%atm_mass(iatom) = mass(i)
969 : user_defined = .TRUE.
970 : EXIT
971 : END IF
972 : END DO
973 : ! If name didn't match let's try with the element
974 278540 : IF (.NOT. user_defined) THEN
975 278540 : upper_sym_1 = TRIM(id2str(atom_info%id_element(iatom)))
976 278540 : CALL get_ptable_info(symbol=upper_sym_1, ielement=ielem_found, amass=atom_info%atm_mass(iatom))
977 : END IF
978 280542 : IF (iw > 0) WRITE (iw, '(7X,A,A5,A,F12.5)') "In topology_set_atm_mass :: element = ", &
979 12599 : id2str(atom_info%id_element(iatom)), " a_mass ", atom_info%atm_mass(iatom)
980 : END DO
981 9791 : DEALLOCATE (keyword)
982 9791 : DEALLOCATE (mass)
983 :
984 9791 : CALL timestop(handle)
985 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
986 9791 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
987 :
988 29373 : END SUBROUTINE topology_set_atm_mass
989 :
990 : ! **************************************************************************************************
991 : !> \brief Check and verify that all molecules of the same kind are bonded the same
992 : !> \param topology ...
993 : !> \param subsys_section ...
994 : ! **************************************************************************************************
995 9736 : SUBROUTINE topology_molecules_check(topology, subsys_section)
996 : TYPE(topology_parameters_type), INTENT(INOUT) :: topology
997 : TYPE(section_vals_type), POINTER :: subsys_section
998 :
999 : CHARACTER(len=*), PARAMETER :: routineN = 'topology_molecules_check'
1000 :
1001 : INTEGER :: counter, first, first_loc, handle, i, &
1002 : iatom, iw, k, loc_counter, mol_num, &
1003 : mol_typ, n, natom
1004 : LOGICAL :: icheck_num, icheck_typ
1005 9736 : TYPE(array1_list_type), ALLOCATABLE, DIMENSION(:) :: atom_bond_list
1006 : TYPE(atom_info_type), POINTER :: atom_info
1007 : TYPE(connectivity_info_type), POINTER :: conn_info
1008 : TYPE(cp_logger_type), POINTER :: logger
1009 :
1010 9736 : CALL timeset(routineN, handle)
1011 :
1012 9736 : NULLIFY (logger)
1013 9736 : logger => cp_get_default_logger()
1014 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
1015 9736 : extension=".subsysLog")
1016 :
1017 9736 : atom_info => topology%atom_info
1018 9736 : conn_info => topology%conn_info
1019 9736 : natom = topology%natoms
1020 :
1021 9736 : IF (iw > 0) THEN
1022 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
1023 26 : "FORCEFIELD| Checking consistency of the generated molecules"
1024 : END IF
1025 :
1026 780065 : ALLOCATE (atom_bond_list(natom))
1027 :
1028 760593 : DO I = 1, natom
1029 760593 : ALLOCATE (atom_bond_list(I)%array1(0))
1030 : END DO
1031 9736 : N = 0
1032 9736 : IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
1033 9736 : CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
1034 :
1035 9736 : mol_typ = atom_info%map_mol_typ(1)
1036 9736 : mol_num = atom_info%map_mol_num(1)
1037 9736 : counter = 1
1038 9736 : loc_counter = 1
1039 9736 : first = 1
1040 9736 : first_loc = 1
1041 750857 : DO iatom = 2, natom
1042 741121 : icheck_num = (atom_info%map_mol_num(iatom) == mol_num)
1043 741121 : icheck_typ = (atom_info%map_mol_typ(iatom) == mol_typ)
1044 741121 : IF ((icheck_typ .AND. (.NOT. icheck_num)) .OR. (.NOT. icheck_typ)) THEN
1045 : !-----------------------------------------------------------------------------
1046 : ! 1. Check each molecule have the same number of atoms
1047 : !-----------------------------------------------------------------------------
1048 288992 : IF (counter /= loc_counter) THEN
1049 : CALL cp_abort(__LOCATION__, &
1050 : "Different number of atoms for the same molecule kind found "// &
1051 : "(molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// &
1052 : ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// &
1053 : ", expected number of atoms = "// &
1054 : TRIM(ADJUSTL(cp_to_string(counter)))// &
1055 : ", number of atoms found = "// &
1056 0 : TRIM(ADJUSTL(cp_to_string(loc_counter)))//").")
1057 : END IF
1058 : END IF
1059 741121 : IF (.NOT. icheck_typ) THEN
1060 119861 : first = iatom
1061 119861 : first_loc = iatom
1062 119861 : counter = 1
1063 119861 : loc_counter = 1
1064 119861 : mol_typ = atom_info%map_mol_typ(iatom)
1065 : END IF
1066 741121 : IF (icheck_num) THEN
1067 571322 : IF (icheck_typ) loc_counter = loc_counter + 1
1068 : !-----------------------------------------------------------------------------
1069 : ! 2. Check that each molecule has the same atom sequences
1070 : !-----------------------------------------------------------------------------
1071 571322 : IF (atom_info%id_atmname(iatom) /= &
1072 : atom_info%id_atmname(first + loc_counter - 1)) THEN
1073 : CALL cp_abort(__LOCATION__, &
1074 : "Different atom names for the same molecule kind found "// &
1075 : "(atom number = "//TRIM(ADJUSTL(cp_to_string(iatom)))// &
1076 : ", molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// &
1077 : ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// &
1078 : ", expected atom name = "// &
1079 : TRIM(ADJUSTL(id2str(atom_info%id_atmname(first + loc_counter - 1))))// &
1080 : ", atom name found = "// &
1081 0 : TRIM(ADJUSTL(id2str(atom_info%id_atmname(iatom))))//").")
1082 : END IF
1083 : !-----------------------------------------------------------------------------
1084 : ! 3. Check that each molecule have the same bond sequences
1085 : !-----------------------------------------------------------------------------
1086 571322 : IF (SIZE(atom_bond_list(iatom)%array1) /= SIZE(atom_bond_list(first + loc_counter - 1)%array1)) THEN
1087 : CALL cp_abort(__LOCATION__, &
1088 : "Different number of bonds for the same molecule kind found "// &
1089 : "(molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// &
1090 : ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// &
1091 : ", expected number of bonds = "// &
1092 : TRIM(ADJUSTL(cp_to_string(SIZE(atom_bond_list(first + loc_counter - 1)%array1))))// &
1093 : ", number of bonds found = "// &
1094 : TRIM(ADJUSTL(cp_to_string(SIZE(atom_bond_list(iatom)%array1))))// &
1095 0 : "). Check the connectivity of your system.")
1096 : END IF
1097 1279072 : DO k = 1, SIZE(atom_bond_list(iatom)%array1)
1098 1023103 : IF (ALL(atom_bond_list(first + loc_counter - 1)%array1 - first /= &
1099 571322 : atom_bond_list(iatom)%array1(k) - first_loc)) THEN
1100 : CALL cp_abort(__LOCATION__, &
1101 : "Different sequence of bonds for the same molecule kind found "// &
1102 : "(molecule type = "//TRIM(ADJUSTL(cp_to_string(mol_typ)))// &
1103 : ", molecule number = "//TRIM(ADJUSTL(cp_to_string(mol_num)))// &
1104 0 : "). Check the connectivity of your system.")
1105 : END IF
1106 : END DO
1107 : ELSE
1108 169799 : mol_num = atom_info%map_mol_num(iatom)
1109 169799 : loc_counter = 1
1110 169799 : first_loc = iatom
1111 : END IF
1112 750857 : IF (mol_num == 1 .AND. icheck_typ) counter = counter + 1
1113 : END DO
1114 :
1115 9736 : IF (iw > 0) THEN
1116 : WRITE (UNIT=iw, FMT="(T2,A)") &
1117 26 : "FORCEFIELD| Consistency check for molecules completed"
1118 : END IF
1119 :
1120 760593 : DO I = 1, natom
1121 760593 : DEALLOCATE (atom_bond_list(I)%array1)
1122 : END DO
1123 9736 : DEALLOCATE (atom_bond_list)
1124 :
1125 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1126 9736 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
1127 :
1128 9736 : CALL timestop(handle)
1129 :
1130 19472 : END SUBROUTINE topology_molecules_check
1131 :
1132 : ! **************************************************************************************************
1133 : !> \brief Check and returns the ELEMENT label
1134 : !> \param element_in ...
1135 : !> \param atom_name_in ...
1136 : !> \param element_out ...
1137 : !> \param subsys_section ...
1138 : !> \param use_mm_map_first ...
1139 : !> \par History
1140 : !> 12.2005 created [teo]
1141 : !> \author Teodoro Laino
1142 : ! **************************************************************************************************
1143 54934 : SUBROUTINE check_subsys_element(element_in, atom_name_in, element_out, subsys_section, use_mm_map_first)
1144 : CHARACTER(len=*), INTENT(IN) :: element_in, atom_name_in
1145 : CHARACTER(len=default_string_length), INTENT(OUT) :: element_out
1146 : TYPE(section_vals_type), POINTER :: subsys_section
1147 : LOGICAL :: use_mm_map_first
1148 :
1149 : CHARACTER(len=default_string_length) :: atom_name, element_symbol, keyword
1150 : INTEGER :: i, i_rep, n_rep
1151 : LOGICAL :: defined_kind_section, found, in_ptable
1152 : TYPE(section_vals_type), POINTER :: kind_section
1153 :
1154 27467 : found = .FALSE.
1155 27467 : element_symbol = element_in
1156 27467 : atom_name = atom_name_in
1157 27467 : element_out = ""
1158 27467 : defined_kind_section = .FALSE.
1159 :
1160 : ! First check if a KIND section is overriding the element
1161 : ! definition
1162 27467 : CALL uppercase(atom_name)
1163 27467 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
1164 27467 : CALL section_vals_get(kind_section, n_repetition=n_rep)
1165 81266 : DO i_rep = 1, n_rep
1166 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1167 55045 : c_val=keyword, i_rep_section=i_rep)
1168 55045 : CALL uppercase(keyword)
1169 81266 : IF (TRIM(keyword) == TRIM(atom_name)) THEN
1170 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
1171 14733 : keyword_name="ELEMENT", n_rep_val=i)
1172 14733 : IF (i > 0) THEN
1173 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
1174 1246 : keyword_name="ELEMENT", c_val=element_symbol)
1175 : defined_kind_section = .TRUE.
1176 : EXIT
1177 : ELSE
1178 13487 : element_symbol = element_in
1179 : defined_kind_section = .TRUE.
1180 : END IF
1181 : END IF
1182 : END DO
1183 :
1184 : ! Let's check the validity of the element so far stored..
1185 : ! if we are not having a connectivity file, we must first match against the ptable.
1186 : ! this helps to resolve Ca/CA (calcium and Calpha) or Cn/CN7 (Coppernicum (112) CN) conflicts
1187 : ! so, in the presence of connectivity CA will be 'C', while in the absence of connectivity CA will be 'Ca'
1188 26221 : IF (defined_kind_section .OR. .NOT. use_mm_map_first) THEN
1189 : ! lengths larger than 2 should not match, because 'trailing' characters are ignored.
1190 22987 : in_ptable = .FALSE.
1191 22987 : IF (LEN_TRIM(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
1192 22987 : IF (in_ptable) THEN
1193 22895 : element_out = TRIM(element_symbol)
1194 22895 : found = .TRUE.
1195 : END IF
1196 : END IF
1197 :
1198 : ! This is clearly a user error
1199 27467 : IF (.NOT. found .AND. defined_kind_section) &
1200 : CALL cp_abort(__LOCATION__, "Element <"//TRIM(element_symbol)// &
1201 : "> provided for KIND <"//TRIM(atom_name_in)//"> "// &
1202 : "which cannot be mapped with any standard element label. Please correct your "// &
1203 0 : "input file!")
1204 :
1205 : ! Last chance.. are these atom_kinds of AMBER or CHARMM or GROMOS FF ?
1206 27467 : CALL uppercase(element_symbol)
1207 27467 : IF ((.NOT. found) .AND. (ASSOCIATED(amber_map))) THEN
1208 : ! First we go through the AMBER library
1209 145742 : DO i = 1, SIZE(amber_map%kind)
1210 145742 : IF (element_symbol == amber_map%kind(i)) THEN
1211 : found = .TRUE.
1212 : EXIT
1213 : END IF
1214 : END DO
1215 4572 : IF (found) THEN
1216 4136 : element_out = amber_map%element(i)
1217 : END IF
1218 : END IF
1219 4572 : IF ((.NOT. found) .AND. (ASSOCIATED(charmm_map))) THEN
1220 : ! Then we go through the CHARMM library
1221 39242 : DO i = 1, SIZE(charmm_map%kind)
1222 39242 : IF (element_symbol == charmm_map%kind(i)) THEN
1223 : found = .TRUE.
1224 : EXIT
1225 : END IF
1226 : END DO
1227 436 : IF (found) THEN
1228 196 : element_out = charmm_map%element(i)
1229 : END IF
1230 : END IF
1231 436 : IF ((.NOT. found) .AND. (ASSOCIATED(gromos_map))) THEN
1232 : ! Then we go through the GROMOS library
1233 5520 : DO i = 1, SIZE(gromos_map%kind)
1234 5520 : IF (element_symbol == gromos_map%kind(i)) THEN
1235 : found = .TRUE.
1236 : EXIT
1237 : END IF
1238 : END DO
1239 240 : IF (found) THEN
1240 0 : element_out = gromos_map%element(i)
1241 : END IF
1242 : END IF
1243 :
1244 : ! final check. We has a connectivity, so we first tried to match against the ff_maps, but the element was not there.
1245 : ! Try again the periodic table.
1246 0 : IF (.NOT. found) THEN
1247 240 : in_ptable = .FALSE.
1248 240 : IF (LEN_TRIM(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
1249 240 : IF (in_ptable) THEN
1250 240 : element_out = TRIM(element_symbol)
1251 : found = .TRUE.
1252 : END IF
1253 : END IF
1254 :
1255 : ! If no element is found the job stops here.
1256 0 : IF (.NOT. found) THEN
1257 : CALL cp_abort(__LOCATION__, &
1258 : " Unknown element for KIND <"//TRIM(atom_name_in)//">."// &
1259 : " This problem can be fixed specifying properly elements in PDB"// &
1260 : " or specifying a KIND section or getting in touch with one of"// &
1261 0 : " the developers!")
1262 : END IF
1263 :
1264 27467 : END SUBROUTINE check_subsys_element
1265 :
1266 0 : END MODULE topology_util
|
