Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Handles XTL (Molecular Simulations, Inc (MSI)) files
10 : !> \author Teodoro Laino [tlaino]
11 : !> \date 05.2009
12 : ! **************************************************************************************************
13 : MODULE topology_xtl
14 : USE cell_methods, ONLY: cell_create,&
15 : set_cell_param,&
16 : write_cell
17 : USE cell_types, ONLY: cell_release,&
18 : cell_type,&
19 : pbc,&
20 : scaled_to_real
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type,&
23 : cp_to_string
24 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
25 : cp_print_key_unit_nr
26 : USE cp_parser_methods, ONLY: parser_get_next_line,&
27 : parser_get_object,&
28 : parser_search_string
29 : USE cp_parser_types, ONLY: cp_parser_type,&
30 : parser_create,&
31 : parser_release
32 : USE cp_units, ONLY: cp_unit_to_cp2k
33 : USE input_section_types, ONLY: section_get_rval,&
34 : section_vals_type
35 : USE kinds, ONLY: default_string_length,&
36 : dp
37 : USE memory_utilities, ONLY: reallocate
38 : USE message_passing, ONLY: mp_para_env_type
39 : USE string_table, ONLY: id2str,&
40 : s2s,&
41 : str2id
42 : USE topology_types, ONLY: atom_info_type,&
43 : topology_parameters_type
44 : #include "./base/base_uses.f90"
45 :
46 : IMPLICIT NONE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'topology_xtl'
49 :
50 : PRIVATE
51 : PUBLIC :: read_coordinate_xtl
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief Performs the real task of reading the proper information from the XTL
57 : !> file
58 : !> \param topology ...
59 : !> \param para_env ...
60 : !> \param subsys_section ...
61 : !> \date 05.2009
62 : !> \par Format Information implemented:
63 : !> TITLE
64 : !> DIMENSION
65 : !> CELL
66 : !> SYMMETRY
67 : !> SYM MAT
68 : !> ATOMS
69 : !> EOF
70 : !>
71 : !> \author Teodoro Laino [tlaino]
72 : ! **************************************************************************************************
73 6 : SUBROUTINE read_coordinate_xtl(topology, para_env, subsys_section)
74 : TYPE(topology_parameters_type) :: topology
75 : TYPE(mp_para_env_type), POINTER :: para_env
76 : TYPE(section_vals_type), POINTER :: subsys_section
77 :
78 : CHARACTER(len=*), PARAMETER :: routineN = 'read_coordinate_xtl'
79 : INTEGER, PARAMETER :: nblock = 1000
80 : REAL(KIND=dp), PARAMETER :: threshold = 1.0E-6_dp
81 :
82 : CHARACTER(LEN=default_string_length) :: strtmp
83 : INTEGER :: dimensions, handle, icol, ii, isym, iw, &
84 : jj, natom, natom_orig, newsize
85 : INTEGER, DIMENSION(3) :: periodic
86 : LOGICAL :: check, found, my_end
87 : REAL(KIND=dp) :: pfactor, threshold2
88 : REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths, r, r1, r2, s, &
89 : transl_vec
90 : REAL(KIND=dp), DIMENSION(3, 3) :: rot_mat
91 : TYPE(atom_info_type), POINTER :: atom_info
92 : TYPE(cell_type), POINTER :: cell
93 : TYPE(cp_logger_type), POINTER :: logger
94 : TYPE(cp_parser_type) :: parser
95 :
96 2 : NULLIFY (logger)
97 4 : logger => cp_get_default_logger()
98 : iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XTL_INFO", &
99 2 : extension=".subsysLog")
100 2 : CALL timeset(routineN, handle)
101 :
102 2 : pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
103 : ! Element is assigned on the basis of the atm_name
104 2 : topology%aa_element = .TRUE.
105 :
106 2 : atom_info => topology%atom_info
107 2 : CALL reallocate(atom_info%id_molname, 1, nblock)
108 2 : CALL reallocate(atom_info%id_resname, 1, nblock)
109 2 : CALL reallocate(atom_info%resid, 1, nblock)
110 2 : CALL reallocate(atom_info%id_atmname, 1, nblock)
111 2 : CALL reallocate(atom_info%r, 1, 3, 1, nblock)
112 2 : CALL reallocate(atom_info%atm_mass, 1, nblock)
113 2 : CALL reallocate(atom_info%atm_charge, 1, nblock)
114 2 : CALL reallocate(atom_info%occup, 1, nblock)
115 2 : CALL reallocate(atom_info%beta, 1, nblock)
116 2 : CALL reallocate(atom_info%id_element, 1, nblock)
117 :
118 2 : IF (iw > 0) WRITE (iw, *) " Reading in XTL file ", TRIM(topology%coord_file_name)
119 2 : CALL parser_create(parser, topology%coord_file_name, para_env=para_env)
120 :
121 : ! Check for TITLE
122 : CALL parser_search_string(parser, "TITLE", ignore_case=.FALSE., found=found, &
123 2 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
124 2 : IF (found) THEN
125 2 : IF (iw > 0) WRITE (iw, '(/,A)') " XTL_INFO| TITLE :: "//TRIM(parser%input_line(parser%icol:))
126 : END IF
127 :
128 : ! Check for _chemical_formula_sum
129 : CALL parser_search_string(parser, "DIMENSION", ignore_case=.FALSE., found=found, &
130 2 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
131 2 : IF (found) THEN
132 2 : IF (iw > 0) WRITE (iw, '(A)') " XTL_INFO| DIMENSION :: "//TRIM(parser%input_line(parser%icol:))
133 2 : CALL parser_get_object(parser, dimensions)
134 2 : IF (dimensions /= 3) THEN
135 0 : CPABORT("XTL file with working DIMENSION different from 3 cannot be parsed!")
136 : END IF
137 : ELSE
138 : ! Assuming by default we work in 3D-periodic systems
139 0 : dimensions = 3
140 : END IF
141 :
142 : ! Parsing cell infos
143 8 : periodic = 1
144 : ! Check for _cell_length_a
145 : CALL parser_search_string(parser, "CELL", ignore_case=.FALSE., found=found, &
146 2 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
147 2 : IF (.NOT. found) THEN
148 0 : CPABORT("The field CELL was not found in XTL file! ")
149 : END IF
150 2 : CALL parser_get_next_line(parser, 1)
151 : ! CELL LENGTH A
152 2 : CALL parser_get_object(parser, cell_lengths(1))
153 2 : cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
154 : ! CELL LENGTH B
155 2 : CALL parser_get_object(parser, cell_lengths(2))
156 2 : cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
157 : ! CELL LENGTH C
158 2 : CALL parser_get_object(parser, cell_lengths(3))
159 2 : cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
160 :
161 : ! CELL ANGLE ALPHA
162 2 : CALL parser_get_object(parser, cell_angles(1))
163 2 : cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
164 : ! CELL ANGLE BETA
165 2 : CALL parser_get_object(parser, cell_angles(2))
166 2 : cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
167 : ! CELL ANGLE GAMMA
168 2 : CALL parser_get_object(parser, cell_angles(3))
169 2 : cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
170 :
171 : ! Create cell
172 2 : NULLIFY (cell)
173 2 : CALL cell_create(cell, tag="CELL_XTL")
174 : CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
175 2 : do_init_cell=.TRUE.)
176 2 : CALL write_cell(cell, subsys_section)
177 :
178 : ! Parse atoms info and fractional coordinates
179 : ! Check for _atom_site_label
180 : CALL parser_search_string(parser, "ATOMS", ignore_case=.FALSE., found=found, &
181 2 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
182 2 : IF (.NOT. found) THEN
183 0 : CPABORT("The field ATOMS was not found in XTL file! ")
184 : END IF
185 2 : CALL parser_get_next_line(parser, 1)
186 : ! Paranoic syntax check.. if this fails one should improve the description of XTL files
187 2 : found = (INDEX(parser%input_line, "NAME X Y Z") /= 0)
188 2 : IF (.NOT. found) THEN
189 0 : CPABORT("The field ATOMS in XTL file, is not followed by name and coordinates tags! ")
190 : END IF
191 2 : CALL parser_get_next_line(parser, 1)
192 : ! Parse real info
193 2 : natom = 0
194 62 : DO WHILE (INDEX(parser%input_line, "EOF") == 0)
195 60 : natom = natom + 1
196 : ! Resize in case needed
197 60 : IF (natom > SIZE(atom_info%id_molname)) THEN
198 0 : newsize = INT(pfactor*natom)
199 0 : CALL reallocate(atom_info%id_molname, 1, newsize)
200 0 : CALL reallocate(atom_info%id_resname, 1, newsize)
201 0 : CALL reallocate(atom_info%resid, 1, newsize)
202 0 : CALL reallocate(atom_info%id_atmname, 1, newsize)
203 0 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
204 0 : CALL reallocate(atom_info%atm_mass, 1, newsize)
205 0 : CALL reallocate(atom_info%atm_charge, 1, newsize)
206 0 : CALL reallocate(atom_info%occup, 1, newsize)
207 0 : CALL reallocate(atom_info%beta, 1, newsize)
208 0 : CALL reallocate(atom_info%id_element, 1, newsize)
209 : END IF
210 : ! NAME
211 60 : CALL parser_get_object(parser, strtmp)
212 60 : atom_info%id_atmname(natom) = str2id(strtmp)
213 60 : atom_info%id_molname(natom) = str2id(s2s("MOL"//TRIM(ADJUSTL(cp_to_string(natom)))))
214 60 : atom_info%id_resname(natom) = atom_info%id_molname(natom)
215 60 : atom_info%resid(natom) = 1
216 60 : atom_info%id_element(natom) = atom_info%id_atmname(natom)
217 : ! X
218 60 : CALL parser_get_object(parser, atom_info%r(1, natom))
219 : ! Y
220 60 : CALL parser_get_object(parser, atom_info%r(2, natom))
221 : ! Z
222 60 : CALL parser_get_object(parser, atom_info%r(3, natom))
223 240 : s = atom_info%r(1:3, natom)
224 60 : CALL scaled_to_real(atom_info%r(1:3, natom), s, cell)
225 60 : CALL parser_get_next_line(parser, 1, at_end=my_end)
226 62 : IF (my_end) EXIT
227 : END DO
228 : !
229 2 : threshold2 = threshold*threshold
230 : ! Preliminary check: check if atoms provided are really unique.. this is a paranoic
231 : ! check since they should be REALLY unique.. anyway..
232 62 : DO ii = 1, natom
233 240 : r1 = atom_info%r(1:3, ii)
234 932 : DO jj = ii + 1, natom
235 3480 : r2 = atom_info%r(1:3, jj)
236 3480 : r = pbc(r1 - r2, cell)
237 : ! SQRT(DOT_PRODUCT(r, r)) >= threshold
238 3480 : check = (DOT_PRODUCT(r, r) >= threshold2)
239 930 : CPASSERT(check)
240 : END DO
241 : END DO
242 : ! Parse Symmetry Group and generation elements..
243 : ! Check for SYMMETRY
244 : CALL parser_search_string(parser, "SYMMETRY", ignore_case=.FALSE., found=found, &
245 2 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
246 2 : IF (found) THEN
247 2 : IF (iw > 0) WRITE (iw, '(A)') " XTL_INFO| Symmetry Infos :: "//TRIM(parser%input_line(parser%icol:))
248 : END IF
249 :
250 : ! Check for SYM MAT
251 : CALL parser_search_string(parser, "SYM MAT", ignore_case=.FALSE., found=found, &
252 2 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
253 2 : CPWARN_IF(.NOT. found, "The field SYM MAT was not found in XTL file! ")
254 2 : IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of atoms before applying symmetry operations :: ", natom
255 122 : 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)
256 2 : IF (found) THEN
257 : ! Apply symmetry elements and generate the whole set of atoms in the unit cell
258 2 : isym = 0
259 2 : natom_orig = natom
260 4 : DO WHILE (found)
261 2 : isym = isym + 1
262 2 : icol = INDEX(parser%input_line, "SYM MAT") + 8
263 8 : READ (parser%input_line(icol:), *) ((rot_mat(ii, jj), jj=1, 3), ii=1, 3), transl_vec(1:3)
264 62 : Loop_over_unique_atoms: DO ii = 1, natom_orig
265 : ! Rotate and apply translation
266 960 : r1 = MATMUL(rot_mat, atom_info%r(1:3, ii)) + transl_vec
267 : ! Verify if this atom is really unique..
268 60 : check = .TRUE.
269 930 : DO jj = 1, natom
270 3720 : r2 = atom_info%r(1:3, jj)
271 3720 : r = pbc(r1 - r2, cell)
272 : ! SQRT(DOT_PRODUCT(r, r)) <= threshold
273 3720 : IF (DOT_PRODUCT(r, r) <= threshold2) THEN
274 : check = .FALSE.
275 : EXIT
276 : END IF
277 : END DO
278 : ! If the atom generated is unique let's add to the atom set..
279 62 : IF (check) THEN
280 0 : natom = natom + 1
281 : ! Resize in case needed
282 0 : IF (natom > SIZE(atom_info%id_molname)) THEN
283 0 : newsize = INT(pfactor*natom)
284 0 : CALL reallocate(atom_info%id_molname, 1, newsize)
285 0 : CALL reallocate(atom_info%id_resname, 1, newsize)
286 0 : CALL reallocate(atom_info%resid, 1, newsize)
287 0 : CALL reallocate(atom_info%id_atmname, 1, newsize)
288 0 : CALL reallocate(atom_info%r, 1, 3, 1, newsize)
289 0 : CALL reallocate(atom_info%atm_mass, 1, newsize)
290 0 : CALL reallocate(atom_info%atm_charge, 1, newsize)
291 0 : CALL reallocate(atom_info%occup, 1, newsize)
292 0 : CALL reallocate(atom_info%beta, 1, newsize)
293 0 : CALL reallocate(atom_info%id_element, 1, newsize)
294 : END IF
295 0 : atom_info%id_atmname(natom) = atom_info%id_atmname(ii)
296 0 : atom_info%id_molname(natom) = atom_info%id_molname(ii)
297 0 : atom_info%id_resname(natom) = atom_info%id_resname(ii)
298 0 : atom_info%resid(natom) = atom_info%resid(ii)
299 0 : atom_info%id_element(natom) = atom_info%id_element(ii)
300 0 : atom_info%r(1:3, natom) = r1
301 : END IF
302 : END DO Loop_over_unique_atoms
303 : CALL parser_search_string(parser, "SYM MAT", ignore_case=.FALSE., found=found, &
304 4 : begin_line=.FALSE., search_from_begin_of_file=.FALSE.)
305 : END DO
306 : END IF
307 2 : IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of symmetry operations :: ", isym
308 2 : IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of total atoms :: ", natom
309 122 : 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)
310 :
311 : ! Releasse local cell type and parser
312 2 : CALL cell_release(cell)
313 2 : CALL parser_release(parser)
314 :
315 : ! Reallocate all structures with the exact NATOM size
316 2 : CALL reallocate(atom_info%id_molname, 1, natom)
317 2 : CALL reallocate(atom_info%id_resname, 1, natom)
318 2 : CALL reallocate(atom_info%resid, 1, natom)
319 2 : CALL reallocate(atom_info%id_atmname, 1, natom)
320 2 : CALL reallocate(atom_info%r, 1, 3, 1, natom)
321 2 : CALL reallocate(atom_info%atm_mass, 1, natom)
322 2 : CALL reallocate(atom_info%atm_charge, 1, natom)
323 2 : CALL reallocate(atom_info%occup, 1, natom)
324 2 : CALL reallocate(atom_info%beta, 1, natom)
325 2 : CALL reallocate(atom_info%id_element, 1, natom)
326 :
327 2 : topology%natoms = natom
328 2 : topology%molname_generated = .TRUE.
329 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
330 2 : "PRINT%TOPOLOGY_INFO/XTL_INFO")
331 2 : CALL timestop(handle)
332 6 : END SUBROUTINE read_coordinate_xtl
333 :
334 : END MODULE topology_xtl
|
