(git:b279b6b)
cell_methods.F
Go to the documentation of this file.
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 ! **************************************************************************************************
16  USE cell_types, ONLY: &
24  cp_logger_type
28  parser_get_object,&
30  USE cp_parser_types, ONLY: cp_parser_type,&
33  USE cp_units, ONLY: cp_unit_from_cp2k,&
35  USE input_constants, ONLY: do_cell_cif,&
36  do_cell_cp2k,&
40  enumeration_type
41  USE input_keyword_types, ONLY: keyword_get,&
42  keyword_type
43  USE input_section_types, ONLY: &
47  USE kinds, ONLY: default_path_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, &
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  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  cpassert(.NOT. ASSOCIATED(cell))
93  ALLOCATE (cell)
94  cell%ref_count = 1
95  IF (PRESENT(periodic)) THEN
96  cell%perd = periodic
97  ELSE
98  cell%perd = 1
99  END IF
100  cell%orthorhombic = .false.
101  cell%symmetry_id = cell_sym_none
102  IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
103  IF (PRESENT(tag)) cell%tag = tag
104 
105  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  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  cpassert(ASSOCIATED(cell))
132 
133  IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
134  IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
135 
136  cell%deth = abs(det_3x3(cell%hmat))
137 
138  IF (cell%deth < 1.0e-10_dp) THEN
139  CALL write_cell_low(cell, "angstrom", default_output_unit)
140  CALL cp_abort(__location__, &
141  "An invalid set of cell vectors was specified. "// &
142  "The cell volume is too small")
143  END IF
144 
145  SELECT CASE (cell%symmetry_id)
146  CASE (cell_sym_cubic, &
151  CALL get_cell(cell=cell, abc=abc)
152  abc(2) = plane_distance(0, 1, 0, cell=cell)
153  abc(3) = plane_distance(0, 0, 1, cell=cell)
154  SELECT CASE (cell%symmetry_id)
155  CASE (cell_sym_cubic)
156  abc(1:3) = sum(abc(1:3))/3.0_dp
157  CASE (cell_sym_tetragonal_ab, &
160  SELECT CASE (cell%symmetry_id)
162  a = 0.5_dp*(abc(1) + abc(2))
163  abc(1) = a
164  abc(2) = a
166  a = 0.5_dp*(abc(1) + abc(3))
167  abc(1) = a
168  abc(3) = a
170  a = 0.5_dp*(abc(2) + abc(3))
171  abc(2) = a
172  abc(3) = a
173  END SELECT
174  END SELECT
175  cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
176  cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
177  cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
179  CALL get_cell(cell=cell, abc=abc)
180  a = 0.5_dp*(abc(1) + abc(2))
181  acosg = 0.5_dp*a
182  asing = sqrt3*acosg
183  IF (cell%symmetry_id == cell_sym_hexagonal_gamma_120) acosg = -acosg
184  cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
185  cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
186  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  CALL get_cell(cell=cell, abc=abc)
189  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  angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
193  acosa = a*cos(alpha)
194  asina = a*sin(alpha)
195  acosah = a*cos(0.5_dp*alpha)
196  asinah = a*sin(0.5_dp*alpha)
197  norm = acosa/acosah
198  norm_c = sqrt(1.0_dp - norm*norm)
199  cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
200  cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
201  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  CALL get_cell(cell=cell, abc=abc)
204  beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
205  cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*cos(beta)
206  cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
207  cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*sin(beta)
209  ! Cell symmetry with a = b, alpha = beta = 90 degree and gammma not equal 90 degree
210  CALL get_cell(cell=cell, abc=abc)
211  a = 0.5_dp*(abc(1) + abc(2))
212  gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
213  acosg = a*cos(gamma)
214  asing = a*sin(gamma)
215  cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosg; cell%hmat(1, 3) = 0.0_dp
216  cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asing; cell%hmat(2, 3) = 0.0_dp
217  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  (abs(cell%hmat(3, 1)) < eps_hmat) .AND. (abs(cell%hmat(3, 2)) < eps_hmat)) THEN
226  cell%orthorhombic = .true.
227  ELSE
228  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  IF (cell%orthorhombic) THEN
234  cell%hmat(1, 2) = 0.0_dp
235  cell%hmat(1, 3) = 0.0_dp
236  cell%hmat(2, 1) = 0.0_dp
237  cell%hmat(2, 3) = 0.0_dp
238  cell%hmat(3, 1) = 0.0_dp
239  cell%hmat(3, 2) = 0.0_dp
240  END IF
241 
242  dim = count(cell%perd == 1)
243  IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
244  cpabort("Non-orthorhombic and not periodic")
245  END IF
246 
247  ! Update deth and hmat_inv with enforced symmetry
248  cell%deth = abs(det_3x3(cell%hmat))
249  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  "the requested cell symmetry. The cell volume is too small")
253  END IF
254  cell%h_inv = inv_3x3(cell%hmat)
255 
256  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  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: cell_angles, cell_par
285  TYPE(section_vals_type), POINTER :: cell_ref_section
286 
287  my_check = .true.
288  NULLIFY (cell_ref_section, cell_par, multiple_unit_cell)
289  IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell, tag="CELL")
290  IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref, tag="CELL_REF")
291  IF (PRESENT(check_for_ref)) my_check = check_for_ref
292 
293  cell%deth = 0.0_dp
294  cell%orthorhombic = .false.
295  cell%perd(:) = 1
296  cell%symmetry_id = cell_sym_none
297  cell%hmat(:, :) = 0.0_dp
298  cell%h_inv(:, :) = 0.0_dp
299  cell_read_file = .false.
300  cell_read_a = .false.
301  cell_read_b = .false.
302  cell_read_c = .false.
303  ! Trying to read cell info from file
304  CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
305  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  CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
311  IF (cell_read_a) THEN
312  CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
313  cell%hmat(:, 1) = cell_par(:)
314  END IF
315  CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
316  IF (cell_read_b) THEN
317  CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
318  cell%hmat(:, 2) = cell_par(:)
319  END IF
320  CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
321  IF (cell_read_c) THEN
322  CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
323  cell%hmat(:, 3) = cell_par(:)
324  END IF
325  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  "vectors were provided! Cell setup may be incomplete!")
330 
331  ! Very last option.. Trying to read cell info from ABC keyword
332  CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
333  IF (cell_read_abc) THEN
334  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  " The definition of the ABC keyword will override the one provided by A,B and C.")
339  cell%hmat = 0.0_dp
340  CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
341  CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_angles)
342  CALL set_cell_param(cell, cell_par, cell_angles, do_init_cell=.false.)
343  END IF
344 
345  ! Multiple unit cell
346  CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
347  IF (any(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
348 
349  CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
350  SELECT CASE (my_per)
351  CASE (use_perd_x)
352  cell%perd = (/1, 0, 0/)
353  CASE (use_perd_y)
354  cell%perd = (/0, 1, 0/)
355  CASE (use_perd_z)
356  cell%perd = (/0, 0, 1/)
357  CASE (use_perd_xy)
358  cell%perd = (/1, 1, 0/)
359  CASE (use_perd_xz)
360  cell%perd = (/1, 0, 1/)
361  CASE (use_perd_yz)
362  cell%perd = (/0, 1, 1/)
363  CASE (use_perd_xyz)
364  cell%perd = (/1, 1, 1/)
365  CASE (use_perd_none)
366  cell%perd = (/0, 0, 0/)
367  CASE DEFAULT
368  cpabort("")
369  END SELECT
370 
371  ! Load requested cell symmetry
372  CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
373 
374  ! Initialize cell
375  CALL init_cell(cell)
376 
377  IF (my_check) THEN
378  ! Recursive check for reference cell requested
379  cell_ref_section => section_vals_get_subs_vals(cell_section, "CELL_REF")
380  IF (parsed_cp2k_input(cell_ref_section, check_this_section=.true.)) THEN
381  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  para_env=para_env)
385  ELSE
386  CALL cell_clone(cell, cell_ref, tag="CELL_REF")
387  IF (PRESENT(use_ref_cell)) use_ref_cell = .false.
388  END IF
389  END IF
390 
391  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  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  my_check = .false.
411  IF (PRESENT(check_this_section)) my_check = check_this_section
412  res = ASSOCIATED(input_file)
413  IF (res) THEN
414  cpassert(input_file%ref_count > 0)
415  IF (.NOT. my_check) THEN
416  glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
417  CALL section_vals_get(glob_section, explicit=res)
418  ELSE
419  CALL section_vals_get(input_file, explicit=res)
420  END IF
421  END IF
422 
423  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  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  cpassert(ASSOCIATED(cell))
450  cpassert(all(cell_angle /= 0.0_dp))
451 
452  cos_gamma = cos(cell_angle(3)); IF (abs(cos_gamma) < eps) cos_gamma = 0.0_dp
453  IF (abs(abs(cos_gamma) - 1.0_dp) < eps) cos_gamma = sign(1.0_dp, cos_gamma)
454  sin_gamma = sin(cell_angle(3)); IF (abs(sin_gamma) < eps) sin_gamma = 0.0_dp
455  IF (abs(abs(sin_gamma) - 1.0_dp) < eps) sin_gamma = sign(1.0_dp, sin_gamma)
456  cos_beta = cos(cell_angle(2)); IF (abs(cos_beta) < eps) cos_beta = 0.0_dp
457  IF (abs(abs(cos_beta) - 1.0_dp) < eps) cos_beta = sign(1.0_dp, cos_beta)
458  cos_alpha = cos(cell_angle(1)); IF (abs(cos_alpha) < eps) cos_alpha = 0.0_dp
459  IF (abs(abs(cos_alpha) - 1.0_dp) < eps) cos_alpha = sign(1.0_dp, cos_alpha)
460 
461  cell%hmat(:, 1) = (/1.0_dp, 0.0_dp, 0.0_dp/)
462  cell%hmat(:, 2) = (/cos_gamma, sin_gamma, 0.0_dp/)
463  cell%hmat(:, 3) = (/cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp/)
464  cell%hmat(3, 3) = sqrt(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
465 
466  cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
467  cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
468  cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
469 
470  IF (do_init_cell) THEN
471  IF (PRESENT(periodic)) THEN
472  CALL init_cell(cell=cell, periodic=periodic)
473  ELSE
474  CALL init_cell(cell=cell)
475  END IF
476  END IF
477 
478  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  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  cpassert(ASSOCIATED(cell))
494 
495  ! Abort, if one of the value is set to zero
496  IF (any(multiple_unit_cell <= 0)) &
497  CALL cp_abort(__location__, &
498  "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
499  "A value of 0 or negative is meaningless!")
500 
501  ! Scale abc according to user request
502  cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
503  cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
504  cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
505 
506  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: cell_par
527  TYPE(cell_type), POINTER :: cell
528  TYPE(cp_parser_type) :: parser
529 
530  CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
531  CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=my_format)
532  SELECT CASE (my_format)
533  CASE (do_cell_cp2k)
534  CALL parser_create(parser, cell_file_name, para_env=para_env)
535  CALL parser_get_next_line(parser, 1)
536  my_end = .false.
537  DO WHILE (.NOT. my_end)
538  READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
539  CALL parser_get_next_line(parser, 1, at_end=my_end)
540  END DO
541  CALL parser_release(parser)
542  ! Convert to CP2K units
543  DO i = 1, 3
544  DO j = 1, 3
545  hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
546  END DO
547  END DO
548  CASE (do_cell_xsc)
549  CALL parser_create(parser, cell_file_name, para_env=para_env)
550  CALL parser_get_next_line(parser, 1)
551  READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
552  CALL parser_release(parser)
553  ! Convert to CP2K units
554  DO i = 1, 3
555  DO j = 1, 3
556  hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
557  END DO
558  END DO
559  CASE (do_cell_cif)
560  NULLIFY (cell)
561  CALL cell_create(cell)
562  CALL read_cell_cif(cell_file_name, cell, para_env)
563  hmat = cell%hmat
564  CALL cell_release(cell)
565  END SELECT
566  CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
567  CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
568  ! Check if the cell was already defined
569  explicit = .false.
570  CALL section_vals_val_get(cell_section, "A", n_rep_val=n_rep)
571  explicit = explicit .OR. (n_rep == 1)
572  CALL section_vals_val_get(cell_section, "B", n_rep_val=n_rep)
573  explicit = explicit .OR. (n_rep == 1)
574  CALL section_vals_val_get(cell_section, "C", n_rep_val=n_rep)
575  explicit = explicit .OR. (n_rep == 1)
576  CALL section_vals_val_get(cell_section, "ABC", n_rep_val=n_rep)
577  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  "external file! ")
584  ! Copy cell information in the A, B, C fields..(we may need them later on..)
585  ALLOCATE (cell_par(3))
586  cell_par = hmat(:, 1)
587  CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
588  ALLOCATE (cell_par(3))
589  cell_par = hmat(:, 2)
590  CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
591  ALLOCATE (cell_par(3))
592  cell_par = hmat(:, 3)
593  CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
594  ! Unset possible keywords
595  CALL section_vals_val_unset(cell_section, "ABC")
596  CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
597 
598  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  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  CALL timeset(routinen, handle)
632 
633  CALL parser_create(parser, cif_file_name, &
634  para_env=para_env, apply_preprocessing=.false.)
635 
636  ! Parsing cell infos
637  periodic = 1
638  ! Check for _cell_length_a
639  CALL parser_search_string(parser, "_cell_length_a", ignore_case=.false., found=found, &
640  begin_line=.false., search_from_begin_of_file=.true.)
641  IF (.NOT. found) &
642  cpabort("The field (_cell_length_a) was not found in CIF file! ")
643  CALL cif_get_real(parser, cell_lengths(1))
644  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  begin_line=.false., search_from_begin_of_file=.true.)
649  IF (.NOT. found) &
650  cpabort("The field (_cell_length_b) was not found in CIF file! ")
651  CALL cif_get_real(parser, cell_lengths(2))
652  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  begin_line=.false., search_from_begin_of_file=.true.)
657  IF (.NOT. found) &
658  cpabort("The field (_cell_length_c) was not found in CIF file! ")
659  CALL cif_get_real(parser, cell_lengths(3))
660  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  begin_line=.false., search_from_begin_of_file=.true.)
665  IF (.NOT. found) &
666  cpabort("The field (_cell_angle_alpha) was not found in CIF file! ")
667  CALL cif_get_real(parser, cell_angles(1))
668  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  begin_line=.false., search_from_begin_of_file=.true.)
673  IF (.NOT. found) &
674  cpabort("The field (_cell_angle_beta) was not found in CIF file! ")
675  CALL cif_get_real(parser, cell_angles(2))
676  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  begin_line=.false., search_from_begin_of_file=.true.)
681  IF (.NOT. found) &
682  cpabort("The field (_cell_angle_gamma) was not found in CIF file! ")
683  CALL cif_get_real(parser, cell_angles(3))
684  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  do_init_cell=.true.)
689 
690  CALL parser_release(parser)
691 
692  CALL timestop(handle)
693 
694  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  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  CALL parser_get_object(parser, s_tag)
713  iln = len_trim(s_tag)
714  IF (index(s_tag, "(") /= 0) iln = index(s_tag, "(") - 1
715  READ (s_tag(1:iln), *) r
716 
717  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  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  NULLIFY (logger)
741  logger => cp_get_default_logger()
742  IF (PRESENT(tag)) THEN
743  label = trim(tag)//"|"
744  ELSE
745  label = trim(cell%tag)//"|"
746  END IF
747 
748  output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", extension=".Log")
749  CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
750  CALL write_cell_low(cell, unit_str, output_unit, label)
751  CALL cp_print_key_finished_output(output_unit, logger, subsys_section, "PRINT%CELL")
752 
753  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  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  NULLIFY (enum)
783  NULLIFY (keyword)
784  NULLIFY (section)
785 
786  IF (output_unit > 0) THEN
787  CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma, tag=tag)
788  IF (PRESENT(label)) THEN
789  my_label = label
790  ELSE
791  my_label = trim(tag)//"|"
792  END IF
793  val = cp_unit_from_cp2k(cell%deth, trim(unit_str)//"^3")
794  WRITE (unit=output_unit, fmt="(/,T2,A,T61,F20.6)") &
795  trim(my_label)//" Volume ["//trim(unit_str)//"^3]:", val
796  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  trim(my_label)//" Vector a ["//trim(unit_str)//"]:", cell%hmat(:, 1)*val, &
799  "|a| = ", abc(1)*val, &
800  trim(my_label)//" Vector b ["//trim(unit_str)//"]:", cell%hmat(:, 2)*val, &
801  "|b| = ", abc(2)*val, &
802  trim(my_label)//" Vector c ["//trim(unit_str)//"]:", cell%hmat(:, 3)*val, &
803  "|c| = ", abc(3)*val
804  WRITE (unit=output_unit, fmt="(T2,A,T69,F12.6)") &
805  trim(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
806  trim(my_label)//" Angle (a,c), beta [degree]: ", beta, &
807  trim(my_label)//" Angle (a,b), gamma [degree]: ", gamma
808  IF (cell%symmetry_id /= cell_sym_none) THEN
809  CALL create_cell_section(section)
810  keyword => section_get_keyword(section, "SYMMETRY")
811  CALL keyword_get(keyword, enum=enum)
812  WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
813  trim(my_label)//" Requested initial symmetry: ", &
814  adjustr(trim(enum_i2c(enum, cell%symmetry_id)))
815  CALL section_release(section)
816  END IF
817  IF (cell%orthorhombic) THEN
818  WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
819  trim(my_label)//" Numerically orthorhombic: ", "YES"
820  ELSE
821  WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
822  trim(my_label)//" Numerically orthorhombic: ", " NO"
823  END IF
824  IF (sum(cell%perd(1:3)) == 0) THEN
825  WRITE (unit=output_unit, fmt="(T2,A,T77,A4)") &
826  trim(my_label)//" Periodicity", "NONE"
827  ELSE
828  string = ""
829  IF (cell%perd(1) == 1) string = trim(string)//"X"
830  IF (cell%perd(2) == 1) string = trim(string)//"Y"
831  IF (cell%perd(3) == 1) string = trim(string)//"Z"
832  WRITE (unit=output_unit, fmt="(T2,A,T78,A3)") &
833  trim(my_label)//" Periodicity", adjustr(string)
834  END IF
835  END IF
836 
837  END SUBROUTINE write_cell_low
838 
839 END MODULE cell_methods
real(kind=dp) function det_3x3(a)
...
Definition: dumpdcd.F:1091
Handles all functions related to the CELL.
Definition: cell_methods.F:15
recursive subroutine, public read_cell(cell, cell_ref, use_ref_cell, cell_section, check_for_ref, para_env)
...
Definition: cell_methods.F:272
subroutine, public write_cell(cell, subsys_section, tag)
Write the cell parameters to the output unit.
Definition: cell_methods.F:731
subroutine, public read_cell_cif(cif_file_name, cell, para_env)
Reads cell information from CIF file.
Definition: cell_methods.F:618
subroutine, public set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma) using the convention: a parall...
Definition: cell_methods.F:439
subroutine, public init_cell(cell, hmat, periodic)
Initialise/readjust a simulation cell after hmat has been changed.
Definition: cell_methods.F:117
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition: cell_methods.F:85
Handles all functions related to the CELL.
Definition: cell_types.F:15
integer, parameter, public use_perd_xyz
Definition: cell_types.F:42
integer, parameter, public cell_sym_monoclinic
Definition: cell_types.F:29
integer, parameter, public use_perd_y
Definition: cell_types.F:42
integer, parameter, public cell_sym_triclinic
Definition: cell_types.F:29
integer, parameter, public cell_sym_tetragonal_ab
Definition: cell_types.F:29
integer, parameter, public use_perd_xz
Definition: cell_types.F:42
integer, parameter, public cell_sym_rhombohedral
Definition: cell_types.F:29
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition: cell_types.F:559
integer, parameter, public use_perd_x
Definition: cell_types.F:42
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition: cell_types.F:107
integer, parameter, public cell_sym_tetragonal_ac
Definition: cell_types.F:29
integer, parameter, public use_perd_z
Definition: cell_types.F:42
integer, parameter, public use_perd_yz
Definition: cell_types.F:42
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition: cell_types.F:195
integer, parameter, public use_perd_none
Definition: cell_types.F:42
integer, parameter, public cell_sym_hexagonal_gamma_60
Definition: cell_types.F:29
integer, parameter, public cell_sym_orthorhombic
Definition: cell_types.F:29
integer, parameter, public cell_sym_none
Definition: cell_types.F:29
integer, parameter, public cell_sym_hexagonal_gamma_120
Definition: cell_types.F:29
integer, parameter, public cell_sym_monoclinic_gamma_ab
Definition: cell_types.F:29
integer, parameter, public cell_sym_cubic
Definition: cell_types.F:29
integer, parameter, public use_perd_xy
Definition: cell_types.F:42
integer, parameter, public cell_sym_tetragonal_bc
Definition: cell_types.F:29
real(kind=dp) function, public plane_distance(h, k, l, cell)
Calculate the distance between two lattice planes as defined by a triple of Miller indices (hkl).
Definition: cell_types.F:252
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
subroutine, public parser_search_string(parser, string, ignore_case, found, line, begin_line, search_from_begin_of_file)
Search a string pattern in a file defined by its logical unit number "unit". A case sensitive search ...
Utility routines to read data from files. Kept as close as possible to the old parser because.
subroutine, public parser_release(parser)
releases the parser
subroutine, public parser_create(parser, file_name, unit_nr, para_env, end_section_label, separator_chars, comment_char, continuation_char, quote_char, section_char, parse_white_lines, initial_variables, apply_preprocessing)
Start a parser run. Initial variables allow to @SET stuff before opening the file.
unit conversion facility
Definition: cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition: cp_units.F:1150
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition: gamma.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_cell_cif
integer, parameter, public do_cell_xsc
integer, parameter, public do_cell_cp2k
builds the subsystem section of the input
subroutine, public create_cell_section(section, periodic)
creates the cell section
represents an enumeration, i.e. a mapping between integers and strings
character(len=default_string_length) function, public enum_i2c(enum, i)
maps an integer to a string
represents keywords in an input
subroutine, public keyword_get(keyword, names, usage, description, type_of_var, n_var, default_value, lone_keyword_value, repeats, enum, citations)
...
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_unset(section_vals, keyword_name, i_rep_section, i_rep_val)
unsets (removes) the requested value (if it is a keyword repetitions removes the repetition,...
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
recursive type(keyword_type) function, pointer, public section_get_keyword(section, keyword_name)
returns the requested keyword
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
integer, parameter, public default_output_unit
Definition: machine.F:45
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public degree
real(kind=dp), parameter, public sqrt3
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
pure real(kind=dp) function, public angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition: mathlib.F:175
pure real(kind=dp) function, dimension(3, 3), public inv_3x3(a)
Returns the inverse of the 3 x 3 matrix a.
Definition: mathlib.F:516
Interface to the message passing library MPI.