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 Handles CIF (Crystallographic Information File) files
10 : !> \author Teodoro Laino [tlaino]
11 : !> \date 12.2008
12 : ! **************************************************************************************************
13 : MODULE topology_cif
14 : USE cell_methods, ONLY: cell_create,&
15 : read_cell_cif,&
16 : write_cell
17 : USE cell_types, ONLY: cell_release,&
18 : cell_type,&
19 : pbc,&
20 : real_to_scaled,&
21 : scaled_to_real
22 : USE cp_log_handling, ONLY: cp_get_default_logger,&
23 : cp_logger_type,&
24 : cp_to_string
25 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
26 : cp_print_key_unit_nr
27 : USE cp_parser_methods, ONLY: parser_get_next_line,&
28 : parser_get_object,&
29 : parser_search_string
30 : USE cp_parser_types, ONLY: cp_parser_type,&
31 : parser_create,&
32 : parser_release
33 : USE fparser, ONLY: evalf,&
34 : finalizef,&
35 : initf,&
36 : parsef
37 : USE input_section_types, ONLY: section_get_rval,&
38 : section_vals_type
39 : USE kinds, ONLY: default_string_length,&
40 : dp
41 : USE memory_utilities, ONLY: reallocate
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE string_table, ONLY: id2str,&
44 : s2s,&
45 : str2id
46 : USE string_utilities, ONLY: s2a
47 : USE topology_types, ONLY: atom_info_type,&
48 : topology_parameters_type
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 :
53 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_cif'
54 :
55 : PRIVATE
56 : PUBLIC :: read_coordinate_cif
57 :
58 : CONTAINS
59 :
60 : ! **************************************************************************************************
61 : !> \brief Performs the real task of reading the proper information from the CIF
62 : !> file
63 : !> \param topology ...
64 : !> \param para_env ...
65 : !> \param subsys_section ...
66 : !> \date 12.2008
67 : !> \par Format Information implemented:
68 : !> _chemical_name
69 : !> _chemical_formula_sum
70 : !> _cell_length_a
71 : !> _cell_length_b
72 : !> _cell_length_c
73 : !> _cell_angle_alpha
74 : !> _cell_angle_beta
75 : !> _cell_angle_gamma
76 : !> _symmetry_space_group_name_h-m
77 : !> _symmetry_equiv_pos_as_xyz
78 : !> _space_group_symop_operation_xyz
79 : !> _atom_site_label
80 : !> _atom_site_type_symbol
81 : !> _atom_site_fract_x
82 : !> _atom_site_fract_y
83 : !> _atom_site_fract_z
84 : !>
85 : !> \author Teodoro Laino [tlaino]
86 : ! **************************************************************************************************
87 28 : SUBROUTINE read_coordinate_cif(topology, para_env, subsys_section)
88 : TYPE(topology_parameters_type) :: topology
89 : TYPE(mp_para_env_type), POINTER :: para_env
90 : TYPE(section_vals_type), POINTER :: subsys_section
91 :
92 : CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_cif'
93 : INTEGER, PARAMETER :: nblock = 1000
94 : REAL(KIND=dp), PARAMETER :: threshold = 1.0E-3_dp
95 :
96 : CHARACTER(LEN=1) :: sep
97 : CHARACTER(LEN=default_string_length) :: s_tag, strtmp
98 : INTEGER :: field_kind, field_label, field_symbol, field_x, field_y, field_z, handle, ii, &
99 : iln0, iln1, iln2, iln3, isym, iw, jj, natom, natom_orig, newsize, num_fields
100 : LOGICAL :: check, found, my_end
101 : REAL(KIND=dp) :: pfactor
102 : REAL(KIND=dp), DIMENSION(3) :: r, r1, r2, s, s_tmp
103 : TYPE(atom_info_type), POINTER :: atom_info
104 : TYPE(cell_type), POINTER :: cell
105 : TYPE(cp_logger_type), POINTER :: logger
106 : TYPE(cp_parser_type) :: parser
107 :
108 14 : NULLIFY (logger)
109 28 : logger => cp_get_default_logger()
110 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/CIF_INFO", &
111 14 : extension=".subsysLog")
112 14 : CALL timeset(routineN, handle)
113 14 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
114 :
115 : ! Element is assigned on the basis of the atm_name
116 14 : topology%aa_element = .TRUE.
117 :
118 14 : atom_info => topology%atom_info
119 14 : CALL reallocate(atom_info%id_molname, 1, nblock)
120 14 : CALL reallocate(atom_info%id_resname, 1, nblock)
121 14 : CALL reallocate(atom_info%resid, 1, nblock)
122 14 : CALL reallocate(atom_info%id_atmname, 1, nblock)
123 14 : CALL reallocate(atom_info%r, 1, 3, 1, nblock)
124 14 : CALL reallocate(atom_info%atm_mass, 1, nblock)
125 14 : CALL reallocate(atom_info%atm_charge, 1, nblock)
126 14 : CALL reallocate(atom_info%occup, 1, nblock)
127 14 : CALL reallocate(atom_info%beta, 1, nblock)
128 14 : CALL reallocate(atom_info%id_element, 1, nblock)
129 :
130 14 : IF (iw > 0) WRITE (iw, "(/,A,A)") " Reading in CIF file ", TRIM(topology%coord_file_name)
131 :
132 : ! Create cell
133 14 : NULLIFY (cell)
134 14 : CALL cell_create(cell, tag="CELL_CIF")
135 14 : CALL read_cell_cif(topology%coord_file_name, cell, para_env)
136 14 : CALL write_cell(cell, subsys_section)
137 :
138 : CALL parser_create(parser, topology%coord_file_name, &
139 14 : para_env=para_env, apply_preprocessing=.FALSE.)
140 :
141 : ! Check for _chemical_name
142 : CALL parser_search_string(parser, "_chemical_name", ignore_case=.FALSE., found=found, &
143 14 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
144 14 : IF (found) THEN
145 6 : IF (iw > 0) WRITE (iw, '(/,A)') " CIF_INFO| _chemical_name :: "//TRIM(parser%input_line(parser%icol:))
146 : END IF
147 :
148 : ! Check for _chemical_formula_sum
149 : CALL parser_search_string(parser, "_chemical_formula_sum", ignore_case=.FALSE., found=found, &
150 14 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
151 14 : IF (found) THEN
152 14 : IF (iw > 0) WRITE (iw, '(A)') " CIF_INFO| _chemical_formula_sum :: "//TRIM(parser%input_line(parser%icol:))
153 : END IF
154 :
155 : ! Parse atoms info and fractional coordinates
156 14 : num_fields = 0
157 14 : field_label = -1; field_symbol = -1; field_x = -1; field_y = -1; field_z = -1
158 : CALL parser_search_string(parser, "_atom_site_", ignore_case=.FALSE., found=found, &
159 14 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
160 110 : DO WHILE (INDEX(parser%input_line, "_atom_site_") /= 0)
161 96 : num_fields = num_fields + 1
162 96 : IF (INDEX(parser%input_line, "_atom_site_label") /= 0) field_label = num_fields
163 96 : IF (INDEX(parser%input_line, "_atom_site_type_symbol") /= 0) field_symbol = num_fields
164 96 : IF (INDEX(parser%input_line, "_atom_site_fract_x") /= 0) field_x = num_fields
165 96 : IF (INDEX(parser%input_line, "_atom_site_fract_y") /= 0) field_y = num_fields
166 96 : IF (INDEX(parser%input_line, "_atom_site_fract_z") /= 0) field_z = num_fields
167 110 : CALL parser_get_next_line(parser, 1)
168 : END DO
169 :
170 : ! Check for missing fields.
171 14 : IF (field_label < 0) CPABORT("Field _atom_site_label not found in CIF file.")
172 14 : IF (field_x < 0) CPABORT("Field _atom_site_fract_x not found in CIF file.")
173 14 : IF (field_y < 0) CPABORT("Field _atom_site_fract_y not found in CIF file.")
174 14 : IF (field_z < 0) CPABORT("Field _atom_site_fract_z not found in CIF file.")
175 :
176 : ! Read atom kind from the symbol field if it's present, otherwise fall back to the label field.
177 14 : IF (field_symbol > 0) THEN
178 : field_kind = field_symbol
179 : ELSE
180 6 : field_kind = field_label
181 : END IF
182 :
183 : ! Parse real info
184 14 : natom = 0
185 112 : DO WHILE ((INDEX(parser%input_line, "loop_") == 0) .AND. (parser%input_line(1:1) /= "_"))
186 110 : natom = natom + 1
187 : ! Resize in case needed
188 110 : IF (natom > SIZE(atom_info%id_molname)) THEN
189 0 : newsize = INT(pfactor*natom)
190 0 : CALL reallocate(atom_info%id_molname, 1, newsize)
191 0 : CALL reallocate(atom_info%id_resname, 1, newsize)
192 0 : CALL reallocate(atom_info%resid, 1, newsize)
193 0 : CALL reallocate(atom_info%id_atmname, 1, newsize)
194 0 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
195 0 : CALL reallocate(atom_info%atm_mass, 1, newsize)
196 0 : CALL reallocate(atom_info%atm_charge, 1, newsize)
197 0 : CALL reallocate(atom_info%occup, 1, newsize)
198 0 : CALL reallocate(atom_info%beta, 1, newsize)
199 0 : CALL reallocate(atom_info%id_element, 1, newsize)
200 : END IF
201 1150 : DO ii = 1, num_fields
202 1150 : IF (ii == field_kind) THEN
203 110 : CALL parser_get_object(parser, strtmp)
204 110 : atom_info%id_atmname(natom) = str2id(strtmp)
205 110 : atom_info%id_molname(natom) = str2id(s2s("MOL"//TRIM(ADJUSTL(cp_to_string(natom)))))
206 110 : atom_info%id_resname(natom) = atom_info%id_molname(natom)
207 110 : atom_info%resid(natom) = 1
208 110 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
209 930 : ELSE IF (ii == field_x) THEN
210 110 : CALL cif_get_real(parser, atom_info%r(1, natom))
211 820 : ELSE IF (ii == field_y) THEN
212 110 : CALL cif_get_real(parser, atom_info%r(2, natom))
213 710 : ELSE IF (ii == field_z) THEN
214 110 : CALL cif_get_real(parser, atom_info%r(3, natom))
215 : ELSE
216 600 : CALL parser_get_object(parser, s_tag) ! Skip this field
217 : END IF
218 : END DO
219 440 : s = atom_info%r(1:3, natom)
220 110 : CALL scaled_to_real(atom_info%r(1:3, natom), s, cell)
221 110 : CALL parser_get_next_line(parser, 1, at_end=my_end)
222 112 : IF (my_end) EXIT
223 : END DO
224 : ! Preliminary check: check if atoms provided are really unique.. this is a paranoic
225 : ! check since they should be REALLY unique.. anyway..
226 124 : DO ii = 1, natom
227 440 : r1 = atom_info%r(1:3, ii)
228 1018 : DO jj = ii + 1, natom
229 3576 : r2 = atom_info%r(1:3, jj)
230 3576 : r = pbc(r1 - r2, cell)
231 : ! check = (SQRT(DOT_PRODUCT(r, r)) >= threshold)
232 3576 : check = (DOT_PRODUCT(r, r) >= (threshold*threshold))
233 1004 : CPASSERT(check)
234 : END DO
235 : END DO
236 : ! Parse Symmetry Group and generation elements..
237 : ! Check for _symmetry_space_group_name_h-m
238 : CALL parser_search_string(parser, "_symmetry_space_group_name_h-m", ignore_case=.FALSE., found=found, &
239 14 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
240 14 : IF (found) THEN
241 4 : IF (iw > 0) WRITE (iw, '(A)') " CIF_INFO| _symmetry_space_group_name_h-m :: "//TRIM(parser%input_line(parser%icol:))
242 : END IF
243 :
244 : ! Check for _symmetry_equiv_pos_as_xyz
245 : ! Check for _space_group_symop_operation_xyz
246 : CALL parser_search_string(parser, "_symmetry_equiv_pos_as_xyz", ignore_case=.FALSE., found=found, &
247 14 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
248 14 : IF (.NOT. found) THEN
249 : CALL parser_search_string(parser, "_space_group_symop_operation_xyz", ignore_case=.FALSE., found=found, &
250 8 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
251 : END IF
252 14 : IF (.NOT. found) &
253 : CALL cp_warn(__LOCATION__, "The fields (_symmetry_equiv_pos_as_xyz) or "// &
254 0 : "(_space_group_symop_operation_xyz) were not found in CIF file!")
255 14 : IF (iw > 0) WRITE (iw, '(A,I0)') " CIF_INFO| Number of atoms before applying symmetry operations :: ", natom
256 43 : IF (iw > 0) WRITE (iw, '(A10,1X,3F12.6)') (TRIM(id2str(atom_info%id_atmname(ii))), atom_info%r(1:3, ii), ii=1, natom)
257 14 : isym = 0
258 14 : IF (found) THEN
259 : ! Apply symmetry elements and generate the whole set of atoms in the unit cell
260 14 : CALL parser_get_next_line(parser, 1)
261 : isym = 0
262 14 : natom_orig = natom
263 794 : DO WHILE ((INDEX(parser%input_line, "loop_") == 0) .AND. (parser%input_line(1:1) /= "_"))
264 780 : isym = isym + 1
265 : ! find seprator ' or "
266 780 : sep = "'"
267 780 : IF (INDEX(parser%input_line(1:), '"') > 0) sep = '"'
268 780 : iln0 = INDEX(parser%input_line(1:), sep)
269 780 : iln1 = INDEX(parser%input_line(iln0 + 1:), ",") + iln0
270 780 : iln2 = INDEX(parser%input_line(iln1 + 1:), ",") + iln1
271 780 : IF (iln0 == 0) THEN
272 768 : iln3 = LEN_TRIM(parser%input_line) + 1
273 : ELSE
274 12 : iln3 = INDEX(parser%input_line(iln2 + 1:), sep) + iln2
275 : END IF
276 780 : CPASSERT(iln1 /= 0)
277 780 : CPASSERT(iln2 /= iln1)
278 780 : CPASSERT(iln3 /= iln2)
279 780 : CALL initf(3)
280 780 : CALL parsef(1, TRIM(parser%input_line(iln0 + 1:iln1 - 1)), s2a("x", "y", "z"))
281 780 : CALL parsef(2, TRIM(parser%input_line(iln1 + 1:iln2 - 1)), s2a("x", "y", "z"))
282 780 : CALL parsef(3, TRIM(parser%input_line(iln2 + 1:iln3 - 1)), s2a("x", "y", "z"))
283 2468 : Loop_over_unique_atoms: DO ii = 1, natom_orig
284 1688 : CALL real_to_scaled(s_tmp, atom_info%r(1:3, ii), cell)
285 6752 : s(1) = evalf(1, (/s_tmp(1), s_tmp(2), s_tmp(3)/))
286 6752 : s(2) = evalf(2, (/s_tmp(1), s_tmp(2), s_tmp(3)/))
287 6752 : s(3) = evalf(3, (/s_tmp(1), s_tmp(2), s_tmp(3)/))
288 1688 : CALL scaled_to_real(r1, s, cell)
289 1688 : check = .TRUE.
290 9706 : DO jj = 1, natom
291 38536 : r2 = atom_info%r(1:3, jj)
292 38536 : r = pbc(r1 - r2, cell)
293 : ! SQRT(DOT_PRODUCT(r, r)) <= threshold
294 38608 : IF (DOT_PRODUCT(r, r) <= (threshold*threshold)) THEN
295 : check = .FALSE.
296 : EXIT
297 : END IF
298 : END DO
299 : ! If the atom generated is unique let's add to the atom set..
300 2468 : IF (check) THEN
301 72 : natom = natom + 1
302 : ! Resize in case needed
303 72 : IF (natom > SIZE(atom_info%id_molname)) THEN
304 0 : newsize = INT(pfactor*natom)
305 0 : CALL reallocate(atom_info%id_molname, 1, newsize)
306 0 : CALL reallocate(atom_info%id_resname, 1, newsize)
307 0 : CALL reallocate(atom_info%resid, 1, newsize)
308 0 : CALL reallocate(atom_info%id_atmname, 1, newsize)
309 0 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
310 0 : CALL reallocate(atom_info%atm_mass, 1, newsize)
311 0 : CALL reallocate(atom_info%atm_charge, 1, newsize)
312 0 : CALL reallocate(atom_info%occup, 1, newsize)
313 0 : CALL reallocate(atom_info%beta, 1, newsize)
314 0 : CALL reallocate(atom_info%id_element, 1, newsize)
315 : END IF
316 72 : atom_info%id_atmname(natom) = atom_info%id_atmname(ii)
317 72 : atom_info%id_molname(natom) = atom_info%id_molname(ii)
318 72 : atom_info%id_resname(natom) = atom_info%id_resname(ii)
319 72 : atom_info%id_element(natom) = atom_info%id_element(ii)
320 72 : atom_info%resid(natom) = atom_info%resid(ii)
321 288 : atom_info%r(1:3, natom) = r1
322 : END IF
323 : END DO Loop_over_unique_atoms
324 780 : CALL finalizef()
325 780 : CALL parser_get_next_line(parser, 1, at_end=my_end)
326 794 : IF (my_end) EXIT
327 : END DO
328 : END IF
329 14 : IF (iw > 0) WRITE (iw, '(A,I0)') " CIF_INFO| Number of symmetry operations :: ", isym
330 14 : IF (iw > 0) WRITE (iw, '(A,I0)') " CIF_INFO| Number of total atoms :: ", natom
331 79 : IF (iw > 0) WRITE (iw, '(A10,1X,3F12.6)') (TRIM(id2str(atom_info%id_atmname(ii))), atom_info%r(1:3, ii), ii=1, natom)
332 :
333 : ! Releasse local cell type and parser
334 14 : CALL cell_release(cell)
335 14 : CALL parser_release(parser)
336 :
337 : ! Reallocate all structures with the exact NATOM size
338 14 : CALL reallocate(atom_info%id_molname, 1, natom)
339 14 : CALL reallocate(atom_info%id_resname, 1, natom)
340 14 : CALL reallocate(atom_info%resid, 1, natom)
341 14 : CALL reallocate(atom_info%id_atmname, 1, natom)
342 14 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
343 14 : CALL reallocate(atom_info%atm_mass, 1, natom)
344 14 : CALL reallocate(atom_info%atm_charge, 1, natom)
345 14 : CALL reallocate(atom_info%occup, 1, natom)
346 14 : CALL reallocate(atom_info%beta, 1, natom)
347 14 : CALL reallocate(atom_info%id_element, 1, natom)
348 :
349 14 : topology%natoms = natom
350 14 : topology%molname_generated = .TRUE.
351 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
352 14 : "PRINT%TOPOLOGY_INFO/CIF_INFO")
353 14 : CALL timestop(handle)
354 42 : END SUBROUTINE read_coordinate_cif
355 :
356 : ! **************************************************************************************************
357 : !> \brief Reads REAL from the CIF file.. This wrapper is needed in order to
358 : !> treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
359 : !> \param parser ...
360 : !> \param r ...
361 : !> \date 12.2008
362 : !> \author Teodoro Laino [tlaino]
363 : ! **************************************************************************************************
364 330 : SUBROUTINE cif_get_real(parser, r)
365 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
366 : REAL(KIND=dp), INTENT(OUT) :: r
367 :
368 : CHARACTER(LEN=default_string_length) :: s_tag
369 : INTEGER :: iln
370 :
371 330 : CALL parser_get_object(parser, s_tag)
372 330 : iln = LEN_TRIM(s_tag)
373 330 : IF (INDEX(s_tag, "(") /= 0) iln = INDEX(s_tag, "(") - 1
374 330 : READ (s_tag(1:iln), *) r
375 330 : END SUBROUTINE cif_get_real
376 :
377 : END MODULE topology_cif
|
