Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Collection of subroutine needed for topology related things
10 : !> \par History
11 : !> jgh (23-05-2004) Last atom of molecule information added
12 : ! **************************************************************************************************
13 : MODULE topology_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 :
77 : CONTAINS
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 314 : 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 314 : CHARACTER(LEN=default_string_length), POINTER :: molname(:), telement(:), &
99 314 : tlabel_atmname(:), tlabel_molname(:), &
100 314 : 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 314 : INTEGER, DIMENSION(:), POINTER :: mm_link_atoms, qm_atom_index
105 314 : INTEGER, POINTER :: atm_map1(:), atm_map2(:), atm_map3(:), map_atom_type(:), &
106 314 : map_mol_hash(:), mm_indexes_n(:), mm_indexes_v(:), mol_bnd(:, :), mol_hash(:), &
107 314 : mol_num(:), new_position(:), order(:), tmp_n(:), tmp_v(:), wrk(:)
108 : LOGICAL :: explicit, matches, my_qmmm
109 314 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: tr
110 314 : REAL(KIND=dp), POINTER :: tatm_charge(:), tatm_mass(:)
111 314 : 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 314 : 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 314 : TYPE(vertex), DIMENSION(:), POINTER :: reference, unordered
120 :
121 314 : NULLIFY (logger, generate_section, isolated_section, tmp_v, tmp_n)
122 628 : logger => cp_get_default_logger()
123 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
124 314 : extension=".subsysLog")
125 314 : CALL timeset(routineN, handle)
126 314 : output_unit = cp_logger_get_default_io_unit(logger)
127 314 : IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER | ")')
128 :
129 314 : my_qmmm = .FALSE.
130 314 : IF (PRESENT(qmmm) .AND. PRESENT(qmmm_env_mm)) my_qmmm = qmmm
131 :
132 314 : atom_info => topology%atom_info
133 314 : conn_info => topology%conn_info
134 314 : natom = topology%natoms
135 :
136 314 : NULLIFY (new_position, reference_set)
137 314 : NULLIFY (tlabel_atmname, telement, mol_num, tlabel_molname, tlabel_resname)
138 314 : 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 314 : CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_num))
142 314 : CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_typ))
143 314 : CPASSERT(.NOT. ASSOCIATED(atom_info%map_mol_res))
144 : !ALLOCATE all the temporary arrays needed for this routine
145 942 : ALLOCATE (new_position(natom))
146 628 : ALLOCATE (mol_num(natom))
147 942 : ALLOCATE (molname(natom))
148 628 : ALLOCATE (tlabel_atmname(natom))
149 628 : ALLOCATE (tlabel_molname(natom))
150 628 : ALLOCATE (tlabel_resname(natom))
151 942 : ALLOCATE (tr(3, natom))
152 942 : ALLOCATE (tatm_charge(natom))
153 628 : ALLOCATE (tatm_mass(natom))
154 628 : ALLOCATE (telement(natom))
155 628 : ALLOCATE (atm_map1(natom))
156 628 : 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 628 : ALLOCATE (map_atom_type(natom))
161 : ! 1. Build a list of bonds
162 9500 : ALLOCATE (atom_bond_list(natom))
163 8872 : DO I = 1, natom
164 8558 : map_atom_type(I) = atom_info%id_atmname(i)
165 8872 : ALLOCATE (atom_bond_list(I)%array1(0))
166 : END DO
167 314 : N = 0
168 314 : IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
169 314 : CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
170 628 : ALLOCATE (atom_info%map_mol_num(natom))
171 8872 : atom_info%map_mol_num = -1
172 314 : CALL find_molecule(atom_bond_list, atom_info%map_mol_num, atom_info%id_molname)
173 8872 : max_mol_num = MAXVAL(atom_info%map_mol_num)
174 : ! In atom_info%map_mol_num have already been mapped molecules
175 942 : ALLOCATE (mol_bnd(2, max_mol_num))
176 942 : ALLOCATE (mol_hash(max_mol_num))
177 628 : 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 314 : CALL sort(atom_info%map_mol_num, natom, atm_map1)
181 314 : old_mol = 0
182 314 : iindex = 0
183 314 : imol = 0
184 8872 : DO i = 1, natom
185 8558 : IF (old_mol .NE. atom_info%map_mol_num(I)) THEN
186 1420 : old_mol = atom_info%map_mol_num(I)
187 1420 : iindex = 0
188 1420 : IF (imol > 0) THEN
189 1106 : mol_bnd(2, imol) = i - 1
190 : END IF
191 1420 : imol = imol + 1
192 1420 : mol_bnd(1, imol) = i
193 : END IF
194 8558 : iindex = iindex + 1
195 8872 : atm_map2(atm_map1(i)) = iindex
196 : END DO
197 314 : mol_bnd(2, imol) = natom
198 : ! Indexes of the two molecules to check
199 314 : iref = 1
200 314 : iund = max_mol_num/2 + 1
201 : ! Allocate reference and unordered
202 314 : NULLIFY (reference, unordered)
203 : ! This is the real matching of graphs
204 1734 : DO j = 1, max_mol_num
205 : CALL setup_graph(j, reference, map_atom_type, &
206 1420 : atom_bond_list, mol_bnd, atm_map1, atm_map2)
207 :
208 4260 : ALLOCATE (order(SIZE(reference)))
209 1420 : CALL hash_molecule(reference, order, mol_hash(j))
210 :
211 1420 : DEALLOCATE (order)
212 9978 : DO I = 1, SIZE(reference)
213 9978 : DEALLOCATE (reference(I)%bonds)
214 : END DO
215 1734 : DEALLOCATE (reference)
216 : END DO
217 : ! Reorder molecules hashes
218 314 : CALL sort(mol_hash, max_mol_num, map_mol_hash)
219 : ! Now find unique molecules and reorder atoms too (if molecules match)
220 314 : old_hash = -1
221 314 : unique_mol = 0
222 314 : natom_loc = 0
223 314 : IF (output_unit > 0) THEN
224 : WRITE (output_unit, '(T2,"REORDER | ",A)') &
225 157 : "Reordering Molecules. The Reordering of molecules will override all", &
226 157 : "information regarding molecule names and residue names.", &
227 314 : "New ones will be provided on the basis of the connectivity!"
228 : END IF
229 1734 : DO j = 1, max_mol_num
230 1734 : IF (mol_hash(j) .NE. old_hash) THEN
231 326 : unique_mol = unique_mol + 1
232 326 : 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 326 : atm_map3)
236 : ! Reorder Last added reference
237 326 : mol_id = TRIM(ADJUSTL(cp_to_string(unique_mol)))
238 5636 : DO i = 1, SIZE(atm_map3)
239 5310 : natom_loc = natom_loc + 1
240 5310 : new_position(natom_loc) = atm_map3(i)
241 5310 : molname(natom_loc) = mol_id
242 5636 : mol_num(natom_loc) = unique_mol
243 : END DO
244 326 : DEALLOCATE (atm_map3)
245 : ELSE
246 : CALL setup_graph(map_mol_hash(j), unordered, map_atom_type, &
247 1094 : atom_bond_list, mol_bnd, atm_map1, atm_map2, atm_map3)
248 1150 : DO imol_ref = 1, unique_mol
249 : !
250 1150 : reference => reference_set(imol_ref)%graph
251 3450 : ALLOCATE (order(SIZE(reference)))
252 1150 : CALL reorder_graph(reference, unordered, order, matches)
253 1150 : IF (matches) EXIT
254 2300 : DEALLOCATE (order)
255 : END DO
256 1094 : IF (matches) THEN
257 : ! Reorder according the correct index
258 3282 : ALLOCATE (wrk(SIZE(order)))
259 1094 : CALL sort(order, SIZE(order), wrk)
260 4342 : DO i = 1, SIZE(order)
261 3248 : natom_loc = natom_loc + 1
262 3248 : new_position(natom_loc) = atm_map3(wrk(i))
263 3248 : molname(natom_loc) = mol_id
264 4342 : mol_num(natom_loc) = unique_mol
265 : END DO
266 : !
267 1094 : DEALLOCATE (order)
268 1094 : DEALLOCATE (wrk)
269 : ELSE
270 0 : 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 0 : atm_map3)
274 : ! Reorder Last added reference
275 0 : mol_id = TRIM(ADJUSTL(cp_to_string(unique_mol)))
276 0 : DO i = 1, SIZE(atm_map3)
277 0 : natom_loc = natom_loc + 1
278 0 : new_position(natom_loc) = atm_map3(i)
279 0 : molname(natom_loc) = mol_id
280 0 : mol_num(natom_loc) = unique_mol
281 : END DO
282 0 : DEALLOCATE (atm_map3)
283 : END IF
284 4342 : DO I = 1, SIZE(unordered)
285 4342 : DEALLOCATE (unordered(I)%bonds)
286 : END DO
287 1094 : DEALLOCATE (unordered)
288 1094 : DEALLOCATE (atm_map3)
289 : END IF
290 : END DO
291 314 : IF (output_unit > 0) THEN
292 157 : WRITE (output_unit, '(T2,"REORDER | ",A,I7,A)') "Number of unique molecules found:", unique_mol, "."
293 : END IF
294 314 : CPASSERT(natom_loc == natom)
295 314 : DEALLOCATE (map_atom_type)
296 314 : DEALLOCATE (atm_map1)
297 314 : DEALLOCATE (atm_map2)
298 314 : DEALLOCATE (mol_bnd)
299 314 : DEALLOCATE (mol_hash)
300 314 : DEALLOCATE (map_mol_hash)
301 : ! Deallocate working arrays
302 8872 : DO I = 1, natom
303 8872 : DEALLOCATE (atom_bond_list(I)%array1)
304 : END DO
305 314 : DEALLOCATE (atom_bond_list)
306 314 : DEALLOCATE (atom_info%map_mol_num)
307 : ! Deallocate reference_set
308 640 : DO i = 1, SIZE(reference_set)
309 5636 : DO j = 1, SIZE(reference_set(i)%graph)
310 5636 : DEALLOCATE (reference_set(i)%graph(j)%bonds)
311 : END DO
312 640 : DEALLOCATE (reference_set(i)%graph)
313 : END DO
314 314 : DEALLOCATE (reference_set)
315 : !Lets swap the atoms now
316 8872 : DO iatm = 1, natom
317 8558 : location = new_position(iatm)
318 8558 : tlabel_molname(iatm) = id2str(atom_info%id_molname(location))
319 8558 : tlabel_resname(iatm) = id2str(atom_info%id_resname(location))
320 8558 : tlabel_atmname(iatm) = id2str(atom_info%id_atmname(location))
321 8558 : telement(iatm) = id2str(atom_info%id_element(location))
322 8558 : tr(1, iatm) = atom_info%r(1, location)
323 8558 : tr(2, iatm) = atom_info%r(2, location)
324 8558 : tr(3, iatm) = atom_info%r(3, location)
325 8558 : tatm_charge(iatm) = atom_info%atm_charge(location)
326 8872 : tatm_mass(iatm) = atom_info%atm_mass(location)
327 : END DO
328 314 : IF (topology%create_molecules) THEN
329 8858 : DO iatm = 1, natom
330 8546 : tlabel_molname(iatm) = "MOL"//TRIM(molname(iatm))
331 8858 : tlabel_resname(iatm) = "R"//TRIM(molname(iatm))
332 : END DO
333 312 : topology%create_molecules = .FALSE.
334 : END IF
335 8872 : DO iatm = 1, natom
336 8558 : atom_info%id_molname(iatm) = str2id(tlabel_molname(iatm))
337 8558 : atom_info%id_resname(iatm) = str2id(tlabel_resname(iatm))
338 8558 : atom_info%id_atmname(iatm) = str2id(tlabel_atmname(iatm))
339 8558 : atom_info%id_element(iatm) = str2id(telement(iatm))
340 8558 : atom_info%resid(iatm) = mol_num(iatm)
341 8558 : atom_info%r(1, iatm) = tr(1, iatm)
342 8558 : atom_info%r(2, iatm) = tr(2, iatm)
343 8558 : atom_info%r(3, iatm) = tr(3, iatm)
344 8558 : atom_info%atm_charge(iatm) = tatm_charge(iatm)
345 8872 : atom_info%atm_mass(iatm) = tatm_mass(iatm)
346 : END DO
347 :
348 : ! Let's reorder all the list provided in the input..
349 942 : ALLOCATE (wrk(SIZE(new_position)))
350 314 : 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 314 : IF (my_qmmm) THEN
356 : ! Update the qmmm_env_mm
357 2 : CPASSERT(SIZE(qmmm_env_mm%qm_atom_index) /= 0)
358 8 : CPASSERT(ALL(qmmm_env_mm%qm_atom_index /= 0))
359 6 : ALLOCATE (qm_atom_index(SIZE(qmmm_env_mm%qm_atom_index)))
360 8 : DO iatm = 1, SIZE(qmmm_env_mm%qm_atom_index)
361 8 : qm_atom_index(iatm) = wrk(qmmm_env_mm%qm_atom_index(iatm))
362 : END DO
363 16 : qmmm_env_mm%qm_atom_index = qm_atom_index
364 8 : CPASSERT(ALL(qmmm_env_mm%qm_atom_index /= 0))
365 2 : DEALLOCATE (qm_atom_index)
366 : ! Update the qmmm_section: MM_INDEX of the QM_KIND
367 2 : qmmm_section => section_vals_get_subs_vals(force_env_section, "QMMM")
368 2 : qm_kinds => section_vals_get_subs_vals(qmmm_section, "QM_KIND")
369 2 : CALL section_vals_get(qm_kinds, n_repetition=nkind)
370 6 : DO ikind = 1, nkind
371 4 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, n_rep_val=n_var)
372 10 : DO k = 1, n_var
373 : CALL section_vals_val_get(qm_kinds, "MM_INDEX", i_rep_section=ikind, i_rep_val=k, &
374 4 : i_vals=mm_indexes_v)
375 12 : ALLOCATE (mm_indexes_n(SIZE(mm_indexes_v)))
376 8 : DO j = 1, SIZE(mm_indexes_v)
377 8 : 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 8 : i_vals_ptr=mm_indexes_n, i_rep_val=k)
381 : END DO
382 : END DO
383 : ! Handle the link atoms
384 4 : IF (qmmm_env_mm%qmmm_link) THEN
385 : ! Update the qmmm_env_mm
386 2 : CPASSERT(SIZE(qmmm_env_mm%mm_link_atoms) > 0)
387 6 : ALLOCATE (mm_link_atoms(SIZE(qmmm_env_mm%mm_link_atoms)))
388 4 : DO iatm = 1, SIZE(qmmm_env_mm%mm_link_atoms)
389 4 : mm_link_atoms(iatm) = wrk(qmmm_env_mm%mm_link_atoms(iatm))
390 : END DO
391 8 : qmmm_env_mm%mm_link_atoms = mm_link_atoms
392 4 : CPASSERT(ALL(qmmm_env_mm%mm_link_atoms /= 0))
393 2 : DEALLOCATE (mm_link_atoms)
394 : ! Update the qmmm_section: MM_INDEX,QM_INDEX of the LINK atom list
395 2 : qmmm_link_section => section_vals_get_subs_vals(qmmm_section, "LINK")
396 2 : CALL section_vals_get(qmmm_link_section, n_repetition=nlinks)
397 2 : CPASSERT(nlinks /= 0)
398 6 : DO ikind = 1, nlinks
399 2 : CALL section_vals_val_get(qmmm_link_section, "QM_INDEX", i_rep_section=ikind, i_val=qm_index)
400 2 : CALL section_vals_val_get(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
401 2 : mm_index = wrk(mm_index)
402 2 : qm_index = wrk(qm_index)
403 2 : CALL section_vals_val_set(qmmm_link_section, "MM_INDEX", i_rep_section=ikind, i_val=mm_index)
404 6 : 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 314 : generate_section => section_vals_get_subs_vals(subsys_section, "TOPOLOGY%GENERATE")
412 314 : isolated_section => section_vals_get_subs_vals(generate_section, "ISOLATED_ATOMS")
413 314 : CALL section_vals_get(isolated_section, explicit=explicit)
414 314 : IF (explicit) THEN
415 4 : CALL section_vals_val_get(isolated_section, "LIST", n_rep_val=n_rep)
416 10 : DO i = 1, n_rep
417 6 : CALL section_vals_val_get(isolated_section, "LIST", i_vals=tmp_v, i_rep_val=i)
418 18 : ALLOCATE (tmp_n(SIZE(tmp_v)))
419 50 : DO j = 1, SIZE(tmp_v)
420 50 : tmp_n(j) = wrk(tmp_v(j))
421 : END DO
422 10 : CALL section_vals_val_set(isolated_section, "LIST", i_vals_ptr=tmp_n, i_rep_val=i)
423 : END DO
424 : END IF
425 314 : DEALLOCATE (wrk)
426 : !DEALLOCATE all the temporary arrays needed for this routine
427 314 : DEALLOCATE (new_position)
428 314 : DEALLOCATE (tlabel_atmname)
429 314 : DEALLOCATE (tlabel_molname)
430 314 : DEALLOCATE (tlabel_resname)
431 314 : DEALLOCATE (telement)
432 314 : DEALLOCATE (tr)
433 314 : DEALLOCATE (tatm_charge)
434 314 : DEALLOCATE (tatm_mass)
435 314 : DEALLOCATE (molname)
436 314 : DEALLOCATE (mol_num)
437 :
438 : ! DEALLOCATE the bond structures in the connectivity
439 314 : DEALLOCATE (conn_info%bond_a)
440 314 : DEALLOCATE (conn_info%bond_b)
441 314 : IF (output_unit > 0) WRITE (output_unit, '(T2,"REORDER | ")')
442 314 : CALL timestop(handle)
443 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
444 314 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
445 1256 : 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 326 : 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 326 : TYPE(graph_type), DIMENSION(:), POINTER :: tmp_graph_set
471 :
472 326 : ldim = 0
473 326 : NULLIFY (tmp_graph_set)
474 326 : IF (ASSOCIATED(graph_set)) THEN
475 12 : ldim = SIZE(graph_set)
476 12 : CPASSERT(ldim + 1 == idim)
477 : MARK_USED(idim)
478 : NULLIFY (tmp_graph_set)
479 12 : CALL allocate_graph_set(graph_set, tmp_graph_set)
480 : END IF
481 326 : CALL allocate_graph_set(tmp_graph_set, graph_set, ldim, ldim + 1)
482 : CALL setup_graph(ind, graph_set(ldim + 1)%graph, array2, &
483 326 : atom_bond_list, map_mol, atm_map1, atm_map2, atm_map3)
484 :
485 326 : 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 338 : 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 338 : CPASSERT(.NOT. ASSOCIATED(graph_set_out))
502 338 : mydim_in = 0
503 338 : mydim_out = 0
504 338 : IF (ASSOCIATED(graph_set_in)) THEN
505 24 : mydim_in = SIZE(graph_set_in)
506 24 : mydim_out = SIZE(graph_set_in)
507 : END IF
508 338 : IF (PRESENT(ldim_in)) mydim_in = ldim_in
509 338 : IF (PRESENT(ldim_out)) mydim_out = ldim_out
510 1372 : ALLOCATE (graph_set_out(mydim_out))
511 696 : DO i = 1, mydim_out
512 696 : NULLIFY (graph_set_out(i)%graph)
513 : END DO
514 : ! Copy graph structure into the temporary array
515 370 : DO i = 1, mydim_in
516 32 : v_dim = SIZE(graph_set_in(i)%graph)
517 332 : ALLOCATE (graph_set_out(i)%graph(v_dim))
518 268 : DO j = 1, v_dim
519 236 : graph_set_out(i)%graph(j)%kind = graph_set_in(i)%graph(j)%kind
520 236 : b_dim = SIZE(graph_set_in(i)%graph(j)%bonds)
521 700 : ALLOCATE (graph_set_out(i)%graph(j)%bonds(b_dim))
522 1304 : graph_set_out(i)%graph(j)%bonds = graph_set_in(i)%graph(j)%bonds
523 268 : DEALLOCATE (graph_set_in(i)%graph(j)%bonds)
524 : END DO
525 370 : DEALLOCATE (graph_set_in(i)%graph)
526 : END DO
527 338 : IF (ASSOCIATED(graph_set_in)) THEN
528 24 : DEALLOCATE (graph_set_in)
529 : END IF
530 :
531 338 : 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 2840 : 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 2840 : IF (PRESENT(atm_map3)) THEN
559 1420 : CPASSERT(.NOT. ASSOCIATED(atm_map3))
560 : END IF
561 2840 : CPASSERT(.NOT. ASSOCIATED(graph))
562 : ! Setup reference graph
563 2840 : idim = 0
564 2840 : ifirst = map_mol(1, ind)
565 2840 : ilast = map_mol(2, ind)
566 2840 : nelement = ilast - ifirst + 1
567 25636 : ALLOCATE (graph(nelement))
568 2840 : IF (PRESENT(atm_map3)) THEN
569 4260 : ALLOCATE (atm_map3(nelement))
570 : END IF
571 19956 : DO i = ifirst, ilast
572 17116 : idim = idim + 1
573 17116 : graph(idim)%kind = array2(atm_map1(i))
574 17116 : nbonds = SIZE(atom_bond_list(atm_map1(i))%array1)
575 51260 : ALLOCATE (graph(idim)%bonds(nbonds))
576 49220 : DO j = 1, nbonds
577 49220 : graph(idim)%bonds(j) = atm_map2(atom_bond_list(atm_map1(i))%array1(j))
578 : END DO
579 19956 : IF (PRESENT(atm_map3)) THEN
580 8558 : atm_map3(idim) = atm_map1(i)
581 : END IF
582 : END DO
583 :
584 2840 : 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 260692 : 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 260692 : INTEGER, DIMENSION(:), POINTER :: tmp, wrk
603 :
604 0 : CPASSERT(nsize > 0)
605 757522 : ALLOCATE (wrk(Ndim))
606 260692 : CALL sort(Ilist1, Ndim, wrk)
607 260692 : IF (nsize /= 1) THEN
608 211671 : ALLOCATE (tmp(Ndim))
609 1053466 : tmp = Ilist2(1:Ndim)
610 526733 : DO i = 1, Ndim
611 526733 : Ilist2(i) = tmp(wrk(i))
612 : END DO
613 19622 : SELECT CASE (nsize)
614 : CASE (3)
615 292760 : tmp = Ilist3(1:Ndim)
616 146380 : DO i = 1, Ndim
617 146380 : Ilist3(i) = tmp(wrk(i))
618 : END DO
619 : CASE (4)
620 96196 : tmp = Ilist3(1:Ndim)
621 48098 : DO i = 1, Ndim
622 48098 : Ilist3(i) = tmp(wrk(i))
623 : END DO
624 96196 : tmp = Ilist4(1:Ndim)
625 161655 : DO i = 1, Ndim
626 48098 : Ilist4(i) = tmp(wrk(i))
627 : END DO
628 : END SELECT
629 113557 : DEALLOCATE (tmp)
630 113557 : istart = 1
631 526733 : DO i = 1, Ndim
632 526733 : IF (Ilist1(i) /= Ilist1(istart)) THEN
633 133632 : iend = i - 1
634 133632 : ldim = iend - istart + 1
635 : CALL reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
636 133632 : ldim, istart, iend)
637 133632 : istart = i
638 : END IF
639 : END DO
640 : ! Last term to sort
641 113557 : iend = Ndim
642 113557 : ldim = iend - istart + 1
643 : CALL reorder_list_array_low(Ilist2, Ilist3, Ilist4, nsize, &
644 113557 : ldim, istart, iend)
645 : END IF
646 260692 : DEALLOCATE (wrk)
647 260692 : 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 247189 : 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 247189 : INTEGER, DIMENSION(:), POINTER :: tmp_2, tmp_3, tmp_4
666 :
667 394324 : SELECT CASE (nsize)
668 : CASE (2)
669 432294 : ALLOCATE (tmp_2(ldim))
670 778198 : tmp_2(:) = Ilist2(istart:iend)
671 147135 : CALL reorder_list_array(tmp_2, nsize=nsize - 1, ndim=ldim)
672 778198 : Ilist2(istart:iend) = tmp_2(:)
673 147135 : DEALLOCATE (tmp_2)
674 : CASE (3)
675 243552 : ALLOCATE (tmp_2(ldim))
676 240342 : ALLOCATE (tmp_3(ldim))
677 418024 : tmp_2(:) = Ilist2(istart:iend)
678 497068 : tmp_3(:) = Ilist3(istart:iend)
679 82254 : CALL reorder_list_array(tmp_2, tmp_3, nsize=nsize - 1, ndim=ldim)
680 418024 : Ilist2(istart:iend) = tmp_2(:)
681 418024 : Ilist3(istart:iend) = tmp_3(:)
682 82254 : DEALLOCATE (tmp_2)
683 82254 : DEALLOCATE (tmp_3)
684 : CASE (4)
685 50278 : ALLOCATE (tmp_2(ldim))
686 47156 : ALLOCATE (tmp_3(ldim))
687 47156 : ALLOCATE (tmp_4(ldim))
688 124508 : tmp_2(:) = Ilist2(istart:iend)
689 139186 : tmp_3(:) = Ilist3(istart:iend)
690 139186 : tmp_4(:) = Ilist4(istart:iend)
691 17800 : CALL reorder_list_array(tmp_2, tmp_3, tmp_4, nsize=nsize - 1, ndim=ldim)
692 124508 : Ilist2(istart:iend) = tmp_2(:)
693 124508 : Ilist3(istart:iend) = tmp_3(:)
694 124508 : Ilist4(istart:iend) = tmp_4(:)
695 17800 : DEALLOCATE (tmp_2)
696 17800 : DEALLOCATE (tmp_3)
697 17800 : DEALLOCATE (tmp_4)
698 : END SELECT
699 :
700 247189 : 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 207718 : 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 207718 : IF (mol_map(i) == my_mol) THEN
723 207710 : icheck(i) = .TRUE.
724 422110 : DO j = 1, SIZE(bond_list(i)%array1)
725 214400 : k = bond_list(i)%array1(j)
726 214400 : IF (icheck(k)) CYCLE
727 106200 : mol_natom = mol_natom + 1
728 422110 : 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 207718 : 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 6706 : 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 6706 : icheck(i) = my_mol
753 17898 : DO j = 1, SIZE(bond_list(i)%array1)
754 11192 : k = bond_list(i)%array1(j)
755 11192 : IF (k <= i) CYCLE
756 17898 : CALL tag_molecule(icheck, bond_list, k, my_mol)
757 : END DO
758 :
759 6706 : 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 71215 : 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 71215 : INTEGER, DIMENSION(:), POINTER :: wrk_tmp
778 :
779 3372765 : DO I = 1, N
780 3301550 : index1 = list1(I)
781 3301550 : index2 = list2(I)
782 :
783 3301550 : wrk_tmp => work(index1)%array1
784 3301550 : Nsize = SIZE(wrk_tmp)
785 9904650 : ALLOCATE (work(index1)%array1(Nsize + 1))
786 6453595 : work(index1)%array1(1:Nsize) = wrk_tmp
787 3301550 : work(index1)%array1(Nsize + 1) = index2
788 3301550 : DEALLOCATE (wrk_tmp)
789 :
790 3301550 : wrk_tmp => work(index2)%array1
791 3301550 : Nsize = SIZE(wrk_tmp)
792 9904650 : ALLOCATE (work(index2)%array1(Nsize + 1))
793 4519031 : work(index2)%array1(1:Nsize) = wrk_tmp
794 3301550 : work(index2)%array1(Nsize + 1) = index1
795 3372765 : DEALLOCATE (wrk_tmp)
796 : END DO
797 :
798 71215 : 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 8035 : 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 8035 : INTEGER, DIMENSION(:), POINTER :: wrk_tmp
818 :
819 1141082 : DO I = 1, N
820 1133047 : index1 = list1(I)
821 1133047 : index2 = list2(I)
822 1133047 : index3 = list3(I)
823 :
824 1133047 : wrk_tmp => work(index1)%array1
825 1133047 : Nsize = SIZE(wrk_tmp)
826 3399141 : ALLOCATE (work(index1)%array1(Nsize + 1))
827 23824901 : work(index1)%array1(1:Nsize) = wrk_tmp
828 1133047 : work(index1)%array1(Nsize + 1) = index2
829 1133047 : DEALLOCATE (wrk_tmp)
830 :
831 1133047 : wrk_tmp => work(index2)%array1
832 1133047 : Nsize = SIZE(wrk_tmp)
833 3399141 : ALLOCATE (work(index2)%array1(Nsize + 1))
834 15722189 : work(index2)%array1(1:Nsize) = wrk_tmp
835 1133047 : work(index2)%array1(Nsize + 1) = index1
836 1133047 : DEALLOCATE (wrk_tmp)
837 :
838 1133047 : wrk_tmp => work(index1)%array2
839 1133047 : Nsize = SIZE(wrk_tmp)
840 3399141 : ALLOCATE (work(index1)%array2(Nsize + 1))
841 23824901 : work(index1)%array2(1:Nsize) = wrk_tmp
842 1133047 : work(index1)%array2(Nsize + 1) = index3
843 1133047 : DEALLOCATE (wrk_tmp)
844 :
845 1133047 : wrk_tmp => work(index2)%array2
846 1133047 : Nsize = SIZE(wrk_tmp)
847 3399141 : ALLOCATE (work(index2)%array2(Nsize + 1))
848 15722189 : work(index2)%array2(1:Nsize) = wrk_tmp
849 1133047 : work(index2)%array2(Nsize + 1) = -index3
850 1141082 : DEALLOCATE (wrk_tmp)
851 : END DO
852 :
853 8035 : 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 10470 : 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 10470 : N = SIZE(atom_bond_list)
871 10470 : nmol = 0
872 772169 : DO I = 1, N
873 772169 : IF (mol_info(I) == -1) THEN
874 329474 : nmol = nmol + 1
875 329474 : my_mol_name = mol_name(I)
876 : CALL spread_mol(atom_bond_list, mol_info, i, nmol, my_mol_name, &
877 329474 : mol_name)
878 : END IF
879 : END DO
880 10470 : 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 761699 : 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 761699 : mol_info(iatom) = imol
902 1794129 : DO I = 1, SIZE(atom_bond_list(iatom)%array1)
903 1032430 : 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 1794129 : IF (mol_info(atom_b) == -1 .AND. my_mol_name == mol_name(atom_b)) THEN
908 432225 : CALL spread_mol(atom_bond_list, mol_info, atom_b, imol, my_mol_name, mol_name)
909 432225 : IF (mol_info(atom_b) /= imol) CPABORT("internal error")
910 : END IF
911 : END DO
912 761699 : 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 9675 : 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 9675 : DIMENSION(:), POINTER :: keyword
929 : INTEGER :: handle, i, i_rep, iatom, ielem_found, &
930 : iw, n_rep, natom
931 : LOGICAL :: user_defined
932 9675 : 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 9675 : NULLIFY (logger)
938 19350 : logger => cp_get_default_logger()
939 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
940 9675 : extension=".subsysLog")
941 9675 : CALL timeset(routineN, handle)
942 :
943 9675 : atom_info => topology%atom_info
944 9675 : natom = topology%natoms
945 :
946 : ! Available external info
947 9675 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
948 9675 : CALL section_vals_get(kind_section, n_repetition=n_rep)
949 25345 : ALLOCATE (keyword(n_rep))
950 25345 : ALLOCATE (mass(n_rep))
951 21702 : mass = HUGE(0.0_dp)
952 21702 : DO i_rep = 1, n_rep
953 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
954 12027 : c_val=keyword(i_rep), i_rep_section=i_rep)
955 12027 : CALL uppercase(keyword(i_rep))
956 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
957 12027 : keyword_name="MASS", n_rep_val=i)
958 12027 : IF (i > 0) CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
959 21798 : keyword_name="MASS", r_val=mass(i_rep))
960 : END DO
961 : !
962 287393 : 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 277718 : user_defined = .FALSE.
966 390024 : DO i = 1, SIZE(keyword)
967 112900 : atmname_upper = id2str(atom_info%id_atmname(iatom))
968 112900 : CALL uppercase(atmname_upper)
969 390024 : IF (TRIM(atmname_upper) == TRIM(keyword(i)) .AND. mass(i) /= HUGE(0.0_dp)) THEN
970 594 : 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 277124 : IF (.NOT. user_defined) THEN
977 277124 : upper_sym_1 = id2str(atom_info%id_element(iatom))
978 277124 : CALL get_ptable_info(symbol=upper_sym_1, ielement=ielem_found, amass=atom_info%atm_mass(iatom))
979 : END IF
980 279122 : IF (iw > 0) WRITE (iw, '(7X,A,A5,A,F12.5)') "In topology_set_atm_mass :: element = ", &
981 12483 : id2str(atom_info%id_element(iatom)), " a_mass ", atom_info%atm_mass(iatom)
982 : END DO
983 9675 : DEALLOCATE (keyword)
984 9675 : DEALLOCATE (mass)
985 :
986 9675 : CALL timestop(handle)
987 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
988 9675 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
989 :
990 29025 : 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 9628 : 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 9628 : 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 9628 : NULLIFY (logger)
1013 9628 : logger => cp_get_default_logger()
1014 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/UTIL_INFO", &
1015 9628 : extension=".subsysLog")
1016 9628 : CALL timeset(routineN, handle)
1017 :
1018 9628 : atom_info => topology%atom_info
1019 9628 : conn_info => topology%conn_info
1020 9628 : natom = topology%natoms
1021 :
1022 9654 : IF (iw > 0) WRITE (iw, '(A)') "Start of Molecule_Check", &
1023 52 : " Checking consistency between the generated molecules"
1024 :
1025 778361 : ALLOCATE (atom_bond_list(natom))
1026 759105 : DO I = 1, natom
1027 759105 : ALLOCATE (atom_bond_list(I)%array1(0))
1028 : END DO
1029 9628 : N = 0
1030 9628 : IF (ASSOCIATED(conn_info%bond_a)) N = SIZE(conn_info%bond_a)
1031 9628 : CALL reorder_structure(atom_bond_list, conn_info%bond_a, conn_info%bond_b, N)
1032 :
1033 9628 : mol_typ = atom_info%map_mol_typ(1)
1034 9628 : mol_num = atom_info%map_mol_num(1)
1035 9628 : counter = 1
1036 9628 : loc_counter = 1
1037 9628 : first = 1
1038 9628 : first_loc = 1
1039 749477 : DO iatom = 2, natom
1040 739849 : icheck_num = (atom_info%map_mol_num(iatom) == mol_num)
1041 739849 : icheck_typ = (atom_info%map_mol_typ(iatom) == mol_typ)
1042 739849 : 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 286840 : 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 0 : cp_to_string(loc_counter))
1054 : END IF
1055 : END IF
1056 739849 : IF (.NOT. icheck_typ) THEN
1057 118599 : first = iatom
1058 118599 : first_loc = iatom
1059 118599 : counter = 1
1060 118599 : loc_counter = 1
1061 118599 : mol_typ = atom_info%map_mol_typ(iatom)
1062 : END IF
1063 739849 : IF (icheck_num) THEN
1064 570916 : IF (icheck_typ) loc_counter = loc_counter + 1
1065 : !-----------------------------------------------------------------------------
1066 : !-----------------------------------------------------------------------------
1067 : ! 2. Check that each molecule has the same atom sequences
1068 : !-----------------------------------------------------------------------------
1069 570916 : 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 0 : " 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 570916 : 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 0 : " NOT FOUND! Check the connectivity of your system.")
1092 : END IF
1093 :
1094 1278456 : DO k = 1, SIZE(atom_bond_list(iatom)%array1)
1095 1022827 : IF (ALL(atom_bond_list(first + loc_counter - 1)%array1 - first /= &
1096 570916 : 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 0 : " NOT FOUND! Check the connectivity of your system.")
1102 : END IF
1103 : END DO
1104 : ELSE
1105 168933 : mol_num = atom_info%map_mol_num(iatom)
1106 168933 : loc_counter = 1
1107 168933 : first_loc = iatom
1108 : END IF
1109 749477 : IF (mol_num == 1 .AND. icheck_typ) counter = counter + 1
1110 : END DO
1111 9628 : IF (iw > 0) WRITE (iw, '(A)') "End of Molecule_Check"
1112 :
1113 759105 : DO I = 1, natom
1114 759105 : DEALLOCATE (atom_bond_list(I)%array1)
1115 : END DO
1116 9628 : DEALLOCATE (atom_bond_list)
1117 9628 : CALL timestop(handle)
1118 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
1119 9628 : "PRINT%TOPOLOGY_INFO/UTIL_INFO")
1120 :
1121 19256 : 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 54326 : 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 27163 : found = .FALSE.
1146 27163 : element_symbol = element_in
1147 27163 : atom_name = atom_name_in
1148 27163 : element_out = ""
1149 27163 : defined_kind_section = .FALSE.
1150 :
1151 : ! First check if a KIND section is overriding the element
1152 : ! definition
1153 27163 : CALL uppercase(atom_name)
1154 27163 : kind_section => section_vals_get_subs_vals(subsys_section, "KIND")
1155 27163 : CALL section_vals_get(kind_section, n_repetition=n_rep)
1156 80432 : DO i_rep = 1, n_rep
1157 : CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
1158 54509 : c_val=keyword, i_rep_section=i_rep)
1159 54509 : CALL uppercase(keyword)
1160 80432 : IF (TRIM(keyword) == TRIM(atom_name)) THEN
1161 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
1162 14541 : keyword_name="ELEMENT", n_rep_val=i)
1163 14541 : IF (i > 0) THEN
1164 : CALL section_vals_val_get(kind_section, i_rep_section=i_rep, &
1165 1240 : keyword_name="ELEMENT", c_val=element_symbol)
1166 : defined_kind_section = .TRUE.
1167 : EXIT
1168 : ELSE
1169 13301 : 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 25923 : 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 22687 : in_ptable = .FALSE.
1182 22687 : IF (LEN_TRIM(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
1183 22687 : IF (in_ptable) THEN
1184 22593 : element_out = TRIM(element_symbol)
1185 22593 : found = .TRUE.
1186 : END IF
1187 : END IF
1188 :
1189 : ! This is clearly a user error
1190 27163 : 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 0 : "input file!")
1195 :
1196 : ! Last chance.. are these atom_kinds of AMBER or CHARMM or GROMOS FF ?
1197 27163 : CALL uppercase(element_symbol)
1198 27163 : IF ((.NOT. found) .AND. (ASSOCIATED(amber_map))) THEN
1199 : ! First we go through the AMBER library
1200 145670 : DO i = 1, SIZE(amber_map%kind)
1201 145670 : IF (element_symbol == amber_map%kind(i)) THEN
1202 : found = .TRUE.
1203 : EXIT
1204 : END IF
1205 : END DO
1206 4570 : IF (found) THEN
1207 4134 : element_out = amber_map%element(i)
1208 : END IF
1209 : END IF
1210 4570 : IF ((.NOT. found) .AND. (ASSOCIATED(charmm_map))) THEN
1211 : ! Then we go through the CHARMM library
1212 39242 : DO i = 1, SIZE(charmm_map%kind)
1213 39242 : IF (element_symbol == charmm_map%kind(i)) THEN
1214 : found = .TRUE.
1215 : EXIT
1216 : END IF
1217 : END DO
1218 436 : IF (found) THEN
1219 196 : element_out = charmm_map%element(i)
1220 : END IF
1221 : END IF
1222 436 : IF ((.NOT. found) .AND. (ASSOCIATED(gromos_map))) THEN
1223 : ! Then we go through the GROMOS library
1224 5520 : DO i = 1, SIZE(gromos_map%kind)
1225 5520 : IF (element_symbol == gromos_map%kind(i)) THEN
1226 : found = .TRUE.
1227 : EXIT
1228 : END IF
1229 : END DO
1230 240 : IF (found) THEN
1231 0 : 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 0 : IF (.NOT. found) THEN
1238 240 : in_ptable = .FALSE.
1239 240 : IF (LEN_TRIM(element_symbol) <= 2) CALL get_ptable_info(element_symbol, found=in_ptable)
1240 240 : IF (in_ptable) THEN
1241 240 : 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 0 : 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 0 : " the developers!")
1253 : END IF
1254 :
1255 27163 : END SUBROUTINE check_subsys_element
1256 :
1257 0 : END MODULE topology_util
|
