(git:34ef472)
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 ! **************************************************************************************************
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  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  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 
332 END MODULE topology_xtl
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
cell_methods
Handles all functions related to the CELL.
Definition: cell_methods.F:15
cell_methods::write_cell
subroutine, public write_cell(cell, subsys_section, tag)
Write the cell parameters to the output unit.
Definition: cell_methods.F:731
cell_methods::set_cell_param
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...
Definition: cell_methods.F:439
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition: cell_methods.F:85
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition: cell_types.F:516
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition: cell_types.F:559
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
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)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
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,...
Definition: cp_output_handling.F:1013
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
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 ...
Definition: cp_parser_methods.F:421
cp_parser_methods::parser_search_string
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 ...
Definition: cp_parser_methods.F:786
cp_parser_types
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition: cp_parser_types.F:21
cp_parser_types::parser_release
subroutine, public parser_release(parser)
releases the parser
Definition: cp_parser_types.F:109
cp_parser_types::parser_create
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.
Definition: cp_parser_types.F:147
cp_units
unit conversion facility
Definition: cp_units.F:30
cp_units::cp_unit_to_cp2k
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
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_get_rval
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Definition: input_section_types.F:941
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
string_table
generates a unique id number for a string (str2id) that can be used two compare two strings....
Definition: string_table.F:22
string_table::s2s
character(len=default_string_length) function, public s2s(str)
converts a string in a string of default_string_length
Definition: string_table.F:141
string_table::str2id
integer function, public str2id(str)
returns a unique id for a given string, and stores the string for later retrieval using the id.
Definition: string_table.F:72
string_table::id2str
character(len=default_string_length) function, public id2str(id)
returns the string associated with a given id
Definition: string_table.F:115
topology_types
Definition: topology_types.F:18
topology_xtl
Handles XTL (Molecular Simulations, Inc (MSI)) files.
Definition: topology_xtl.F:13
topology_xtl::read_coordinate_xtl
subroutine, public read_coordinate_xtl(topology, para_env, subsys_section)
Performs the real task of reading the proper information from the XTL file.
Definition: topology_xtl.F:74
topology
Control for reading in different topologies and coordinates.
Definition: topology.F:13
topology_types::atom_info_type
Definition: topology_types.F:35
topology_types::topology_parameters_type
Definition: topology_types.F:142