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