92 CHARACTER(len=*),
PARAMETER :: routinen =
'read_coordinate_cif'
93 INTEGER,
PARAMETER :: nblock = 1000
94 REAL(kind=
dp),
PARAMETER :: threshold = 1.0e-3_dp
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
111 extension=
".subsysLog")
112 CALL timeset(routinen, handle)
113 pfactor =
section_get_rval(subsys_section,
"TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
119 CALL reallocate(atom_info%id_molname, 1, nblock)
120 CALL reallocate(atom_info%id_resname, 1, nblock)
122 CALL reallocate(atom_info%id_atmname, 1, nblock)
124 CALL reallocate(atom_info%atm_mass, 1, nblock)
125 CALL reallocate(atom_info%atm_charge, 1, nblock)
128 CALL reallocate(atom_info%id_element, 1, nblock)
130 IF (iw > 0)
WRITE (iw,
"(/,A,A)")
" Reading in CIF file ", trim(
topology%coord_file_name)
139 para_env=para_env, apply_preprocessing=.false.)
143 begin_line=.false., search_from_begin_of_file=.true.)
145 IF (iw > 0)
WRITE (iw,
'(/,A)')
" CIF_INFO| _chemical_name :: "//trim(parser%input_line(parser%icol:))
150 begin_line=.false., search_from_begin_of_file=.true.)
152 IF (iw > 0)
WRITE (iw,
'(A)')
" CIF_INFO| _chemical_formula_sum :: "//trim(parser%input_line(parser%icol:))
157 field_label = -1; field_symbol = -1; field_x = -1; field_y = -1; field_z = -1
159 begin_line=.false., search_from_begin_of_file=.true.)
160 DO WHILE (index(parser%input_line,
"_atom_site_") /= 0)
161 num_fields = num_fields + 1
162 IF (index(parser%input_line,
"_atom_site_label") /= 0) field_label = num_fields
163 IF (index(parser%input_line,
"_atom_site_type_symbol") /= 0) field_symbol = num_fields
164 IF (index(parser%input_line,
"_atom_site_fract_x") /= 0) field_x = num_fields
165 IF (index(parser%input_line,
"_atom_site_fract_y") /= 0) field_y = num_fields
166 IF (index(parser%input_line,
"_atom_site_fract_z") /= 0) field_z = num_fields
171 IF (field_label < 0) cpabort(
"Field _atom_site_label not found in CIF file.")
172 IF (field_x < 0) cpabort(
"Field _atom_site_fract_x not found in CIF file.")
173 IF (field_y < 0) cpabort(
"Field _atom_site_fract_y not found in CIF file.")
174 IF (field_z < 0) cpabort(
"Field _atom_site_fract_z not found in CIF file.")
177 IF (field_symbol > 0)
THEN
178 field_kind = field_symbol
180 field_kind = field_label
185 DO WHILE ((index(parser%input_line,
"loop_") == 0) .AND. (parser%input_line(1:1) /=
"_"))
188 IF (natom >
SIZE(atom_info%id_molname))
THEN
189 newsize = int(pfactor*natom)
190 CALL reallocate(atom_info%id_molname, 1, newsize)
191 CALL reallocate(atom_info%id_resname, 1, newsize)
193 CALL reallocate(atom_info%id_atmname, 1, newsize)
194 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
195 CALL reallocate(atom_info%atm_mass, 1, newsize)
196 CALL reallocate(atom_info%atm_charge, 1, newsize)
199 CALL reallocate(atom_info%id_element, 1, newsize)
201 DO ii = 1, num_fields
202 IF (ii == field_kind)
THEN
204 atom_info%id_atmname(natom) =
str2id(strtmp)
206 atom_info%id_resname(natom) = atom_info%id_molname(natom)
207 atom_info%resid(natom) = 1
208 atom_info%id_element(natom) = atom_info%id_atmname(natom)
209 ELSE IF (ii == field_x)
THEN
210 CALL cif_get_real(parser, atom_info%r(1, natom))
211 ELSE IF (ii == field_y)
THEN
212 CALL cif_get_real(parser, atom_info%r(2, natom))
213 ELSE IF (ii == field_z)
THEN
214 CALL cif_get_real(parser, atom_info%r(3, natom))
219 s = atom_info%r(1:3, natom)
227 r1 = atom_info%r(1:3, ii)
228 DO jj = ii + 1, natom
229 r2 = atom_info%r(1:3, jj)
230 r =
pbc(r1 - r2, cell)
232 check = (dot_product(r, r) >= (threshold*threshold))
238 CALL parser_search_string(parser,
"_symmetry_space_group_name_h-m", ignore_case=.false., found=found, &
239 begin_line=.false., search_from_begin_of_file=.true.)
241 IF (iw > 0)
WRITE (iw,
'(A)')
" CIF_INFO| _symmetry_space_group_name_h-m :: "//trim(parser%input_line(parser%icol:))
246 CALL parser_search_string(parser,
"_symmetry_equiv_pos_as_xyz", ignore_case=.false., found=found, &
247 begin_line=.false., search_from_begin_of_file=.true.)
248 IF (.NOT. found)
THEN
249 CALL parser_search_string(parser,
"_space_group_symop_operation_xyz", ignore_case=.false., found=found, &
250 begin_line=.false., search_from_begin_of_file=.true.)
253 CALL cp_warn(__location__,
"The fields (_symmetry_equiv_pos_as_xyz) or "// &
254 "(_space_group_symop_operation_xyz) were not found in CIF file!")
255 IF (iw > 0)
WRITE (iw,
'(A,I0)')
" CIF_INFO| Number of atoms before applying symmetry operations :: ", natom
256 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)
263 DO WHILE ((index(parser%input_line,
"loop_") == 0) .AND. (parser%input_line(1:1) /=
"_"))
267 IF (index(parser%input_line(1:),
'"') > 0) sep =
'"'
268 iln0 = index(parser%input_line(1:), sep)
269 iln1 = index(parser%input_line(iln0 + 1:),
",") + iln0
270 iln2 = index(parser%input_line(iln1 + 1:),
",") + iln1
272 iln3 = len_trim(parser%input_line) + 1
274 iln3 = index(parser%input_line(iln2 + 1:), sep) + iln2
277 cpassert(iln2 /= iln1)
278 cpassert(iln3 /= iln2)
280 CALL parsef(1, trim(parser%input_line(iln0 + 1:iln1 - 1)),
s2a(
"x",
"y",
"z"))
281 CALL parsef(2, trim(parser%input_line(iln1 + 1:iln2 - 1)),
s2a(
"x",
"y",
"z"))
282 CALL parsef(3, trim(parser%input_line(iln2 + 1:iln3 - 1)),
s2a(
"x",
"y",
"z"))
283 loop_over_unique_atoms:
DO ii = 1, natom_orig
285 s(1) =
evalf(1, (/s_tmp(1), s_tmp(2), s_tmp(3)/))
286 s(2) =
evalf(2, (/s_tmp(1), s_tmp(2), s_tmp(3)/))
287 s(3) =
evalf(3, (/s_tmp(1), s_tmp(2), s_tmp(3)/))
291 r2 = atom_info%r(1:3, jj)
292 r =
pbc(r1 - r2, cell)
294 IF (dot_product(r, r) <= (threshold*threshold))
THEN
303 IF (natom >
SIZE(atom_info%id_molname))
THEN
304 newsize = int(pfactor*natom)
305 CALL reallocate(atom_info%id_molname, 1, newsize)
306 CALL reallocate(atom_info%id_resname, 1, newsize)
308 CALL reallocate(atom_info%id_atmname, 1, newsize)
309 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
310 CALL reallocate(atom_info%atm_mass, 1, newsize)
311 CALL reallocate(atom_info%atm_charge, 1, newsize)
314 CALL reallocate(atom_info%id_element, 1, newsize)
316 atom_info%id_atmname(natom) = atom_info%id_atmname(ii)
317 atom_info%id_molname(natom) = atom_info%id_molname(ii)
318 atom_info%id_resname(natom) = atom_info%id_resname(ii)
319 atom_info%id_element(natom) = atom_info%id_element(ii)
320 atom_info%resid(natom) = atom_info%resid(ii)
321 atom_info%r(1:3, natom) = r1
323 END DO loop_over_unique_atoms
329 IF (iw > 0)
WRITE (iw,
'(A,I0)')
" CIF_INFO| Number of symmetry operations :: ", isym
330 IF (iw > 0)
WRITE (iw,
'(A,I0)')
" CIF_INFO| Number of total atoms :: ", natom
331 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)
338 CALL reallocate(atom_info%id_molname, 1, natom)
339 CALL reallocate(atom_info%id_resname, 1, natom)
341 CALL reallocate(atom_info%id_atmname, 1, natom)
344 CALL reallocate(atom_info%atm_charge, 1, natom)
347 CALL reallocate(atom_info%id_element, 1, natom)
352 "PRINT%TOPOLOGY_INFO/CIF_INFO")
353 CALL timestop(handle)