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