(git:374b731)
Loading...
Searching...
No Matches
topology_xtl.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Handles XTL (Molecular Simulations, Inc (MSI)) files
10!> \author Teodoro Laino [tlaino]
11!> \date 05.2009
12! **************************************************************************************************
14 USE cell_methods, ONLY: cell_create,&
17 USE cell_types, ONLY: cell_release,&
18 cell_type,&
19 pbc,&
32 USE cp_units, ONLY: cp_unit_to_cp2k
35 USE kinds, ONLY: default_string_length,&
36 dp
39 USE string_table, ONLY: id2str,&
40 s2s,&
41 str2id
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
53CONTAINS
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 SUBROUTINE read_coordinate_xtl(topology, para_env, subsys_section)
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 NULLIFY (logger)
97 logger => cp_get_default_logger()
98 iw = cp_print_key_unit_nr(logger, subsys_section, "PRINT%TOPOLOGY_INFO/XTL_INFO", &
99 extension=".subsysLog")
100 CALL timeset(routinen, handle)
101
102 pfactor = section_get_rval(subsys_section, "TOPOLOGY%MEMORY_PROGRESSION_FACTOR")
103 ! Element is assigned on the basis of the atm_name
104 topology%aa_element = .true.
105
106 atom_info => topology%atom_info
107 CALL reallocate(atom_info%id_molname, 1, nblock)
108 CALL reallocate(atom_info%id_resname, 1, nblock)
109 CALL reallocate(atom_info%resid, 1, nblock)
110 CALL reallocate(atom_info%id_atmname, 1, nblock)
111 CALL reallocate(atom_info%r, 1, 3, 1, nblock)
112 CALL reallocate(atom_info%atm_mass, 1, nblock)
113 CALL reallocate(atom_info%atm_charge, 1, nblock)
114 CALL reallocate(atom_info%occup, 1, nblock)
115 CALL reallocate(atom_info%beta, 1, nblock)
116 CALL reallocate(atom_info%id_element, 1, nblock)
117
118 IF (iw > 0) WRITE (iw, *) " Reading in XTL file ", trim(topology%coord_file_name)
119 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 begin_line=.false., search_from_begin_of_file=.true.)
124 IF (found) THEN
125 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 begin_line=.false., search_from_begin_of_file=.true.)
131 IF (found) THEN
132 IF (iw > 0) WRITE (iw, '(A)') " XTL_INFO| DIMENSION :: "//trim(parser%input_line(parser%icol:))
133 CALL parser_get_object(parser, dimensions)
134 IF (dimensions /= 3) THEN
135 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 dimensions = 3
140 END IF
141
142 ! Parsing cell infos
143 periodic = 1
144 ! Check for _cell_length_a
145 CALL parser_search_string(parser, "CELL", ignore_case=.false., found=found, &
146 begin_line=.false., search_from_begin_of_file=.true.)
147 IF (.NOT. found) &
148 cpabort("The field CELL was not found in XTL file! ")
149 CALL parser_get_next_line(parser, 1)
150 ! CELL LENGTH A
151 CALL parser_get_object(parser, cell_lengths(1))
152 cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
153 ! CELL LENGTH B
154 CALL parser_get_object(parser, cell_lengths(2))
155 cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
156 ! CELL LENGTH C
157 CALL parser_get_object(parser, cell_lengths(3))
158 cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
159
160 ! CELL ANGLE ALPHA
161 CALL parser_get_object(parser, cell_angles(1))
162 cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
163 ! CELL ANGLE BETA
164 CALL parser_get_object(parser, cell_angles(2))
165 cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
166 ! CELL ANGLE GAMMA
167 CALL parser_get_object(parser, cell_angles(3))
168 cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
169
170 ! Create cell
171 NULLIFY (cell)
172 CALL cell_create(cell, tag="CELL_XTL")
173 CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
174 do_init_cell=.true.)
175 CALL write_cell(cell, subsys_section)
176
177 ! Parse atoms info and fractional coordinates
178 ! Check for _atom_site_label
179 CALL parser_search_string(parser, "ATOMS", ignore_case=.false., found=found, &
180 begin_line=.false., search_from_begin_of_file=.true.)
181 IF (.NOT. found) &
182 cpabort("The field ATOMS was not found in XTL file! ")
183 CALL parser_get_next_line(parser, 1)
184 ! Paranoic syntax check.. if this fails one should improve the description of XTL files
185 found = (index(parser%input_line, "NAME X Y Z") /= 0)
186 IF (.NOT. found) &
187 cpabort("The field ATOMS in XTL file, is not followed by name and coordinates tags! ")
188 CALL parser_get_next_line(parser, 1)
189 ! Parse real info
190 natom = 0
191 DO WHILE (index(parser%input_line, "EOF") == 0)
192 natom = natom + 1
193 ! Resize in case needed
194 IF (natom > SIZE(atom_info%id_molname)) THEN
195 newsize = int(pfactor*natom)
196 CALL reallocate(atom_info%id_molname, 1, newsize)
197 CALL reallocate(atom_info%id_resname, 1, newsize)
198 CALL reallocate(atom_info%resid, 1, newsize)
199 CALL reallocate(atom_info%id_atmname, 1, newsize)
200 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
201 CALL reallocate(atom_info%atm_mass, 1, newsize)
202 CALL reallocate(atom_info%atm_charge, 1, newsize)
203 CALL reallocate(atom_info%occup, 1, newsize)
204 CALL reallocate(atom_info%beta, 1, newsize)
205 CALL reallocate(atom_info%id_element, 1, newsize)
206 END IF
207 ! NAME
208 CALL parser_get_object(parser, strtmp)
209 atom_info%id_atmname(natom) = str2id(strtmp)
210 atom_info%id_molname(natom) = str2id(s2s("MOL"//trim(adjustl(cp_to_string(natom)))))
211 atom_info%id_resname(natom) = atom_info%id_molname(natom)
212 atom_info%resid(natom) = 1
213 atom_info%id_element(natom) = atom_info%id_atmname(natom)
214 ! X
215 CALL parser_get_object(parser, atom_info%r(1, natom))
216 ! Y
217 CALL parser_get_object(parser, atom_info%r(2, natom))
218 ! Z
219 CALL parser_get_object(parser, atom_info%r(3, natom))
220 s = atom_info%r(1:3, natom)
221 CALL scaled_to_real(atom_info%r(1:3, natom), s, cell)
222 CALL parser_get_next_line(parser, 1, at_end=my_end)
223 IF (my_end) EXIT
224 END DO
225 !
226 threshold2 = threshold*threshold
227 ! Preliminary check: check if atoms provided are really unique.. this is a paranoic
228 ! check since they should be REALLY unique.. anyway..
229 DO ii = 1, natom
230 r1 = atom_info%r(1:3, ii)
231 DO jj = ii + 1, natom
232 r2 = atom_info%r(1:3, jj)
233 r = pbc(r1 - r2, cell)
234 ! SQRT(DOT_PRODUCT(r, r)) >= threshold
235 check = (dot_product(r, r) >= threshold2)
236 cpassert(check)
237 END DO
238 END DO
239 ! Parse Symmetry Group and generation elements..
240 ! Check for SYMMETRY
241 CALL parser_search_string(parser, "SYMMETRY", ignore_case=.false., found=found, &
242 begin_line=.false., search_from_begin_of_file=.true.)
243 IF (found) THEN
244 IF (iw > 0) WRITE (iw, '(A)') " XTL_INFO| Symmetry Infos :: "//trim(parser%input_line(parser%icol:))
245 END IF
246
247 ! Check for SYM MAT
248 CALL parser_search_string(parser, "SYM MAT", ignore_case=.false., found=found, &
249 begin_line=.false., search_from_begin_of_file=.true.)
250 IF (.NOT. found) &
251 cpwarn("The field SYM MAT was not found in XTL file! ")
252 IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of atoms before applying symmetry operations :: ", natom
253 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)
254 IF (found) THEN
255 ! Apply symmetry elements and generate the whole set of atoms in the unit cell
256 isym = 0
257 natom_orig = natom
258 DO WHILE (found)
259 isym = isym + 1
260 icol = index(parser%input_line, "SYM MAT") + 8
261 READ (parser%input_line(icol:), *) ((rot_mat(ii, jj), jj=1, 3), ii=1, 3), transl_vec(1:3)
262 loop_over_unique_atoms: DO ii = 1, natom_orig
263 ! Rotate and apply translation
264 r1 = matmul(rot_mat, atom_info%r(1:3, ii)) + transl_vec
265 ! Verify if this atom is really unique..
266 check = .true.
267 DO jj = 1, natom
268 r2 = atom_info%r(1:3, jj)
269 r = pbc(r1 - r2, cell)
270 ! SQRT(DOT_PRODUCT(r, r)) <= threshold
271 IF (dot_product(r, r) <= threshold2) THEN
272 check = .false.
273 EXIT
274 END IF
275 END DO
276 ! If the atom generated is unique let's add to the atom set..
277 IF (check) THEN
278 natom = natom + 1
279 ! Resize in case needed
280 IF (natom > SIZE(atom_info%id_molname)) THEN
281 newsize = int(pfactor*natom)
282 CALL reallocate(atom_info%id_molname, 1, newsize)
283 CALL reallocate(atom_info%id_resname, 1, newsize)
284 CALL reallocate(atom_info%resid, 1, newsize)
285 CALL reallocate(atom_info%id_atmname, 1, newsize)
286 CALL reallocate(atom_info%r, 1, 3, 1, newsize)
287 CALL reallocate(atom_info%atm_mass, 1, newsize)
288 CALL reallocate(atom_info%atm_charge, 1, newsize)
289 CALL reallocate(atom_info%occup, 1, newsize)
290 CALL reallocate(atom_info%beta, 1, newsize)
291 CALL reallocate(atom_info%id_element, 1, newsize)
292 END IF
293 atom_info%id_atmname(natom) = atom_info%id_atmname(ii)
294 atom_info%id_molname(natom) = atom_info%id_molname(ii)
295 atom_info%id_resname(natom) = atom_info%id_resname(ii)
296 atom_info%resid(natom) = atom_info%resid(ii)
297 atom_info%id_element(natom) = atom_info%id_element(ii)
298 atom_info%r(1:3, natom) = r1
299 END IF
300 END DO loop_over_unique_atoms
301 CALL parser_search_string(parser, "SYM MAT", ignore_case=.false., found=found, &
302 begin_line=.false., search_from_begin_of_file=.false.)
303 END DO
304 END IF
305 IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of symmetry operations :: ", isym
306 IF (iw > 0) WRITE (iw, '(A,I0)') " XTL_INFO| Number of total atoms :: ", natom
307 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)
308
309 ! Releasse local cell type and parser
310 CALL cell_release(cell)
311 CALL parser_release(parser)
312
313 ! Reallocate all structures with the exact NATOM size
314 CALL reallocate(atom_info%id_molname, 1, natom)
315 CALL reallocate(atom_info%id_resname, 1, natom)
316 CALL reallocate(atom_info%resid, 1, natom)
317 CALL reallocate(atom_info%id_atmname, 1, natom)
318 CALL reallocate(atom_info%r, 1, 3, 1, natom)
319 CALL reallocate(atom_info%atm_mass, 1, natom)
320 CALL reallocate(atom_info%atm_charge, 1, natom)
321 CALL reallocate(atom_info%occup, 1, natom)
322 CALL reallocate(atom_info%beta, 1, natom)
323 CALL reallocate(atom_info%id_element, 1, natom)
324
325 topology%natoms = natom
326 topology%molname_generated = .true.
327 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
328 "PRINT%TOPOLOGY_INFO/XTL_INFO")
329 CALL timestop(handle)
330 END SUBROUTINE read_coordinate_xtl
331
332END MODULE topology_xtl
Handles all functions related to the CELL.
subroutine, public write_cell(cell, subsys_section, tag)
Write the cell parameters to the output unit.
subroutine, public set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma) using the convention: a parall...
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:516
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:559
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1150
objects that represent the structure of input sections and the data contained in an input section
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Utility routines for the memory handling.
Interface to the message passing library MPI.
generates a unique id number for a string (str2id) that can be used two compare two strings....
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Handles XTL (Molecular Simulations, Inc (MSI)) files.
subroutine, public read_coordinate_xtl(topology, para_env, subsys_section)
Performs the real task of reading the proper information from the XTL file.
Control for reading in different topologies and coordinates.
Definition topology.F:13
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores all the informations relevant to an mpi environment