Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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: do_cell_cif,&
36 : do_cell_cp2k,&
37 : do_cell_xsc
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 : #include "./base/base_uses.f90"
58 :
59 : IMPLICIT NONE
60 :
61 : PRIVATE
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_methods'
64 :
65 : PUBLIC :: cell_create, &
66 : init_cell, &
67 : read_cell, &
68 : read_cell_cif, &
69 : set_cell_param, &
70 : write_cell
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief allocates and initializes a cell
76 : !> \param cell the cell to initialize
77 : !> \param hmat the h matrix that defines the cell
78 : !> \param periodic periodicity of the cell
79 : !> \param tag ...
80 : !> \par History
81 : !> 09.2003 created [fawzi]
82 : !> \author Fawzi Mohamed
83 : ! **************************************************************************************************
84 57736 : SUBROUTINE cell_create(cell, hmat, periodic, tag)
85 :
86 : TYPE(cell_type), POINTER :: cell
87 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
88 : OPTIONAL :: hmat
89 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
90 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
91 :
92 0 : CPASSERT(.NOT. ASSOCIATED(cell))
93 1732080 : ALLOCATE (cell)
94 57736 : cell%ref_count = 1
95 57736 : IF (PRESENT(periodic)) THEN
96 1332 : cell%perd = periodic
97 : ELSE
98 229612 : cell%perd = 1
99 : END IF
100 : cell%orthorhombic = .FALSE.
101 57736 : cell%symmetry_id = cell_sym_none
102 57736 : IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
103 57736 : IF (PRESENT(tag)) cell%tag = tag
104 :
105 57736 : END SUBROUTINE cell_create
106 :
107 : ! **************************************************************************************************
108 : !> \brief Initialise/readjust a simulation cell after hmat has been changed
109 : !> \param cell ...
110 : !> \param hmat ...
111 : !> \param periodic ...
112 : !> \date 16.01.2002
113 : !> \author Matthias Krack
114 : !> \version 1.0
115 : ! **************************************************************************************************
116 342047 : SUBROUTINE init_cell(cell, hmat, periodic)
117 :
118 : TYPE(cell_type), POINTER :: cell
119 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
120 : OPTIONAL :: hmat
121 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
122 :
123 : REAL(KIND=dp), PARAMETER :: eps_hmat = 1.0E-14_dp
124 :
125 : INTEGER :: dim
126 : REAL(KIND=dp) :: a, acosa, acosah, acosg, alpha, asina, &
127 : asinah, asing, beta, gamma, norm, &
128 : norm_c
129 : REAL(KIND=dp), DIMENSION(3) :: abc
130 :
131 342047 : CPASSERT(ASSOCIATED(cell))
132 :
133 441995 : IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
134 342125 : IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
135 :
136 342047 : cell%deth = ABS(det_3x3(cell%hmat))
137 :
138 342047 : IF (cell%deth < 1.0E-10_dp) THEN
139 0 : CALL write_cell_low(cell, "angstrom", default_output_unit)
140 : CALL cp_abort(__LOCATION__, &
141 : "An invalid set of cell vectors was specified. "// &
142 0 : "The cell volume is too small")
143 : END IF
144 :
145 346091 : SELECT CASE (cell%symmetry_id)
146 : CASE (cell_sym_cubic, &
147 : cell_sym_tetragonal_ab, &
148 : cell_sym_tetragonal_ac, &
149 : cell_sym_tetragonal_bc, &
150 : cell_sym_orthorhombic)
151 4044 : CALL get_cell(cell=cell, abc=abc)
152 4044 : abc(2) = plane_distance(0, 1, 0, cell=cell)
153 4044 : abc(3) = plane_distance(0, 0, 1, cell=cell)
154 4044 : SELECT CASE (cell%symmetry_id)
155 : CASE (cell_sym_cubic)
156 7621 : abc(1:3) = SUM(abc(1:3))/3.0_dp
157 : CASE (cell_sym_tetragonal_ab, &
158 : cell_sym_tetragonal_ac, &
159 : cell_sym_tetragonal_bc)
160 1322 : SELECT CASE (cell%symmetry_id)
161 : CASE (cell_sym_tetragonal_ab)
162 1322 : a = 0.5_dp*(abc(1) + abc(2))
163 1322 : abc(1) = a
164 1322 : abc(2) = a
165 : CASE (cell_sym_tetragonal_ac)
166 661 : a = 0.5_dp*(abc(1) + abc(3))
167 661 : abc(1) = a
168 661 : abc(3) = a
169 : CASE (cell_sym_tetragonal_bc)
170 612 : a = 0.5_dp*(abc(2) + abc(3))
171 612 : abc(2) = a
172 612 : abc(3) = a
173 : END SELECT
174 : END SELECT
175 4044 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
176 4044 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
177 4044 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
178 : CASE (cell_sym_hexagonal_gamma_60, cell_sym_hexagonal_gamma_120)
179 2634 : CALL get_cell(cell=cell, abc=abc)
180 2634 : a = 0.5_dp*(abc(1) + abc(2))
181 2634 : acosg = 0.5_dp*a
182 2634 : asing = sqrt3*acosg
183 2634 : IF (cell%symmetry_id == cell_sym_hexagonal_gamma_120) acosg = -acosg
184 2634 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
185 2634 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
186 2634 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
187 : CASE (cell_sym_rhombohedral)
188 665 : CALL get_cell(cell=cell, abc=abc)
189 2660 : a = SUM(abc(1:3))/3.0_dp
190 : alpha = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
191 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
192 665 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
193 665 : acosa = a*COS(alpha)
194 665 : asina = a*SIN(alpha)
195 665 : acosah = a*COS(0.5_dp*alpha)
196 665 : asinah = a*SIN(0.5_dp*alpha)
197 665 : norm = acosa/acosah
198 665 : norm_c = SQRT(1.0_dp - norm*norm)
199 665 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
200 665 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
201 665 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = a*norm_c
202 : CASE (cell_sym_monoclinic)
203 880 : CALL get_cell(cell=cell, abc=abc)
204 880 : beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
205 880 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)
206 880 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
207 880 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)
208 : CASE (cell_sym_monoclinic_gamma_ab)
209 : ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
210 742 : CALL get_cell(cell=cell, abc=abc)
211 742 : a = 0.5_dp*(abc(1) + abc(2))
212 742 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
213 742 : acosg = a*COS(gamma)
214 742 : asing = a*SIN(gamma)
215 742 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
216 742 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
217 342789 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
218 : CASE (cell_sym_triclinic)
219 : ! Nothing to do
220 : END SELECT
221 :
222 : ! Do we have an (almost) orthorhombic cell?
223 : IF ((ABS(cell%hmat(1, 2)) < eps_hmat) .AND. (ABS(cell%hmat(1, 3)) < eps_hmat) .AND. &
224 : (ABS(cell%hmat(2, 1)) < eps_hmat) .AND. (ABS(cell%hmat(2, 3)) < eps_hmat) .AND. &
225 342047 : (ABS(cell%hmat(3, 1)) < eps_hmat) .AND. (ABS(cell%hmat(3, 2)) < eps_hmat)) THEN
226 317361 : cell%orthorhombic = .TRUE.
227 : ELSE
228 24686 : cell%orthorhombic = .FALSE.
229 : END IF
230 :
231 : ! Retain an exact orthorhombic cell
232 : ! (off-diagonal elements must remain zero identically to keep QS fast)
233 342047 : IF (cell%orthorhombic) THEN
234 317361 : cell%hmat(1, 2) = 0.0_dp
235 317361 : cell%hmat(1, 3) = 0.0_dp
236 317361 : cell%hmat(2, 1) = 0.0_dp
237 317361 : cell%hmat(2, 3) = 0.0_dp
238 317361 : cell%hmat(3, 1) = 0.0_dp
239 317361 : cell%hmat(3, 2) = 0.0_dp
240 : END IF
241 :
242 1368188 : dim = COUNT(cell%perd == 1)
243 342047 : IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
244 0 : CPABORT("Non-orthorhombic and not periodic")
245 : END IF
246 :
247 : ! Update deth and hmat_inv with enforced symmetry
248 342047 : cell%deth = ABS(det_3x3(cell%hmat))
249 342047 : IF (cell%deth < 1.0E-10_dp) THEN
250 : CALL cp_abort(__LOCATION__, &
251 : "An invalid set of cell vectors was obtained after applying "// &
252 0 : "the requested cell symmetry. The cell volume is too small")
253 : END IF
254 4446611 : cell%h_inv = inv_3x3(cell%hmat)
255 :
256 342047 : END SUBROUTINE init_cell
257 :
258 : ! **************************************************************************************************
259 : !> \brief ...
260 : !> \param cell ...
261 : !> \param cell_ref ...
262 : !> \param use_ref_cell ...
263 : !> \param cell_section ...
264 : !> \param check_for_ref ...
265 : !> \param para_env ...
266 : !> \par History
267 : !> 03.2005 created [teo]
268 : !> \author Teodoro Laino
269 : ! **************************************************************************************************
270 129885 : RECURSIVE SUBROUTINE read_cell(cell, cell_ref, use_ref_cell, cell_section, &
271 : check_for_ref, para_env)
272 :
273 : TYPE(cell_type), POINTER :: cell, cell_ref
274 : LOGICAL, INTENT(INOUT), OPTIONAL :: use_ref_cell
275 : TYPE(section_vals_type), OPTIONAL, POINTER :: cell_section
276 : LOGICAL, INTENT(IN), OPTIONAL :: check_for_ref
277 : TYPE(mp_para_env_type), POINTER :: para_env
278 :
279 : INTEGER :: my_per
280 18555 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
281 : LOGICAL :: cell_read_a, cell_read_abc, cell_read_b, &
282 : cell_read_c, cell_read_file, check, &
283 : my_check
284 18555 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_angles, cell_par
285 : TYPE(section_vals_type), POINTER :: cell_ref_section
286 :
287 18555 : my_check = .TRUE.
288 18555 : NULLIFY (cell_ref_section, cell_par, multiple_unit_cell)
289 18511 : IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
290 18555 : IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
291 18555 : IF (PRESENT(check_for_ref)) my_check = check_for_ref
292 :
293 18555 : cell%deth = 0.0_dp
294 18555 : cell%orthorhombic = .FALSE.
295 74220 : cell%perd(:) = 1
296 18555 : cell%symmetry_id = cell_sym_none
297 241215 : cell%hmat(:, :) = 0.0_dp
298 241215 : cell%h_inv(:, :) = 0.0_dp
299 : cell_read_file = .FALSE.
300 : cell_read_a = .FALSE.
301 18555 : cell_read_b = .FALSE.
302 18555 : cell_read_c = .FALSE.
303 : ! Trying to read cell info from file
304 18555 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
305 18555 : IF (cell_read_file) CALL read_cell_from_external_file(cell_section, para_env)
306 :
307 : ! Trying to read cell info from the separate A, B, C vectors
308 : ! If cell information is provided through file A,B,C contain the file information..
309 : ! a print warning is shown on screen..
310 18555 : CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
311 18555 : IF (cell_read_a) THEN
312 2422 : CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
313 19376 : cell%hmat(:, 1) = cell_par(:)
314 : END IF
315 18555 : CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
316 18555 : IF (cell_read_b) THEN
317 2422 : CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
318 19376 : cell%hmat(:, 2) = cell_par(:)
319 : END IF
320 18555 : CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
321 18555 : IF (cell_read_c) THEN
322 2422 : CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
323 19376 : cell%hmat(:, 3) = cell_par(:)
324 : END IF
325 18555 : check = ((cell_read_a .EQV. cell_read_b) .AND. (cell_read_b .EQV. cell_read_c))
326 : IF (.NOT. check) &
327 : CALL cp_warn(__LOCATION__, &
328 : "Cell Information provided through vectors A, B and C. Not all three "// &
329 0 : "vectors were provided! Cell setup may be incomplete!")
330 :
331 : ! Very last option.. Trying to read cell info from ABC keyword
332 18555 : CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
333 18555 : IF (cell_read_abc) THEN
334 16133 : check = (cell_read_a .OR. cell_read_b .OR. cell_read_c)
335 : IF (check) &
336 : CALL cp_warn(__LOCATION__, &
337 : "Cell Information provided through vectors A, B and C in conjunction with ABC."// &
338 0 : " The definition of the ABC keyword will override the one provided by A,B and C.")
339 209729 : cell%hmat = 0.0_dp
340 16133 : CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
341 16133 : CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_angles)
342 16133 : CALL set_cell_param(cell, cell_par, cell_angles, do_init_cell=.FALSE.)
343 : END IF
344 :
345 : ! Multiple unit cell
346 18555 : CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
347 73340 : IF (ANY(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
348 :
349 18555 : CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
350 0 : SELECT CASE (my_per)
351 : CASE (use_perd_x)
352 0 : cell%perd = (/1, 0, 0/)
353 : CASE (use_perd_y)
354 0 : cell%perd = (/0, 1, 0/)
355 : CASE (use_perd_z)
356 32 : cell%perd = (/0, 0, 1/)
357 : CASE (use_perd_xy)
358 112 : cell%perd = (/1, 1, 0/)
359 : CASE (use_perd_xz)
360 48 : cell%perd = (/1, 0, 1/)
361 : CASE (use_perd_yz)
362 208 : cell%perd = (/0, 1, 1/)
363 : CASE (use_perd_xyz)
364 51788 : cell%perd = (/1, 1, 1/)
365 : CASE (use_perd_none)
366 22032 : cell%perd = (/0, 0, 0/)
367 : CASE DEFAULT
368 18555 : CPABORT("")
369 : END SELECT
370 :
371 : ! Load requested cell symmetry
372 18555 : CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
373 :
374 : ! Initialize cell
375 18555 : CALL init_cell(cell)
376 :
377 18555 : IF (my_check) THEN
378 : ! Recursive check for reference cell requested
379 18117 : cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
380 18117 : IF (parsed_cp2k_input(cell_ref_section, check_this_section=.TRUE.)) THEN
381 44 : IF (PRESENT(use_ref_cell)) use_ref_cell = .TRUE.
382 : CALL read_cell(cell_ref, cell_ref, use_ref_cell=use_ref_cell, &
383 : cell_section=cell_ref_section, check_for_ref=.FALSE., &
384 44 : para_env=para_env)
385 : ELSE
386 18073 : CALL cell_clone(cell, cell_ref, tag="CELL_REF")
387 18073 : IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE.
388 : END IF
389 : END IF
390 :
391 18555 : END SUBROUTINE read_cell
392 :
393 : ! **************************************************************************************************
394 : !> \brief utility function to ease the transition to the new input.
395 : !> returns true if the new input was parsed
396 : !> \param input_file the parsed input file
397 : !> \param check_this_section ...
398 : !> \return ...
399 : !> \author fawzi
400 : ! **************************************************************************************************
401 18117 : FUNCTION parsed_cp2k_input(input_file, check_this_section) RESULT(res)
402 :
403 : TYPE(section_vals_type), POINTER :: input_file
404 : LOGICAL, INTENT(IN), OPTIONAL :: check_this_section
405 : LOGICAL :: res
406 :
407 : LOGICAL :: my_check
408 : TYPE(section_vals_type), POINTER :: glob_section
409 :
410 18117 : my_check = .FALSE.
411 18117 : IF (PRESENT(check_this_section)) my_check = check_this_section
412 18117 : res = ASSOCIATED(input_file)
413 18117 : IF (res) THEN
414 18117 : CPASSERT(input_file%ref_count > 0)
415 18117 : IF (.NOT. my_check) THEN
416 0 : glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
417 0 : CALL section_vals_get(glob_section, explicit=res)
418 : ELSE
419 18117 : CALL section_vals_get(input_file, explicit=res)
420 : END IF
421 : END IF
422 :
423 18117 : END FUNCTION parsed_cp2k_input
424 :
425 : ! **************************************************************************************************
426 : !> \brief Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma)
427 : !> using the convention: a parallel to the x axis, b in the x-y plane and
428 : !> and c univoquely determined; gamma is the angle between a and b; beta
429 : !> is the angle between c and a and alpha is the angle between c and b
430 : !> \param cell ...
431 : !> \param cell_length ...
432 : !> \param cell_angle ...
433 : !> \param periodic ...
434 : !> \param do_init_cell ...
435 : !> \date 03.2008
436 : !> \author Teodoro Laino
437 : ! **************************************************************************************************
438 16159 : SUBROUTINE set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
439 :
440 : TYPE(cell_type), POINTER :: cell
441 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: cell_length, cell_angle
442 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
443 : LOGICAL, INTENT(IN) :: do_init_cell
444 :
445 : REAL(KIND=dp), PARAMETER :: eps = EPSILON(0.0_dp)
446 :
447 : REAL(KIND=dp) :: cos_alpha, cos_beta, cos_gamma, sin_gamma
448 :
449 16159 : CPASSERT(ASSOCIATED(cell))
450 64636 : CPASSERT(ALL(cell_angle /= 0.0_dp))
451 :
452 16159 : cos_gamma = COS(cell_angle(3)); IF (ABS(cos_gamma) < eps) cos_gamma = 0.0_dp
453 16159 : IF (ABS(ABS(cos_gamma) - 1.0_dp) < eps) cos_gamma = SIGN(1.0_dp, cos_gamma)
454 16159 : sin_gamma = SIN(cell_angle(3)); IF (ABS(sin_gamma) < eps) sin_gamma = 0.0_dp
455 16159 : IF (ABS(ABS(sin_gamma) - 1.0_dp) < eps) sin_gamma = SIGN(1.0_dp, sin_gamma)
456 16159 : cos_beta = COS(cell_angle(2)); IF (ABS(cos_beta) < eps) cos_beta = 0.0_dp
457 16159 : IF (ABS(ABS(cos_beta) - 1.0_dp) < eps) cos_beta = SIGN(1.0_dp, cos_beta)
458 16159 : cos_alpha = COS(cell_angle(1)); IF (ABS(cos_alpha) < eps) cos_alpha = 0.0_dp
459 16159 : IF (ABS(ABS(cos_alpha) - 1.0_dp) < eps) cos_alpha = SIGN(1.0_dp, cos_alpha)
460 :
461 64636 : cell%hmat(:, 1) = (/1.0_dp, 0.0_dp, 0.0_dp/)
462 64636 : cell%hmat(:, 2) = (/cos_gamma, sin_gamma, 0.0_dp/)
463 64636 : cell%hmat(:, 3) = (/cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp/)
464 16159 : cell%hmat(3, 3) = SQRT(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
465 :
466 64636 : cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
467 64636 : cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
468 64636 : cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
469 :
470 16159 : IF (do_init_cell) THEN
471 26 : IF (PRESENT(periodic)) THEN
472 26 : CALL init_cell(cell=cell, periodic=periodic)
473 : ELSE
474 0 : CALL init_cell(cell=cell)
475 : END IF
476 : END IF
477 :
478 16159 : END SUBROUTINE set_cell_param
479 :
480 : ! **************************************************************************************************
481 : !> \brief Setup of the multiple unit_cell
482 : !> \param cell ...
483 : !> \param multiple_unit_cell ...
484 : !> \date 05.2009
485 : !> \author Teodoro Laino [tlaino]
486 : !> \version 1.0
487 : ! **************************************************************************************************
488 300 : SUBROUTINE set_multiple_unit_cell(cell, multiple_unit_cell)
489 :
490 : TYPE(cell_type), POINTER :: cell
491 : INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell
492 :
493 300 : CPASSERT(ASSOCIATED(cell))
494 :
495 : ! Abort, if one of the value is set to zero
496 1200 : IF (ANY(multiple_unit_cell <= 0)) &
497 : CALL cp_abort(__LOCATION__, &
498 : "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
499 0 : "A value of 0 or negative is meaningless!")
500 :
501 : ! Scale abc according to user request
502 1200 : cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
503 1200 : cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
504 1200 : cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
505 :
506 300 : END SUBROUTINE set_multiple_unit_cell
507 :
508 : ! **************************************************************************************************
509 : !> \brief Read cell information from an external file
510 : !> \param cell_section ...
511 : !> \param para_env ...
512 : !> \date 02.2008
513 : !> \author Teodoro Laino [tlaino] - University of Zurich
514 : !> \version 1.0
515 : ! **************************************************************************************************
516 12 : SUBROUTINE read_cell_from_external_file(cell_section, para_env)
517 :
518 : TYPE(section_vals_type), POINTER :: cell_section
519 : TYPE(mp_para_env_type), POINTER :: para_env
520 :
521 : CHARACTER(LEN=default_path_length) :: cell_file_name
522 : INTEGER :: i, idum, j, my_format, n_rep
523 : LOGICAL :: explicit, my_end
524 : REAL(KIND=dp) :: xdum
525 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
526 12 : REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par
527 : TYPE(cell_type), POINTER :: cell
528 : TYPE(cp_parser_type) :: parser
529 :
530 12 : CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
531 12 : CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=my_format)
532 2 : SELECT CASE (my_format)
533 : CASE (do_cell_cp2k)
534 2 : CALL parser_create(parser, cell_file_name, para_env=para_env)
535 2 : CALL parser_get_next_line(parser, 1)
536 2 : my_end = .FALSE.
537 24 : DO WHILE (.NOT. my_end)
538 22 : READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
539 24 : CALL parser_get_next_line(parser, 1, at_end=my_end)
540 : END DO
541 2 : CALL parser_release(parser)
542 : ! Convert to CP2K units
543 8 : DO i = 1, 3
544 26 : DO j = 1, 3
545 24 : hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
546 : END DO
547 : END DO
548 : CASE (do_cell_xsc)
549 0 : CALL parser_create(parser, cell_file_name, para_env=para_env)
550 0 : CALL parser_get_next_line(parser, 1)
551 0 : READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
552 0 : CALL parser_release(parser)
553 : ! Convert to CP2K units
554 0 : DO i = 1, 3
555 0 : DO j = 1, 3
556 0 : hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
557 : END DO
558 : END DO
559 : CASE (do_cell_cif)
560 10 : NULLIFY (cell)
561 10 : CALL cell_create(cell)
562 10 : CALL read_cell_cif(cell_file_name, cell, para_env)
563 130 : hmat = cell%hmat
564 22 : CALL cell_release(cell)
565 : END SELECT
566 12 : CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
567 12 : CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
568 : ! Check if the cell was already defined
569 12 : explicit = .FALSE.
570 12 : CALL section_vals_val_get(cell_section, "A", n_rep_val=n_rep)
571 12 : explicit = explicit .OR. (n_rep == 1)
572 12 : CALL section_vals_val_get(cell_section, "B", n_rep_val=n_rep)
573 12 : explicit = explicit .OR. (n_rep == 1)
574 12 : CALL section_vals_val_get(cell_section, "C", n_rep_val=n_rep)
575 12 : explicit = explicit .OR. (n_rep == 1)
576 12 : CALL section_vals_val_get(cell_section, "ABC", n_rep_val=n_rep)
577 12 : explicit = explicit .OR. (n_rep == 1)
578 : ! Possibly print a warning
579 : IF (explicit) &
580 : CALL cp_warn(__LOCATION__, &
581 : "Cell specification (A,B,C or ABC) provided together with the external "// &
582 : "cell setup! Ignoring (A,B,C or ABC) and proceeding with info read from the "// &
583 0 : "external file! ")
584 : ! Copy cell information in the A, B, C fields..(we may need them later on..)
585 12 : ALLOCATE (cell_par(3))
586 48 : cell_par = hmat(:, 1)
587 12 : CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
588 12 : ALLOCATE (cell_par(3))
589 48 : cell_par = hmat(:, 2)
590 12 : CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
591 12 : ALLOCATE (cell_par(3))
592 48 : cell_par = hmat(:, 3)
593 12 : CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
594 : ! Unset possible keywords
595 12 : CALL section_vals_val_unset(cell_section, "ABC")
596 12 : CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
597 :
598 84 : END SUBROUTINE read_cell_from_external_file
599 :
600 : ! **************************************************************************************************
601 : !> \brief Reads cell information from CIF file
602 : !> \param cif_file_name ...
603 : !> \param cell ...
604 : !> \param para_env ...
605 : !> \date 12.2008
606 : !> \par Format Information implemented:
607 : !> _cell_length_a
608 : !> _cell_length_b
609 : !> _cell_length_c
610 : !> _cell_angle_alpha
611 : !> _cell_angle_beta
612 : !> _cell_angle_gamma
613 : !>
614 : !> \author Teodoro Laino [tlaino]
615 : !> moved from topology_cif (1/2019 JHU)
616 : ! **************************************************************************************************
617 48 : SUBROUTINE read_cell_cif(cif_file_name, cell, para_env)
618 :
619 : CHARACTER(len=*) :: cif_file_name
620 : TYPE(cell_type), POINTER :: cell
621 : TYPE(mp_para_env_type), POINTER :: para_env
622 :
623 : CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_cif'
624 :
625 : INTEGER :: handle
626 : INTEGER, DIMENSION(3) :: periodic
627 : LOGICAL :: found
628 : REAL(KIND=dp), DIMENSION(3) :: cell_angles, cell_lengths
629 : TYPE(cp_parser_type) :: parser
630 :
631 24 : CALL timeset(routineN, handle)
632 :
633 : CALL parser_create(parser, cif_file_name, &
634 24 : para_env=para_env, apply_preprocessing=.FALSE.)
635 :
636 : ! Parsing cell infos
637 96 : periodic = 1
638 : ! Check for _cell_length_a
639 : CALL parser_search_string(parser, "_cell_length_a", ignore_case=.FALSE., found=found, &
640 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
641 24 : IF (.NOT. found) &
642 0 : CPABORT("The field (_cell_length_a) was not found in CIF file! ")
643 24 : CALL cif_get_real(parser, cell_lengths(1))
644 24 : cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
645 :
646 : ! Check for _cell_length_b
647 : CALL parser_search_string(parser, "_cell_length_b", ignore_case=.FALSE., found=found, &
648 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
649 24 : IF (.NOT. found) &
650 0 : CPABORT("The field (_cell_length_b) was not found in CIF file! ")
651 24 : CALL cif_get_real(parser, cell_lengths(2))
652 24 : cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
653 :
654 : ! Check for _cell_length_c
655 : CALL parser_search_string(parser, "_cell_length_c", ignore_case=.FALSE., found=found, &
656 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
657 24 : IF (.NOT. found) &
658 0 : CPABORT("The field (_cell_length_c) was not found in CIF file! ")
659 24 : CALL cif_get_real(parser, cell_lengths(3))
660 24 : cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
661 :
662 : ! Check for _cell_angle_alpha
663 : CALL parser_search_string(parser, "_cell_angle_alpha", ignore_case=.FALSE., found=found, &
664 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
665 24 : IF (.NOT. found) &
666 0 : CPABORT("The field (_cell_angle_alpha) was not found in CIF file! ")
667 24 : CALL cif_get_real(parser, cell_angles(1))
668 24 : cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
669 :
670 : ! Check for _cell_angle_beta
671 : CALL parser_search_string(parser, "_cell_angle_beta", ignore_case=.FALSE., found=found, &
672 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
673 24 : IF (.NOT. found) &
674 0 : CPABORT("The field (_cell_angle_beta) was not found in CIF file! ")
675 24 : CALL cif_get_real(parser, cell_angles(2))
676 24 : cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
677 :
678 : ! Check for _cell_angle_gamma
679 : CALL parser_search_string(parser, "_cell_angle_gamma", ignore_case=.FALSE., found=found, &
680 24 : begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
681 24 : IF (.NOT. found) &
682 0 : CPABORT("The field (_cell_angle_gamma) was not found in CIF file! ")
683 24 : CALL cif_get_real(parser, cell_angles(3))
684 24 : cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
685 :
686 : ! Create cell
687 : CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
688 24 : do_init_cell=.TRUE.)
689 :
690 24 : CALL parser_release(parser)
691 :
692 24 : CALL timestop(handle)
693 :
694 72 : END SUBROUTINE read_cell_cif
695 :
696 : ! **************************************************************************************************
697 : !> \brief Reads REAL from the CIF file.. This wrapper is needed in order to
698 : !> treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
699 : !> \param parser ...
700 : !> \param r ...
701 : !> \date 12.2008
702 : !> \author Teodoro Laino [tlaino]
703 : ! **************************************************************************************************
704 144 : SUBROUTINE cif_get_real(parser, r)
705 :
706 : TYPE(cp_parser_type), INTENT(INOUT) :: parser
707 : REAL(KIND=dp), INTENT(OUT) :: r
708 :
709 : CHARACTER(LEN=default_string_length) :: s_tag
710 : INTEGER :: iln
711 :
712 144 : CALL parser_get_object(parser, s_tag)
713 144 : iln = LEN_TRIM(s_tag)
714 144 : IF (INDEX(s_tag, "(") /= 0) iln = INDEX(s_tag, "(") - 1
715 144 : READ (s_tag(1:iln), *) r
716 :
717 144 : END SUBROUTINE cif_get_real
718 :
719 : ! **************************************************************************************************
720 : !> \brief Write the cell parameters to the output unit.
721 : !> \param cell ...
722 : !> \param subsys_section ...
723 : !> \param tag ...
724 : !> \date 02.06.2000
725 : !> \par History
726 : !> - 11.2008 Teodoro Laino [tlaino] - rewrite and enabling user driven units
727 : !> \author Matthias Krack
728 : !> \version 1.0
729 : ! **************************************************************************************************
730 34201 : SUBROUTINE write_cell(cell, subsys_section, tag)
731 :
732 : TYPE(cell_type), POINTER :: cell
733 : TYPE(section_vals_type), POINTER :: subsys_section
734 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
735 :
736 : CHARACTER(LEN=default_string_length) :: label, unit_str
737 : INTEGER :: output_unit
738 : TYPE(cp_logger_type), POINTER :: logger
739 :
740 34201 : NULLIFY (logger)
741 34201 : logger => cp_get_default_logger()
742 34201 : IF (PRESENT(tag)) THEN
743 8984 : label = TRIM(tag)//"|"
744 : ELSE
745 25217 : label = TRIM(cell%tag)//"|"
746 : END IF
747 :
748 34201 : output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", extension=".Log")
749 34201 : CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
750 34201 : CALL write_cell_low(cell, unit_str, output_unit, label)
751 34201 : CALL cp_print_key_finished_output(output_unit, logger, subsys_section, "PRINT%CELL")
752 :
753 34201 : END SUBROUTINE write_cell
754 :
755 : ! **************************************************************************************************
756 : !> \brief Write the cell parameters to the output unit
757 : !> \param cell ...
758 : !> \param unit_str ...
759 : !> \param output_unit ...
760 : !> \param label ...
761 : !> \date 17.05.2023
762 : !> \par History
763 : !> - Extracted from write_cell (17.05.2023, MK)
764 : !> \version 1.0
765 : ! **************************************************************************************************
766 34201 : SUBROUTINE write_cell_low(cell, unit_str, output_unit, label)
767 :
768 : TYPE(cell_type), POINTER :: cell
769 : CHARACTER(LEN=*), INTENT(IN) :: unit_str
770 : INTEGER, INTENT(IN) :: output_unit
771 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: label
772 :
773 : CHARACTER(LEN=12) :: tag
774 : CHARACTER(LEN=3) :: string
775 : CHARACTER(LEN=default_string_length) :: my_label
776 : REAL(KIND=dp) :: alpha, beta, gamma, val
777 : REAL(KIND=dp), DIMENSION(3) :: abc
778 : TYPE(enumeration_type), POINTER :: enum
779 : TYPE(keyword_type), POINTER :: keyword
780 : TYPE(section_type), POINTER :: section
781 :
782 34201 : NULLIFY (enum)
783 34201 : NULLIFY (keyword)
784 34201 : NULLIFY (section)
785 :
786 45373 : IF (output_unit > 0) THEN
787 11172 : CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma, tag=tag)
788 11172 : IF (PRESENT(label)) THEN
789 11172 : my_label = label
790 : ELSE
791 0 : my_label = TRIM(tag)//"|"
792 : END IF
793 11172 : val = cp_unit_from_cp2k(cell%deth, TRIM(unit_str)//"^3")
794 : WRITE (UNIT=output_unit, FMT="(/,T2,A,T61,F20.6)") &
795 11172 : TRIM(my_label)//" Volume ["//TRIM(unit_str)//"^3]:", val
796 11172 : val = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
797 : WRITE (UNIT=output_unit, FMT="(T2,A,T30,3F10.3,3X,A6,F12.6)") &
798 44688 : TRIM(my_label)//" Vector a ["//TRIM(unit_str)//"]:", cell%hmat(:, 1)*val, &
799 11172 : "|a| = ", abc(1)*val, &
800 44688 : TRIM(my_label)//" Vector b ["//TRIM(unit_str)//"]:", cell%hmat(:, 2)*val, &
801 11172 : "|b| = ", abc(2)*val, &
802 44688 : TRIM(my_label)//" Vector c ["//TRIM(unit_str)//"]:", cell%hmat(:, 3)*val, &
803 22344 : "|c| = ", abc(3)*val
804 : WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.6)") &
805 11172 : TRIM(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
806 11172 : TRIM(my_label)//" Angle (a,c), beta [degree]: ", beta, &
807 22344 : TRIM(my_label)//" Angle (a,b), gamma [degree]: ", gamma
808 11172 : IF (cell%symmetry_id /= cell_sym_none) THEN
809 1453 : CALL create_cell_section(section)
810 1453 : keyword => section_get_keyword(section, "SYMMETRY")
811 1453 : CALL keyword_get(keyword, enum=enum)
812 : WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
813 1453 : TRIM(my_label)//" Requested initial symmetry: ", &
814 2906 : ADJUSTR(TRIM(enum_i2c(enum, cell%symmetry_id)))
815 1453 : CALL section_release(section)
816 : END IF
817 11172 : IF (cell%orthorhombic) THEN
818 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
819 7942 : TRIM(my_label)//" Numerically orthorhombic: ", "YES"
820 : ELSE
821 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
822 3230 : TRIM(my_label)//" Numerically orthorhombic: ", " NO"
823 : END IF
824 44688 : IF (SUM(cell%perd(1:3)) == 0) THEN
825 : WRITE (UNIT=output_unit, FMT="(T2,A,T77,A4)") &
826 2312 : TRIM(my_label)//" Periodicity", "NONE"
827 : ELSE
828 8860 : string = ""
829 8860 : IF (cell%perd(1) == 1) string = TRIM(string)//"X"
830 8860 : IF (cell%perd(2) == 1) string = TRIM(string)//"Y"
831 8860 : IF (cell%perd(3) == 1) string = TRIM(string)//"Z"
832 : WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
833 8860 : TRIM(my_label)//" Periodicity", ADJUSTR(string)
834 : END IF
835 : END IF
836 :
837 34201 : END SUBROUTINE write_cell_low
838 :
839 : END MODULE cell_methods
|
