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 : USE machine, ONLY: default_output_unit
51 : USE mathconstants, ONLY: degree,&
52 : sqrt3
53 : USE mathlib, ONLY: angle,&
54 : det_3x3,&
55 : inv_3x3
56 : USE message_passing, ONLY: mp_para_env_type
57 : USE string_utilities, ONLY: uppercase
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, &
67 : cell_finalize_canonical_input, &
68 : init_cell, &
69 : read_cell, &
70 : read_cell_cif, &
71 : read_cell_cp2k, &
72 : read_cell_extxyz, &
73 : read_cell_pdb, &
74 : read_cell_xsc, &
75 : set_cell_param, &
76 : write_cell, &
77 : write_cell_low
78 :
79 : CONTAINS
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 79608 : 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 0 : CPASSERT(.NOT. ASSOCIATED(cell))
100 5492952 : ALLOCATE (cell)
101 79608 : cell%ref_count = 1
102 79608 : IF (PRESENT(periodic)) THEN
103 1372 : cell%perd = periodic
104 : ELSE
105 317060 : cell%perd = 1
106 : END IF
107 79608 : cell%orthorhombic = .FALSE.
108 79608 : cell%input_cell_canonicalized = .FALSE.
109 1034904 : cell%input_hmat(:, :) = 0.0_dp
110 1034904 : cell%input_to_canonical(:, :) = 0.0_dp
111 1034904 : cell%input_recip_to_canonical(:, :) = 0.0_dp
112 79608 : cell%symmetry_id = cell_sym_none
113 79608 : IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
114 79608 : IF (PRESENT(tag)) cell%tag = tag
115 :
116 79608 : 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 106 : 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 106 : CPASSERT(ASSOCIATED(cell))
134 :
135 1378 : IF (MAXVAL(ABS(hmat_canonical - hmat_input)) <= eps_hmat) THEN
136 14 : cell%input_cell_canonicalized = .FALSE.
137 182 : cell%input_hmat(:, :) = 0.0_dp
138 182 : cell%input_to_canonical(:, :) = 0.0_dp
139 182 : cell%input_recip_to_canonical(:, :) = 0.0_dp
140 : ELSE
141 6164 : tmat = MATMUL(hmat_canonical, inv_3x3(hmat_input))
142 92 : cell%input_cell_canonicalized = .TRUE.
143 1196 : cell%input_hmat(:, :) = hmat_input(:, :)
144 1196 : cell%input_to_canonical(:, :) = tmat(:, :)
145 1196 : cell%input_recip_to_canonical(:, :) = TRANSPOSE(inv_3x3(tmat))
146 : END IF
147 :
148 106 : 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 106 : SUBROUTINE canonicalize_cell_matrix(cell)
155 :
156 : TYPE(cell_type), POINTER :: cell
157 :
158 : REAL(KIND=dp), DIMENSION(3) :: abc, cell_angle
159 :
160 106 : CPASSERT(ASSOCIATED(cell))
161 :
162 106 : CALL get_cell(cell=cell, abc=abc)
163 106 : cell_angle(1) = angle(cell%hmat(:, 2), cell%hmat(:, 3))
164 106 : cell_angle(2) = angle(cell%hmat(:, 1), cell%hmat(:, 3))
165 106 : 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 106 : periodic=cell%perd, do_init_cell=.TRUE.)
169 :
170 106 : 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 340158 : 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 340158 : CPASSERT(ASSOCIATED(cell))
197 :
198 500610 : IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
199 340560 : IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
200 :
201 340158 : cell%deth = ABS(det_3x3(cell%hmat))
202 :
203 340158 : IF (cell%deth < 1.0E-10_dp) THEN
204 0 : CALL write_cell_low(cell, "angstrom", default_output_unit)
205 : CALL cp_abort(__LOCATION__, &
206 : "An invalid set of cell vectors was specified. "// &
207 0 : "The cell volume is too small")
208 : END IF
209 :
210 344241 : SELECT CASE (cell%symmetry_id)
211 : CASE (cell_sym_cubic, &
212 : cell_sym_tetragonal_ab, &
213 : cell_sym_tetragonal_ac, &
214 : cell_sym_tetragonal_bc, &
215 : cell_sym_orthorhombic)
216 4083 : CALL get_cell(cell=cell, abc=abc)
217 4083 : abc(2) = plane_distance(0, 1, 0, cell=cell)
218 4083 : abc(3) = plane_distance(0, 0, 1, cell=cell)
219 4083 : SELECT CASE (cell%symmetry_id)
220 : CASE (cell_sym_cubic)
221 5565 : abc(1:3) = SUM(abc(1:3))/3.0_dp
222 : CASE (cell_sym_tetragonal_ab, &
223 : cell_sym_tetragonal_ac, &
224 : cell_sym_tetragonal_bc)
225 4083 : SELECT CASE (cell%symmetry_id)
226 : CASE (cell_sym_tetragonal_ab)
227 1318 : a = 0.5_dp*(abc(1) + abc(2))
228 1318 : abc(1) = a
229 1318 : abc(2) = a
230 : CASE (cell_sym_tetragonal_ac)
231 631 : a = 0.5_dp*(abc(1) + abc(3))
232 631 : abc(1) = a
233 631 : abc(3) = a
234 : CASE (cell_sym_tetragonal_bc)
235 610 : a = 0.5_dp*(abc(2) + abc(3))
236 610 : abc(2) = a
237 2559 : abc(3) = a
238 : END SELECT
239 : END SELECT
240 4083 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
241 4083 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
242 4083 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
243 : CASE (cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120)
244 2630 : CALL get_cell(cell=cell, abc=abc)
245 2630 : a = 0.5_dp*(abc(1) + abc(2))
246 2630 : acosg = 0.5_dp*a
247 2630 : asing = sqrt3*acosg
248 2630 : IF (cell%symmetry_id == cell_sym_hexagonal_gamma_120) acosg = -acosg
249 2630 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
250 2630 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
251 2630 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
252 : CASE (cell_sym_rhombohedral)
253 649 : CALL get_cell(cell=cell, abc=abc)
254 2596 : 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 649 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
258 649 : acosa = a*COS(alpha)
259 649 : asina = a*SIN(alpha)
260 649 : acosah = a*COS(0.5_dp*alpha)
261 649 : asinah = a*SIN(0.5_dp*alpha)
262 649 : norm = acosa/acosah
263 649 : norm_c = SQRT(1.0_dp - norm*norm)
264 649 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
265 649 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
266 649 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = a*norm_c
267 : CASE (cell_sym_monoclinic)
268 6486 : CALL get_cell(cell=cell, abc=abc)
269 6486 : beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
270 6486 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)
271 6486 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
272 6486 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)
273 : CASE (cell_sym_monoclinic_gamma_ab)
274 : ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
275 738 : CALL get_cell(cell=cell, abc=abc)
276 738 : a = 0.5_dp*(abc(1) + abc(2))
277 738 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
278 738 : acosg = a*COS(gamma)
279 738 : asing = a*SIN(gamma)
280 738 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
281 738 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
282 340896 : 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 340158 : (ABS(cell%hmat(3, 1)) < eps_hmat) .AND. (ABS(cell%hmat(3, 2)) < eps_hmat)) THEN
291 311456 : cell%orthorhombic = .TRUE.
292 : ELSE
293 28702 : 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 340158 : IF (cell%orthorhombic) THEN
299 311456 : cell%hmat(1, 2) = 0.0_dp
300 311456 : cell%hmat(1, 3) = 0.0_dp
301 311456 : cell%hmat(2, 1) = 0.0_dp
302 311456 : cell%hmat(2, 3) = 0.0_dp
303 311456 : cell%hmat(3, 1) = 0.0_dp
304 311456 : cell%hmat(3, 2) = 0.0_dp
305 : END IF
306 :
307 1360632 : dim = COUNT(cell%perd == 1)
308 340158 : IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
309 0 : CPABORT("Non-orthorhombic and not periodic")
310 : END IF
311 :
312 : ! Update deth and hmat_inv with enforced symmetry
313 340158 : cell%deth = ABS(det_3x3(cell%hmat))
314 340158 : 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 0 : "the requested cell symmetry. The cell volume is too small")
318 : END IF
319 4422054 : cell%h_inv = inv_3x3(cell%hmat)
320 :
321 340158 : 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 112400 : 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 11240 : 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 11240 : 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 11240 : my_check_ref = .TRUE.
363 11240 : 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 11240 : CALL cell_create(cell_tmp)
371 11240 : IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
372 11240 : IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
373 11240 : IF (PRESENT(check_for_ref)) my_check_ref = check_for_ref
374 :
375 11240 : cell%deth = 0.0_dp
376 11240 : cell%orthorhombic = .FALSE.
377 44960 : cell%perd(:) = 1
378 11240 : cell%symmetry_id = cell_sym_none
379 146120 : cell%hmat(:, :) = 0.0_dp
380 146120 : cell%h_inv(:, :) = 0.0_dp
381 11240 : cell%input_cell_canonicalized = .FALSE.
382 146120 : cell%input_hmat(:, :) = 0.0_dp
383 146120 : cell%input_to_canonical(:, :) = 0.0_dp
384 146120 : 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 11240 : hmat_input(:, :) = 0.0_dp
392 11240 : read_mat(:, :) = 0.0_dp
393 11240 : read_ang(:) = 0.0_dp
394 11240 : 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 11240 : CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
406 11240 : CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
407 11240 : CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
408 11240 : CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
409 11240 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", explicit=cell_read_alpha_beta_gamma)
410 11240 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
411 11240 : CALL section_vals_val_get(cell_section, "CANONICALIZE", i_val=canonicalize_mode)
412 11240 : canonicalize_cell = (canonicalize_mode == canonicalize_cell_true)
413 :
414 : ! Case 4
415 11240 : tmp_comb_top = (.NOT. (cell_read_file .OR. cell_read_abc))
416 1560 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_a))
417 0 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_b))
418 0 : 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 0 : "cell information.")
425 0 : IF (ASSOCIATED(topology_section)) THEN
426 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", explicit=topo_read_coord)
427 0 : IF (topo_read_coord) THEN
428 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", c_val=coord_file_name)
429 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_FORMAT", i_val=coord_file_format)
430 0 : SELECT CASE (coord_file_format) ! Add formats with both cell and coord parser manually
431 : CASE (do_coord_cif)
432 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
433 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cif)
434 : CASE (do_coord_cp2k)
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_cp2k)
437 : CASE (do_coord_pdb)
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_pdb)
440 : CASE (do_coord_xyz)
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_extxyz)
443 : CASE DEFAULT
444 : CALL cp_abort(__LOCATION__, &
445 : "COORD_FILE_FORMAT is not set to one of the implemented "// &
446 0 : "CELL_FILE_FORMAT options and cannot be parsed for cell information!")
447 : END SELECT
448 : ELSE
449 : CALL cp_abort(__LOCATION__, &
450 0 : "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 0 : "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 11240 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
460 11240 : IF (cell_read_file) THEN ! Case 1
461 16 : 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 0 : "other ones.")
467 16 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
468 16 : CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=cell_file_format)
469 2 : SELECT CASE (cell_file_format)
470 : CASE (do_cell_cp2k)
471 2 : CALL read_cell_cp2k(cell_file_name, cell_tmp, para_env)
472 : CASE (do_cell_xsc)
473 0 : CALL read_cell_xsc(cell_file_name, cell_tmp, para_env)
474 : CASE (do_cell_extxyz)
475 2 : CALL read_cell_extxyz(cell_file_name, cell_tmp, para_env)
476 : CASE (do_cell_pdb)
477 2 : CALL read_cell_pdb(cell_file_name, cell_tmp, para_env)
478 : CASE (do_cell_cif)
479 10 : 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 16 : "options and cannot be parsed for cell information!")
484 : END SELECT
485 208 : read_mat = cell_tmp%hmat
486 : ELSE
487 11224 : IF (cell_read_abc) THEN ! Case 2
488 9664 : CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
489 38656 : read_len = cell_par
490 9664 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_par)
491 38656 : read_ang = cell_par
492 9664 : 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 0 : "The definition of the ABC keyword will override the one provided by A, B and C.")
496 : ELSE ! Case 3
497 1560 : 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 1560 : CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
500 6240 : read_mat(:, 1) = cell_par(:)
501 1560 : CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
502 6240 : read_mat(:, 2) = cell_par(:)
503 1560 : CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
504 6240 : read_mat(:, 3) = cell_par(:)
505 1560 : 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 0 : "keyword ABC.")
509 : ELSE
510 : CALL cp_abort(__LOCATION__, &
511 : "Neither of the keywords CELL_FILE_NAME or ABC are specified, "// &
512 0 : "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 127340 : IF (ANY(read_mat(:, :) > eps)) THEN
519 : ! Make a warning before storing cell vectors that
520 : ! do not form a triangular matrix.
521 1576 : 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 348 : 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 346 : "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 2 : "B is in the XY plane.")
538 : END IF
539 : END IF
540 20488 : cell%hmat = read_mat
541 : ELSE
542 19328 : 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 9664 : do_init_cell=.FALSE.)
545 : ELSE
546 : CALL cp_abort(__LOCATION__, &
547 0 : "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 11240 : CALL reset_cell_section_by_cell_mat(cell, cell_section)
552 :
553 : ! Multiple unit cell
554 11240 : CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
555 44526 : IF (ANY(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
556 :
557 11240 : CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
558 16 : SELECT CASE (my_per)
559 : CASE (use_perd_x)
560 64 : cell%perd = [1, 0, 0]
561 : CASE (use_perd_y)
562 16 : cell%perd = [0, 1, 0]
563 : CASE (use_perd_z)
564 24 : cell%perd = [0, 0, 1]
565 : CASE (use_perd_xy)
566 144 : cell%perd = [1, 1, 0]
567 : CASE (use_perd_xz)
568 32 : cell%perd = [1, 0, 1]
569 : CASE (use_perd_yz)
570 160 : cell%perd = [0, 1, 1]
571 : CASE (use_perd_xyz)
572 28312 : cell%perd = [1, 1, 1]
573 : CASE (use_perd_none)
574 16208 : cell%perd = [0, 0, 0]
575 : CASE DEFAULT
576 11240 : CPABORT("Invalid or not yet implemented cell periodicity")
577 : END SELECT
578 :
579 : ! Load requested cell symmetry
580 11240 : 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 146120 : hmat_input(:, :) = cell%hmat(:, :)
584 11240 : CALL cell_clone(cell, cell_tmp)
585 11240 : CALL init_cell(cell_tmp)
586 145834 : 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 26 : "anyway with the input atomic coordinates."
594 26 : CALL cp_warn(__LOCATION__, error_msg)
595 : END IF
596 11240 : IF (canonicalize_cell) THEN
597 106 : CALL canonicalize_cell_matrix(cell_tmp)
598 106 : CALL cell_finalize_canonical_input(cell_tmp, hmat_input, cell_tmp%hmat)
599 : END IF
600 11240 : CALL cell_clone(cell_tmp, cell)
601 11240 : CALL cell_release(cell_tmp)
602 11240 : CALL reset_cell_section_by_cell_mat(cell, cell_section)
603 :
604 11240 : IF (my_check_ref) THEN
605 : ! Recursive check for reference cell requested
606 10820 : cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
607 10820 : IF (parsed_cp2k_input(cell_ref_section, check_this_section=.TRUE.)) THEN
608 26 : 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 26 : para_env=para_env)
612 : ELSE
613 10794 : CALL cell_clone(cell, cell_ref, tag="CELL_REF")
614 10794 : IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE.
615 : END IF
616 : END IF
617 :
618 11240 : 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 10820 : 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 10820 : my_check = .FALSE.
638 10820 : IF (PRESENT(check_this_section)) my_check = check_this_section
639 10820 : res = ASSOCIATED(input_file)
640 10820 : IF (res) THEN
641 10820 : CPASSERT(input_file%ref_count > 0)
642 10820 : IF (.NOT. my_check) THEN
643 0 : glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
644 0 : CALL section_vals_get(glob_section, explicit=res)
645 : ELSE
646 10820 : CALL section_vals_get(input_file, explicit=res)
647 : END IF
648 : END IF
649 :
650 10820 : 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 9798 : 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 9798 : CPASSERT(ASSOCIATED(cell))
677 39192 : CPASSERT(ALL(cell_angle /= 0.0_dp))
678 :
679 9798 : cos_gamma = COS(cell_angle(3)); IF (ABS(cos_gamma) < eps) cos_gamma = 0.0_dp
680 9798 : IF (ABS(ABS(cos_gamma) - 1.0_dp) < eps) cos_gamma = SIGN(1.0_dp, cos_gamma)
681 9798 : sin_gamma = SIN(cell_angle(3)); IF (ABS(sin_gamma) < eps) sin_gamma = 0.0_dp
682 9798 : IF (ABS(ABS(sin_gamma) - 1.0_dp) < eps) sin_gamma = SIGN(1.0_dp, sin_gamma)
683 9798 : cos_beta = COS(cell_angle(2)); IF (ABS(cos_beta) < eps) cos_beta = 0.0_dp
684 9798 : IF (ABS(ABS(cos_beta) - 1.0_dp) < eps) cos_beta = SIGN(1.0_dp, cos_beta)
685 9798 : cos_alpha = COS(cell_angle(1)); IF (ABS(cos_alpha) < eps) cos_alpha = 0.0_dp
686 9798 : IF (ABS(ABS(cos_alpha) - 1.0_dp) < eps) cos_alpha = SIGN(1.0_dp, cos_alpha)
687 :
688 39192 : cell%hmat(:, 1) = [1.0_dp, 0.0_dp, 0.0_dp]
689 39192 : cell%hmat(:, 2) = [cos_gamma, sin_gamma, 0.0_dp]
690 39192 : cell%hmat(:, 3) = [cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp]
691 9798 : cell%hmat(3, 3) = SQRT(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
692 :
693 39192 : cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
694 39192 : cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
695 39192 : cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
696 :
697 9798 : IF (do_init_cell) THEN
698 134 : IF (PRESENT(periodic)) THEN
699 134 : CALL init_cell(cell=cell, periodic=periodic)
700 : ELSE
701 0 : CALL init_cell(cell=cell)
702 : END IF
703 : END IF
704 :
705 9798 : 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 148 : 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 148 : CPASSERT(ASSOCIATED(cell))
721 :
722 : ! Abort, if one of the value is set to zero
723 592 : IF (ANY(multiple_unit_cell <= 0)) &
724 : CALL cp_abort(__LOCATION__, &
725 : "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
726 0 : "A value of 0 or negative is meaningless!")
727 :
728 : ! Scale abc according to user request
729 592 : cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
730 592 : cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
731 592 : cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
732 :
733 148 : 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 48 : 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 24 : CALL timeset(routineN, handle)
767 :
768 : CALL parser_create(parser, cif_file_name, &
769 24 : para_env=para_env, apply_preprocessing=.FALSE.)
770 :
771 : ! Parsing cell infos
772 96 : 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 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
776 24 : IF (.NOT. found) THEN
777 : CALL parser_search_string(parser, "_cell.length_a", ignore_case=.FALSE., found=found, &
778 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
779 0 : IF (.NOT. found) &
780 0 : CPABORT("The field _cell_length_a or _cell.length_a was not found in CIF file! ")
781 : END IF
782 24 : CALL cif_get_real(parser, cell_lengths(1))
783 24 : 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 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
788 24 : IF (.NOT. found) THEN
789 : CALL parser_search_string(parser, "_cell.length_b", ignore_case=.FALSE., found=found, &
790 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
791 0 : IF (.NOT. found) &
792 0 : CPABORT("The field _cell_length_b or _cell.length_b was not found in CIF file! ")
793 : END IF
794 24 : CALL cif_get_real(parser, cell_lengths(2))
795 24 : 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 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
800 24 : IF (.NOT. found) THEN
801 : CALL parser_search_string(parser, "_cell.length_c", ignore_case=.FALSE., found=found, &
802 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
803 0 : IF (.NOT. found) &
804 0 : CPABORT("The field _cell_length_c or _cell.length_c was not found in CIF file! ")
805 : END IF
806 24 : CALL cif_get_real(parser, cell_lengths(3))
807 24 : 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 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
812 24 : IF (.NOT. found) THEN
813 : CALL parser_search_string(parser, "_cell.angle_alpha", ignore_case=.FALSE., found=found, &
814 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
815 0 : IF (.NOT. found) &
816 0 : CPABORT("The field _cell_angle_alpha or _cell.angle_alpha was not found in CIF file! ")
817 : END IF
818 24 : CALL cif_get_real(parser, cell_angles(1))
819 24 : 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 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
824 24 : IF (.NOT. found) THEN
825 : CALL parser_search_string(parser, "_cell.angle_beta", ignore_case=.FALSE., found=found, &
826 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
827 0 : IF (.NOT. found) &
828 0 : CPABORT("The field _cell_angle_beta or _cell.angle_beta was not found in CIF file! ")
829 : END IF
830 24 : CALL cif_get_real(parser, cell_angles(2))
831 24 : 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 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
836 24 : IF (.NOT. found) THEN
837 : CALL parser_search_string(parser, "_cell.angle_gamma", ignore_case=.FALSE., found=found, &
838 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
839 0 : IF (.NOT. found) &
840 0 : CPABORT("The field _cell_angle_gamma or _cell.angle_gamma was not found in CIF file! ")
841 : END IF
842 24 : CALL cif_get_real(parser, cell_angles(3))
843 24 : 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 24 : do_init_cell=.TRUE.)
848 :
849 24 : CALL parser_release(parser)
850 :
851 24 : CALL timestop(handle)
852 :
853 72 : 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 144 : 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 144 : CALL parser_get_object(parser, s_tag)
872 144 : iln = LEN_TRIM(s_tag)
873 144 : IF (INDEX(s_tag, "(") /= 0) iln = INDEX(s_tag, "(") - 1
874 144 : READ (s_tag(1:iln), *) r
875 :
876 144 : 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 4 : 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 2 : CALL timeset(routineN, handle)
904 :
905 : CALL parser_create(parser, extxyz_file_name, &
906 2 : para_env=para_env, apply_preprocessing=.FALSE.)
907 2 : CALL parser_get_next_line(parser, 2)
908 2 : CALL uppercase(parser%input_line)
909 2 : id1 = INDEX(parser%input_line, "LATTICE=")
910 2 : IF (id1 > 0) THEN
911 2 : id2 = INDEX(parser%input_line(id1 + 9:), '"')
912 2 : READ (parser%input_line(id1 + 9:id1 + id2 + 7), '(A)') raw_cell_str
913 2 : READ (raw_cell_str, *, IOSTAT=ios) hmat(:, 1), hmat(:, 2), hmat(:, 3)
914 2 : IF (ios /= 0) THEN
915 : CALL cp_abort(__LOCATION__, "Error while parsing extended XYZ file "// &
916 : "<"//TRIM(extxyz_file_name)//"> for cell vectors: "// &
917 0 : "found <lattice=> field as <"//TRIM(raw_cell_str)//">")
918 : END IF
919 8 : DO i = 1, 3
920 26 : DO j = 1, 3
921 24 : 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 0 : "extended XYZ specification! ")
930 : END IF
931 :
932 2 : CALL parser_release(parser)
933 :
934 2 : CALL timestop(handle)
935 :
936 6 : 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 4 : 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 2 : CALL timeset(routineN, handle)
972 :
973 : CALL parser_create(parser, pdb_file_name, &
974 2 : para_env=para_env, apply_preprocessing=.FALSE.)
975 :
976 : CALL parser_search_string(parser, "CRYST1", ignore_case=.FALSE., found=found, &
977 2 : begin_line=.TRUE., search_from_begin_of_file=.TRUE.)
978 2 : IF (.NOT. found) &
979 0 : CPABORT("The line <CRYST1> was not found in PDB file! ")
980 :
981 8 : periodic = 1
982 2 : READ (parser%input_line, *, IOSTAT=ios) cryst, cell_lengths(:), cell_angles(:)
983 2 : 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 0 : "found CRYST1 line as <"//TRIM(parser%input_line)//">")
987 : END IF
988 8 : DO i = 1, 3
989 6 : cell_lengths(i) = cp_unit_to_cp2k(cell_lengths(i), "angstrom")
990 8 : 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 2 : do_init_cell=.TRUE.)
994 :
995 2 : CALL parser_release(parser)
996 :
997 2 : CALL timestop(handle)
998 :
999 6 : 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 4 : 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 2 : CALL timeset(routineN, handle)
1024 :
1025 : CALL parser_create(parser, cp2k_file_name, &
1026 2 : para_env=para_env, apply_preprocessing=.FALSE.)
1027 :
1028 2 : CALL parser_get_next_line(parser, 1)
1029 2 : my_end = .FALSE.
1030 24 : DO WHILE (.NOT. my_end)
1031 22 : READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
1032 24 : CALL parser_get_next_line(parser, 1, at_end=my_end)
1033 : END DO
1034 8 : DO i = 1, 3
1035 26 : DO j = 1, 3
1036 24 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
1037 : END DO
1038 : END DO
1039 :
1040 2 : CALL parser_release(parser)
1041 :
1042 2 : CALL timestop(handle)
1043 :
1044 6 : 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 0 : 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 0 : CALL timeset(routineN, handle)
1067 :
1068 : CALL parser_create(parser, xsc_file_name, &
1069 0 : para_env=para_env, apply_preprocessing=.FALSE.)
1070 :
1071 0 : CALL parser_get_next_line(parser, 1)
1072 0 : READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
1073 0 : DO i = 1, 3
1074 0 : DO j = 1, 3
1075 0 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
1076 : END DO
1077 : END DO
1078 :
1079 0 : CALL parser_release(parser)
1080 :
1081 0 : CALL timestop(handle)
1082 :
1083 0 : 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 22480 : 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 22480 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
1099 :
1100 22480 : CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
1101 22480 : CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
1102 22480 : CALL section_vals_val_unset(cell_section, "ABC")
1103 22480 : CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
1104 22480 : CALL section_vals_val_unset(cell_section, "A")
1105 22480 : CALL section_vals_val_unset(cell_section, "B")
1106 22480 : CALL section_vals_val_unset(cell_section, "C")
1107 22480 : ALLOCATE (cell_par(3))
1108 179840 : cell_par = cell%hmat(:, 1)
1109 22480 : CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
1110 22480 : ALLOCATE (cell_par(3))
1111 179840 : cell_par = cell%hmat(:, 2)
1112 22480 : CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
1113 22480 : ALLOCATE (cell_par(3))
1114 179840 : cell_par = cell%hmat(:, 3)
1115 22480 : CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
1116 :
1117 22480 : 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 30957 : 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 30957 : NULLIFY (logger)
1141 30957 : logger => cp_get_default_logger()
1142 30957 : IF (PRESENT(tag)) THEN
1143 22708 : label = TRIM(tag)//"|"
1144 : ELSE
1145 8249 : label = TRIM(cell%tag)//"|"
1146 : END IF
1147 :
1148 30957 : output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", extension=".Log")
1149 30957 : CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
1150 30957 : CALL write_cell_low(cell, unit_str, output_unit, label)
1151 30957 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, "PRINT%CELL")
1152 :
1153 30957 : 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 30957 : 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 30957 : NULLIFY (enum)
1183 30957 : NULLIFY (keyword)
1184 30957 : NULLIFY (section)
1185 :
1186 41083 : IF (output_unit > 0) THEN
1187 10126 : CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma, tag=tag)
1188 10126 : IF (PRESENT(label)) THEN
1189 10126 : my_label = label
1190 : ELSE
1191 0 : my_label = TRIM(tag)//"|"
1192 : END IF
1193 10126 : val = cp_unit_from_cp2k(cell%deth, TRIM(unit_str)//"^3")
1194 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T61,F20.6)") &
1195 10126 : TRIM(my_label)//" Volume ["//TRIM(unit_str)//"^3]:", val
1196 10126 : 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 40504 : TRIM(my_label)//" Vector a ["//TRIM(unit_str)//"]:", cell%hmat(:, 1)*val, &
1199 10126 : "|a| = ", abc(1)*val, &
1200 40504 : TRIM(my_label)//" Vector b ["//TRIM(unit_str)//"]:", cell%hmat(:, 2)*val, &
1201 10126 : "|b| = ", abc(2)*val, &
1202 40504 : TRIM(my_label)//" Vector c ["//TRIM(unit_str)//"]:", cell%hmat(:, 3)*val, &
1203 20252 : "|c| = ", abc(3)*val
1204 : WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.6)") &
1205 10126 : TRIM(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
1206 10126 : TRIM(my_label)//" Angle (a,c), beta [degree]: ", beta, &
1207 20252 : TRIM(my_label)//" Angle (a,b), gamma [degree]: ", gamma
1208 10126 : IF (cell%symmetry_id /= cell_sym_none) THEN
1209 2208 : CALL create_cell_section(section)
1210 2208 : keyword => section_get_keyword(section, "SYMMETRY")
1211 2208 : CALL keyword_get(keyword, enum=enum)
1212 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
1213 2208 : TRIM(my_label)//" Requested initial symmetry: ", &
1214 4416 : ADJUSTR(TRIM(enum_i2c(enum, cell%symmetry_id)))
1215 2208 : CALL section_release(section)
1216 : END IF
1217 10126 : IF (cell%orthorhombic) THEN
1218 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1219 6284 : TRIM(my_label)//" Numerically orthorhombic: ", "YES"
1220 : ELSE
1221 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1222 3842 : TRIM(my_label)//" Numerically orthorhombic: ", " NO"
1223 : END IF
1224 40504 : IF (SUM(cell%perd(1:3)) == 0) THEN
1225 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
1226 1914 : TRIM(my_label)//" Periodicity", "NONE"
1227 : ELSE
1228 8212 : string = ""
1229 8212 : IF (cell%perd(1) == 1) string = TRIM(string)//"X"
1230 8212 : IF (cell%perd(2) == 1) string = TRIM(string)//"Y"
1231 8212 : IF (cell%perd(3) == 1) string = TRIM(string)//"Z"
1232 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1233 8212 : TRIM(my_label)//" Periodicity", ADJUSTR(string)
1234 : END IF
1235 : END IF
1236 :
1237 30957 : END SUBROUTINE write_cell_low
1238 :
1239 : END MODULE cell_methods
|
