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 : do_cell_cif, do_cell_cp2k, do_cell_extxyz, do_cell_pdb, do_cell_xsc, do_coord_cif, &
37 : 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 : init_cell, &
68 : read_cell, &
69 : read_cell_cif, &
70 : read_cell_cp2k, &
71 : read_cell_extxyz, &
72 : read_cell_pdb, &
73 : read_cell_xsc, &
74 : set_cell_param, &
75 : write_cell, &
76 : write_cell_low
77 :
78 : CONTAINS
79 :
80 : ! **************************************************************************************************
81 : !> \brief allocates and initializes a cell
82 : !> \param cell the cell to initialize
83 : !> \param hmat the h matrix that defines the cell
84 : !> \param periodic periodicity of the cell
85 : !> \param tag ...
86 : !> \par History
87 : !> 09.2003 created [fawzi]
88 : !> \author Fawzi Mohamed
89 : ! **************************************************************************************************
90 82792 : SUBROUTINE cell_create(cell, hmat, periodic, tag)
91 :
92 : TYPE(cell_type), POINTER :: cell
93 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
94 : OPTIONAL :: hmat
95 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
96 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
97 :
98 0 : CPASSERT(.NOT. ASSOCIATED(cell))
99 2483760 : ALLOCATE (cell)
100 82792 : cell%ref_count = 1
101 82792 : IF (PRESENT(periodic)) THEN
102 1340 : cell%perd = periodic
103 : ELSE
104 329828 : cell%perd = 1
105 : END IF
106 : cell%orthorhombic = .FALSE.
107 : cell%symmetry_id = cell_sym_none
108 82792 : IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
109 82792 : IF (PRESENT(tag)) cell%tag = tag
110 :
111 82792 : END SUBROUTINE cell_create
112 :
113 : ! **************************************************************************************************
114 : !> \brief Initialise/readjust a simulation cell after hmat has been changed
115 : !> \param cell ...
116 : !> \param hmat ...
117 : !> \param periodic ...
118 : !> \date 16.01.2002
119 : !> \author Matthias Krack
120 : !> \version 1.0
121 : ! **************************************************************************************************
122 346134 : SUBROUTINE init_cell(cell, hmat, periodic)
123 :
124 : TYPE(cell_type), POINTER :: cell
125 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
126 : OPTIONAL :: hmat
127 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
128 :
129 : REAL(KIND=dp), PARAMETER :: eps_hmat = 1.0E-14_dp
130 :
131 : INTEGER :: dim
132 : REAL(KIND=dp) :: a, acosa, acosah, acosg, alpha, asina, &
133 : asinah, asing, beta, gamma, norm, &
134 : norm_c
135 : REAL(KIND=dp), DIMENSION(3) :: abc
136 :
137 346134 : CPASSERT(ASSOCIATED(cell))
138 :
139 450786 : IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
140 346218 : IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
141 :
142 346134 : cell%deth = ABS(det_3x3(cell%hmat))
143 :
144 346134 : IF (cell%deth < 1.0E-10_dp) THEN
145 0 : CALL write_cell_low(cell, "angstrom", default_output_unit)
146 : CALL cp_abort(__LOCATION__, &
147 : "An invalid set of cell vectors was specified. "// &
148 0 : "The cell volume is too small")
149 : END IF
150 :
151 350200 : SELECT CASE (cell%symmetry_id)
152 : CASE (cell_sym_cubic, &
153 : cell_sym_tetragonal_ab, &
154 : cell_sym_tetragonal_ac, &
155 : cell_sym_tetragonal_bc, &
156 : cell_sym_orthorhombic)
157 4066 : CALL get_cell(cell=cell, abc=abc)
158 4066 : abc(2) = plane_distance(0, 1, 0, cell=cell)
159 4066 : abc(3) = plane_distance(0, 0, 1, cell=cell)
160 4066 : SELECT CASE (cell%symmetry_id)
161 : CASE (cell_sym_cubic)
162 5376 : abc(1:3) = SUM(abc(1:3))/3.0_dp
163 : CASE (cell_sym_tetragonal_ab, &
164 : cell_sym_tetragonal_ac, &
165 : cell_sym_tetragonal_bc)
166 4066 : SELECT CASE (cell%symmetry_id)
167 : CASE (cell_sym_tetragonal_ab)
168 1322 : a = 0.5_dp*(abc(1) + abc(2))
169 1322 : abc(1) = a
170 1322 : abc(2) = a
171 : CASE (cell_sym_tetragonal_ac)
172 633 : a = 0.5_dp*(abc(1) + abc(3))
173 633 : abc(1) = a
174 633 : abc(3) = a
175 : CASE (cell_sym_tetragonal_bc)
176 612 : a = 0.5_dp*(abc(2) + abc(3))
177 612 : abc(2) = a
178 2567 : abc(3) = a
179 : END SELECT
180 : END SELECT
181 4066 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
182 4066 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
183 4066 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
184 : CASE (cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120)
185 2638 : CALL get_cell(cell=cell, abc=abc)
186 2638 : a = 0.5_dp*(abc(1) + abc(2))
187 2638 : acosg = 0.5_dp*a
188 2638 : asing = sqrt3*acosg
189 2638 : IF (cell%symmetry_id == cell_sym_hexagonal_gamma_120) acosg = -acosg
190 2638 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
191 2638 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
192 2638 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
193 : CASE (cell_sym_rhombohedral)
194 665 : CALL get_cell(cell=cell, abc=abc)
195 2660 : a = SUM(abc(1:3))/3.0_dp
196 : alpha = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
197 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
198 665 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
199 665 : acosa = a*COS(alpha)
200 665 : asina = a*SIN(alpha)
201 665 : acosah = a*COS(0.5_dp*alpha)
202 665 : asinah = a*SIN(0.5_dp*alpha)
203 665 : norm = acosa/acosah
204 665 : norm_c = SQRT(1.0_dp - norm*norm)
205 665 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
206 665 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
207 665 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = a*norm_c
208 : CASE (cell_sym_monoclinic)
209 873 : CALL get_cell(cell=cell, abc=abc)
210 873 : beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
211 873 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)
212 873 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
213 873 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)
214 : CASE (cell_sym_monoclinic_gamma_ab)
215 : ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
216 742 : CALL get_cell(cell=cell, abc=abc)
217 742 : a = 0.5_dp*(abc(1) + abc(2))
218 742 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
219 742 : acosg = a*COS(gamma)
220 742 : asing = a*SIN(gamma)
221 742 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
222 742 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
223 346876 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
224 : CASE (cell_sym_triclinic)
225 : ! Nothing to do
226 : END SELECT
227 :
228 : ! Do we have an (almost) orthorhombic cell?
229 : IF ((ABS(cell%hmat(1, 2)) < eps_hmat) .AND. (ABS(cell%hmat(1, 3)) < eps_hmat) .AND. &
230 : (ABS(cell%hmat(2, 1)) < eps_hmat) .AND. (ABS(cell%hmat(2, 3)) < eps_hmat) .AND. &
231 346134 : (ABS(cell%hmat(3, 1)) < eps_hmat) .AND. (ABS(cell%hmat(3, 2)) < eps_hmat)) THEN
232 320471 : cell%orthorhombic = .TRUE.
233 : ELSE
234 25663 : cell%orthorhombic = .FALSE.
235 : END IF
236 :
237 : ! Retain an exact orthorhombic cell
238 : ! (off-diagonal elements must remain zero identically to keep QS fast)
239 346134 : IF (cell%orthorhombic) THEN
240 320471 : cell%hmat(1, 2) = 0.0_dp
241 320471 : cell%hmat(1, 3) = 0.0_dp
242 320471 : cell%hmat(2, 1) = 0.0_dp
243 320471 : cell%hmat(2, 3) = 0.0_dp
244 320471 : cell%hmat(3, 1) = 0.0_dp
245 320471 : cell%hmat(3, 2) = 0.0_dp
246 : END IF
247 :
248 1384536 : dim = COUNT(cell%perd == 1)
249 346134 : IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
250 0 : CPABORT("Non-orthorhombic and not periodic")
251 : END IF
252 :
253 : ! Update deth and hmat_inv with enforced symmetry
254 346134 : cell%deth = ABS(det_3x3(cell%hmat))
255 346134 : IF (cell%deth < 1.0E-10_dp) THEN
256 : CALL cp_abort(__LOCATION__, &
257 : "An invalid set of cell vectors was obtained after applying "// &
258 0 : "the requested cell symmetry. The cell volume is too small")
259 : END IF
260 4499742 : cell%h_inv = inv_3x3(cell%hmat)
261 :
262 346134 : END SUBROUTINE init_cell
263 :
264 : ! **************************************************************************************************
265 : !> \brief ...
266 : !> \param cell ...
267 : !> \param cell_ref ...
268 : !> \param use_ref_cell ...
269 : !> \param cell_section ...
270 : !> \param topology_section ...
271 : !> \param check_for_ref ...
272 : !> \param para_env ...
273 : !> \par History
274 : !> 03.2005 created [teo]
275 : !> 03.2026 revamped logic with pdb and extxyz parsers
276 : !> \author Teodoro Laino
277 : ! **************************************************************************************************
278 185697 : RECURSIVE SUBROUTINE read_cell(cell, cell_ref, use_ref_cell, cell_section, &
279 : topology_section, check_for_ref, para_env)
280 :
281 : TYPE(cell_type), POINTER :: cell, cell_ref
282 : LOGICAL, INTENT(INOUT), OPTIONAL :: use_ref_cell
283 : TYPE(section_vals_type), OPTIONAL, POINTER :: cell_section, topology_section
284 : LOGICAL, INTENT(IN), OPTIONAL :: check_for_ref
285 : TYPE(mp_para_env_type), POINTER :: para_env
286 :
287 : REAL(KIND=dp), PARAMETER :: eps = 1.0E-14_dp
288 :
289 : CHARACTER(LEN=default_path_length) :: cell_file_name, coord_file_name, &
290 : error_msg
291 : INTEGER :: cell_file_format, coord_file_format, &
292 : my_per
293 20633 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
294 : LOGICAL :: cell_read_a, cell_read_abc, cell_read_alpha_beta_gamma, cell_read_b, cell_read_c, &
295 : cell_read_file, my_check_ref, tmp_comb_abc, tmp_comb_cell, tmp_comb_top, topo_read_coord
296 : REAL(KIND=dp), DIMENSION(3) :: read_ang, read_len
297 : REAL(KIND=dp), DIMENSION(3, 3) :: read_mat
298 20633 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
299 : TYPE(cell_type), POINTER :: cell_tmp
300 : TYPE(section_vals_type), POINTER :: cell_ref_section
301 :
302 20633 : my_check_ref = .TRUE.
303 20633 : NULLIFY (cell_ref_section, cell_par, cell_tmp, multiple_unit_cell)
304 : ! cell_tmp has two purposes:
305 : ! 1. for transferring matrix of cell vectors from individual
306 : ! file parser subroutines to read_mat here, assuming that
307 : ! unit conversion has been done in those subroutines;
308 : ! 2. for testing whether enforcing symmetry makes a new set
309 : ! of cell vectors significantly different from parsed input
310 20633 : CALL cell_create(cell_tmp)
311 20633 : IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
312 20633 : IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
313 20633 : IF (PRESENT(check_for_ref)) my_check_ref = check_for_ref
314 :
315 20633 : cell%deth = 0.0_dp
316 20633 : cell%orthorhombic = .FALSE.
317 82532 : cell%perd(:) = 1
318 20633 : cell%symmetry_id = cell_sym_none
319 268229 : cell%hmat(:, :) = 0.0_dp
320 268229 : cell%h_inv(:, :) = 0.0_dp
321 : cell_read_file = .FALSE.
322 : cell_read_a = .FALSE.
323 : cell_read_b = .FALSE.
324 : cell_read_c = .FALSE.
325 : cell_read_abc = .FALSE.
326 : cell_read_alpha_beta_gamma = .FALSE.
327 20633 : read_mat(:, :) = 0.0_dp
328 20633 : read_ang(:) = 0.0_dp
329 20633 : read_len(:) = 0.0_dp
330 :
331 : ! Precedence of retrieving cell information from input:
332 : ! 1. CELL/CELL_FILE_NAME
333 : ! 2. CELL/ABC and optionally CELL/ALPHA_BETA_GAMMA
334 : ! 3. CELL/A, CELL/B, CELL/C
335 : ! 4. TOPOLOGY/COORD_FILE_NAME, if topology_section is present
336 : ! The actual order of processing is 4 -> 1 -> 2 -> 3, with
337 : ! case 4 merged to case 1 (if file format permits) first.
338 : ! Store data into either read_mat or read_ang and read_len
339 : ! in CP2K units, which will be converted to cell%hmat and A, B, C.
340 20633 : CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
341 20633 : CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
342 20633 : CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
343 20633 : CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
344 20633 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", explicit=cell_read_alpha_beta_gamma)
345 20633 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
346 :
347 : ! Case 4
348 20633 : tmp_comb_top = (.NOT. (cell_read_file .OR. cell_read_abc))
349 11709 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_a))
350 0 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_b))
351 0 : tmp_comb_top = (tmp_comb_top .AND. (.NOT. cell_read_c))
352 : IF (tmp_comb_top) THEN
353 : CALL cp_warn(__LOCATION__, &
354 : "None of the keywords CELL_FILE_NAME, ABC, or A, B, C "// &
355 : "are specified in CELL section. CP2K will now attempt to read "// &
356 : "TOPOLOGY/COORD_FILE_NAME if its format can be parsed for "// &
357 0 : "cell information.")
358 0 : IF (ASSOCIATED(topology_section)) THEN
359 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", explicit=topo_read_coord)
360 0 : IF (topo_read_coord) THEN
361 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_NAME", c_val=coord_file_name)
362 0 : CALL section_vals_val_get(topology_section, "COORD_FILE_FORMAT", i_val=coord_file_format)
363 0 : SELECT CASE (coord_file_format) ! Add formats with both cell and coord parser manually
364 : CASE (do_coord_cif)
365 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
366 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cif)
367 : CASE (do_coord_cp2k)
368 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
369 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_cp2k)
370 : CASE (do_coord_pdb)
371 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
372 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_pdb)
373 : CASE (do_coord_xyz)
374 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_NAME", c_val=coord_file_name)
375 0 : CALL section_vals_val_set(cell_section, "CELL_FILE_FORMAT", i_val=do_cell_extxyz)
376 : CASE DEFAULT
377 : CALL cp_abort(__LOCATION__, &
378 : "COORD_FILE_FORMAT is not set to one of the implemented "// &
379 0 : "CELL_FILE_FORMAT options and cannot be parsed for cell information!")
380 : END SELECT
381 : ELSE
382 : CALL cp_abort(__LOCATION__, &
383 0 : "COORD_FILE_NAME is not set, so no cell information is available!")
384 : END IF
385 : ELSE
386 : CALL cp_warn(__LOCATION__, &
387 : "TOPOLOGY section is not available, so COORD_FILE_NAME cannot "// &
388 0 : "be parsed for cell information in lieu of missing CELL settings.")
389 : END IF
390 : END IF
391 : ! Former logic in SUBROUTINE read_cell_from_external_file is moved here
392 20633 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
393 20633 : IF (cell_read_file) THEN ! Case 1
394 16 : tmp_comb_cell = (cell_read_abc .OR. (cell_read_a .OR. (cell_read_b .OR. cell_read_c)))
395 : IF (tmp_comb_cell) &
396 : CALL cp_warn(__LOCATION__, &
397 : "Cell Information provided through A, B, C, or ABC in conjunction "// &
398 : "with CELL_FILE_NAME. The definition in external file will override "// &
399 0 : "other ones.")
400 16 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
401 16 : CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=cell_file_format)
402 2 : SELECT CASE (cell_file_format)
403 : CASE (do_cell_cp2k)
404 2 : CALL read_cell_cp2k(cell_file_name, cell_tmp, para_env)
405 : CASE (do_cell_xsc)
406 0 : CALL read_cell_xsc(cell_file_name, cell_tmp, para_env)
407 : CASE (do_cell_extxyz)
408 2 : CALL read_cell_extxyz(cell_file_name, cell_tmp, para_env)
409 : CASE (do_cell_pdb)
410 2 : CALL read_cell_pdb(cell_file_name, cell_tmp, para_env)
411 : CASE (do_cell_cif)
412 10 : CALL read_cell_cif(cell_file_name, cell_tmp, para_env)
413 : CASE DEFAULT
414 : CALL cp_abort(__LOCATION__, &
415 : "CELL_FILE_FORMAT is not set to one of the implemented "// &
416 16 : "options and cannot be parsed for cell information!")
417 : END SELECT
418 208 : read_mat = cell_tmp%hmat
419 : ELSE
420 20617 : IF (cell_read_abc) THEN ! Case 2
421 8908 : CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
422 35632 : read_len = cell_par
423 8908 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_par)
424 35632 : read_ang = cell_par
425 8908 : IF (cell_read_a .OR. cell_read_b .OR. cell_read_c) &
426 : CALL cp_warn(__LOCATION__, &
427 : "Cell information provided through vectors A, B or C in conjunction with ABC. "// &
428 0 : "The definition of the ABC keyword will override the one provided by A, B and C.")
429 : ELSE ! Case 3
430 11709 : tmp_comb_abc = ((cell_read_a .EQV. cell_read_b) .AND. (cell_read_b .EQV. cell_read_c))
431 : IF (tmp_comb_abc) THEN
432 11709 : CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
433 46836 : read_mat(:, 1) = cell_par(:)
434 11709 : CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
435 46836 : read_mat(:, 2) = cell_par(:)
436 11709 : CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
437 46836 : read_mat(:, 3) = cell_par(:)
438 11709 : IF (cell_read_alpha_beta_gamma) &
439 : CALL cp_warn(__LOCATION__, &
440 : "The keyword ALPHA_BETA_GAMMA is ignored because it was used without the "// &
441 0 : "keyword ABC.")
442 : ELSE
443 : CALL cp_abort(__LOCATION__, &
444 : "Neither of the keywords CELL_FILE_NAME or ABC are specified, "// &
445 0 : "and cell vector settings in A, B, C are incomplete!")
446 : END IF
447 : END IF
448 : END IF
449 :
450 : ! Convert read_mat or read_len and read_ang to actual cell%hmat
451 127777 : IF (ANY(read_mat(:, :) > eps)) THEN
452 : ! Make a warning before storing cell vectors that
453 : ! do not form a triangular matrix.
454 : IF ((ABS(read_mat(2, 1)) > eps) .OR. &
455 11725 : (ABS(read_mat(3, 1)) > eps) .OR. &
456 : (ABS(read_mat(3, 2)) > eps)) &
457 : CALL cp_warn(__LOCATION__, &
458 : "Cell vectors are read but cell matrix is not "// &
459 : "a lower triangle, not conforming to the program "// &
460 : "convention that A lies along the X-axis and "// &
461 868 : "B is in the XY plane.")
462 152425 : cell%hmat = read_mat
463 : ELSE
464 17816 : IF (ANY(read_ang(:) > eps) .AND. ANY(read_len(:) > eps)) THEN
465 : CALL set_cell_param(cell, cell_length=read_len, cell_angle=read_ang, &
466 8908 : do_init_cell=.FALSE.)
467 : ELSE
468 : CALL cp_abort(__LOCATION__, &
469 0 : "No meaningful cell information is read from parser!")
470 : END IF
471 : END IF
472 : ! Reset cell section so that only A, B, C are kept
473 20633 : CALL reset_cell_section_by_cell_mat(cell, cell_section)
474 :
475 : ! Multiple unit cell
476 20633 : CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
477 81652 : IF (ANY(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
478 :
479 20633 : CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
480 0 : SELECT CASE (my_per)
481 : CASE (use_perd_x)
482 0 : cell%perd = [1, 0, 0]
483 : CASE (use_perd_y)
484 0 : cell%perd = [0, 1, 0]
485 : CASE (use_perd_z)
486 32 : cell%perd = [0, 0, 1]
487 : CASE (use_perd_xy)
488 112 : cell%perd = [1, 1, 0]
489 : CASE (use_perd_xz)
490 48 : cell%perd = [1, 0, 1]
491 : CASE (use_perd_yz)
492 192 : cell%perd = [0, 1, 1]
493 : CASE (use_perd_xyz)
494 53084 : cell%perd = [1, 1, 1]
495 : CASE (use_perd_none)
496 29064 : cell%perd = [0, 0, 0]
497 : CASE DEFAULT
498 20633 : CPABORT("")
499 : END SELECT
500 :
501 : ! Load requested cell symmetry
502 20633 : CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
503 : ! Try enforcing symmetry by initializing a temporary copy of cell
504 : ! and see if the resulting cell matrix differ significantly
505 20633 : CALL cell_clone(cell, cell_tmp)
506 20633 : CALL init_cell(cell_tmp)
507 267657 : IF (ANY(ABS(cell_tmp%hmat - cell%hmat) > eps)) THEN
508 : WRITE (UNIT=error_msg, FMT="(A)") &
509 : "When initializing cell vectors with requested symmetry, one "// &
510 : "or more elements of the cell matrix has varied significantly. "// &
511 : "The input parameters are either deviating from the symmetry, "// &
512 : "or not conforming to the program convention that cell matrix "// &
513 : "is a lower triangle. The symmetrized cell vectors will be used "// &
514 52 : "anyway with the input atomic coordinates."
515 52 : CALL cp_warn(__LOCATION__, error_msg)
516 : END IF
517 20633 : CALL cell_clone(cell_tmp, cell)
518 20633 : CALL cell_release(cell_tmp)
519 :
520 20633 : IF (my_check_ref) THEN
521 : ! Recursive check for reference cell requested
522 20187 : cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
523 20187 : IF (parsed_cp2k_input(cell_ref_section, check_this_section=.TRUE.)) THEN
524 52 : IF (PRESENT(use_ref_cell)) use_ref_cell = .TRUE.
525 : CALL read_cell(cell_ref, cell_ref, use_ref_cell=use_ref_cell, &
526 : cell_section=cell_ref_section, check_for_ref=.FALSE., &
527 52 : para_env=para_env)
528 : ELSE
529 20135 : CALL cell_clone(cell, cell_ref, tag="CELL_REF")
530 20135 : IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE.
531 : END IF
532 : END IF
533 :
534 20633 : END SUBROUTINE read_cell
535 :
536 : ! **************************************************************************************************
537 : !> \brief utility function to ease the transition to the new input.
538 : !> returns true if the new input was parsed
539 : !> \param input_file the parsed input file
540 : !> \param check_this_section ...
541 : !> \return ...
542 : !> \author fawzi
543 : ! **************************************************************************************************
544 20187 : FUNCTION parsed_cp2k_input(input_file, check_this_section) RESULT(res)
545 :
546 : TYPE(section_vals_type), POINTER :: input_file
547 : LOGICAL, INTENT(IN), OPTIONAL :: check_this_section
548 : LOGICAL :: res
549 :
550 : LOGICAL :: my_check
551 : TYPE(section_vals_type), POINTER :: glob_section
552 :
553 20187 : my_check = .FALSE.
554 20187 : IF (PRESENT(check_this_section)) my_check = check_this_section
555 20187 : res = ASSOCIATED(input_file)
556 20187 : IF (res) THEN
557 20187 : CPASSERT(input_file%ref_count > 0)
558 20187 : IF (.NOT. my_check) THEN
559 0 : glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
560 0 : CALL section_vals_get(glob_section, explicit=res)
561 : ELSE
562 20187 : CALL section_vals_get(input_file, explicit=res)
563 : END IF
564 : END IF
565 :
566 20187 : END FUNCTION parsed_cp2k_input
567 :
568 : ! **************************************************************************************************
569 : !> \brief Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma)
570 : !> using the convention: a parallel to the x axis, b in the x-y plane and
571 : !> and c univoquely determined; gamma is the angle between a and b; beta
572 : !> is the angle between c and a and alpha is the angle between c and b
573 : !> \param cell ...
574 : !> \param cell_length ...
575 : !> \param cell_angle ...
576 : !> \param periodic ...
577 : !> \param do_init_cell ...
578 : !> \date 03.2008
579 : !> \author Teodoro Laino
580 : ! **************************************************************************************************
581 8936 : SUBROUTINE set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
582 :
583 : TYPE(cell_type), POINTER :: cell
584 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: cell_length, cell_angle
585 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
586 : LOGICAL, INTENT(IN) :: do_init_cell
587 :
588 : REAL(KIND=dp), PARAMETER :: eps = EPSILON(0.0_dp)
589 :
590 : REAL(KIND=dp) :: cos_alpha, cos_beta, cos_gamma, sin_gamma
591 :
592 8936 : CPASSERT(ASSOCIATED(cell))
593 35744 : CPASSERT(ALL(cell_angle /= 0.0_dp))
594 :
595 8936 : cos_gamma = COS(cell_angle(3)); IF (ABS(cos_gamma) < eps) cos_gamma = 0.0_dp
596 8936 : IF (ABS(ABS(cos_gamma) - 1.0_dp) < eps) cos_gamma = SIGN(1.0_dp, cos_gamma)
597 8936 : sin_gamma = SIN(cell_angle(3)); IF (ABS(sin_gamma) < eps) sin_gamma = 0.0_dp
598 8936 : IF (ABS(ABS(sin_gamma) - 1.0_dp) < eps) sin_gamma = SIGN(1.0_dp, sin_gamma)
599 8936 : cos_beta = COS(cell_angle(2)); IF (ABS(cos_beta) < eps) cos_beta = 0.0_dp
600 8936 : IF (ABS(ABS(cos_beta) - 1.0_dp) < eps) cos_beta = SIGN(1.0_dp, cos_beta)
601 8936 : cos_alpha = COS(cell_angle(1)); IF (ABS(cos_alpha) < eps) cos_alpha = 0.0_dp
602 8936 : IF (ABS(ABS(cos_alpha) - 1.0_dp) < eps) cos_alpha = SIGN(1.0_dp, cos_alpha)
603 :
604 35744 : cell%hmat(:, 1) = [1.0_dp, 0.0_dp, 0.0_dp]
605 35744 : cell%hmat(:, 2) = [cos_gamma, sin_gamma, 0.0_dp]
606 35744 : cell%hmat(:, 3) = [cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp]
607 8936 : cell%hmat(3, 3) = SQRT(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
608 :
609 35744 : cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
610 35744 : cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
611 35744 : cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
612 :
613 8936 : IF (do_init_cell) THEN
614 28 : IF (PRESENT(periodic)) THEN
615 28 : CALL init_cell(cell=cell, periodic=periodic)
616 : ELSE
617 0 : CALL init_cell(cell=cell)
618 : END IF
619 : END IF
620 :
621 8936 : END SUBROUTINE set_cell_param
622 :
623 : ! **************************************************************************************************
624 : !> \brief Setup of the multiple unit_cell
625 : !> \param cell ...
626 : !> \param multiple_unit_cell ...
627 : !> \date 05.2009
628 : !> \author Teodoro Laino [tlaino]
629 : !> \version 1.0
630 : ! **************************************************************************************************
631 300 : SUBROUTINE set_multiple_unit_cell(cell, multiple_unit_cell)
632 :
633 : TYPE(cell_type), POINTER :: cell
634 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
635 :
636 300 : CPASSERT(ASSOCIATED(cell))
637 :
638 : ! Abort, if one of the value is set to zero
639 1200 : IF (ANY(multiple_unit_cell <= 0)) &
640 : CALL cp_abort(__LOCATION__, &
641 : "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
642 0 : "A value of 0 or negative is meaningless!")
643 :
644 : ! Scale abc according to user request
645 1200 : cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
646 1200 : cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
647 1200 : cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
648 :
649 300 : END SUBROUTINE set_multiple_unit_cell
650 :
651 : ! **************************************************************************************************
652 : !> \brief Reads cell information from CIF file
653 : !> \param cif_file_name ...
654 : !> \param cell ...
655 : !> \param para_env ...
656 : !> \date 12.2008
657 : !> \par Format Information implemented:
658 : !> _cell_length_a (_cell.length_a)
659 : !> _cell_length_b (_cell.length_b)
660 : !> _cell_length_c (_cell.length_c)
661 : !> _cell_angle_alpha (_cell.length_alpha)
662 : !> _cell_angle_beta (_cell.length_beta)
663 : !> _cell_angle_gamma (_cell.length_gamma)
664 : !>
665 : !> \author Teodoro Laino [tlaino]
666 : !> moved from topology_cif (1/2019 JHU)
667 : ! **************************************************************************************************
668 48 : SUBROUTINE read_cell_cif(cif_file_name, cell, para_env)
669 :
670 : CHARACTER(len=*) :: cif_file_name
671 : TYPE(cell_type), POINTER :: cell
672 : TYPE(mp_para_env_type), POINTER :: para_env
673 :
674 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_cif'
675 :
676 : INTEGER :: handle
677 : INTEGER, DIMENSION(3) :: periodic
678 : LOGICAL :: found
679 : REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths
680 : TYPE(cp_parser_type) :: parser
681 :
682 24 : CALL timeset(routineN, handle)
683 :
684 : CALL parser_create(parser, cif_file_name, &
685 24 : para_env=para_env, apply_preprocessing=.FALSE.)
686 :
687 : ! Parsing cell infos
688 96 : periodic = 1
689 : ! Check for _cell_length_a or _cell.length_a
690 : CALL parser_search_string(parser, "_cell_length_a", ignore_case=.FALSE., found=found, &
691 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
692 24 : IF (.NOT. found) THEN
693 : CALL parser_search_string(parser, "_cell.length_a", ignore_case=.FALSE., found=found, &
694 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
695 0 : IF (.NOT. found) &
696 0 : CPABORT("The field _cell_length_a or _cell.length_a was not found in CIF file! ")
697 : END IF
698 24 : CALL cif_get_real(parser, cell_lengths(1))
699 24 : cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
700 :
701 : ! Check for _cell_length_b or _cell.length_b
702 : CALL parser_search_string(parser, "_cell_length_b", ignore_case=.FALSE., found=found, &
703 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
704 24 : IF (.NOT. found) THEN
705 : CALL parser_search_string(parser, "_cell.length_b", ignore_case=.FALSE., found=found, &
706 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
707 0 : IF (.NOT. found) &
708 0 : CPABORT("The field _cell_length_b or _cell.length_b was not found in CIF file! ")
709 : END IF
710 24 : CALL cif_get_real(parser, cell_lengths(2))
711 24 : cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
712 :
713 : ! Check for _cell_length_c or _cell.length_c
714 : CALL parser_search_string(parser, "_cell_length_c", ignore_case=.FALSE., found=found, &
715 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
716 24 : IF (.NOT. found) THEN
717 : CALL parser_search_string(parser, "_cell.length_c", ignore_case=.FALSE., found=found, &
718 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
719 0 : IF (.NOT. found) &
720 0 : CPABORT("The field _cell_length_c or _cell.length_c was not found in CIF file! ")
721 : END IF
722 24 : CALL cif_get_real(parser, cell_lengths(3))
723 24 : cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
724 :
725 : ! Check for _cell_angle_alpha or _cell.angle_alpha
726 : CALL parser_search_string(parser, "_cell_angle_alpha", ignore_case=.FALSE., found=found, &
727 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
728 24 : IF (.NOT. found) THEN
729 : CALL parser_search_string(parser, "_cell.angle_alpha", ignore_case=.FALSE., found=found, &
730 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
731 0 : IF (.NOT. found) &
732 0 : CPABORT("The field _cell_angle_alpha or _cell.angle_alpha was not found in CIF file! ")
733 : END IF
734 24 : CALL cif_get_real(parser, cell_angles(1))
735 24 : cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
736 :
737 : ! Check for _cell_angle_beta or _cell.angle_beta
738 : CALL parser_search_string(parser, "_cell_angle_beta", ignore_case=.FALSE., found=found, &
739 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
740 24 : IF (.NOT. found) THEN
741 : CALL parser_search_string(parser, "_cell.angle_beta", ignore_case=.FALSE., found=found, &
742 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
743 0 : IF (.NOT. found) &
744 0 : CPABORT("The field _cell_angle_beta or _cell.angle_beta was not found in CIF file! ")
745 : END IF
746 24 : CALL cif_get_real(parser, cell_angles(2))
747 24 : cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
748 :
749 : ! Check for _cell_angle_gamma or _cell.angle_gamma
750 : CALL parser_search_string(parser, "_cell_angle_gamma", ignore_case=.FALSE., found=found, &
751 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
752 24 : IF (.NOT. found) THEN
753 : CALL parser_search_string(parser, "_cell.angle_gamma", ignore_case=.FALSE., found=found, &
754 0 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
755 0 : IF (.NOT. found) &
756 0 : CPABORT("The field _cell_angle_gamma or _cell.angle_gamma was not found in CIF file! ")
757 : END IF
758 24 : CALL cif_get_real(parser, cell_angles(3))
759 24 : cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
760 :
761 : ! Create cell
762 : CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
763 24 : do_init_cell=.TRUE.)
764 :
765 24 : CALL parser_release(parser)
766 :
767 24 : CALL timestop(handle)
768 :
769 72 : END SUBROUTINE read_cell_cif
770 :
771 : ! **************************************************************************************************
772 : !> \brief Reads REAL from the CIF file.. This wrapper is needed in order to
773 : !> treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
774 : !> \param parser ...
775 : !> \param r ...
776 : !> \date 12.2008
777 : !> \author Teodoro Laino [tlaino]
778 : ! **************************************************************************************************
779 144 : SUBROUTINE cif_get_real(parser, r)
780 :
781 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
782 : REAL(KIND=dp), INTENT(OUT) :: r
783 :
784 : CHARACTER(LEN=default_string_length) :: s_tag
785 : INTEGER :: iln
786 :
787 144 : CALL parser_get_object(parser, s_tag)
788 144 : iln = LEN_TRIM(s_tag)
789 144 : IF (INDEX(s_tag, "(") /= 0) iln = INDEX(s_tag, "(") - 1
790 144 : READ (s_tag(1:iln), *) r
791 :
792 144 : END SUBROUTINE cif_get_real
793 :
794 : ! **************************************************************************************************
795 : !> \brief Reads cell information from comment line of extended xyz file
796 : !> \param extxyz_file_name ...
797 : !> \param cell ...
798 : !> \param para_env ...
799 : !> \date 03.2026
800 : !> \par Extended xyz format has comment on the second line in the form as
801 : !> Lattice="Ax Ay Az Bx By Bz Cx Cy Cz"
802 : !> where Ax, Ay, Az are three Cartesian components of cell vector A,
803 : !> Bx, By, Bz are components of B, Cx, Cy, Cz are components of C,
804 : !> all in the unit of angstrom.
805 : ! **************************************************************************************************
806 4 : SUBROUTINE read_cell_extxyz(extxyz_file_name, cell, para_env)
807 :
808 : CHARACTER(len=*) :: extxyz_file_name
809 : TYPE(cell_type), POINTER :: cell
810 : TYPE(mp_para_env_type), POINTER :: para_env
811 :
812 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_extxyz'
813 :
814 : INTEGER :: handle, i, id1, id2, j
815 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
816 : TYPE(cp_parser_type) :: parser
817 :
818 2 : CALL timeset(routineN, handle)
819 :
820 : CALL parser_create(parser, extxyz_file_name, &
821 2 : para_env=para_env, apply_preprocessing=.FALSE.)
822 2 : CALL parser_get_next_line(parser, 2)
823 2 : CALL uppercase(parser%input_line)
824 2 : id1 = INDEX(parser%input_line, "LATTICE=")
825 2 : IF (id1 > 0) THEN
826 2 : id2 = INDEX(parser%input_line(id1 + 9:), '"')
827 2 : READ (parser%input_line(id1 + 9:id1 + id2 + 7), *) hmat(:, 1), hmat(:, 2), hmat(:, 3)
828 8 : DO i = 1, 3
829 26 : DO j = 1, 3
830 24 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
831 : END DO
832 : END DO
833 :
834 : ELSE
835 : CALL cp_abort(__LOCATION__, &
836 : "The field (lattice=) was not found on comment line "// &
837 : "of XYZ file, so cell information cannot be set via "// &
838 0 : "extended XYZ specification! ")
839 : END IF
840 :
841 2 : CALL parser_release(parser)
842 :
843 2 : CALL timestop(handle)
844 :
845 6 : END SUBROUTINE read_cell_extxyz
846 :
847 : ! **************************************************************************************************
848 : !> \brief Reads cell information from CRYST1 record of PDB file
849 : !> \param pdb_file_name ...
850 : !> \param cell ...
851 : !> \param para_env ...
852 : !> \date 03.2026
853 : !> \par CRYST1 record may contain space group and Z value at the end,
854 : !> but here only the first entries are read:
855 : !> COLUMNS DATA TYPE FIELD DEFINITION
856 : !> -------------------------------------------------------------
857 : !> 1 - 6 Record name "CRYST1"
858 : !> 7 - 15 Real(9.3) a a (Angstroms).
859 : !> 16 - 24 Real(9.3) b b (Angstroms).
860 : !> 25 - 33 Real(9.3) c c (Angstroms).
861 : !> 34 - 40 Real(7.2) alpha alpha (degrees).
862 : !> 41 - 47 Real(7.2) beta beta (degrees).
863 : !> 48 - 54 Real(7.2) gamma gamma (degrees).
864 : ! **************************************************************************************************
865 4 : SUBROUTINE read_cell_pdb(pdb_file_name, cell, para_env)
866 :
867 : CHARACTER(len=*) :: pdb_file_name
868 : TYPE(cell_type), POINTER :: cell
869 : TYPE(mp_para_env_type), POINTER :: para_env
870 :
871 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_pdb'
872 :
873 : CHARACTER(LEN=default_string_length) :: cryst
874 : INTEGER :: handle, i
875 : INTEGER, DIMENSION(3) :: periodic
876 : LOGICAL :: found
877 : REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths
878 : TYPE(cp_parser_type) :: parser
879 :
880 2 : CALL timeset(routineN, handle)
881 :
882 : CALL parser_create(parser, pdb_file_name, &
883 2 : para_env=para_env, apply_preprocessing=.FALSE.)
884 :
885 : CALL parser_search_string(parser, "CRYST1", ignore_case=.FALSE., found=found, &
886 2 : begin_line=.TRUE., search_from_begin_of_file=.TRUE.)
887 2 : IF (.NOT. found) &
888 0 : CPABORT("The field (CRYST1) was not found in PDB file! ")
889 :
890 8 : periodic = 1
891 2 : READ (parser%input_line, *) cryst, cell_lengths(:), cell_angles(:)
892 8 : DO i = 1, 3
893 6 : cell_lengths(i) = cp_unit_to_cp2k(cell_lengths(i), "angstrom")
894 8 : cell_angles(i) = cp_unit_to_cp2k(cell_angles(i), "deg")
895 : END DO
896 : CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
897 2 : do_init_cell=.TRUE.)
898 :
899 2 : CALL parser_release(parser)
900 :
901 2 : CALL timestop(handle)
902 :
903 6 : END SUBROUTINE read_cell_pdb
904 :
905 : ! **************************************************************************************************
906 : !> \brief Reads cell information from cp2k file
907 : !> \param cp2k_file_name ...
908 : !> \param cell ...
909 : !> \param para_env ...
910 : !> \date 03.2026
911 : !> \par Isolated from former read_cell_from_external_file
912 : ! **************************************************************************************************
913 4 : SUBROUTINE read_cell_cp2k(cp2k_file_name, cell, para_env)
914 :
915 : CHARACTER(len=*) :: cp2k_file_name
916 : TYPE(cell_type), POINTER :: cell
917 : TYPE(mp_para_env_type), POINTER :: para_env
918 :
919 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_cp2k'
920 :
921 : INTEGER :: handle, i, idum, j
922 : LOGICAL :: my_end
923 : REAL(KIND=dp) :: xdum
924 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
925 : TYPE(cp_parser_type) :: parser
926 :
927 2 : CALL timeset(routineN, handle)
928 :
929 : CALL parser_create(parser, cp2k_file_name, &
930 2 : para_env=para_env, apply_preprocessing=.FALSE.)
931 :
932 2 : CALL parser_get_next_line(parser, 1)
933 2 : my_end = .FALSE.
934 24 : DO WHILE (.NOT. my_end)
935 22 : READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
936 24 : CALL parser_get_next_line(parser, 1, at_end=my_end)
937 : END DO
938 8 : DO i = 1, 3
939 26 : DO j = 1, 3
940 24 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
941 : END DO
942 : END DO
943 :
944 2 : CALL parser_release(parser)
945 :
946 2 : CALL timestop(handle)
947 :
948 6 : END SUBROUTINE read_cell_cp2k
949 :
950 : ! **************************************************************************************************
951 : !> \brief Reads cell information from xsc file
952 : !> \param xsc_file_name ...
953 : !> \param cell ...
954 : !> \param para_env ...
955 : !> \date 03.2026
956 : !> \par Isolated from former read_cell_from_external_file
957 : ! **************************************************************************************************
958 0 : SUBROUTINE read_cell_xsc(xsc_file_name, cell, para_env)
959 :
960 : CHARACTER(len=*) :: xsc_file_name
961 : TYPE(cell_type), POINTER :: cell
962 : TYPE(mp_para_env_type), POINTER :: para_env
963 :
964 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_xsc'
965 :
966 : INTEGER :: handle, i, idum, j
967 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
968 : TYPE(cp_parser_type) :: parser
969 :
970 0 : CALL timeset(routineN, handle)
971 :
972 : CALL parser_create(parser, xsc_file_name, &
973 0 : para_env=para_env, apply_preprocessing=.FALSE.)
974 :
975 0 : CALL parser_get_next_line(parser, 1)
976 0 : READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
977 0 : DO i = 1, 3
978 0 : DO j = 1, 3
979 0 : cell%hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
980 : END DO
981 : END DO
982 :
983 0 : CALL parser_release(parser)
984 :
985 0 : CALL timestop(handle)
986 :
987 0 : END SUBROUTINE read_cell_xsc
988 :
989 : ! **************************************************************************************************
990 : !> \brief Reset cell section by matrix in cell-type pointer
991 : !> \param cell ...
992 : !> \param cell_section ...
993 : !> \date 03.2026
994 : !> \par Alternative keywords for cell settings will be unset
995 : !> except MULTIPLE_UNIT_CELL, PERIODIC and SYMMETRY.
996 : ! **************************************************************************************************
997 20633 : SUBROUTINE reset_cell_section_by_cell_mat(cell, cell_section)
998 :
999 : TYPE(cell_type), POINTER :: cell
1000 : TYPE(section_vals_type), POINTER :: cell_section
1001 :
1002 20633 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
1003 :
1004 20633 : CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
1005 20633 : CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
1006 20633 : CALL section_vals_val_unset(cell_section, "ABC")
1007 20633 : CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
1008 20633 : CALL section_vals_val_unset(cell_section, "A")
1009 20633 : CALL section_vals_val_unset(cell_section, "B")
1010 20633 : CALL section_vals_val_unset(cell_section, "C")
1011 20633 : ALLOCATE (cell_par(3))
1012 165064 : cell_par = cell%hmat(:, 1)
1013 20633 : CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
1014 20633 : ALLOCATE (cell_par(3))
1015 165064 : cell_par = cell%hmat(:, 2)
1016 20633 : CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
1017 20633 : ALLOCATE (cell_par(3))
1018 165064 : cell_par = cell%hmat(:, 3)
1019 20633 : CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
1020 :
1021 20633 : END SUBROUTINE reset_cell_section_by_cell_mat
1022 :
1023 : ! **************************************************************************************************
1024 : !> \brief Write the cell parameters to the output unit.
1025 : !> \param cell ...
1026 : !> \param subsys_section ...
1027 : !> \param tag ...
1028 : !> \date 02.06.2000
1029 : !> \par History
1030 : !> - 11.2008 Teodoro Laino [tlaino] - rewrite and enabling user driven units
1031 : !> \author Matthias Krack
1032 : !> \version 1.0
1033 : ! **************************************************************************************************
1034 37561 : SUBROUTINE write_cell(cell, subsys_section, tag)
1035 :
1036 : TYPE(cell_type), POINTER :: cell
1037 : TYPE(section_vals_type), POINTER :: subsys_section
1038 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
1039 :
1040 : CHARACTER(LEN=default_string_length) :: label, unit_str
1041 : INTEGER :: output_unit
1042 : TYPE(cp_logger_type), POINTER :: logger
1043 :
1044 37561 : NULLIFY (logger)
1045 37561 : logger => cp_get_default_logger()
1046 37561 : IF (PRESENT(tag)) THEN
1047 10012 : label = TRIM(tag)//"|"
1048 : ELSE
1049 27549 : label = TRIM(cell%tag)//"|"
1050 : END IF
1051 :
1052 37561 : output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", extension=".Log")
1053 37561 : CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
1054 37561 : CALL write_cell_low(cell, unit_str, output_unit, label)
1055 37561 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, "PRINT%CELL")
1056 :
1057 37561 : END SUBROUTINE write_cell
1058 :
1059 : ! **************************************************************************************************
1060 : !> \brief Write the cell parameters to the output unit
1061 : !> \param cell ...
1062 : !> \param unit_str ...
1063 : !> \param output_unit ...
1064 : !> \param label ...
1065 : !> \date 17.05.2023
1066 : !> \par History
1067 : !> - Extracted from write_cell (17.05.2023, MK)
1068 : !> \version 1.0
1069 : ! **************************************************************************************************
1070 37561 : SUBROUTINE write_cell_low(cell, unit_str, output_unit, label)
1071 :
1072 : TYPE(cell_type), POINTER :: cell
1073 : CHARACTER(LEN=*), INTENT(IN) :: unit_str
1074 : INTEGER, INTENT(IN) :: output_unit
1075 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: label
1076 :
1077 : CHARACTER(LEN=12) :: tag
1078 : CHARACTER(LEN=3) :: string
1079 : CHARACTER(LEN=default_string_length) :: my_label
1080 : REAL(KIND=dp) :: alpha, beta, gamma, val
1081 : REAL(KIND=dp), DIMENSION(3) :: abc
1082 : TYPE(enumeration_type), POINTER :: enum
1083 : TYPE(keyword_type), POINTER :: keyword
1084 : TYPE(section_type), POINTER :: section
1085 :
1086 37561 : NULLIFY (enum)
1087 37561 : NULLIFY (keyword)
1088 37561 : NULLIFY (section)
1089 :
1090 49107 : IF (output_unit > 0) THEN
1091 11546 : CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma, tag=tag)
1092 11546 : IF (PRESENT(label)) THEN
1093 11546 : my_label = label
1094 : ELSE
1095 0 : my_label = TRIM(tag)//"|"
1096 : END IF
1097 11546 : val = cp_unit_from_cp2k(cell%deth, TRIM(unit_str)//"^3")
1098 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T61,F20.6)") &
1099 11546 : TRIM(my_label)//" Volume ["//TRIM(unit_str)//"^3]:", val
1100 11546 : val = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
1101 : WRITE (UNIT=output_unit, FMT="(T2,A,T30,3F10.3,3X,A6,F12.6)") &
1102 46184 : TRIM(my_label)//" Vector a ["//TRIM(unit_str)//"]:", cell%hmat(:, 1)*val, &
1103 11546 : "|a| = ", abc(1)*val, &
1104 46184 : TRIM(my_label)//" Vector b ["//TRIM(unit_str)//"]:", cell%hmat(:, 2)*val, &
1105 11546 : "|b| = ", abc(2)*val, &
1106 46184 : TRIM(my_label)//" Vector c ["//TRIM(unit_str)//"]:", cell%hmat(:, 3)*val, &
1107 23092 : "|c| = ", abc(3)*val
1108 : WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.6)") &
1109 11546 : TRIM(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
1110 11546 : TRIM(my_label)//" Angle (a,c), beta [degree]: ", beta, &
1111 23092 : TRIM(my_label)//" Angle (a,b), gamma [degree]: ", gamma
1112 11546 : IF (cell%symmetry_id /= cell_sym_none) THEN
1113 1464 : CALL create_cell_section(section)
1114 1464 : keyword => section_get_keyword(section, "SYMMETRY")
1115 1464 : CALL keyword_get(keyword, enum=enum)
1116 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
1117 1464 : TRIM(my_label)//" Requested initial symmetry: ", &
1118 2928 : ADJUSTR(TRIM(enum_i2c(enum, cell%symmetry_id)))
1119 1464 : CALL section_release(section)
1120 : END IF
1121 11546 : IF (cell%orthorhombic) THEN
1122 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1123 8185 : TRIM(my_label)//" Numerically orthorhombic: ", "YES"
1124 : ELSE
1125 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1126 3361 : TRIM(my_label)//" Numerically orthorhombic: ", " NO"
1127 : END IF
1128 46184 : IF (SUM(cell%perd(1:3)) == 0) THEN
1129 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
1130 2475 : TRIM(my_label)//" Periodicity", "NONE"
1131 : ELSE
1132 9071 : string = ""
1133 9071 : IF (cell%perd(1) == 1) string = TRIM(string)//"X"
1134 9071 : IF (cell%perd(2) == 1) string = TRIM(string)//"Y"
1135 9071 : IF (cell%perd(3) == 1) string = TRIM(string)//"Z"
1136 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
1137 9071 : TRIM(my_label)//" Periodicity", ADJUSTR(string)
1138 : END IF
1139 : END IF
1140 :
1141 37561 : END SUBROUTINE write_cell_low
1142 :
1143 : END MODULE cell_methods
|
