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 all functions related to the CELL
10 : !> \par History
11 : !> 11.2008 Teodoro Laino [tlaino] - deeply cleaning cell_type from units
12 : !> 10.2014 Moved many routines to cell_types.F.
13 : !> \author Matthias KracK (16.01.2002, based on a earlier version of CJM, JGH)
14 : ! **************************************************************************************************
15 : MODULE cell_methods
16 : USE cell_types, ONLY: &
17 : cell_clone, cell_release, cell_sym_cubic, cell_sym_hexagonal_gamma_120, &
18 : cell_sym_hexagonal_gamma_60, cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, &
19 : cell_sym_none, cell_sym_orthorhombic, cell_sym_rhombohedral, cell_sym_tetragonal_ab, &
20 : cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type, get_cell, &
21 : plane_distance, use_perd_none, use_perd_x, use_perd_xy, use_perd_xyz, use_perd_xz, &
22 : use_perd_y, use_perd_yz, use_perd_z
23 : USE cp_log_handling, ONLY: cp_get_default_logger,&
24 : cp_logger_type
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 cp_units, ONLY: cp_unit_from_cp2k,&
34 : cp_unit_to_cp2k
35 : USE input_constants, ONLY: &
36 : canonicalize_cell_auto, canonicalize_cell_true, do_cell_cif, do_cell_cp2k, do_cell_extxyz, &
37 : do_cell_pdb, do_cell_xsc, do_coord_cif, do_coord_cp2k, do_coord_pdb, do_coord_xyz
38 : USE input_cp2k_subsys, ONLY: create_cell_section
39 : USE input_enumeration_types, ONLY: enum_i2c,&
40 : enumeration_type
41 : USE input_keyword_types, ONLY: keyword_get,&
42 : keyword_type
43 : USE input_section_types, ONLY: &
44 : section_get_keyword, section_release, section_type, section_vals_get, &
45 : section_vals_get_subs_vals, section_vals_type, section_vals_val_get, section_vals_val_set, &
46 : section_vals_val_unset
47 : USE kinds, ONLY: default_path_length,&
48 : default_string_length,&
49 : dp,&
50 : max_line_length
51 : USE machine, ONLY: default_output_unit
52 : USE mathconstants, ONLY: degree,&
53 : sqrt3
54 : USE mathlib, ONLY: angle,&
55 : det_3x3,&
56 : inv_3x3
57 : USE message_passing, ONLY: mp_para_env_type
58 : USE string_utilities, ONLY: uppercase
59 : #include "./base/base_uses.f90"
60 :
61 : IMPLICIT NONE
62 :
63 : PRIVATE
64 :
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_methods'
66 :
67 : PUBLIC :: cell_create, &
68 : cell_finalize_canonical_input, &
69 : init_cell, &
70 : read_cell, &
71 : read_cell_cif, &
72 : read_cell_cp2k, &
73 : read_cell_xyz, &
74 : read_cell_pdb, &
75 : read_cell_xsc, &
76 : read_xyz_comment, &
77 : set_cell_param, &
78 : write_cell, &
79 : write_cell_low
80 :
81 : CONTAINS
82 :
83 : ! **************************************************************************************************
84 : !> \brief allocates and initializes a cell
85 : !> \param cell the cell to initialize
86 : !> \param hmat the h matrix that defines the cell
87 : !> \param periodic periodicity of the cell
88 : !> \param tag ...
89 : !> \par History
90 : !> 09.2003 created [fawzi]
91 : !> \author Fawzi Mohamed
92 : ! **************************************************************************************************
93 79754 : SUBROUTINE cell_create(cell, hmat, periodic, tag)
94 :
95 : TYPE(cell_type), POINTER :: cell
96 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
97 : OPTIONAL :: hmat
98 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
99 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
100 :
101 0 : CPASSERT(.NOT. ASSOCIATED(cell))
102 5503026 : ALLOCATE (cell)
103 79754 : cell%ref_count = 1
104 79754 : IF (PRESENT(periodic)) THEN
105 1372 : cell%perd = periodic
106 : ELSE
107 317644 : cell%perd = 1
108 : END IF
109 79754 : cell%orthorhombic = .FALSE.
110 79754 : cell%input_cell_canonicalized = .FALSE.
111 1036802 : cell%input_hmat(:, :) = 0.0_dp
112 1036802 : cell%input_to_canonical(:, :) = 0.0_dp
113 1036802 : cell%input_recip_to_canonical(:, :) = 0.0_dp
114 79754 : cell%symmetry_id = cell_sym_none
115 79754 : IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
116 79754 : IF (PRESENT(tag)) cell%tag = tag
117 :
118 79754 : END SUBROUTINE cell_create
119 :
120 : ! **************************************************************************************************
121 : !> \brief Store the transform between the user input cell and the canonical cell.
122 : !> \param cell ...
123 : !> \param hmat_input ...
124 : !> \param hmat_canonical ...
125 : ! **************************************************************************************************
126 106 : SUBROUTINE cell_finalize_canonical_input(cell, hmat_input, hmat_canonical)
127 :
128 : TYPE(cell_type), POINTER :: cell
129 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: hmat_input, hmat_canonical
130 :
131 : REAL(KIND=dp), PARAMETER :: eps_hmat = 1.0E-12_dp
132 :
133 : REAL(KIND=dp), DIMENSION(3, 3) :: tmat
134 :
135 106 : CPASSERT(ASSOCIATED(cell))
136 :
137 1378 : IF (MAXVAL(ABS(hmat_canonical - hmat_input)) <= eps_hmat) THEN
138 14 : cell%input_cell_canonicalized = .FALSE.
139 182 : cell%input_hmat(:, :) = 0.0_dp
140 182 : cell%input_to_canonical(:, :) = 0.0_dp
141 182 : cell%input_recip_to_canonical(:, :) = 0.0_dp
142 : ELSE
143 6164 : tmat = MATMUL(hmat_canonical, inv_3x3(hmat_input))
144 92 : cell%input_cell_canonicalized = .TRUE.
145 1196 : cell%input_hmat(:, :) = hmat_input(:, :)
146 1196 : cell%input_to_canonical(:, :) = tmat(:, :)
147 1196 : cell%input_recip_to_canonical(:, :) = TRANSPOSE(inv_3x3(tmat))
148 : END IF
149 :
150 106 : END SUBROUTINE cell_finalize_canonical_input
151 :
152 : ! **************************************************************************************************
153 : !> \brief Canonicalize a general cell matrix without changing lengths and angles.
154 : !> \param cell ...
155 : ! **************************************************************************************************
156 106 : SUBROUTINE canonicalize_cell_matrix(cell)
157 :
158 : TYPE(cell_type), POINTER :: cell
159 :
160 : REAL(KIND=dp), DIMENSION(3) :: abc, cell_angle
161 :
162 106 : CPASSERT(ASSOCIATED(cell))
163 :
164 106 : CALL get_cell(cell=cell, abc=abc)
165 106 : cell_angle(1) = angle(cell%hmat(:, 2), cell%hmat(:, 3))
166 106 : cell_angle(2) = angle(cell%hmat(:, 1), cell%hmat(:, 3))
167 106 : cell_angle(3) = angle(cell%hmat(:, 1), cell%hmat(:, 2))
168 :
169 : CALL set_cell_param(cell, cell_length=abc, cell_angle=cell_angle, &
170 106 : periodic=cell%perd, do_init_cell=.TRUE.)
171 :
172 106 : END SUBROUTINE canonicalize_cell_matrix
173 :
174 : ! **************************************************************************************************
175 : !> \brief Initialise/readjust a simulation cell after hmat has been changed
176 : !> \param cell ...
177 : !> \param hmat ...
178 : !> \param periodic ...
179 : !> \date 16.01.2002
180 : !> \author Matthias Krack
181 : !> \version 1.0
182 : ! **************************************************************************************************
183 340137 : SUBROUTINE init_cell(cell, hmat, periodic)
184 :
185 : TYPE(cell_type), POINTER :: cell
186 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
187 : OPTIONAL :: hmat
188 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
189 :
190 : REAL(KIND=dp), PARAMETER :: eps_hmat = 1.0E-14_dp
191 :
192 : INTEGER :: dim
193 : REAL(KIND=dp) :: a, acosa, acosah, acosg, alpha, asina, &
194 : asinah, asing, beta, gamma, norm, &
195 : norm_c
196 : REAL(KIND=dp), DIMENSION(3) :: abc
197 :
198 340137 : CPASSERT(ASSOCIATED(cell))
199 :
200 500589 : IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
201 340539 : IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
202 :
203 340137 : cell%deth = ABS(det_3x3(cell%hmat))
204 :
205 340137 : IF (cell%deth < 1.0E-10_dp) THEN
206 0 : CALL write_cell_low(cell, "angstrom", default_output_unit)
207 : CALL cp_abort(__LOCATION__, &
208 : "An invalid set of cell vectors was specified. "// &
209 0 : "The cell volume is too small")
210 : END IF
211 :
212 344220 : SELECT CASE (cell%symmetry_id)
213 : CASE (cell_sym_cubic, &
214 : cell_sym_tetragonal_ab, &
215 : cell_sym_tetragonal_ac, &
216 : cell_sym_tetragonal_bc, &
217 : cell_sym_orthorhombic)
218 4083 : CALL get_cell(cell=cell, abc=abc)
219 4083 : abc(2) = plane_distance(0, 1, 0, cell=cell)
220 4083 : abc(3) = plane_distance(0, 0, 1, cell=cell)
221 4083 : SELECT CASE (cell%symmetry_id)
222 : CASE (cell_sym_cubic)
223 5565 : abc(1:3) = SUM(abc(1:3))/3.0_dp
224 : CASE (cell_sym_tetragonal_ab, &
225 : cell_sym_tetragonal_ac, &
226 : cell_sym_tetragonal_bc)
227 4083 : SELECT CASE (cell%symmetry_id)
228 : CASE (cell_sym_tetragonal_ab)
229 1318 : a = 0.5_dp*(abc(1) + abc(2))
230 1318 : abc(1) = a
231 1318 : abc(2) = a
232 : CASE (cell_sym_tetragonal_ac)
233 631 : a = 0.5_dp*(abc(1) + abc(3))
234 631 : abc(1) = a
235 631 : abc(3) = a
236 : CASE (cell_sym_tetragonal_bc)
237 610 : a = 0.5_dp*(abc(2) + abc(3))
238 610 : abc(2) = a
239 2559 : abc(3) = a
240 : END SELECT
241 : END SELECT
242 4083 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
243 4083 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
244 4083 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
245 : CASE (cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120)
246 2630 : CALL get_cell(cell=cell, abc=abc)
247 2630 : a = 0.5_dp*(abc(1) + abc(2))
248 2630 : acosg = 0.5_dp*a
249 2630 : asing = sqrt3*acosg
250 2630 : IF (cell%symmetry_id == cell_sym_hexagonal_gamma_120) acosg = -acosg
251 2630 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
252 2630 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
253 2630 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
254 : CASE (cell_sym_rhombohedral)
255 649 : CALL get_cell(cell=cell, abc=abc)
256 2596 : a = SUM(abc(1:3))/3.0_dp
257 : alpha = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
258 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
259 649 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
260 649 : acosa = a*COS(alpha)
261 649 : asina = a*SIN(alpha)
262 649 : acosah = a*COS(0.5_dp*alpha)
263 649 : asinah = a*SIN(0.5_dp*alpha)
264 649 : norm = acosa/acosah
265 649 : norm_c = SQRT(1.0_dp - norm*norm)
266 649 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
267 649 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
268 649 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = a*norm_c
269 : CASE (cell_sym_monoclinic)
270 6486 : CALL get_cell(cell=cell, abc=abc)
271 6486 : beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
272 6486 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)
273 6486 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
274 6486 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)
275 : CASE (cell_sym_monoclinic_gamma_ab)
276 : ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
277 738 : CALL get_cell(cell=cell, abc=abc)
278 738 : a = 0.5_dp*(abc(1) + abc(2))
279 738 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
280 738 : acosg = a*COS(gamma)
281 738 : asing = a*SIN(gamma)
282 738 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
283 738 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
284 340875 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
285 : CASE (cell_sym_triclinic)
286 : ! Nothing to do
287 : END SELECT
288 :
289 : ! Do we have an (almost) orthorhombic cell?
290 : IF ((ABS(cell%hmat(1, 2)) < eps_hmat) .AND. (ABS(cell%hmat(1, 3)) < eps_hmat) .AND. &
291 : (ABS(cell%hmat(2, 1)) < eps_hmat) .AND. (ABS(cell%hmat(2, 3)) < eps_hmat) .AND. &
292 340137 : (ABS(cell%hmat(3, 1)) < eps_hmat) .AND. (ABS(cell%hmat(3, 2)) < eps_hmat)) THEN
293 311407 : cell%orthorhombic = .TRUE.
294 : ELSE
295 28730 : cell%orthorhombic = .FALSE.
296 : END IF
297 :
298 : ! Retain an exact orthorhombic cell
299 : ! (off-diagonal elements must remain zero identically to keep QS fast)
300 340137 : IF (cell%orthorhombic) THEN
301 311407 : cell%hmat(1, 2) = 0.0_dp
302 311407 : cell%hmat(1, 3) = 0.0_dp
303 311407 : cell%hmat(2, 1) = 0.0_dp
304 311407 : cell%hmat(2, 3) = 0.0_dp
305 311407 : cell%hmat(3, 1) = 0.0_dp
306 311407 : cell%hmat(3, 2) = 0.0_dp
307 : END IF
308 :
309 1360548 : dim = COUNT(cell%perd == 1)
310 340137 : IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
311 0 : CPABORT("Non-orthorhombic and not periodic")
312 : END IF
313 :
314 : ! Update deth and hmat_inv with enforced symmetry
315 340137 : cell%deth = ABS(det_3x3(cell%hmat))
316 340137 : IF (cell%deth < 1.0E-10_dp) THEN
317 : CALL cp_abort(__LOCATION__, &
318 : "An invalid set of cell vectors was obtained after applying "// &
319 0 : "the requested cell symmetry. The cell volume is too small")
320 : END IF
321 4421781 : cell%h_inv = inv_3x3(cell%hmat)
322 :
323 340137 : END SUBROUTINE init_cell
324 :
325 : ! **************************************************************************************************
326 : !> \brief ...
327 : !> \param cell ...
328 : !> \param cell_ref ...
329 : !> \param use_ref_cell ...
330 : !> \param cell_section ...
331 : !> \param topology_section ...
332 : !> \param check_for_ref ...
333 : !> \param para_env ...
334 : !> \par History
335 : !> 03.2005 created [teo]
336 : !> 03.2026 revamped logic with pdb and extxyz parsers
337 : !> \author Teodoro Laino
338 : ! **************************************************************************************************
339 112900 : RECURSIVE SUBROUTINE read_cell(cell, cell_ref, use_ref_cell, cell_section, &
340 : topology_section, check_for_ref, para_env)
341 :
342 : TYPE(cell_type), POINTER :: cell, cell_ref
343 : LOGICAL, INTENT(INOUT), OPTIONAL :: use_ref_cell
344 : TYPE(section_vals_type), OPTIONAL, POINTER :: cell_section, topology_section
345 : LOGICAL, INTENT(IN), OPTIONAL :: check_for_ref
346 : TYPE(mp_para_env_type), POINTER :: para_env
347 :
348 : REAL(KIND=dp), PARAMETER :: eps = 1.0E-14_dp
349 :
350 : CHARACTER(LEN=default_path_length) :: cell_file_name, coord_file_name, &
351 : error_msg
352 : INTEGER :: canonicalize_mode, cell_file_format, &
353 : coord_file_format, my_per
354 11290 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
355 : LOGICAL :: canonicalize_cell, cell_read_a, cell_read_abc, cell_read_alpha_beta_gamma, &
356 : cell_read_b, cell_read_c, cell_read_file, my_check_ref, tmp_comb_abc, tmp_comb_cell, &
357 : tmp_comb_top, topo_read_coord
358 : REAL(KIND=dp), DIMENSION(3) :: read_ang, read_len
359 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat_input, read_mat
360 11290 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
361 : TYPE(cell_type), POINTER :: cell_tmp
362 : TYPE(section_vals_type), POINTER :: cell_ref_section
363 :
364 11290 : my_check_ref = .TRUE.
365 11290 : NULLIFY (cell_ref_section, cell_par, cell_tmp, multiple_unit_cell)
366 : ! cell_tmp has two purposes:
367 : ! 1. for transferring matrix of cell vectors from individual
368 : ! file parser subroutines to read_mat here, assuming that
369 : ! unit conversion has been done in those subroutines;
370 : ! 2. for testing whether enforcing symmetry makes a new set
371 : ! of cell vectors significantly different from parsed input
372 11290 : CALL cell_create(cell_tmp)
373 11290 : IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
374 11290 : IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
375 11290 : IF (PRESENT(check_for_ref)) my_check_ref = check_for_ref
376 :
377 11290 : cell%deth = 0.0_dp
378 11290 : cell%orthorhombic = .FALSE.
379 45160 : cell%perd(:) = 1
380 11290 : cell%symmetry_id = cell_sym_none
381 146770 : cell%hmat(:, :) = 0.0_dp
382 146770 : cell%h_inv(:, :) = 0.0_dp
383 11290 : cell%input_cell_canonicalized = .FALSE.
384 146770 : cell%input_hmat(:, :) = 0.0_dp
385 146770 : cell%input_to_canonical(:, :) = 0.0_dp
386 146770 : cell%input_recip_to_canonical(:, :) = 0.0_dp
387 : cell_read_file = .FALSE.
388 : cell_read_a = .FALSE.
389 : cell_read_b = .FALSE.
390 : cell_read_c = .FALSE.
391 : cell_read_abc = .FALSE.
392 : cell_read_alpha_beta_gamma = .FALSE.
393 11290 : hmat_input(:, :) = 0.0_dp
394 11290 : read_mat(:, :) = 0.0_dp
395 11290 : read_ang(:) = 0.0_dp
396 11290 : read_len(:) = 0.0_dp
397 :
398 : ! Precedence of retrieving cell information from input:
399 : ! 1. CELL/CELL_FILE_NAME
400 : ! 2. CELL/ABC and optionally CELL/ALPHA_BETA_GAMMA
401 : ! 3. CELL/A, CELL/B, CELL/C
402 : ! 4. TOPOLOGY/COORD_FILE_NAME, if topology_section is present
403 : ! The actual order of processing is 4 -> 1 -> 2 -> 3, with
404 : ! case 4 merged to case 1 (if file format permits) first.
405 : ! Store data into either read_mat or read_ang and read_len
406 : ! in CP2K units, which will be converted to cell%hmat and A, B, C.
407 11290 : CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
408 11290 : CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
409 11290 : CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
410 11290 : CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
411 11290 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", explicit=cell_read_alpha_beta_gamma)
412 11290 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
413 11290 : CALL section_vals_val_get(cell_section, "CANONICALIZE", i_val=canonicalize_mode)
414 11290 : canonicalize_cell = (canonicalize_mode == canonicalize_cell_true)
415 :
416 : ! Case 4
417 11290 : tmp_comb_top = (.NOT. (cell_read_file .OR. cell_read_abc))
418 1560 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_a))
419 0 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_b))
420 0 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_c))
421 : IF (tmp_comb_top) THEN
422 : CALL cp_warn(__LOCATION__, &
423 : "None of the keywords CELL_FILE_NAME, ABC, or A, B, C "// &
424 : "are specified in CELL section. CP2K will now attempt to read "// &
425 : "TOPOLOGY/COORD_FILE_NAME if its format can be parsed for "// &
426 0 : "cell information.")
427 0 : IF (ASSOCIATED(topology_section)) THEN
428 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", explicit=topo_read_coord)
429 0 : IF (topo_read_coord) THEN
430 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", c_val=coord_file_name)
431 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_FORMAT", i_val=coord_file_format)
432 0 : SELECT CASE (coord_file_format) ! Add formats with both cell and coord parser manually
433 : CASE (do_coord_cif)
434 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
435 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cif)
436 : CASE (do_coord_cp2k)
437 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
438 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cp2k)
439 : CASE (do_coord_pdb)
440 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
441 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_pdb)
442 : CASE (do_coord_xyz)
443 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
444 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_extxyz)
445 : CASE DEFAULT
446 : CALL cp_abort(__LOCATION__, &
447 : "COORD_FILE_FORMAT is not set to one of the implemented "// &
448 0 : "CELL_FILE_FORMAT options and cannot be parsed for cell information!")
449 : END SELECT
450 : ELSE
451 : CALL cp_abort(__LOCATION__, &
452 0 : "COORD_FILE_NAME is not set, so no cell information is available!")
453 : END IF
454 : ELSE
455 : CALL cp_warn(__LOCATION__, &
456 : "TOPOLOGY section is not available, so COORD_FILE_NAME cannot "// &
457 0 : "be parsed for cell information in lieu of missing CELL settings.")
458 : END IF
459 : END IF
460 : ! Former logic in SUBROUTINE read_cell_from_external_file is moved here
461 11290 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
462 11290 : IF (cell_read_file) THEN ! Case 1
463 16 : tmp_comb_cell = (cell_read_abc .OR. (cell_read_a .OR. (cell_read_b .OR. cell_read_c)))
464 : IF (tmp_comb_cell) &
465 : CALL cp_warn(__LOCATION__, &
466 : "Cell Information provided through A, B, C, or ABC in conjunction "// &
467 : "with CELL_FILE_NAME. The definition in external file will override "// &
468 0 : "other ones.")
469 16 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
470 16 : CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=cell_file_format)
471 2 : SELECT CASE (cell_file_format)
472 : CASE (do_cell_cp2k)
473 2 : CALL read_cell_cp2k(cell_file_name, cell_tmp, para_env)
474 : CASE (do_cell_xsc)
475 0 : CALL read_cell_xsc(cell_file_name, cell_tmp, para_env)
476 : CASE (do_cell_extxyz)
477 2 : CALL read_cell_xyz(cell_file_name, cell_tmp, para_env)
478 : CASE (do_cell_pdb)
479 2 : CALL read_cell_pdb(cell_file_name, cell_tmp, para_env)
480 : CASE (do_cell_cif)
481 10 : CALL read_cell_cif(cell_file_name, cell_tmp, para_env)
482 : CASE DEFAULT
483 : CALL cp_abort(__LOCATION__, &
484 : "CELL_FILE_FORMAT is not set to one of the implemented "// &
485 16 : "options and cannot be parsed for cell information!")
486 : END SELECT
487 208 : read_mat = cell_tmp%hmat
488 : ELSE
489 11274 : IF (cell_read_abc) THEN ! Case 2
490 9714 : CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
491 38856 : read_len = cell_par
492 9714 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_par)
493 38856 : read_ang = cell_par
494 9714 : IF (cell_read_a .OR. cell_read_b .OR. cell_read_c) &
495 : CALL cp_warn(__LOCATION__, &
496 : "Cell information provided through vectors A, B or C in conjunction with ABC. "// &
497 0 : "The definition of the ABC keyword will override the one provided by A, B and C.")
498 : ELSE ! Case 3
499 1560 : tmp_comb_abc = ((cell_read_a .EQV. cell_read_b) .AND. (cell_read_b .EQV. cell_read_c))
500 : IF (tmp_comb_abc) THEN
501 1560 : CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
502 6240 : read_mat(:, 1) = cell_par(:)
503 1560 : CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
504 6240 : read_mat(:, 2) = cell_par(:)
505 1560 : CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
506 6240 : read_mat(:, 3) = cell_par(:)
507 1560 : IF (cell_read_alpha_beta_gamma) &
508 : CALL cp_warn(__LOCATION__, &
509 : "The keyword ALPHA_BETA_GAMMA is ignored because it was used without the "// &
510 0 : "keyword ABC.")
511 : ELSE
512 : CALL cp_abort(__LOCATION__, &
513 : "Neither of the keywords CELL_FILE_NAME or ABC are specified, "// &
514 0 : "and cell vector settings in A, B, C are incomplete!")
515 : END IF
516 : END IF
517 : END IF
518 :
519 : ! Convert read_mat or read_len and read_ang to actual cell%hmat
520 127990 : IF (ANY(read_mat(:, :) > eps)) THEN
521 : ! Make a warning before storing cell vectors that
522 : ! do not form a triangular matrix.
523 1576 : IF (.NOT. canonicalize_cell .AND. &
524 : ((ABS(read_mat(2, 1)) > eps) .OR. &
525 : (ABS(read_mat(3, 1)) > eps) .OR. &
526 : (ABS(read_mat(3, 2)) > eps))) THEN
527 348 : IF (canonicalize_mode == canonicalize_cell_auto) THEN
528 : CALL cp_warn(__LOCATION__, &
529 : "CELL%CANONICALIZE AUTO keeps the general input cell orientation. "// &
530 : "The cell matrix is not a lower triangle and does not conform to the "// &
531 : "program convention that A lies along the X-axis and B is in the XY plane. "// &
532 : "Set CELL%CANONICALIZE TRUE to explicitly transform the cell and supported "// &
533 346 : "cell-dependent input to the canonical internal frame.")
534 : ELSE
535 : CALL cp_warn(__LOCATION__, &
536 : "Cell vectors are read but cell matrix is not "// &
537 : "a lower triangle, not conforming to the program "// &
538 : "convention that A lies along the X-axis and "// &
539 2 : "B is in the XY plane.")
540 : END IF
541 : END IF
542 20488 : cell%hmat = read_mat
543 : ELSE
544 19428 : IF (ANY(read_ang(:) > eps) .AND. ANY(read_len(:) > eps)) THEN
545 : CALL set_cell_param(cell, cell_length=read_len, cell_angle=read_ang, &
546 9714 : do_init_cell=.FALSE.)
547 : ELSE
548 : CALL cp_abort(__LOCATION__, &
549 0 : "No meaningful cell information is read from parser!")
550 : END IF
551 : END IF
552 : ! Reset cell section so that only A, B, C are kept
553 11290 : CALL reset_cell_section_by_cell_mat(cell, cell_section)
554 :
555 : ! Multiple unit cell
556 11290 : CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
557 44726 : IF (ANY(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
558 :
559 11290 : CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
560 16 : SELECT CASE (my_per)
561 : CASE (use_perd_x)
562 64 : cell%perd = [1, 0, 0]
563 : CASE (use_perd_y)
564 16 : cell%perd = [0, 1, 0]
565 : CASE (use_perd_z)
566 24 : cell%perd = [0, 0, 1]
567 : CASE (use_perd_xy)
568 144 : cell%perd = [1, 1, 0]
569 : CASE (use_perd_xz)
570 32 : cell%perd = [1, 0, 1]
571 : CASE (use_perd_yz)
572 160 : cell%perd = [0, 1, 1]
573 : CASE (use_perd_xyz)
574 28504 : cell%perd = [1, 1, 1]
575 : CASE (use_perd_none)
576 16216 : cell%perd = [0, 0, 0]
577 : CASE DEFAULT
578 11290 : CPABORT("Invalid or not yet implemented cell periodicity")
579 : END SELECT
580 :
581 : ! Load requested cell symmetry
582 11290 : CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
583 : ! Try enforcing symmetry by initializing a temporary copy of cell
584 : ! and see if the resulting cell matrix differ significantly
585 146770 : hmat_input(:, :) = cell%hmat(:, :)
586 11290 : CALL cell_clone(cell, cell_tmp)
587 11290 : CALL init_cell(cell_tmp)
588 146484 : IF (.NOT. canonicalize_cell .AND. ANY(ABS(cell_tmp%hmat - cell%hmat) > eps)) THEN
589 : WRITE (UNIT=error_msg, FMT="(A)") &
590 : "When initializing cell vectors with requested symmetry, one "// &
591 : "or more elements of the cell matrix has varied significantly. "// &
592 : "The input parameters are either deviating from the symmetry, "// &
593 : "or not conforming to the program convention that cell matrix "// &
594 : "is a lower triangle. The symmetrized cell vectors will be used "// &
595 26 : "anyway with the input atomic coordinates."
596 26 : CALL cp_warn(__LOCATION__, error_msg)
597 : END IF
598 11290 : IF (canonicalize_cell) THEN
599 106 : CALL canonicalize_cell_matrix(cell_tmp)
600 106 : CALL cell_finalize_canonical_input(cell_tmp, hmat_input, cell_tmp%hmat)
601 : END IF
602 11290 : CALL cell_clone(cell_tmp, cell)
603 11290 : CALL cell_release(cell_tmp)
604 11290 : CALL reset_cell_section_by_cell_mat(cell, cell_section)
605 :
606 11290 : IF (my_check_ref) THEN
607 : ! Recursive check for reference cell requested
608 10870 : cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
609 10870 : IF (parsed_cp2k_input(cell_ref_section, check_this_section=.TRUE.)) THEN
610 26 : IF (PRESENT(use_ref_cell)) use_ref_cell = .TRUE.
611 : CALL read_cell(cell_ref, cell_ref, use_ref_cell=use_ref_cell, &
612 : cell_section=cell_ref_section, check_for_ref=.FALSE., &
613 26 : para_env=para_env)
614 : ELSE
615 10844 : CALL cell_clone(cell, cell_ref, tag="CELL_REF")
616 10844 : IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE.
617 : END IF
618 : END IF
619 :
620 11290 : END SUBROUTINE read_cell
621 :
622 : ! **************************************************************************************************
623 : !> \brief utility function to ease the transition to the new input.
624 : !> returns true if the new input was parsed
625 : !> \param input_file the parsed input file
626 : !> \param check_this_section ...
627 : !> \return ...
628 : !> \author fawzi
629 : ! **************************************************************************************************
630 10870 : FUNCTION parsed_cp2k_input(input_file, check_this_section) RESULT(res)
631 :
632 : TYPE(section_vals_type), POINTER :: input_file
633 : LOGICAL, INTENT(IN), OPTIONAL :: check_this_section
634 : LOGICAL :: res
635 :
636 : LOGICAL :: my_check
637 : TYPE(section_vals_type), POINTER :: glob_section
638 :
639 10870 : my_check = .FALSE.
640 10870 : IF (PRESENT(check_this_section)) my_check = check_this_section
641 10870 : res = ASSOCIATED(input_file)
642 10870 : IF (res) THEN
643 10870 : CPASSERT(input_file%ref_count > 0)
644 10870 : IF (.NOT. my_check) THEN
645 0 : glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
646 0 : CALL section_vals_get(glob_section, explicit=res)
647 : ELSE
648 10870 : CALL section_vals_get(input_file, explicit=res)
649 : END IF
650 : END IF
651 :
652 10870 : END FUNCTION parsed_cp2k_input
653 :
654 : ! **************************************************************************************************
655 : !> \brief Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma)
656 : !> using the convention: a parallel to the x axis, b in the x-y plane and
657 : !> and c univoquely determined; gamma is the angle between a and b; beta
658 : !> is the angle between c and a and alpha is the angle between c and b
659 : !> \param cell ...
660 : !> \param cell_length ...
661 : !> \param cell_angle ...
662 : !> \param periodic ...
663 : !> \param do_init_cell ...
664 : !> \date 03.2008
665 : !> \author Teodoro Laino
666 : ! **************************************************************************************************
667 9848 : SUBROUTINE set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
668 :
669 : TYPE(cell_type), POINTER :: cell
670 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: cell_length, cell_angle
671 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
672 : LOGICAL, INTENT(IN) :: do_init_cell
673 :
674 : REAL(KIND=dp), PARAMETER :: eps = EPSILON(0.0_dp)
675 :
676 : REAL(KIND=dp) :: cos_alpha, cos_beta, cos_gamma, sin_gamma
677 :
678 9848 : CPASSERT(ASSOCIATED(cell))
679 39392 : CPASSERT(ALL(cell_angle /= 0.0_dp))
680 :
681 9848 : cos_gamma = COS(cell_angle(3)); IF (ABS(cos_gamma) < eps) cos_gamma = 0.0_dp
682 9848 : IF (ABS(ABS(cos_gamma) - 1.0_dp) < eps) cos_gamma = SIGN(1.0_dp, cos_gamma)
683 9848 : sin_gamma = SIN(cell_angle(3)); IF (ABS(sin_gamma) < eps) sin_gamma = 0.0_dp
684 9848 : IF (ABS(ABS(sin_gamma) - 1.0_dp) < eps) sin_gamma = SIGN(1.0_dp, sin_gamma)
685 9848 : cos_beta = COS(cell_angle(2)); IF (ABS(cos_beta) < eps) cos_beta = 0.0_dp
686 9848 : IF (ABS(ABS(cos_beta) - 1.0_dp) < eps) cos_beta = SIGN(1.0_dp, cos_beta)
687 9848 : cos_alpha = COS(cell_angle(1)); IF (ABS(cos_alpha) < eps) cos_alpha = 0.0_dp
688 9848 : IF (ABS(ABS(cos_alpha) - 1.0_dp) < eps) cos_alpha = SIGN(1.0_dp, cos_alpha)
689 :
690 39392 : cell%hmat(:, 1) = [1.0_dp, 0.0_dp, 0.0_dp]
691 39392 : cell%hmat(:, 2) = [cos_gamma, sin_gamma, 0.0_dp]
692 39392 : cell%hmat(:, 3) = [cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp]
693 9848 : cell%hmat(3, 3) = SQRT(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
694 :
695 39392 : cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
696 39392 : cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
697 39392 : cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
698 :
699 9848 : IF (do_init_cell) THEN
700 134 : IF (PRESENT(periodic)) THEN
701 134 : CALL init_cell(cell=cell, periodic=periodic)
702 : ELSE
703 0 : CALL init_cell(cell=cell)
704 : END IF
705 : END IF
706 :
707 9848 : END SUBROUTINE set_cell_param
708 :
709 : ! **************************************************************************************************
710 : !> \brief Setup of the multiple unit_cell
711 : !> \param cell ...
712 : !> \param multiple_unit_cell ...
713 : !> \date 05.2009
714 : !> \author Teodoro Laino [tlaino]
715 : !> \version 1.0
716 : ! **************************************************************************************************
717 148 : SUBROUTINE set_multiple_unit_cell(cell, multiple_unit_cell)
718 :
719 : TYPE(cell_type), POINTER :: cell
720 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
721 :
722 148 : CPASSERT(ASSOCIATED(cell))
723 :
724 : ! Abort, if one of the value is set to zero
725 592 : IF (ANY(multiple_unit_cell <= 0)) &
726 : CALL cp_abort(__LOCATION__, &
727 : "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
728 0 : "A value of 0 or negative is meaningless!")
729 :
730 : ! Scale abc according to user request
731 592 : cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
732 592 : cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
733 592 : cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
734 :
735 148 : END SUBROUTINE set_multiple_unit_cell
736 :
737 : ! **************************************************************************************************
738 : !> \brief Reads cell information from CIF file
739 : !> \param cif_file_name ...
740 : !> \param cell ...
741 : !> \param para_env ...
742 : !> \date 12.2008
743 : !> \par Format Information implemented:
744 : !> _cell_length_a (_cell.length_a)
745 : !> _cell_length_b (_cell.length_b)
746 : !> _cell_length_c (_cell.length_c)
747 : !> _cell_angle_alpha (_cell.length_alpha)
748 : !> _cell_angle_beta (_cell.length_beta)
749 : !> _cell_angle_gamma (_cell.length_gamma)
750 : !>
751 : !> \author Teodoro Laino [tlaino]
752 : !> moved from topology_cif (1/2019 JHU)
753 : ! **************************************************************************************************
754 48 : SUBROUTINE read_cell_cif(cif_file_name, cell, para_env)
755 :
756 : CHARACTER(len=*) :: cif_file_name
757 : TYPE(cell_type), POINTER :: cell
758 : TYPE(mp_para_env_type), POINTER :: para_env
759 :
760 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_cif'
761 :
762 : INTEGER :: handle
763 : INTEGER, DIMENSION(3) :: periodic
764 : LOGICAL :: found
765 : REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths
766 : TYPE(cp_parser_type) :: parser
767 :
768 24 : CALL timeset(routineN, handle)
769 :
770 : CALL parser_create(parser, cif_file_name, &
771 24 : para_env=para_env, apply_preprocessing=.FALSE.)
772 :
773 : ! Parsing cell infos
774 96 : periodic = 1
775 : ! Check for _cell_length_a or _cell.length_a
776 : CALL parser_search_string(parser, "_cell_length_a", ignore_case=.FALSE., found=found, &
777 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
778 24 : IF (.NOT. found) THEN
779 : CALL parser_search_string(parser, "_cell.length_a", ignore_case=.FALSE., found=found, &
780 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
781 0 : IF (.NOT. found) &
782 0 : CPABORT("The field _cell_length_a or _cell.length_a was not found in CIF file! ")
783 : END IF
784 24 : CALL cif_get_real(parser, cell_lengths(1))
785 24 : cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
786 :
787 : ! Check for _cell_length_b or _cell.length_b
788 : CALL parser_search_string(parser, "_cell_length_b", ignore_case=.FALSE., found=found, &
789 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
790 24 : IF (.NOT. found) THEN
791 : CALL parser_search_string(parser, "_cell.length_b", ignore_case=.FALSE., found=found, &
792 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
793 0 : IF (.NOT. found) &
794 0 : CPABORT("The field _cell_length_b or _cell.length_b was not found in CIF file! ")
795 : END IF
796 24 : CALL cif_get_real(parser, cell_lengths(2))
797 24 : cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
798 :
799 : ! Check for _cell_length_c or _cell.length_c
800 : CALL parser_search_string(parser, "_cell_length_c", ignore_case=.FALSE., found=found, &
801 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
802 24 : IF (.NOT. found) THEN
803 : CALL parser_search_string(parser, "_cell.length_c", ignore_case=.FALSE., found=found, &
804 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
805 0 : IF (.NOT. found) &
806 0 : CPABORT("The field _cell_length_c or _cell.length_c was not found in CIF file! ")
807 : END IF
808 24 : CALL cif_get_real(parser, cell_lengths(3))
809 24 : cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
810 :
811 : ! Check for _cell_angle_alpha or _cell.angle_alpha
812 : CALL parser_search_string(parser, "_cell_angle_alpha", ignore_case=.FALSE., found=found, &
813 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
814 24 : IF (.NOT. found) THEN
815 : CALL parser_search_string(parser, "_cell.angle_alpha", ignore_case=.FALSE., found=found, &
816 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
817 0 : IF (.NOT. found) &
818 0 : CPABORT("The field _cell_angle_alpha or _cell.angle_alpha was not found in CIF file! ")
819 : END IF
820 24 : CALL cif_get_real(parser, cell_angles(1))
821 24 : cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
822 :
823 : ! Check for _cell_angle_beta or _cell.angle_beta
824 : CALL parser_search_string(parser, "_cell_angle_beta", ignore_case=.FALSE., found=found, &
825 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
826 24 : IF (.NOT. found) THEN
827 : CALL parser_search_string(parser, "_cell.angle_beta", ignore_case=.FALSE., found=found, &
828 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
829 0 : IF (.NOT. found) &
830 0 : CPABORT("The field _cell_angle_beta or _cell.angle_beta was not found in CIF file! ")
831 : END IF
832 24 : CALL cif_get_real(parser, cell_angles(2))
833 24 : cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
834 :
835 : ! Check for _cell_angle_gamma or _cell.angle_gamma
836 : CALL parser_search_string(parser, "_cell_angle_gamma", ignore_case=.FALSE., found=found, &
837 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
838 24 : IF (.NOT. found) THEN
839 : CALL parser_search_string(parser, "_cell.angle_gamma", ignore_case=.FALSE., found=found, &
840 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
841 0 : IF (.NOT. found) &
842 0 : CPABORT("The field _cell_angle_gamma or _cell.angle_gamma was not found in CIF file! ")
843 : END IF
844 24 : CALL cif_get_real(parser, cell_angles(3))
845 24 : cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
846 :
847 : ! Create cell
848 : CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
849 24 : do_init_cell=.TRUE.)
850 :
851 24 : CALL parser_release(parser)
852 :
853 24 : CALL timestop(handle)
854 :
855 72 : END SUBROUTINE read_cell_cif
856 :
857 : ! **************************************************************************************************
858 : !> \brief Reads REAL from the CIF file.. This wrapper is needed in order to
859 : !> treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
860 : !> \param parser ...
861 : !> \param r ...
862 : !> \date 12.2008
863 : !> \author Teodoro Laino [tlaino]
864 : ! **************************************************************************************************
865 144 : SUBROUTINE cif_get_real(parser, r)
866 :
867 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
868 : REAL(KIND=dp), INTENT(OUT) :: r
869 :
870 : CHARACTER(LEN=default_string_length) :: s_tag
871 : INTEGER :: iln
872 :
873 144 : CALL parser_get_object(parser, s_tag)
874 144 : iln = LEN_TRIM(s_tag)
875 144 : IF (INDEX(s_tag, "(") /= 0) iln = INDEX(s_tag, "(") - 1
876 144 : READ (s_tag(1:iln), *) r
877 :
878 144 : END SUBROUTINE cif_get_real
879 :
880 : ! **************************************************************************************************
881 : !> \brief Reads xyz file and pass comments on the second line to get cell information
882 : !> \param xyz_file_name ...
883 : !> \param cell ...
884 : !> \param para_env ...
885 : !> \par History
886 : !> 03.2026 - Created as read_cell_extxyz with extended XYZ parser
887 : !> 06.2026 - Refactored the parser to allow for reftraj use
888 : !> \author HE Zilong
889 : ! **************************************************************************************************
890 6 : SUBROUTINE read_cell_xyz(xyz_file_name, cell, para_env)
891 :
892 : CHARACTER(len=*) :: xyz_file_name
893 : TYPE(cell_type), POINTER :: cell
894 : TYPE(mp_para_env_type), POINTER :: para_env
895 :
896 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_xyz'
897 :
898 : INTEGER :: handle
899 : LOGICAL :: has_cell
900 : TYPE(cp_parser_type) :: parser
901 :
902 2 : CALL timeset(routineN, handle)
903 :
904 : CALL parser_create(parser, xyz_file_name, &
905 2 : para_env=para_env, apply_preprocessing=.FALSE.)
906 2 : CALL parser_get_next_line(parser, 2) ! Skip number of atoms
907 2 : CALL read_xyz_comment(parser%input_line, cell, has_cell)
908 2 : IF (.NOT. has_cell) THEN
909 : CALL cp_abort(__LOCATION__, &
910 : "The keyword CELL_FILE_FORMAT requested cell information "// &
911 : "from XYZ file, but it is not available from the file <"// &
912 0 : TRIM(ADJUSTL(xyz_file_name))//"> as CELL_FILE_NAME specified!")
913 : END IF
914 2 : CALL parser_release(parser)
915 2 : CALL timestop(handle)
916 :
917 6 : END SUBROUTINE read_cell_xyz
918 :
919 : ! **************************************************************************************************
920 : !> \brief Reads comment line of XYZ files to get cell, step, time and energy info
921 : !> \param line the single comment line of an XYZ file
922 : !> \param cell the pointer to which cell is written
923 : !> \param has_cell a flag for presence of cell information
924 : !> \param step an integer value for step number, if requested; HUGE(0) if not available
925 : !> \param time a real value for time, if requested; HUGE(0.0_dp) if not available
926 : !> \param ener a real value for energy, if requested; HUGE(0.0_dp) if not available
927 : !> \par Intended for both the FORCE_EVAL/SUBSYS/CELL and MOTION/MD/REFTRAJ.
928 : !> At minimum, should work with outputs from write_trajectory() in
929 : !> src/motion_utils.F, around line 879 of src/motion/dumpdcd.F and
930 : !> write_final_structure() in src/particle_methods.F (for cell only).
931 : !> Recognized formats (case insensitive, no hard restriction on data width):
932 : !> (1) Extended xyz format, whose comment on the second line contains fields:
933 : !> Lattice="Ax Ay Az Bx By Bz Cx Cy Cz"
934 : !> where Ax, Ay, Az are three Cartesian components of cell vector A,
935 : !> Bx, By, Bz are components of B, Cx, Cy, Cz are components of C,
936 : !> all in the unit of angstrom, and must occur;
937 : !> Step=S
938 : !> where S is the integer step number;
939 : !> Time=T
940 : !> where T is the time in femtoseconds;
941 : !> Energy=E
942 : !> where E is the energy.
943 : !> No whitespace around the = sign is present; the whitespace is used as
944 : !> the delimiter between fields; apart from lattice= at the front, other
945 : !> fields are optional and do not have a fixed order.
946 : !> (2) dumpdcd format, whose comment on the second line contains fields:
947 : !> a = A, b = B, c = C, alpha = ALPHA, beta = BETA, gamma = GAMMA
948 : !> where A, B, C are three lengths of cell vectors in angstrom, ALPHA,
949 : !> BETA, GAMMA are three angles between cell vectors in degrees;
950 : !> i = I,
951 : !> where I is the integer step number;
952 : !> time = T,
953 : !> where T is the time in femtoseconds;
954 : !> E = ENER,
955 : !> where ENER is the energy.
956 : !> There is one whitespace before and after each equal sign; the comma
957 : !> is used as the delimiter between fields; the cell information is
958 : !> optional.
959 : !>
960 : !> History
961 : !> 06.2026 - Created by combining the extxyz parser from former read_cell_extxyz
962 : !> and the parser for reftraj in src/motion/integrator.F
963 : !> \author HE Zilong
964 : ! **************************************************************************************************
965 284 : SUBROUTINE read_xyz_comment(line, cell, has_cell, step, time, ener)
966 :
967 : CHARACTER(LEN=*), INTENT(IN) :: line
968 : TYPE(cell_type), INTENT(INOUT), POINTER :: cell
969 : LOGICAL, INTENT(OUT) :: has_cell
970 : INTEGER, INTENT(OUT), OPTIONAL :: step
971 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: time, ener
972 :
973 : CHARACTER(LEN=3) :: abc
974 : CHARACTER(LEN=max_line_length) :: my_line, raw_str
975 : INTEGER :: i, id1, id2, ios, j, my_step
976 : REAL(KIND=dp) :: my_ener, my_time
977 : REAL(KIND=dp), DIMENSION(3) :: my_abc, my_albega
978 : REAL(KIND=dp), DIMENSION(3, 3) :: my_hmat
979 :
980 284 : has_cell = .FALSE.
981 284 : my_step = HUGE(0)
982 284 : my_time = HUGE(0.0_dp)
983 284 : my_ener = HUGE(0.0_dp)
984 284 : my_hmat = 0.0_dp
985 284 : my_abc = 0.0_dp
986 284 : my_albega = 0.0_dp
987 :
988 284 : my_line = line
989 284 : CALL uppercase(my_line)
990 284 : id1 = INDEX(my_line, "LATTICE=")
991 284 : IF (id1 > 0) THEN ! Extended XYZ
992 20 : id2 = INDEX(my_line(id1 + 9:), '"') ! Strip 'LATTICE="' and find the next quote
993 20 : READ (my_line(id1 + 9:id1 + id2 + 7), '(A)') raw_str
994 20 : READ (raw_str, *, IOSTAT=ios) my_hmat(:, 1), my_hmat(:, 2), my_hmat(:, 3)
995 20 : IF (ios /= 0) THEN
996 : CALL cp_abort(__LOCATION__, "Error while parsing input line for cell vectors as "// &
997 : "extended XYZ format: expected 9 real values in the <lattice=> "// &
998 0 : "quoted field, found <"//TRIM(raw_str)//"> which is invalid!")
999 : ELSE
1000 20 : has_cell = .TRUE.
1001 80 : DO i = 1, 3
1002 260 : DO j = 1, 3
1003 240 : cell%hmat(j, i) = cp_unit_to_cp2k(my_hmat(j, i), "angstrom")
1004 : END DO
1005 : END DO
1006 : END IF
1007 20 : IF (PRESENT(step)) THEN
1008 18 : id1 = INDEX(my_line, "STEP=")
1009 18 : IF (id1 > 0) THEN
1010 18 : READ (my_line(id1 + 5:), '(A)') raw_str
1011 18 : READ (raw_str, *, IOSTAT=ios) my_step
1012 18 : IF (ios /= 0) THEN
1013 : CALL cp_abort(__LOCATION__, &
1014 : "Error while parsing input line for step as extended "// &
1015 : "XYZ format: expected 1 integer value in the <step=> "// &
1016 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1017 : END IF
1018 : END IF
1019 18 : step = my_step
1020 : END IF
1021 20 : IF (PRESENT(time)) THEN
1022 18 : id1 = INDEX(my_line, "TIME=")
1023 18 : IF (id1 > 0) THEN
1024 18 : READ (my_line(id1 + 5:), '(A)') raw_str
1025 18 : READ (raw_str, *, IOSTAT=ios) my_time
1026 18 : IF (ios /= 0) THEN
1027 : CALL cp_abort(__LOCATION__, &
1028 : "Error while parsing input line for time as extended "// &
1029 : "XYZ format: expected 1 real value in the <time=> "// &
1030 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1031 : END IF
1032 : END IF
1033 18 : time = my_time
1034 : END IF
1035 20 : IF (PRESENT(ener)) THEN
1036 18 : id1 = INDEX(my_line, "ENERGY=")
1037 18 : IF (id1 > 0) THEN
1038 18 : READ (my_line(id1 + 7:), '(A)') raw_str
1039 18 : READ (raw_str, *, IOSTAT=ios) my_ener
1040 18 : IF (ios /= 0) THEN
1041 : CALL cp_abort(__LOCATION__, &
1042 : "Error while parsing input line for energy as extended "// &
1043 : "XYZ format: expected 1 real value in the <energy=> "// &
1044 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1045 : END IF
1046 : END IF
1047 18 : ener = my_ener
1048 : END IF
1049 : ELSE ! May or may not be dumpdcd format, and may or may not has cell
1050 264 : abc = "ABC"
1051 1056 : DO i = 1, 3
1052 792 : id1 = INDEX(my_line, " "//abc(i:i)//" = ")
1053 1056 : IF (id1 > 0) THEN
1054 0 : READ (my_line(id1 + 5:), '(A)') raw_str
1055 0 : READ (raw_str, *, IOSTAT=ios) my_abc(i)
1056 0 : IF (ios /= 0) THEN
1057 : CALL cp_abort(__LOCATION__, &
1058 : "Error while parsing input line for cell vector as dumpdcd "// &
1059 : "XYZ format: expected 1 real value in the <"//abc(i:i)//" = > "// &
1060 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1061 : ELSE
1062 0 : my_abc(i) = cp_unit_to_cp2k(my_abc(i), "angstrom")
1063 : END IF
1064 : END IF
1065 : END DO
1066 264 : id1 = INDEX(my_line, " ALPHA = ")
1067 264 : IF (id1 > 0) THEN
1068 0 : READ (my_line(id1 + 9:), '(A)') raw_str
1069 0 : READ (raw_str, *, IOSTAT=ios) my_albega(1)
1070 0 : IF (ios /= 0) THEN
1071 : CALL cp_abort(__LOCATION__, &
1072 : "Error while parsing input line for cell angle alpha as dumpdcd "// &
1073 : "XYZ format: expected 1 real value in the <alpha = > "// &
1074 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1075 : ELSE
1076 0 : my_albega(1) = cp_unit_to_cp2k(my_albega(1), "deg")
1077 : END IF
1078 : END IF
1079 264 : id1 = INDEX(my_line, " BETA = ")
1080 264 : IF (id1 > 0) THEN
1081 0 : READ (my_line(id1 + 8:), '(A)') raw_str
1082 0 : READ (raw_str, *, IOSTAT=ios) my_albega(2)
1083 0 : IF (ios /= 0) THEN
1084 : CALL cp_abort(__LOCATION__, &
1085 : "Error while parsing input line for cell angle beta as dumpdcd "// &
1086 : "XYZ format: expected 1 real value in the <beta = > "// &
1087 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1088 : ELSE
1089 0 : my_albega(2) = cp_unit_to_cp2k(my_albega(2), "deg")
1090 : END IF
1091 : END IF
1092 264 : id1 = INDEX(my_line, " GAMMA = ")
1093 264 : IF (id1 > 0) THEN
1094 0 : READ (my_line(id1 + 9:), '(A)') raw_str
1095 0 : READ (raw_str, *, IOSTAT=ios) my_albega(3)
1096 0 : IF (ios /= 0) THEN
1097 : CALL cp_abort(__LOCATION__, &
1098 : "Error while parsing input line for cell angle gamma as dumpdcd "// &
1099 : "XYZ format: expected 1 real value in the <gamma = > "// &
1100 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1101 : ELSE
1102 0 : my_albega(3) = cp_unit_to_cp2k(my_albega(3), "deg")
1103 : END IF
1104 : END IF
1105 528 : IF (ALL(my_abc(1:3) > 0.0_dp) .AND. ALL(my_albega(1:3) > 0.0_dp)) THEN
1106 0 : has_cell = .TRUE.
1107 0 : CALL set_cell_param(cell, my_abc, my_albega, do_init_cell=.FALSE.)
1108 : END IF
1109 264 : IF (PRESENT(step)) THEN
1110 264 : id1 = INDEX(my_line, " I = ")
1111 264 : IF (id1 > 0) THEN
1112 214 : READ (my_line(id1 + 5:), '(A)') raw_str
1113 214 : READ (raw_str, *, IOSTAT=ios) my_step
1114 214 : IF (ios /= 0) THEN
1115 : CALL cp_abort(__LOCATION__, &
1116 : "Error while parsing input line for step as dumpdcd "// &
1117 : "XYZ format: expected 1 integer value in the <i = > "// &
1118 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1119 : END IF
1120 : END IF
1121 264 : step = my_step
1122 : END IF
1123 264 : IF (PRESENT(time)) THEN
1124 264 : id1 = INDEX(my_line, " TIME = ")
1125 264 : IF (id1 > 0) THEN
1126 214 : READ (my_line(id1 + 8:), '(A)') raw_str
1127 214 : READ (raw_str, *, IOSTAT=ios) my_time
1128 214 : IF (ios /= 0) THEN
1129 : CALL cp_abort(__LOCATION__, &
1130 : "Error while parsing input line for time as dumpdcd "// &
1131 : "XYZ format: expected 1 real value in the <time = > "// &
1132 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1133 : END IF
1134 : END IF
1135 264 : time = my_time
1136 : END IF
1137 264 : IF (PRESENT(ener)) THEN
1138 264 : id1 = INDEX(my_line, " E = ")
1139 264 : IF (id1 > 0) THEN
1140 214 : READ (my_line(id1 + 5:), '(A)') raw_str
1141 214 : READ (raw_str, *, IOSTAT=ios) my_ener
1142 214 : IF (ios /= 0) THEN
1143 : CALL cp_abort(__LOCATION__, &
1144 : "Error while parsing input line for energy as dumpdcd "// &
1145 : "XYZ format: expected 1 real value in the <E = > "// &
1146 0 : "field, found <"//TRIM(raw_str)//"> which is invalid!")
1147 : END IF
1148 : END IF
1149 264 : ener = my_ener
1150 : END IF
1151 : END IF
1152 :
1153 284 : END SUBROUTINE read_xyz_comment
1154 :
1155 : ! **************************************************************************************************
1156 : !> \brief Reads cell information from CRYST1 record of PDB file
1157 : !> \param pdb_file_name ...
1158 : !> \param cell ...
1159 : !> \param para_env ...
1160 : !> \date 03.2026
1161 : !> \par CRYST1 record may contain space group and Z value at the end,
1162 : !> but here only the first entries are read:
1163 : !> COLUMNS DATA TYPE FIELD DEFINITION
1164 : !> -------------------------------------------------------------
1165 : !> 1 - 6 Record name "CRYST1"
1166 : !> 7 - 15 Real(9.3) a a (Angstroms).
1167 : !> 16 - 24 Real(9.3) b b (Angstroms).
1168 : !> 25 - 33 Real(9.3) c c (Angstroms).
1169 : !> 34 - 40 Real(7.2) alpha alpha (degrees).
1170 : !> 41 - 47 Real(7.2) beta beta (degrees).
1171 : !> 48 - 54 Real(7.2) gamma gamma (degrees).
1172 : ! **************************************************************************************************
1173 4 : SUBROUTINE read_cell_pdb(pdb_file_name, cell, para_env)
1174 :
1175 : CHARACTER(len=*) :: pdb_file_name
1176 : TYPE(cell_type), POINTER :: cell
1177 : TYPE(mp_para_env_type), POINTER :: para_env
1178 :
1179 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_pdb'
1180 :
1181 : CHARACTER(LEN=default_string_length) :: cryst
1182 : INTEGER :: handle, i, ios
1183 : INTEGER, DIMENSION(3) :: periodic
1184 : LOGICAL :: found
1185 : REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths
1186 : TYPE(cp_parser_type) :: parser
1187 :
1188 2 : CALL timeset(routineN, handle)
1189 :
1190 : CALL parser_create(parser, pdb_file_name, &
1191 2 : para_env=para_env, apply_preprocessing=.FALSE.)
1192 :
1193 : CALL parser_search_string(parser, "CRYST1", ignore_case=.FALSE., found=found, &
1194 2 : begin_line=.TRUE., search_from_begin_of_file=.TRUE.)
1195 2 : IF (.NOT. found) &
1196 0 : CPABORT("The line <CRYST1> was not found in PDB file! ")
1197 :
1198 8 : periodic = 1
1199 2 : READ (parser%input_line, *, IOSTAT=ios) cryst, cell_lengths(:), cell_angles(:)
1200 2 : IF (ios /= 0) THEN
1201 : CALL cp_abort(__LOCATION__, "Error while parsing PDB file "// &
1202 : "<"//TRIM(pdb_file_name)//"> for cell lengths and angles: "// &
1203 0 : "found CRYST1 line as <"//TRIM(parser%input_line)//">")
1204 : END IF
1205 8 : DO i = 1, 3
1206 6 : cell_lengths(i) = cp_unit_to_cp2k(cell_lengths(i), "angstrom")
1207 8 : cell_angles(i) = cp_unit_to_cp2k(cell_angles(i), "deg")
1208 : END DO
1209 : CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
1210 2 : do_init_cell=.TRUE.)
1211 :
1212 2 : CALL parser_release(parser)
1213 :
1214 2 : CALL timestop(handle)
1215 :
1216 6 : END SUBROUTINE read_cell_pdb
1217 :
1218 : ! **************************************************************************************************
1219 : !> \brief Reads cell information from cp2k file
1220 : !> \param cp2k_file_name ...
1221 : !> \param cell ...
1222 : !> \param para_env ...
1223 : !> \date 03.2026
1224 : !> \par Isolated from former read_cell_from_external_file
1225 : ! **************************************************************************************************
1226 4 : SUBROUTINE read_cell_cp2k(cp2k_file_name, cell, para_env)
1227 :
1228 : CHARACTER(len=*) :: cp2k_file_name
1229 : TYPE(cell_type), POINTER :: cell
1230 : TYPE(mp_para_env_type), POINTER :: para_env
1231 :
1232 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_cp2k'
1233 :
1234 : INTEGER :: handle, i, idum, j
1235 : LOGICAL :: my_end
1236 : REAL(KIND=dp) :: xdum
1237 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1238 : TYPE(cp_parser_type) :: parser
1239 :
1240 2 : CALL timeset(routineN, handle)
1241 :
1242 : CALL parser_create(parser, cp2k_file_name, &
1243 2 : para_env=para_env, apply_preprocessing=.FALSE.)
1244 :
1245 2 : CALL parser_get_next_line(parser, 1)
1246 2 : my_end = .FALSE.
1247 24 : DO WHILE (.NOT. my_end)
1248 22 : READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
1249 24 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1250 : END DO
1251 8 : DO i = 1, 3
1252 26 : DO j = 1, 3
1253 24 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
1254 : END DO
1255 : END DO
1256 :
1257 2 : CALL parser_release(parser)
1258 :
1259 2 : CALL timestop(handle)
1260 :
1261 6 : END SUBROUTINE read_cell_cp2k
1262 :
1263 : ! **************************************************************************************************
1264 : !> \brief Reads cell information from xsc file
1265 : !> \param xsc_file_name ...
1266 : !> \param cell ...
1267 : !> \param para_env ...
1268 : !> \date 03.2026
1269 : !> \par Isolated from former read_cell_from_external_file
1270 : ! **************************************************************************************************
1271 0 : SUBROUTINE read_cell_xsc(xsc_file_name, cell, para_env)
1272 :
1273 : CHARACTER(len=*) :: xsc_file_name
1274 : TYPE(cell_type), POINTER :: cell
1275 : TYPE(mp_para_env_type), POINTER :: para_env
1276 :
1277 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_xsc'
1278 :
1279 : INTEGER :: handle, i, idum, j
1280 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1281 : TYPE(cp_parser_type) :: parser
1282 :
1283 0 : CALL timeset(routineN, handle)
1284 :
1285 : CALL parser_create(parser, xsc_file_name, &
1286 0 : para_env=para_env, apply_preprocessing=.FALSE.)
1287 :
1288 0 : CALL parser_get_next_line(parser, 1)
1289 0 : READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
1290 0 : DO i = 1, 3
1291 0 : DO j = 1, 3
1292 0 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
1293 : END DO
1294 : END DO
1295 :
1296 0 : CALL parser_release(parser)
1297 :
1298 0 : CALL timestop(handle)
1299 :
1300 0 : END SUBROUTINE read_cell_xsc
1301 :
1302 : ! **************************************************************************************************
1303 : !> \brief Reset cell section by matrix in cell-type pointer
1304 : !> \param cell ...
1305 : !> \param cell_section ...
1306 : !> \date 03.2026
1307 : !> \par Alternative keywords for cell settings will be unset
1308 : !> except MULTIPLE_UNIT_CELL, PERIODIC and SYMMETRY.
1309 : ! **************************************************************************************************
1310 22580 : SUBROUTINE reset_cell_section_by_cell_mat(cell, cell_section)
1311 :
1312 : TYPE(cell_type), POINTER :: cell
1313 : TYPE(section_vals_type), POINTER :: cell_section
1314 :
1315 22580 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
1316 :
1317 22580 : CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
1318 22580 : CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
1319 22580 : CALL section_vals_val_unset(cell_section, "ABC")
1320 22580 : CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
1321 22580 : CALL section_vals_val_unset(cell_section, "A")
1322 22580 : CALL section_vals_val_unset(cell_section, "B")
1323 22580 : CALL section_vals_val_unset(cell_section, "C")
1324 22580 : ALLOCATE (cell_par(3))
1325 180640 : cell_par = cell%hmat(:, 1)
1326 22580 : CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
1327 22580 : ALLOCATE (cell_par(3))
1328 180640 : cell_par = cell%hmat(:, 2)
1329 22580 : CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
1330 22580 : ALLOCATE (cell_par(3))
1331 180640 : cell_par = cell%hmat(:, 3)
1332 22580 : CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
1333 :
1334 22580 : END SUBROUTINE reset_cell_section_by_cell_mat
1335 :
1336 : ! **************************************************************************************************
1337 : !> \brief Write the cell parameters to the output unit.
1338 : !> \param cell ...
1339 : !> \param subsys_section ...
1340 : !> \param tag ...
1341 : !> \date 02.06.2000
1342 : !> \par History
1343 : !> - 11.2008 Teodoro Laino [tlaino] - rewrite and enabling user driven units
1344 : !> \author Matthias Krack
1345 : !> \version 1.0
1346 : ! **************************************************************************************************
1347 31057 : SUBROUTINE write_cell(cell, subsys_section, tag)
1348 :
1349 : TYPE(cell_type), POINTER :: cell
1350 : TYPE(section_vals_type), POINTER :: subsys_section
1351 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
1352 :
1353 : CHARACTER(LEN=default_string_length) :: label, unit_str
1354 : INTEGER :: output_unit
1355 : TYPE(cp_logger_type), POINTER :: logger
1356 :
1357 31057 : NULLIFY (logger)
1358 31057 : logger => cp_get_default_logger()
1359 31057 : IF (PRESENT(tag)) THEN
1360 22808 : label = TRIM(tag)//"|"
1361 : ELSE
1362 8249 : label = TRIM(cell%tag)//"|"
1363 : END IF
1364 :
1365 31057 : output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", extension=".Log")
1366 31057 : CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
1367 31057 : CALL write_cell_low(cell, unit_str, output_unit, label)
1368 31057 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, "PRINT%CELL")
1369 :
1370 31057 : END SUBROUTINE write_cell
1371 :
1372 : ! **************************************************************************************************
1373 : !> \brief Write the cell parameters to the output unit
1374 : !> \param cell ...
1375 : !> \param unit_str ...
1376 : !> \param output_unit ...
1377 : !> \param label ...
1378 : !> \date 17.05.2023
1379 : !> \par History
1380 : !> - Extracted from write_cell (17.05.2023, MK)
1381 : !> \version 1.0
1382 : ! **************************************************************************************************
1383 31057 : SUBROUTINE write_cell_low(cell, unit_str, output_unit, label)
1384 :
1385 : TYPE(cell_type), POINTER :: cell
1386 : CHARACTER(LEN=*), INTENT(IN) :: unit_str
1387 : INTEGER, INTENT(IN) :: output_unit
1388 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: label
1389 :
1390 : CHARACTER(LEN=12) :: tag
1391 : CHARACTER(LEN=3) :: string
1392 : CHARACTER(LEN=default_string_length) :: my_label
1393 : REAL(KIND=dp) :: alpha, beta, gamma, val
1394 : REAL(KIND=dp), DIMENSION(3) :: abc
1395 : TYPE(enumeration_type), POINTER :: enum
1396 : TYPE(keyword_type), POINTER :: keyword
1397 : TYPE(section_type), POINTER :: section
1398 :
1399 31057 : NULLIFY (enum)
1400 31057 : NULLIFY (keyword)
1401 31057 : NULLIFY (section)
1402 :
1403 41183 : IF (output_unit > 0) THEN
1404 10126 : CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma, tag=tag)
1405 10126 : IF (PRESENT(label)) THEN
1406 10126 : my_label = label
1407 : ELSE
1408 0 : my_label = TRIM(tag)//"|"
1409 : END IF
1410 10126 : val = cp_unit_from_cp2k(cell%deth, TRIM(unit_str)//"^3")
1411 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T61,F20.6)") &
1412 10126 : TRIM(my_label)//" Volume ["//TRIM(unit_str)//"^3]:", val
1413 10126 : val = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
1414 : WRITE (UNIT=output_unit, FMT="(T2,A,T30,3F10.3,3X,A6,F12.6)") &
1415 40504 : TRIM(my_label)//" Vector a ["//TRIM(unit_str)//"]:", cell%hmat(:, 1)*val, &
1416 10126 : "|a| = ", abc(1)*val, &
1417 40504 : TRIM(my_label)//" Vector b ["//TRIM(unit_str)//"]:", cell%hmat(:, 2)*val, &
1418 10126 : "|b| = ", abc(2)*val, &
1419 40504 : TRIM(my_label)//" Vector c ["//TRIM(unit_str)//"]:", cell%hmat(:, 3)*val, &
1420 20252 : "|c| = ", abc(3)*val
1421 : WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.6)") &
1422 10126 : TRIM(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
1423 10126 : TRIM(my_label)//" Angle (a,c), beta [degree]: ", beta, &
1424 20252 : TRIM(my_label)//" Angle (a,b), gamma [degree]: ", gamma
1425 10126 : IF (cell%symmetry_id /= cell_sym_none) THEN
1426 2208 : CALL create_cell_section(section)
1427 2208 : keyword => section_get_keyword(section, "SYMMETRY")
1428 2208 : CALL keyword_get(keyword, enum=enum)
1429 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
1430 2208 : TRIM(my_label)//" Requested initial symmetry: ", &
1431 4416 : ADJUSTR(TRIM(enum_i2c(enum, cell%symmetry_id)))
1432 2208 : CALL section_release(section)
1433 : END IF
1434 10126 : IF (cell%orthorhombic) THEN
1435 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1436 6284 : TRIM(my_label)//" Numerically orthorhombic: ", "YES"
1437 : ELSE
1438 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1439 3842 : TRIM(my_label)//" Numerically orthorhombic: ", " NO"
1440 : END IF
1441 40504 : IF (SUM(cell%perd(1:3)) == 0) THEN
1442 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
1443 1914 : TRIM(my_label)//" Periodicity", "NONE"
1444 : ELSE
1445 8212 : string = ""
1446 8212 : IF (cell%perd(1) == 1) string = TRIM(string)//"X"
1447 8212 : IF (cell%perd(2) == 1) string = TRIM(string)//"Y"
1448 8212 : IF (cell%perd(3) == 1) string = TRIM(string)//"Z"
1449 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1450 8212 : TRIM(my_label)//" Periodicity", ADJUSTR(string)
1451 : END IF
1452 : END IF
1453 :
1454 31057 : END SUBROUTINE write_cell_low
1455 :
1456 : END MODULE cell_methods
|
