(git:561f475)
Loading...
Searching...
No Matches
cell_methods.F
Go to the documentation of this file.
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! **************************************************************************************************
16 USE cell_types, ONLY: &
33 USE cp_units, ONLY: cp_unit_from_cp2k,&
35 USE input_constants, ONLY: &
43 USE input_section_types, ONLY: &
47 USE kinds, ONLY: default_path_length,&
49 dp
51 USE mathconstants, ONLY: degree,&
52 sqrt3
53 USE mathlib, ONLY: angle,&
54 det_3x3,&
58#include "./base/base_uses.f90"
59
60 IMPLICIT NONE
61
62 PRIVATE
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_methods'
65
66 PUBLIC :: cell_create, &
68 init_cell, &
69 read_cell, &
76 write_cell, &
78
79CONTAINS
80
81! **************************************************************************************************
82!> \brief allocates and initializes a cell
83!> \param cell the cell to initialize
84!> \param hmat the h matrix that defines the cell
85!> \param periodic periodicity of the cell
86!> \param tag ...
87!> \par History
88!> 09.2003 created [fawzi]
89!> \author Fawzi Mohamed
90! **************************************************************************************************
91 SUBROUTINE cell_create(cell, hmat, periodic, tag)
92
93 TYPE(cell_type), POINTER :: cell
94 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
95 OPTIONAL :: hmat
96 INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
97 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
98
99 cpassert(.NOT. ASSOCIATED(cell))
100 ALLOCATE (cell)
101 cell%ref_count = 1
102 IF (PRESENT(periodic)) THEN
103 cell%perd = periodic
104 ELSE
105 cell%perd = 1
106 END IF
107 cell%orthorhombic = .false.
108 cell%input_cell_canonicalized = .false.
109 cell%input_hmat(:, :) = 0.0_dp
110 cell%input_to_canonical(:, :) = 0.0_dp
111 cell%input_recip_to_canonical(:, :) = 0.0_dp
112 cell%symmetry_id = cell_sym_none
113 IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
114 IF (PRESENT(tag)) cell%tag = tag
115
116 END SUBROUTINE cell_create
117
118! **************************************************************************************************
119!> \brief Store the transform between the user input cell and the canonical cell.
120!> \param cell ...
121!> \param hmat_input ...
122!> \param hmat_canonical ...
123! **************************************************************************************************
124 SUBROUTINE cell_finalize_canonical_input(cell, hmat_input, hmat_canonical)
125
126 TYPE(cell_type), POINTER :: cell
127 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: hmat_input, hmat_canonical
128
129 REAL(kind=dp), PARAMETER :: eps_hmat = 1.0e-12_dp
130
131 REAL(kind=dp), DIMENSION(3, 3) :: tmat
132
133 cpassert(ASSOCIATED(cell))
134
135 IF (maxval(abs(hmat_canonical - hmat_input)) <= eps_hmat) THEN
136 cell%input_cell_canonicalized = .false.
137 cell%input_hmat(:, :) = 0.0_dp
138 cell%input_to_canonical(:, :) = 0.0_dp
139 cell%input_recip_to_canonical(:, :) = 0.0_dp
140 ELSE
141 tmat = matmul(hmat_canonical, inv_3x3(hmat_input))
142 cell%input_cell_canonicalized = .true.
143 cell%input_hmat(:, :) = hmat_input(:, :)
144 cell%input_to_canonical(:, :) = tmat(:, :)
145 cell%input_recip_to_canonical(:, :) = transpose(inv_3x3(tmat))
146 END IF
147
148 END SUBROUTINE cell_finalize_canonical_input
149
150! **************************************************************************************************
151!> \brief Canonicalize a general cell matrix without changing lengths and angles.
152!> \param cell ...
153! **************************************************************************************************
154 SUBROUTINE canonicalize_cell_matrix(cell)
155
156 TYPE(cell_type), POINTER :: cell
157
158 REAL(kind=dp), DIMENSION(3) :: abc, cell_angle
159
160 cpassert(ASSOCIATED(cell))
161
162 CALL get_cell(cell=cell, abc=abc)
163 cell_angle(1) = angle(cell%hmat(:, 2), cell%hmat(:, 3))
164 cell_angle(2) = angle(cell%hmat(:, 1), cell%hmat(:, 3))
165 cell_angle(3) = angle(cell%hmat(:, 1), cell%hmat(:, 2))
166
167 CALL set_cell_param(cell, cell_length=abc, cell_angle=cell_angle, &
168 periodic=cell%perd, do_init_cell=.true.)
169
170 END SUBROUTINE canonicalize_cell_matrix
171
172! **************************************************************************************************
173!> \brief Initialise/readjust a simulation cell after hmat has been changed
174!> \param cell ...
175!> \param hmat ...
176!> \param periodic ...
177!> \date 16.01.2002
178!> \author Matthias Krack
179!> \version 1.0
180! **************************************************************************************************
181 SUBROUTINE init_cell(cell, hmat, periodic)
182
183 TYPE(cell_type), POINTER :: cell
184 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN), &
185 OPTIONAL :: hmat
186 INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
187
188 REAL(kind=dp), PARAMETER :: eps_hmat = 1.0e-14_dp
189
190 INTEGER :: dim
191 REAL(kind=dp) :: a, acosa, acosah, acosg, alpha, asina, &
192 asinah, asing, beta, gamma, norm, &
193 norm_c
194 REAL(kind=dp), DIMENSION(3) :: abc
195
196 cpassert(ASSOCIATED(cell))
197
198 IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
199 IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
200
201 cell%deth = abs(det_3x3(cell%hmat))
202
203 IF (cell%deth < 1.0e-10_dp) THEN
204 CALL write_cell_low(cell, "angstrom", default_output_unit)
205 CALL cp_abort(__location__, &
206 "An invalid set of cell vectors was specified. "// &
207 "The cell volume is too small")
208 END IF
209
210 SELECT CASE (cell%symmetry_id)
211 CASE (cell_sym_cubic, &
216 CALL get_cell(cell=cell, abc=abc)
217 abc(2) = plane_distance(0, 1, 0, cell=cell)
218 abc(3) = plane_distance(0, 0, 1, cell=cell)
219 SELECT CASE (cell%symmetry_id)
220 CASE (cell_sym_cubic)
221 abc(1:3) = sum(abc(1:3))/3.0_dp
225 SELECT CASE (cell%symmetry_id)
227 a = 0.5_dp*(abc(1) + abc(2))
228 abc(1) = a
229 abc(2) = a
231 a = 0.5_dp*(abc(1) + abc(3))
232 abc(1) = a
233 abc(3) = a
235 a = 0.5_dp*(abc(2) + abc(3))
236 abc(2) = a
237 abc(3) = a
238 END SELECT
239 END SELECT
240 cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
241 cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
242 cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
244 CALL get_cell(cell=cell, abc=abc)
245 a = 0.5_dp*(abc(1) + abc(2))
246 acosg = 0.5_dp*a
247 asing = sqrt3*acosg
248 IF (cell%symmetry_id == cell_sym_hexagonal_gamma_120) acosg = -acosg
249 cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
250 cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
251 cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
253 CALL get_cell(cell=cell, abc=abc)
254 a = sum(abc(1:3))/3.0_dp
255 alpha = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
256 angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
257 angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
258 acosa = a*cos(alpha)
259 asina = a*sin(alpha)
260 acosah = a*cos(0.5_dp*alpha)
261 asinah = a*sin(0.5_dp*alpha)
262 norm = acosa/acosah
263 norm_c = sqrt(1.0_dp - norm*norm)
264 cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
265 cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
266 cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = a*norm_c
268 CALL get_cell(cell=cell, abc=abc)
269 beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
270 cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*cos(beta)
271 cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
272 cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*sin(beta)
274 ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
275 CALL get_cell(cell=cell, abc=abc)
276 a = 0.5_dp*(abc(1) + abc(2))
277 gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
278 acosg = a*cos(gamma)
279 asing = a*sin(gamma)
280 cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
281 cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
282 cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
283 CASE (cell_sym_triclinic)
284 ! Nothing to do
285 END SELECT
286
287 ! Do we have an (almost) orthorhombic cell?
288 IF ((abs(cell%hmat(1, 2)) < eps_hmat) .AND. (abs(cell%hmat(1, 3)) < eps_hmat) .AND. &
289 (abs(cell%hmat(2, 1)) < eps_hmat) .AND. (abs(cell%hmat(2, 3)) < eps_hmat) .AND. &
290 (abs(cell%hmat(3, 1)) < eps_hmat) .AND. (abs(cell%hmat(3, 2)) < eps_hmat)) THEN
291 cell%orthorhombic = .true.
292 ELSE
293 cell%orthorhombic = .false.
294 END IF
295
296 ! Retain an exact orthorhombic cell
297 ! (off-diagonal elements must remain zero identically to keep QS fast)
298 IF (cell%orthorhombic) THEN
299 cell%hmat(1, 2) = 0.0_dp
300 cell%hmat(1, 3) = 0.0_dp
301 cell%hmat(2, 1) = 0.0_dp
302 cell%hmat(2, 3) = 0.0_dp
303 cell%hmat(3, 1) = 0.0_dp
304 cell%hmat(3, 2) = 0.0_dp
305 END IF
306
307 dim = count(cell%perd == 1)
308 IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
309 cpabort("Non-orthorhombic and not periodic")
310 END IF
311
312 ! Update deth and hmat_inv with enforced symmetry
313 cell%deth = abs(det_3x3(cell%hmat))
314 IF (cell%deth < 1.0e-10_dp) THEN
315 CALL cp_abort(__location__, &
316 "An invalid set of cell vectors was obtained after applying "// &
317 "the requested cell symmetry. The cell volume is too small")
318 END IF
319 cell%h_inv = inv_3x3(cell%hmat)
320
321 END SUBROUTINE init_cell
322
323! **************************************************************************************************
324!> \brief ...
325!> \param cell ...
326!> \param cell_ref ...
327!> \param use_ref_cell ...
328!> \param cell_section ...
329!> \param topology_section ...
330!> \param check_for_ref ...
331!> \param para_env ...
332!> \par History
333!> 03.2005 created [teo]
334!> 03.2026 revamped logic with pdb and extxyz parsers
335!> \author Teodoro Laino
336! **************************************************************************************************
337 RECURSIVE SUBROUTINE read_cell(cell, cell_ref, use_ref_cell, cell_section, &
338 topology_section, check_for_ref, para_env)
339
340 TYPE(cell_type), POINTER :: cell, cell_ref
341 LOGICAL, INTENT(INOUT), OPTIONAL :: use_ref_cell
342 TYPE(section_vals_type), OPTIONAL, POINTER :: cell_section, topology_section
343 LOGICAL, INTENT(IN), OPTIONAL :: check_for_ref
344 TYPE(mp_para_env_type), POINTER :: para_env
345
346 REAL(kind=dp), PARAMETER :: eps = 1.0e-14_dp
347
348 CHARACTER(LEN=default_path_length) :: cell_file_name, coord_file_name, &
349 error_msg
350 INTEGER :: canonicalize_mode, cell_file_format, &
351 coord_file_format, my_per
352 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
353 LOGICAL :: canonicalize_cell, cell_read_a, cell_read_abc, cell_read_alpha_beta_gamma, &
354 cell_read_b, cell_read_c, cell_read_file, my_check_ref, tmp_comb_abc, tmp_comb_cell, &
355 tmp_comb_top, topo_read_coord
356 REAL(kind=dp), DIMENSION(3) :: read_ang, read_len
357 REAL(kind=dp), DIMENSION(3, 3) :: hmat_input, read_mat
358 REAL(kind=dp), DIMENSION(:), POINTER :: cell_par
359 TYPE(cell_type), POINTER :: cell_tmp
360 TYPE(section_vals_type), POINTER :: cell_ref_section
361
362 my_check_ref = .true.
363 NULLIFY (cell_ref_section, cell_par, cell_tmp, multiple_unit_cell)
364 ! cell_tmp has two purposes:
365 ! 1. for transferring matrix of cell vectors from individual
366 ! file parser subroutines to read_mat here, assuming that
367 ! unit conversion has been done in those subroutines;
368 ! 2. for testing whether enforcing symmetry makes a new set
369 ! of cell vectors significantly different from parsed input
370 CALL cell_create(cell_tmp)
371 IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
372 IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
373 IF (PRESENT(check_for_ref)) my_check_ref = check_for_ref
374
375 cell%deth = 0.0_dp
376 cell%orthorhombic = .false.
377 cell%perd(:) = 1
378 cell%symmetry_id = cell_sym_none
379 cell%hmat(:, :) = 0.0_dp
380 cell%h_inv(:, :) = 0.0_dp
381 cell%input_cell_canonicalized = .false.
382 cell%input_hmat(:, :) = 0.0_dp
383 cell%input_to_canonical(:, :) = 0.0_dp
384 cell%input_recip_to_canonical(:, :) = 0.0_dp
385 cell_read_file = .false.
386 cell_read_a = .false.
387 cell_read_b = .false.
388 cell_read_c = .false.
389 cell_read_abc = .false.
390 cell_read_alpha_beta_gamma = .false.
391 hmat_input(:, :) = 0.0_dp
392 read_mat(:, :) = 0.0_dp
393 read_ang(:) = 0.0_dp
394 read_len(:) = 0.0_dp
395
396 ! Precedence of retrieving cell information from input:
397 ! 1. CELL/CELL_FILE_NAME
398 ! 2. CELL/ABC and optionally CELL/ALPHA_BETA_GAMMA
399 ! 3. CELL/A, CELL/B, CELL/C
400 ! 4. TOPOLOGY/COORD_FILE_NAME, if topology_section is present
401 ! The actual order of processing is 4 -> 1 -> 2 -> 3, with
402 ! case 4 merged to case 1 (if file format permits) first.
403 ! Store data into either read_mat or read_ang and read_len
404 ! in CP2K units, which will be converted to cell%hmat and A, B, C.
405 CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
406 CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
407 CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
408 CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
409 CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", explicit=cell_read_alpha_beta_gamma)
410 CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
411 CALL section_vals_val_get(cell_section, "CANONICALIZE", i_val=canonicalize_mode)
412 canonicalize_cell = (canonicalize_mode == canonicalize_cell_true)
413
414 ! Case 4
415 tmp_comb_top = (.NOT. (cell_read_file .OR. cell_read_abc))
416 tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_a))
417 tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_b))
418 tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_c))
419 IF (tmp_comb_top) THEN
420 CALL cp_warn(__location__, &
421 "None of the keywords CELL_FILE_NAME, ABC, or A, B, C "// &
422 "are specified in CELL section. CP2K will now attempt to read "// &
423 "TOPOLOGY/COORD_FILE_NAME if its format can be parsed for "// &
424 "cell information.")
425 IF (ASSOCIATED(topology_section)) THEN
426 CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", explicit=topo_read_coord)
427 IF (topo_read_coord) THEN
428 CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", c_val=coord_file_name)
429 CALL section_vals_val_get(topology_section, "COORD_FILE_FORMAT", i_val=coord_file_format)
430 SELECT CASE (coord_file_format) ! Add formats with both cell and coord parser manually
431 CASE (do_coord_cif)
432 CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
433 CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cif)
434 CASE (do_coord_cp2k)
435 CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
436 CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cp2k)
437 CASE (do_coord_pdb)
438 CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
439 CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_pdb)
440 CASE (do_coord_xyz)
441 CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
442 CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_extxyz)
443 CASE DEFAULT
444 CALL cp_abort(__location__, &
445 "COORD_FILE_FORMAT is not set to one of the implemented "// &
446 "CELL_FILE_FORMAT options and cannot be parsed for cell information!")
447 END SELECT
448 ELSE
449 CALL cp_abort(__location__, &
450 "COORD_FILE_NAME is not set, so no cell information is available!")
451 END IF
452 ELSE
453 CALL cp_warn(__location__, &
454 "TOPOLOGY section is not available, so COORD_FILE_NAME cannot "// &
455 "be parsed for cell information in lieu of missing CELL settings.")
456 END IF
457 END IF
458 ! Former logic in SUBROUTINE read_cell_from_external_file is moved here
459 CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
460 IF (cell_read_file) THEN ! Case 1
461 tmp_comb_cell = (cell_read_abc .OR. (cell_read_a .OR. (cell_read_b .OR. cell_read_c)))
462 IF (tmp_comb_cell) &
463 CALL cp_warn(__location__, &
464 "Cell Information provided through A, B, C, or ABC in conjunction "// &
465 "with CELL_FILE_NAME. The definition in external file will override "// &
466 "other ones.")
467 CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
468 CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=cell_file_format)
469 SELECT CASE (cell_file_format)
470 CASE (do_cell_cp2k)
471 CALL read_cell_cp2k(cell_file_name, cell_tmp, para_env)
472 CASE (do_cell_xsc)
473 CALL read_cell_xsc(cell_file_name, cell_tmp, para_env)
474 CASE (do_cell_extxyz)
475 CALL read_cell_extxyz(cell_file_name, cell_tmp, para_env)
476 CASE (do_cell_pdb)
477 CALL read_cell_pdb(cell_file_name, cell_tmp, para_env)
478 CASE (do_cell_cif)
479 CALL read_cell_cif(cell_file_name, cell_tmp, para_env)
480 CASE DEFAULT
481 CALL cp_abort(__location__, &
482 "CELL_FILE_FORMAT is not set to one of the implemented "// &
483 "options and cannot be parsed for cell information!")
484 END SELECT
485 read_mat = cell_tmp%hmat
486 ELSE
487 IF (cell_read_abc) THEN ! Case 2
488 CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
489 read_len = cell_par
490 CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_par)
491 read_ang = cell_par
492 IF (cell_read_a .OR. cell_read_b .OR. cell_read_c) &
493 CALL cp_warn(__location__, &
494 "Cell information provided through vectors A, B or C in conjunction with ABC. "// &
495 "The definition of the ABC keyword will override the one provided by A, B and C.")
496 ELSE ! Case 3
497 tmp_comb_abc = ((cell_read_a .EQV. cell_read_b) .AND. (cell_read_b .EQV. cell_read_c))
498 IF (tmp_comb_abc) THEN
499 CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
500 read_mat(:, 1) = cell_par(:)
501 CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
502 read_mat(:, 2) = cell_par(:)
503 CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
504 read_mat(:, 3) = cell_par(:)
505 IF (cell_read_alpha_beta_gamma) &
506 CALL cp_warn(__location__, &
507 "The keyword ALPHA_BETA_GAMMA is ignored because it was used without the "// &
508 "keyword ABC.")
509 ELSE
510 CALL cp_abort(__location__, &
511 "Neither of the keywords CELL_FILE_NAME or ABC are specified, "// &
512 "and cell vector settings in A, B, C are incomplete!")
513 END IF
514 END IF
515 END IF
516
517 ! Convert read_mat or read_len and read_ang to actual cell%hmat
518 IF (any(read_mat(:, :) > eps)) THEN
519 ! Make a warning before storing cell vectors that
520 ! do not form a triangular matrix.
521 IF (.NOT. canonicalize_cell .AND. &
522 ((abs(read_mat(2, 1)) > eps) .OR. &
523 (abs(read_mat(3, 1)) > eps) .OR. &
524 (abs(read_mat(3, 2)) > eps))) THEN
525 IF (canonicalize_mode == canonicalize_cell_auto) THEN
526 CALL cp_warn(__location__, &
527 "CELL%CANONICALIZE AUTO keeps the general input cell orientation. "// &
528 "The cell matrix is not a lower triangle and does not conform to the "// &
529 "program convention that A lies along the X-axis and B is in the XY plane. "// &
530 "Set CELL%CANONICALIZE TRUE to explicitly transform the cell and supported "// &
531 "cell-dependent input to the canonical internal frame.")
532 ELSE
533 CALL cp_warn(__location__, &
534 "Cell vectors are read but cell matrix is not "// &
535 "a lower triangle, not conforming to the program "// &
536 "convention that A lies along the X-axis and "// &
537 "B is in the XY plane.")
538 END IF
539 END IF
540 cell%hmat = read_mat
541 ELSE
542 IF (any(read_ang(:) > eps) .AND. any(read_len(:) > eps)) THEN
543 CALL set_cell_param(cell, cell_length=read_len, cell_angle=read_ang, &
544 do_init_cell=.false.)
545 ELSE
546 CALL cp_abort(__location__, &
547 "No meaningful cell information is read from parser!")
548 END IF
549 END IF
550 ! Reset cell section so that only A, B, C are kept
551 CALL reset_cell_section_by_cell_mat(cell, cell_section)
552
553 ! Multiple unit cell
554 CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
555 IF (any(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
556
557 CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
558 SELECT CASE (my_per)
559 CASE (use_perd_x)
560 cell%perd = [1, 0, 0]
561 CASE (use_perd_y)
562 cell%perd = [0, 1, 0]
563 CASE (use_perd_z)
564 cell%perd = [0, 0, 1]
565 CASE (use_perd_xy)
566 cell%perd = [1, 1, 0]
567 CASE (use_perd_xz)
568 cell%perd = [1, 0, 1]
569 CASE (use_perd_yz)
570 cell%perd = [0, 1, 1]
571 CASE (use_perd_xyz)
572 cell%perd = [1, 1, 1]
573 CASE (use_perd_none)
574 cell%perd = [0, 0, 0]
575 CASE DEFAULT
576 cpabort("Invalid or not yet implemented cell periodicity")
577 END SELECT
578
579 ! Load requested cell symmetry
580 CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
581 ! Try enforcing symmetry by initializing a temporary copy of cell
582 ! and see if the resulting cell matrix differ significantly
583 hmat_input(:, :) = cell%hmat(:, :)
584 CALL cell_clone(cell, cell_tmp)
585 CALL init_cell(cell_tmp)
586 IF (.NOT. canonicalize_cell .AND. any(abs(cell_tmp%hmat - cell%hmat) > eps)) THEN
587 WRITE (unit=error_msg, fmt="(A)") &
588 "When initializing cell vectors with requested symmetry, one "// &
589 "or more elements of the cell matrix has varied significantly. "// &
590 "The input parameters are either deviating from the symmetry, "// &
591 "or not conforming to the program convention that cell matrix "// &
592 "is a lower triangle. The symmetrized cell vectors will be used "// &
593 "anyway with the input atomic coordinates."
594 CALL cp_warn(__location__, error_msg)
595 END IF
596 IF (canonicalize_cell) THEN
597 CALL canonicalize_cell_matrix(cell_tmp)
598 CALL cell_finalize_canonical_input(cell_tmp, hmat_input, cell_tmp%hmat)
599 END IF
600 CALL cell_clone(cell_tmp, cell)
601 CALL cell_release(cell_tmp)
602 CALL reset_cell_section_by_cell_mat(cell, cell_section)
603
604 IF (my_check_ref) THEN
605 ! Recursive check for reference cell requested
606 cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
607 IF (parsed_cp2k_input(cell_ref_section, check_this_section=.true.)) THEN
608 IF (PRESENT(use_ref_cell)) use_ref_cell = .true.
609 CALL read_cell(cell_ref, cell_ref, use_ref_cell=use_ref_cell, &
610 cell_section=cell_ref_section, check_for_ref=.false., &
611 para_env=para_env)
612 ELSE
613 CALL cell_clone(cell, cell_ref, tag="CELL_REF")
614 IF (PRESENT(use_ref_cell)) use_ref_cell = .false.
615 END IF
616 END IF
617
618 END SUBROUTINE read_cell
619
620! **************************************************************************************************
621!> \brief utility function to ease the transition to the new input.
622!> returns true if the new input was parsed
623!> \param input_file the parsed input file
624!> \param check_this_section ...
625!> \return ...
626!> \author fawzi
627! **************************************************************************************************
628 FUNCTION parsed_cp2k_input(input_file, check_this_section) RESULT(res)
629
630 TYPE(section_vals_type), POINTER :: input_file
631 LOGICAL, INTENT(IN), OPTIONAL :: check_this_section
632 LOGICAL :: res
633
634 LOGICAL :: my_check
635 TYPE(section_vals_type), POINTER :: glob_section
636
637 my_check = .false.
638 IF (PRESENT(check_this_section)) my_check = check_this_section
639 res = ASSOCIATED(input_file)
640 IF (res) THEN
641 cpassert(input_file%ref_count > 0)
642 IF (.NOT. my_check) THEN
643 glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
644 CALL section_vals_get(glob_section, explicit=res)
645 ELSE
646 CALL section_vals_get(input_file, explicit=res)
647 END IF
648 END IF
649
650 END FUNCTION parsed_cp2k_input
651
652! **************************************************************************************************
653!> \brief Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma)
654!> using the convention: a parallel to the x axis, b in the x-y plane and
655!> and c univoquely determined; gamma is the angle between a and b; beta
656!> is the angle between c and a and alpha is the angle between c and b
657!> \param cell ...
658!> \param cell_length ...
659!> \param cell_angle ...
660!> \param periodic ...
661!> \param do_init_cell ...
662!> \date 03.2008
663!> \author Teodoro Laino
664! **************************************************************************************************
665 SUBROUTINE set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
666
667 TYPE(cell_type), POINTER :: cell
668 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: cell_length, cell_angle
669 INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
670 LOGICAL, INTENT(IN) :: do_init_cell
671
672 REAL(kind=dp), PARAMETER :: eps = epsilon(0.0_dp)
673
674 REAL(kind=dp) :: cos_alpha, cos_beta, cos_gamma, sin_gamma
675
676 cpassert(ASSOCIATED(cell))
677 cpassert(all(cell_angle /= 0.0_dp))
678
679 cos_gamma = cos(cell_angle(3)); IF (abs(cos_gamma) < eps) cos_gamma = 0.0_dp
680 IF (abs(abs(cos_gamma) - 1.0_dp) < eps) cos_gamma = sign(1.0_dp, cos_gamma)
681 sin_gamma = sin(cell_angle(3)); IF (abs(sin_gamma) < eps) sin_gamma = 0.0_dp
682 IF (abs(abs(sin_gamma) - 1.0_dp) < eps) sin_gamma = sign(1.0_dp, sin_gamma)
683 cos_beta = cos(cell_angle(2)); IF (abs(cos_beta) < eps) cos_beta = 0.0_dp
684 IF (abs(abs(cos_beta) - 1.0_dp) < eps) cos_beta = sign(1.0_dp, cos_beta)
685 cos_alpha = cos(cell_angle(1)); IF (abs(cos_alpha) < eps) cos_alpha = 0.0_dp
686 IF (abs(abs(cos_alpha) - 1.0_dp) < eps) cos_alpha = sign(1.0_dp, cos_alpha)
687
688 cell%hmat(:, 1) = [1.0_dp, 0.0_dp, 0.0_dp]
689 cell%hmat(:, 2) = [cos_gamma, sin_gamma, 0.0_dp]
690 cell%hmat(:, 3) = [cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp]
691 cell%hmat(3, 3) = sqrt(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
692
693 cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
694 cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
695 cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
696
697 IF (do_init_cell) THEN
698 IF (PRESENT(periodic)) THEN
699 CALL init_cell(cell=cell, periodic=periodic)
700 ELSE
701 CALL init_cell(cell=cell)
702 END IF
703 END IF
704
705 END SUBROUTINE set_cell_param
706
707! **************************************************************************************************
708!> \brief Setup of the multiple unit_cell
709!> \param cell ...
710!> \param multiple_unit_cell ...
711!> \date 05.2009
712!> \author Teodoro Laino [tlaino]
713!> \version 1.0
714! **************************************************************************************************
715 SUBROUTINE set_multiple_unit_cell(cell, multiple_unit_cell)
716
717 TYPE(cell_type), POINTER :: cell
718 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
719
720 cpassert(ASSOCIATED(cell))
721
722 ! Abort, if one of the value is set to zero
723 IF (any(multiple_unit_cell <= 0)) &
724 CALL cp_abort(__location__, &
725 "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
726 "A value of 0 or negative is meaningless!")
727
728 ! Scale abc according to user request
729 cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
730 cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
731 cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
732
733 END SUBROUTINE set_multiple_unit_cell
734
735! **************************************************************************************************
736!> \brief Reads cell information from CIF file
737!> \param cif_file_name ...
738!> \param cell ...
739!> \param para_env ...
740!> \date 12.2008
741!> \par Format Information implemented:
742!> _cell_length_a (_cell.length_a)
743!> _cell_length_b (_cell.length_b)
744!> _cell_length_c (_cell.length_c)
745!> _cell_angle_alpha (_cell.length_alpha)
746!> _cell_angle_beta (_cell.length_beta)
747!> _cell_angle_gamma (_cell.length_gamma)
748!>
749!> \author Teodoro Laino [tlaino]
750!> moved from topology_cif (1/2019 JHU)
751! **************************************************************************************************
752 SUBROUTINE read_cell_cif(cif_file_name, cell, para_env)
753
754 CHARACTER(len=*) :: cif_file_name
755 TYPE(cell_type), POINTER :: cell
756 TYPE(mp_para_env_type), POINTER :: para_env
757
758 CHARACTER(len=*), PARAMETER :: routinen = 'read_cell_cif'
759
760 INTEGER :: handle
761 INTEGER, DIMENSION(3) :: periodic
762 LOGICAL :: found
763 REAL(kind=dp), DIMENSION(3) :: cell_angles, cell_lengths
764 TYPE(cp_parser_type) :: parser
765
766 CALL timeset(routinen, handle)
767
768 CALL parser_create(parser, cif_file_name, &
769 para_env=para_env, apply_preprocessing=.false.)
770
771 ! Parsing cell infos
772 periodic = 1
773 ! Check for _cell_length_a or _cell.length_a
774 CALL parser_search_string(parser, "_cell_length_a", ignore_case=.false., found=found, &
775 begin_line=.false., search_from_begin_of_file=.true.)
776 IF (.NOT. found) THEN
777 CALL parser_search_string(parser, "_cell.length_a", ignore_case=.false., found=found, &
778 begin_line=.false., search_from_begin_of_file=.true.)
779 IF (.NOT. found) &
780 cpabort("The field _cell_length_a or _cell.length_a was not found in CIF file! ")
781 END IF
782 CALL cif_get_real(parser, cell_lengths(1))
783 cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
784
785 ! Check for _cell_length_b or _cell.length_b
786 CALL parser_search_string(parser, "_cell_length_b", ignore_case=.false., found=found, &
787 begin_line=.false., search_from_begin_of_file=.true.)
788 IF (.NOT. found) THEN
789 CALL parser_search_string(parser, "_cell.length_b", ignore_case=.false., found=found, &
790 begin_line=.false., search_from_begin_of_file=.true.)
791 IF (.NOT. found) &
792 cpabort("The field _cell_length_b or _cell.length_b was not found in CIF file! ")
793 END IF
794 CALL cif_get_real(parser, cell_lengths(2))
795 cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
796
797 ! Check for _cell_length_c or _cell.length_c
798 CALL parser_search_string(parser, "_cell_length_c", ignore_case=.false., found=found, &
799 begin_line=.false., search_from_begin_of_file=.true.)
800 IF (.NOT. found) THEN
801 CALL parser_search_string(parser, "_cell.length_c", ignore_case=.false., found=found, &
802 begin_line=.false., search_from_begin_of_file=.true.)
803 IF (.NOT. found) &
804 cpabort("The field _cell_length_c or _cell.length_c was not found in CIF file! ")
805 END IF
806 CALL cif_get_real(parser, cell_lengths(3))
807 cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
808
809 ! Check for _cell_angle_alpha or _cell.angle_alpha
810 CALL parser_search_string(parser, "_cell_angle_alpha", ignore_case=.false., found=found, &
811 begin_line=.false., search_from_begin_of_file=.true.)
812 IF (.NOT. found) THEN
813 CALL parser_search_string(parser, "_cell.angle_alpha", ignore_case=.false., found=found, &
814 begin_line=.false., search_from_begin_of_file=.true.)
815 IF (.NOT. found) &
816 cpabort("The field _cell_angle_alpha or _cell.angle_alpha was not found in CIF file! ")
817 END IF
818 CALL cif_get_real(parser, cell_angles(1))
819 cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
820
821 ! Check for _cell_angle_beta or _cell.angle_beta
822 CALL parser_search_string(parser, "_cell_angle_beta", ignore_case=.false., found=found, &
823 begin_line=.false., search_from_begin_of_file=.true.)
824 IF (.NOT. found) THEN
825 CALL parser_search_string(parser, "_cell.angle_beta", ignore_case=.false., found=found, &
826 begin_line=.false., search_from_begin_of_file=.true.)
827 IF (.NOT. found) &
828 cpabort("The field _cell_angle_beta or _cell.angle_beta was not found in CIF file! ")
829 END IF
830 CALL cif_get_real(parser, cell_angles(2))
831 cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
832
833 ! Check for _cell_angle_gamma or _cell.angle_gamma
834 CALL parser_search_string(parser, "_cell_angle_gamma", ignore_case=.false., found=found, &
835 begin_line=.false., search_from_begin_of_file=.true.)
836 IF (.NOT. found) THEN
837 CALL parser_search_string(parser, "_cell.angle_gamma", ignore_case=.false., found=found, &
838 begin_line=.false., search_from_begin_of_file=.true.)
839 IF (.NOT. found) &
840 cpabort("The field _cell_angle_gamma or _cell.angle_gamma was not found in CIF file! ")
841 END IF
842 CALL cif_get_real(parser, cell_angles(3))
843 cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
844
845 ! Create cell
846 CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
847 do_init_cell=.true.)
848
849 CALL parser_release(parser)
850
851 CALL timestop(handle)
852
853 END SUBROUTINE read_cell_cif
854
855! **************************************************************************************************
856!> \brief Reads REAL from the CIF file.. This wrapper is needed in order to
857!> treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
858!> \param parser ...
859!> \param r ...
860!> \date 12.2008
861!> \author Teodoro Laino [tlaino]
862! **************************************************************************************************
863 SUBROUTINE cif_get_real(parser, r)
864
865 TYPE(cp_parser_type), INTENT(INOUT) :: parser
866 REAL(kind=dp), INTENT(OUT) :: r
867
868 CHARACTER(LEN=default_string_length) :: s_tag
869 INTEGER :: iln
870
871 CALL parser_get_object(parser, s_tag)
872 iln = len_trim(s_tag)
873 IF (index(s_tag, "(") /= 0) iln = index(s_tag, "(") - 1
874 READ (s_tag(1:iln), *) r
875
876 END SUBROUTINE cif_get_real
877
878! **************************************************************************************************
879!> \brief Reads cell information from comment line of extended xyz file
880!> \param extxyz_file_name ...
881!> \param cell ...
882!> \param para_env ...
883!> \date 03.2026
884!> \par Extended xyz format has comment on the second line in the form as
885!> Lattice="Ax Ay Az Bx By Bz Cx Cy Cz"
886!> where Ax, Ay, Az are three Cartesian components of cell vector A,
887!> Bx, By, Bz are components of B, Cx, Cy, Cz are components of C,
888!> all in the unit of angstrom.
889! **************************************************************************************************
890 SUBROUTINE read_cell_extxyz(extxyz_file_name, cell, para_env)
891
892 CHARACTER(len=*) :: extxyz_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_extxyz'
897
898 CHARACTER(len=default_path_length) :: raw_cell_str
899 INTEGER :: handle, i, id1, id2, ios, j
900 REAL(kind=dp), DIMENSION(3, 3) :: hmat
901 TYPE(cp_parser_type) :: parser
902
903 CALL timeset(routinen, handle)
904
905 CALL parser_create(parser, extxyz_file_name, &
906 para_env=para_env, apply_preprocessing=.false.)
907 CALL parser_get_next_line(parser, 2)
908 CALL uppercase(parser%input_line)
909 id1 = index(parser%input_line, "LATTICE=")
910 IF (id1 > 0) THEN
911 id2 = index(parser%input_line(id1 + 9:), '"')
912 READ (parser%input_line(id1 + 9:id1 + id2 + 7), '(A)') raw_cell_str
913 READ (raw_cell_str, *, iostat=ios) hmat(:, 1), hmat(:, 2), hmat(:, 3)
914 IF (ios /= 0) THEN
915 CALL cp_abort(__location__, "Error while parsing extended XYZ file "// &
916 "<"//trim(extxyz_file_name)//"> for cell vectors: "// &
917 "found <lattice=> field as <"//trim(raw_cell_str)//">")
918 END IF
919 DO i = 1, 3
920 DO j = 1, 3
921 cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
922 END DO
923 END DO
924
925 ELSE
926 CALL cp_abort(__location__, &
927 "The field <lattice=> was not found on comment line "// &
928 "of XYZ file, so cell information cannot be set via "// &
929 "extended XYZ specification! ")
930 END IF
931
932 CALL parser_release(parser)
933
934 CALL timestop(handle)
935
936 END SUBROUTINE read_cell_extxyz
937
938! **************************************************************************************************
939!> \brief Reads cell information from CRYST1 record of PDB file
940!> \param pdb_file_name ...
941!> \param cell ...
942!> \param para_env ...
943!> \date 03.2026
944!> \par CRYST1 record may contain space group and Z value at the end,
945!> but here only the first entries are read:
946!> COLUMNS DATA TYPE FIELD DEFINITION
947!> -------------------------------------------------------------
948!> 1 - 6 Record name "CRYST1"
949!> 7 - 15 Real(9.3) a a (Angstroms).
950!> 16 - 24 Real(9.3) b b (Angstroms).
951!> 25 - 33 Real(9.3) c c (Angstroms).
952!> 34 - 40 Real(7.2) alpha alpha (degrees).
953!> 41 - 47 Real(7.2) beta beta (degrees).
954!> 48 - 54 Real(7.2) gamma gamma (degrees).
955! **************************************************************************************************
956 SUBROUTINE read_cell_pdb(pdb_file_name, cell, para_env)
957
958 CHARACTER(len=*) :: pdb_file_name
959 TYPE(cell_type), POINTER :: cell
960 TYPE(mp_para_env_type), POINTER :: para_env
961
962 CHARACTER(len=*), PARAMETER :: routinen = 'read_cell_pdb'
963
964 CHARACTER(LEN=default_string_length) :: cryst
965 INTEGER :: handle, i, ios
966 INTEGER, DIMENSION(3) :: periodic
967 LOGICAL :: found
968 REAL(kind=dp), DIMENSION(3) :: cell_angles, cell_lengths
969 TYPE(cp_parser_type) :: parser
970
971 CALL timeset(routinen, handle)
972
973 CALL parser_create(parser, pdb_file_name, &
974 para_env=para_env, apply_preprocessing=.false.)
975
976 CALL parser_search_string(parser, "CRYST1", ignore_case=.false., found=found, &
977 begin_line=.true., search_from_begin_of_file=.true.)
978 IF (.NOT. found) &
979 cpabort("The line <CRYST1> was not found in PDB file! ")
980
981 periodic = 1
982 READ (parser%input_line, *, iostat=ios) cryst, cell_lengths(:), cell_angles(:)
983 IF (ios /= 0) THEN
984 CALL cp_abort(__location__, "Error while parsing PDB file "// &
985 "<"//trim(pdb_file_name)//"> for cell lengths and angles: "// &
986 "found CRYST1 line as <"//trim(parser%input_line)//">")
987 END IF
988 DO i = 1, 3
989 cell_lengths(i) = cp_unit_to_cp2k(cell_lengths(i), "angstrom")
990 cell_angles(i) = cp_unit_to_cp2k(cell_angles(i), "deg")
991 END DO
992 CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
993 do_init_cell=.true.)
994
995 CALL parser_release(parser)
996
997 CALL timestop(handle)
998
999 END SUBROUTINE read_cell_pdb
1000
1001! **************************************************************************************************
1002!> \brief Reads cell information from cp2k file
1003!> \param cp2k_file_name ...
1004!> \param cell ...
1005!> \param para_env ...
1006!> \date 03.2026
1007!> \par Isolated from former read_cell_from_external_file
1008! **************************************************************************************************
1009 SUBROUTINE read_cell_cp2k(cp2k_file_name, cell, para_env)
1010
1011 CHARACTER(len=*) :: cp2k_file_name
1012 TYPE(cell_type), POINTER :: cell
1013 TYPE(mp_para_env_type), POINTER :: para_env
1014
1015 CHARACTER(len=*), PARAMETER :: routinen = 'read_cell_cp2k'
1016
1017 INTEGER :: handle, i, idum, j
1018 LOGICAL :: my_end
1019 REAL(kind=dp) :: xdum
1020 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1021 TYPE(cp_parser_type) :: parser
1022
1023 CALL timeset(routinen, handle)
1024
1025 CALL parser_create(parser, cp2k_file_name, &
1026 para_env=para_env, apply_preprocessing=.false.)
1027
1028 CALL parser_get_next_line(parser, 1)
1029 my_end = .false.
1030 DO WHILE (.NOT. my_end)
1031 READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
1032 CALL parser_get_next_line(parser, 1, at_end=my_end)
1033 END DO
1034 DO i = 1, 3
1035 DO j = 1, 3
1036 cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
1037 END DO
1038 END DO
1039
1040 CALL parser_release(parser)
1041
1042 CALL timestop(handle)
1043
1044 END SUBROUTINE read_cell_cp2k
1045
1046! **************************************************************************************************
1047!> \brief Reads cell information from xsc file
1048!> \param xsc_file_name ...
1049!> \param cell ...
1050!> \param para_env ...
1051!> \date 03.2026
1052!> \par Isolated from former read_cell_from_external_file
1053! **************************************************************************************************
1054 SUBROUTINE read_cell_xsc(xsc_file_name, cell, para_env)
1055
1056 CHARACTER(len=*) :: xsc_file_name
1057 TYPE(cell_type), POINTER :: cell
1058 TYPE(mp_para_env_type), POINTER :: para_env
1059
1060 CHARACTER(len=*), PARAMETER :: routinen = 'read_cell_xsc'
1061
1062 INTEGER :: handle, i, idum, j
1063 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1064 TYPE(cp_parser_type) :: parser
1065
1066 CALL timeset(routinen, handle)
1067
1068 CALL parser_create(parser, xsc_file_name, &
1069 para_env=para_env, apply_preprocessing=.false.)
1070
1071 CALL parser_get_next_line(parser, 1)
1072 READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
1073 DO i = 1, 3
1074 DO j = 1, 3
1075 cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
1076 END DO
1077 END DO
1078
1079 CALL parser_release(parser)
1080
1081 CALL timestop(handle)
1082
1083 END SUBROUTINE read_cell_xsc
1084
1085! **************************************************************************************************
1086!> \brief Reset cell section by matrix in cell-type pointer
1087!> \param cell ...
1088!> \param cell_section ...
1089!> \date 03.2026
1090!> \par Alternative keywords for cell settings will be unset
1091!> except MULTIPLE_UNIT_CELL, PERIODIC and SYMMETRY.
1092! **************************************************************************************************
1093 SUBROUTINE reset_cell_section_by_cell_mat(cell, cell_section)
1094
1095 TYPE(cell_type), POINTER :: cell
1096 TYPE(section_vals_type), POINTER :: cell_section
1097
1098 REAL(kind=dp), DIMENSION(:), POINTER :: cell_par
1099
1100 CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
1101 CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
1102 CALL section_vals_val_unset(cell_section, "ABC")
1103 CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
1104 CALL section_vals_val_unset(cell_section, "A")
1105 CALL section_vals_val_unset(cell_section, "B")
1106 CALL section_vals_val_unset(cell_section, "C")
1107 ALLOCATE (cell_par(3))
1108 cell_par = cell%hmat(:, 1)
1109 CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
1110 ALLOCATE (cell_par(3))
1111 cell_par = cell%hmat(:, 2)
1112 CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
1113 ALLOCATE (cell_par(3))
1114 cell_par = cell%hmat(:, 3)
1115 CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
1116
1117 END SUBROUTINE reset_cell_section_by_cell_mat
1118
1119! **************************************************************************************************
1120!> \brief Write the cell parameters to the output unit.
1121!> \param cell ...
1122!> \param subsys_section ...
1123!> \param tag ...
1124!> \date 02.06.2000
1125!> \par History
1126!> - 11.2008 Teodoro Laino [tlaino] - rewrite and enabling user driven units
1127!> \author Matthias Krack
1128!> \version 1.0
1129! **************************************************************************************************
1130 SUBROUTINE write_cell(cell, subsys_section, tag)
1131
1132 TYPE(cell_type), POINTER :: cell
1133 TYPE(section_vals_type), POINTER :: subsys_section
1134 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
1135
1136 CHARACTER(LEN=default_string_length) :: label, unit_str
1137 INTEGER :: output_unit
1138 TYPE(cp_logger_type), POINTER :: logger
1139
1140 NULLIFY (logger)
1141 logger => cp_get_default_logger()
1142 IF (PRESENT(tag)) THEN
1143 label = trim(tag)//"|"
1144 ELSE
1145 label = trim(cell%tag)//"|"
1146 END IF
1147
1148 output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", extension=".Log")
1149 CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
1150 CALL write_cell_low(cell, unit_str, output_unit, label)
1151 CALL cp_print_key_finished_output(output_unit, logger, subsys_section, "PRINT%CELL")
1152
1153 END SUBROUTINE write_cell
1154
1155! **************************************************************************************************
1156!> \brief Write the cell parameters to the output unit
1157!> \param cell ...
1158!> \param unit_str ...
1159!> \param output_unit ...
1160!> \param label ...
1161!> \date 17.05.2023
1162!> \par History
1163!> - Extracted from write_cell (17.05.2023, MK)
1164!> \version 1.0
1165! **************************************************************************************************
1166 SUBROUTINE write_cell_low(cell, unit_str, output_unit, label)
1167
1168 TYPE(cell_type), POINTER :: cell
1169 CHARACTER(LEN=*), INTENT(IN) :: unit_str
1170 INTEGER, INTENT(IN) :: output_unit
1171 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: label
1172
1173 CHARACTER(LEN=12) :: tag
1174 CHARACTER(LEN=3) :: string
1175 CHARACTER(LEN=default_string_length) :: my_label
1176 REAL(kind=dp) :: alpha, beta, gamma, val
1177 REAL(kind=dp), DIMENSION(3) :: abc
1178 TYPE(enumeration_type), POINTER :: enum
1179 TYPE(keyword_type), POINTER :: keyword
1180 TYPE(section_type), POINTER :: section
1181
1182 NULLIFY (enum)
1183 NULLIFY (keyword)
1184 NULLIFY (section)
1185
1186 IF (output_unit > 0) THEN
1187 CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma, tag=tag)
1188 IF (PRESENT(label)) THEN
1189 my_label = label
1190 ELSE
1191 my_label = trim(tag)//"|"
1192 END IF
1193 val = cp_unit_from_cp2k(cell%deth, trim(unit_str)//"^3")
1194 WRITE (unit=output_unit, fmt="(/,T2,A,T61,F20.6)") &
1195 trim(my_label)//" Volume ["//trim(unit_str)//"^3]:", val
1196 val = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
1197 WRITE (unit=output_unit, fmt="(T2,A,T30,3F10.3,3X,A6,F12.6)") &
1198 trim(my_label)//" Vector a ["//trim(unit_str)//"]:", cell%hmat(:, 1)*val, &
1199 "|a| = ", abc(1)*val, &
1200 trim(my_label)//" Vector b ["//trim(unit_str)//"]:", cell%hmat(:, 2)*val, &
1201 "|b| = ", abc(2)*val, &
1202 trim(my_label)//" Vector c ["//trim(unit_str)//"]:", cell%hmat(:, 3)*val, &
1203 "|c| = ", abc(3)*val
1204 WRITE (unit=output_unit, fmt="(T2,A,T69,F12.6)") &
1205 trim(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
1206 trim(my_label)//" Angle (a,c), beta [degree]: ", beta, &
1207 trim(my_label)//" Angle (a,b), gamma [degree]: ", gamma
1208 IF (cell%symmetry_id /= cell_sym_none) THEN
1209 CALL create_cell_section(section)
1210 keyword => section_get_keyword(section, "SYMMETRY")
1211 CALL keyword_get(keyword, enum=enum)
1212 WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
1213 trim(my_label)//" Requested initial symmetry: ", &
1214 adjustr(trim(enum_i2c(enum, cell%symmetry_id)))
1215 CALL section_release(section)
1216 END IF
1217 IF (cell%orthorhombic) THEN
1218 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
1219 trim(my_label)//" Numerically orthorhombic: ", "YES"
1220 ELSE
1221 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
1222 trim(my_label)//" Numerically orthorhombic: ", " NO"
1223 END IF
1224 IF (sum(cell%perd(1:3)) == 0) THEN
1225 WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
1226 trim(my_label)//" Periodicity", "NONE"
1227 ELSE
1228 string = ""
1229 IF (cell%perd(1) == 1) string = trim(string)//"X"
1230 IF (cell%perd(2) == 1) string = trim(string)//"Y"
1231 IF (cell%perd(3) == 1) string = trim(string)//"Z"
1232 WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
1233 trim(my_label)//" Periodicity", adjustr(string)
1234 END IF
1235 END IF
1236
1237 END SUBROUTINE write_cell_low
1238
1239END MODULE cell_methods
Handles all functions related to the CELL.
subroutine, public write_cell_low(cell, unit_str, output_unit, label)
Write the cell parameters to the output unit.
subroutine, public write_cell(cell, subsys_section, tag)
Write the cell parameters to the output unit.
subroutine, public read_cell_cp2k(cp2k_file_name, cell, para_env)
Reads cell information from cp2k file.
subroutine, public read_cell_cif(cif_file_name, cell, para_env)
Reads cell information from CIF file.
subroutine, public set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma) using the convention: a parall...
subroutine, public read_cell_pdb(pdb_file_name, cell, para_env)
Reads cell information from CRYST1 record of PDB file.
subroutine, public cell_finalize_canonical_input(cell, hmat_input, hmat_canonical)
Store the transform between the user input cell and the canonical cell.
subroutine, public read_cell_extxyz(extxyz_file_name, cell, para_env)
Reads cell information from comment line of extended xyz file.
recursive subroutine, public read_cell(cell, cell_ref, use_ref_cell, cell_section, topology_section, check_for_ref, para_env)
...
subroutine, public read_cell_xsc(xsc_file_name, cell, para_env)
Reads cell information from xsc file.
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Handles all functions related to the CELL.
Definition cell_types.F:15
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
integer, parameter, public cell_sym_monoclinic
Definition cell_types.F:29
integer, parameter, public use_perd_y
Definition cell_types.F:42
integer, parameter, public cell_sym_triclinic
Definition cell_types.F:29
integer, parameter, public cell_sym_tetragonal_ab
Definition cell_types.F:29
integer, parameter, public use_perd_xz
Definition cell_types.F:42
integer, parameter, public cell_sym_rhombohedral
Definition cell_types.F:29
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:608
integer, parameter, public use_perd_x
Definition cell_types.F:42
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:118
integer, parameter, public cell_sym_tetragonal_ac
Definition cell_types.F:29
integer, parameter, public use_perd_z
Definition cell_types.F:42
integer, parameter, public use_perd_yz
Definition cell_types.F:42
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:210
integer, parameter, public use_perd_none
Definition cell_types.F:42
integer, parameter, public cell_sym_hexagonal_gamma_60
Definition cell_types.F:29
integer, parameter, public cell_sym_orthorhombic
Definition cell_types.F:29
integer, parameter, public cell_sym_none
Definition cell_types.F:29
integer, parameter, public cell_sym_hexagonal_gamma_120
Definition cell_types.F:29
integer, parameter, public cell_sym_monoclinic_gamma_ab
Definition cell_types.F:29
integer, parameter, public cell_sym_cubic
Definition cell_types.F:29
integer, parameter, public use_perd_xy
Definition cell_types.F:42
integer, parameter, public cell_sym_tetragonal_bc
Definition cell_types.F:29
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition cell_types.F:301
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
unit conversion facility
Definition cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1239
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1210
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_cell_extxyz
integer, parameter, public do_cell_pdb
integer, parameter, public canonicalize_cell_true
integer, parameter, public do_coord_cif
integer, parameter, public do_cell_cif
integer, parameter, public do_cell_xsc
integer, parameter, public do_coord_pdb
integer, parameter, public do_cell_cp2k
integer, parameter, public do_coord_xyz
integer, parameter, public canonicalize_cell_auto
integer, parameter, public do_coord_cp2k
builds the subsystem section of the input
subroutine, public create_cell_section(section, periodic)
creates the cell section
represents an enumeration, i.e. a mapping between integers and strings
character(len=default_string_length) function, public enum_i2c(enum, i)
maps an integer to a string
represents keywords in an input
subroutine, public keyword_get(keyword, names, usage, description, type_of_var, n_var, default_value, lone_keyword_value, repeats, enum, citations)
...
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_unset(section_vals, keyword_name, i_rep_section, i_rep_val)
unsets (removes) the requested value (if it is a keyword repetitions removes the repetition,...
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
recursive type(keyword_type) function, pointer, public section_get_keyword(section, keyword_name)
returns the requested keyword
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
integer, parameter, public default_path_length
Definition kinds.F:58
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
integer, parameter, public default_output_unit
Definition machine.F:46
Definition of mathematical constants and functions.
real(kind=dp), parameter, public degree
real(kind=dp), parameter, public sqrt3
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
pure real(kind=dp) function, public angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition mathlib.F:183
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition mathlib.F:523
Interface to the message passing library MPI.
Utilities for string manipulations.
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
type of a logger, at the moment it contains just a print level starting at which level it should be l...
represent a keyword in the input
represent a section of the input file
stores all the informations relevant to an mpi environment