1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Handles all functions related to the CELL
8!> \par History
9!>      11.2008 Teodoro Laino [tlaino] - deeply cleaning cell_type from units
10!>      10.2014 Moved many routines to cell_types.F.
11!> \author Matthias KracK (16.01.2002, based on a earlier version of CJM, JGH)
12! **************************************************************************************************
13MODULE cell_methods
14   USE cell_types,                      ONLY: &
15        cell_clone, cell_create, cell_release, cell_sym_none, cell_type, get_cell, init_cell, &
16        set_cell_param, use_perd_none, use_perd_x, use_perd_xy, use_perd_xyz, use_perd_xz, &
17        use_perd_y, use_perd_yz, use_perd_z
18   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
19                                              cp_logger_type
20   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
21                                              cp_print_key_unit_nr
22   USE cp_para_types,                   ONLY: cp_para_env_type
23   USE cp_parser_methods,               ONLY: parser_get_next_line,&
24                                              parser_get_object,&
25                                              parser_search_string
26   USE cp_parser_types,                 ONLY: cp_parser_type,&
27                                              parser_create,&
28                                              parser_release
29   USE cp_units,                        ONLY: cp_unit_from_cp2k,&
30                                              cp_unit_to_cp2k
31   USE input_constants,                 ONLY: do_cell_cif,&
32                                              do_cell_cp2k,&
33                                              do_cell_xsc
34   USE input_cp2k_subsys,               ONLY: create_cell_section
35   USE input_enumeration_types,         ONLY: enum_i2c,&
36                                              enumeration_type
37   USE input_keyword_types,             ONLY: keyword_get,&
38                                              keyword_type
39   USE input_section_types,             ONLY: &
40        section_get_keyword, section_release, section_type, section_vals_get, &
41        section_vals_get_subs_vals, section_vals_type, section_vals_val_get, section_vals_val_set, &
42        section_vals_val_unset
43   USE kinds,                           ONLY: default_path_length,&
44                                              default_string_length,&
45                                              dp
46   USE mathconstants,                   ONLY: degree
47   USE mathlib,                         ONLY: angle
48#include "./base/base_uses.f90"
49
50   IMPLICIT NONE
51
52   PRIVATE
53
54   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_methods'
55
56   ! Public subroutines
57   PUBLIC :: read_cell, read_cell_cif, write_cell
58
59CONTAINS
60
61! **************************************************************************************************
62!> \brief ...
63!> \param cell ...
64!> \param cell_ref ...
65!> \param use_ref_cell ...
66!> \param cell_section ...
67!> \param check_for_ref ...
68!> \param para_env ...
69!> \par History
70!>      03.2005 created [teo]
71!> \author Teodoro Laino
72! **************************************************************************************************
73   RECURSIVE SUBROUTINE read_cell(cell, cell_ref, use_ref_cell, cell_section, &
74                                  check_for_ref, para_env)
75      TYPE(cell_type), POINTER                           :: cell, cell_ref
76      LOGICAL, INTENT(OUT), OPTIONAL                     :: use_ref_cell
77      TYPE(section_vals_type), OPTIONAL, POINTER         :: cell_section
78      LOGICAL, INTENT(IN), OPTIONAL                      :: check_for_ref
79      TYPE(cp_para_env_type), POINTER                    :: para_env
80
81      CHARACTER(len=*), PARAMETER :: routineN = 'read_cell', routineP = moduleN//':'//routineN
82
83      INTEGER                                            :: my_per
84      INTEGER, DIMENSION(:), POINTER                     :: multiple_unit_cell
85      LOGICAL                                            :: cell_read_a, cell_read_abc, cell_read_b, &
86                                                            cell_read_c, cell_read_file, check, &
87                                                            my_check
88      REAL(KIND=dp), DIMENSION(:), POINTER               :: cell_angles, cell_par
89      TYPE(section_vals_type), POINTER                   :: cell_ref_section
90
91      my_check = .TRUE.
92      NULLIFY (cell_ref_section, cell_par, multiple_unit_cell)
93      IF (.NOT. ASSOCIATED(cell)) CALL cell_create(cell)
94      IF (.NOT. ASSOCIATED(cell_ref)) CALL cell_create(cell_ref)
95      IF (PRESENT(check_for_ref)) my_check = check_for_ref
96
97      cell%deth = 0.0_dp
98      cell%orthorhombic = .FALSE.
99      cell%perd(:) = 1
100      cell%symmetry_id = cell_sym_none
101      cell%hmat(:, :) = 0.0_dp
102      cell%h_inv(:, :) = 0.0_dp
103      cell_read_file = .FALSE.
104      cell_read_a = .FALSE.
105      cell_read_b = .FALSE.
106      cell_read_c = .FALSE.
107      ! Trying to read cell info from file
108      CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", explicit=cell_read_file)
109      IF (cell_read_file) CALL read_cell_from_external_file(cell_section, para_env)
110
111      ! Trying to read cell info from the separate A, B, C vectors
112      ! If cell information is provided through file A,B,C contain the file information..
113      ! a print warning is shown on screen..
114      CALL section_vals_val_get(cell_section, "A", explicit=cell_read_a)
115      IF (cell_read_a) THEN
116         CALL section_vals_val_get(cell_section, "A", r_vals=cell_par)
117         cell%hmat(:, 1) = cell_par(:)
118      END IF
119      CALL section_vals_val_get(cell_section, "B", explicit=cell_read_b)
120      IF (cell_read_b) THEN
121         CALL section_vals_val_get(cell_section, "B", r_vals=cell_par)
122         cell%hmat(:, 2) = cell_par(:)
123      END IF
124      CALL section_vals_val_get(cell_section, "C", explicit=cell_read_c)
125      IF (cell_read_c) THEN
126         CALL section_vals_val_get(cell_section, "C", r_vals=cell_par)
127         cell%hmat(:, 3) = cell_par(:)
128      END IF
129      check = ((cell_read_a .EQV. cell_read_b) .AND. (cell_read_b .EQV. cell_read_c))
130      IF (.NOT. check) &
131         CALL cp_warn(__LOCATION__, &
132                      "Cell Information provided through vectors A, B and C. Not all three "// &
133                      "vectors were provided! Cell setup may be incomplete!")
134
135      ! Very last option.. Trying to read cell info from ABC keyword
136      CALL section_vals_val_get(cell_section, "ABC", explicit=cell_read_abc)
137      IF (cell_read_abc) THEN
138         check = (cell_read_a .OR. cell_read_b .OR. cell_read_c)
139         IF (check) &
140            CALL cp_warn(__LOCATION__, &
141                         "Cell Information provided through vectors A, B and C in conjunction with ABC."// &
142                         " The definition of the ABC keyword will override the one provided by A,B and C.")
143         cell%hmat = 0.0_dp
144         CALL section_vals_val_get(cell_section, "ABC", r_vals=cell_par)
145         CALL section_vals_val_get(cell_section, "ALPHA_BETA_GAMMA", r_vals=cell_angles)
146         CALL set_cell_param(cell, cell_par, cell_angles, do_init_cell=.FALSE.)
147      END IF
148
149      ! Multiple unit cell
150      CALL section_vals_val_get(cell_section, "MULTIPLE_UNIT_CELL", i_vals=multiple_unit_cell)
151      IF (ANY(multiple_unit_cell /= 1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell)
152
153      CALL section_vals_val_get(cell_section, "PERIODIC", i_val=my_per)
154      SELECT CASE (my_per)
155      CASE (use_perd_x)
156         cell%perd = (/1, 0, 0/)
157      CASE (use_perd_y)
158         cell%perd = (/0, 1, 0/)
159      CASE (use_perd_z)
160         cell%perd = (/0, 0, 1/)
161      CASE (use_perd_xy)
162         cell%perd = (/1, 1, 0/)
163      CASE (use_perd_xz)
164         cell%perd = (/1, 0, 1/)
165      CASE (use_perd_yz)
166         cell%perd = (/0, 1, 1/)
167      CASE (use_perd_xyz)
168         cell%perd = (/1, 1, 1/)
169      CASE (use_perd_none)
170         cell%perd = (/0, 0, 0/)
171      CASE DEFAULT
172         CPABORT("")
173      END SELECT
174
175      ! Load requested cell symmetry
176      CALL section_vals_val_get(cell_section, "SYMMETRY", i_val=cell%symmetry_id)
177
178      ! Initialize cell
179      CALL init_cell(cell)
180
181      IF (.NOT. my_check) RETURN
182      cell_ref_section => section_vals_get_subs_vals(cell_section, &
183                                                     "CELL_REF")
184      IF (parsed_cp2k_input(cell_ref_section, check_this_section=.TRUE.)) THEN
185         IF (PRESENT(use_ref_cell)) use_ref_cell = .TRUE.
186         CALL read_cell(cell_ref, cell_ref, use_ref_cell, cell_section=cell_ref_section, &
187                        check_for_ref=.FALSE., para_env=para_env)
188      ELSE
189         CALL cell_clone(cell, cell_ref)
190         IF (PRESENT(use_ref_cell)) use_ref_cell = .FALSE.
191      END IF
192
193   END SUBROUTINE read_cell
194
195! **************************************************************************************************
196!> \brief utility function to ease the transition to the new input.
197!>      returns true if the new input was parsed
198!> \param input_file the parsed input file
199!> \param check_this_section ...
200!> \return ...
201!> \author fawzi
202! **************************************************************************************************
203   FUNCTION parsed_cp2k_input(input_file, check_this_section) RESULT(res)
204      TYPE(section_vals_type), POINTER                   :: input_file
205      LOGICAL, INTENT(IN), OPTIONAL                      :: check_this_section
206      LOGICAL                                            :: res
207
208      CHARACTER(len=*), PARAMETER :: routineN = 'parsed_cp2k_input', &
209         routineP = moduleN//':'//routineN
210
211      LOGICAL                                            :: my_check
212      TYPE(section_vals_type), POINTER                   :: glob_section
213
214      my_check = .FALSE.
215      IF (PRESENT(check_this_section)) my_check = check_this_section
216      res = ASSOCIATED(input_file)
217      IF (res) THEN
218         CPASSERT(input_file%ref_count > 0)
219         IF (.NOT. my_check) THEN
220            glob_section => section_vals_get_subs_vals(input_file, "GLOBAL")
221            CALL section_vals_get(glob_section, explicit=res)
222         ELSE
223            CALL section_vals_get(input_file, explicit=res)
224         END IF
225      END IF
226   END FUNCTION parsed_cp2k_input
227
228! **************************************************************************************************
229!> \brief   Setup of the multiple unit_cell
230!> \param cell ...
231!> \param multiple_unit_cell ...
232!> \date    05.2009
233!> \author  Teodoro Laino [tlaino]
234!> \version 1.0
235! **************************************************************************************************
236   SUBROUTINE set_multiple_unit_cell(cell, multiple_unit_cell)
237
238      TYPE(cell_type), POINTER                           :: cell
239      INTEGER, DIMENSION(:), POINTER                     :: multiple_unit_cell
240
241      CHARACTER(len=*), PARAMETER :: routineN = 'set_multiple_unit_cell', &
242         routineP = moduleN//':'//routineN
243
244! Fail is one of the value is set to zero..
245
246      IF (ANY(multiple_unit_cell <= 0)) &
247         CALL cp_abort(__LOCATION__, &
248                       "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "// &
249                       "A value of 0 or negative is meaningless!")
250
251      ! scale abc accordingly user request
252      cell%hmat(:, 1) = cell%hmat(:, 1)*multiple_unit_cell(1)
253      cell%hmat(:, 2) = cell%hmat(:, 2)*multiple_unit_cell(2)
254      cell%hmat(:, 3) = cell%hmat(:, 3)*multiple_unit_cell(3)
255
256   END SUBROUTINE set_multiple_unit_cell
257
258! **************************************************************************************************
259!> \brief   Read cell information from an external file
260!> \param cell_section ...
261!> \param para_env ...
262!> \date    02.2008
263!> \author  Teodoro Laino [tlaino] - University of Zurich
264!> \version 1.0
265! **************************************************************************************************
266   SUBROUTINE read_cell_from_external_file(cell_section, para_env)
267
268      TYPE(section_vals_type), POINTER                   :: cell_section
269      TYPE(cp_para_env_type), POINTER                    :: para_env
270
271      CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_from_external_file', &
272         routineP = moduleN//':'//routineN
273
274      CHARACTER(LEN=default_path_length)                 :: cell_file_name
275      INTEGER                                            :: i, idum, j, my_format, n_rep
276      LOGICAL                                            :: explicit, my_end
277      REAL(KIND=dp)                                      :: xdum
278      REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
279      REAL(KIND=dp), DIMENSION(:), POINTER               :: cell_par
280      TYPE(cell_type), POINTER                           :: cell
281      TYPE(cp_parser_type), POINTER                      :: parser
282
283      NULLIFY (parser)
284      CALL section_vals_val_get(cell_section, "CELL_FILE_NAME", c_val=cell_file_name)
285      CALL section_vals_val_get(cell_section, "CELL_FILE_FORMAT", i_val=my_format)
286      SELECT CASE (my_format)
287      CASE (do_cell_cp2k)
288         CALL parser_create(parser, cell_file_name, para_env=para_env)
289         CALL parser_get_next_line(parser, 1)
290         my_end = .FALSE.
291         DO WHILE (.NOT. my_end)
292            READ (parser%input_line, *) idum, xdum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
293            CALL parser_get_next_line(parser, 1, at_end=my_end)
294         END DO
295         CALL parser_release(parser)
296         ! Conver to CP2K units
297         DO i = 1, 3
298            DO j = 1, 3
299               hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
300            END DO
301         END DO
302      CASE (do_cell_xsc)
303         CALL parser_create(parser, cell_file_name, para_env=para_env)
304         CALL parser_get_next_line(parser, 1)
305         READ (parser%input_line, *) idum, hmat(:, 1), hmat(:, 2), hmat(:, 3)
306         CALL parser_release(parser)
307         ! Conver to CP2K units
308         DO i = 1, 3
309            DO j = 1, 3
310               hmat(j, i) = cp_unit_to_cp2k(hmat(j, i), "angstrom")
311            END DO
312         END DO
313      CASE (do_cell_cif)
314         NULLIFY (cell)
315         CALL cell_create(cell)
316         CALL read_cell_cif(cell_file_name, cell, para_env)
317         hmat = cell%hmat
318         CALL cell_release(cell)
319      END SELECT
320      CALL section_vals_val_unset(cell_section, "CELL_FILE_NAME")
321      CALL section_vals_val_unset(cell_section, "CELL_FILE_FORMAT")
322      ! Check if the cell was already defined
323      explicit = .FALSE.
324      CALL section_vals_val_get(cell_section, "A", n_rep_val=n_rep)
325      explicit = explicit .OR. (n_rep == 1)
326      CALL section_vals_val_get(cell_section, "B", n_rep_val=n_rep)
327      explicit = explicit .OR. (n_rep == 1)
328      CALL section_vals_val_get(cell_section, "C", n_rep_val=n_rep)
329      explicit = explicit .OR. (n_rep == 1)
330      CALL section_vals_val_get(cell_section, "ABC", n_rep_val=n_rep)
331      explicit = explicit .OR. (n_rep == 1)
332      ! Possibly print a warning
333      IF (explicit) &
334         CALL cp_warn(__LOCATION__, &
335                      "Cell specification (A,B,C or ABC) provided together with the external "// &
336                      "cell setup! Ignoring (A,B,C or ABC) and proceeding with info read from the "// &
337                      "external file! ")
338      ! Copy cell information in the A, B, C fields..(we may need them later on..)
339      ALLOCATE (cell_par(3))
340      cell_par = hmat(:, 1)
341      CALL section_vals_val_set(cell_section, "A", r_vals_ptr=cell_par)
342      ALLOCATE (cell_par(3))
343      cell_par = hmat(:, 2)
344      CALL section_vals_val_set(cell_section, "B", r_vals_ptr=cell_par)
345      ALLOCATE (cell_par(3))
346      cell_par = hmat(:, 3)
347      CALL section_vals_val_set(cell_section, "C", r_vals_ptr=cell_par)
348      ! Unset possible keywords
349      CALL section_vals_val_unset(cell_section, "ABC")
350      CALL section_vals_val_unset(cell_section, "ALPHA_BETA_GAMMA")
351
352   END SUBROUTINE read_cell_from_external_file
353
354! **************************************************************************************************
355!> \brief  Reads cell information from CIF file
356!> \param cif_file_name ...
357!> \param cell ...
358!> \param para_env ...
359!> \date   12.2008
360!> \par    Format Information implemented:
361!>            _cell_length_a
362!>            _cell_length_b
363!>            _cell_length_c
364!>            _cell_angle_alpha
365!>            _cell_angle_beta
366!>            _cell_angle_gamma
367!>
368!> \author Teodoro Laino [tlaino]
369!>         moved from topology_cif (1/2019 JHU)
370! **************************************************************************************************
371   SUBROUTINE read_cell_cif(cif_file_name, cell, para_env)
372      CHARACTER(len=*)                                   :: cif_file_name
373      TYPE(cell_type), POINTER                           :: cell
374      TYPE(cp_para_env_type), POINTER                    :: para_env
375
376      CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_cif', routineP = moduleN//':'//routineN
377
378      INTEGER                                            :: handle
379      INTEGER, DIMENSION(3)                              :: periodic
380      LOGICAL                                            :: found
381      REAL(KIND=dp), DIMENSION(3)                        :: cell_angles, cell_lengths
382      TYPE(cp_parser_type), POINTER                      :: parser
383
384      CALL timeset(routineN, handle)
385
386      NULLIFY (parser)
387      CALL parser_create(parser, cif_file_name, &
388                         para_env=para_env, apply_preprocessing=.FALSE.)
389
390      ! Parsing cell infos
391      periodic = 1
392      ! Check for   _cell_length_a
393      CALL parser_search_string(parser, "_cell_length_a", ignore_case=.FALSE., found=found, &
394                                begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
395      IF (.NOT. found) &
396         CPABORT("The field (_cell_length_a) was not found in CIF file! ")
397      CALL cif_get_real(parser, cell_lengths(1))
398      cell_lengths(1) = cp_unit_to_cp2k(cell_lengths(1), "angstrom")
399
400      ! Check for   _cell_length_b
401      CALL parser_search_string(parser, "_cell_length_b", ignore_case=.FALSE., found=found, &
402                                begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
403      IF (.NOT. found) &
404         CPABORT("The field (_cell_length_b) was not found in CIF file! ")
405      CALL cif_get_real(parser, cell_lengths(2))
406      cell_lengths(2) = cp_unit_to_cp2k(cell_lengths(2), "angstrom")
407
408      ! Check for   _cell_length_c
409      CALL parser_search_string(parser, "_cell_length_c", ignore_case=.FALSE., found=found, &
410                                begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
411      IF (.NOT. found) &
412         CPABORT("The field (_cell_length_c) was not found in CIF file! ")
413      CALL cif_get_real(parser, cell_lengths(3))
414      cell_lengths(3) = cp_unit_to_cp2k(cell_lengths(3), "angstrom")
415
416      ! Check for   _cell_angle_alpha
417      CALL parser_search_string(parser, "_cell_angle_alpha", ignore_case=.FALSE., found=found, &
418                                begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
419      IF (.NOT. found) &
420         CPABORT("The field (_cell_angle_alpha) was not found in CIF file! ")
421      CALL cif_get_real(parser, cell_angles(1))
422      cell_angles(1) = cp_unit_to_cp2k(cell_angles(1), "deg")
423
424      ! Check for   _cell_angle_beta
425      CALL parser_search_string(parser, "_cell_angle_beta", ignore_case=.FALSE., found=found, &
426                                begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
427      IF (.NOT. found) &
428         CPABORT("The field (_cell_angle_beta) was not found in CIF file! ")
429      CALL cif_get_real(parser, cell_angles(2))
430      cell_angles(2) = cp_unit_to_cp2k(cell_angles(2), "deg")
431
432      ! Check for   _cell_angle_gamma
433      CALL parser_search_string(parser, "_cell_angle_gamma", ignore_case=.FALSE., found=found, &
434                                begin_line=.FALSE., search_from_begin_of_file=.TRUE.)
435      IF (.NOT. found) &
436         CPABORT("The field (_cell_angle_gamma) was not found in CIF file! ")
437      CALL cif_get_real(parser, cell_angles(3))
438      cell_angles(3) = cp_unit_to_cp2k(cell_angles(3), "deg")
439
440      ! Create cell
441      CALL set_cell_param(cell, cell_lengths, cell_angles, periodic=periodic, &
442                          do_init_cell=.TRUE.)
443
444      CALL parser_release(parser)
445
446      CALL timestop(handle)
447
448   END SUBROUTINE read_cell_cif
449
450! **************************************************************************************************
451!> \brief  Reads REAL from the CIF file.. This wrapper is needed in order to
452!>         treat properly the accuracy specified in the CIF file, i.e. 3.45(6)
453!> \param parser ...
454!> \param r ...
455!> \date   12.2008
456!> \author Teodoro Laino [tlaino]
457! **************************************************************************************************
458   SUBROUTINE cif_get_real(parser, r)
459      TYPE(cp_parser_type), POINTER                      :: parser
460      REAL(KIND=dp), INTENT(OUT)                         :: r
461
462      CHARACTER(len=*), PARAMETER :: routineN = 'cif_get_real', routineP = moduleN//':'//routineN
463
464      CHARACTER(LEN=default_string_length)               :: s_tag
465      INTEGER                                            :: iln
466
467      CALL parser_get_object(parser, s_tag)
468      iln = LEN_TRIM(s_tag)
469      IF (INDEX(s_tag, "(") /= 0) iln = INDEX(s_tag, "(") - 1
470      READ (s_tag(1:iln), *) r
471   END SUBROUTINE cif_get_real
472
473! **************************************************************************************************
474!> \brief   Write the cell parameters to the output unit.
475!> \param cell ...
476!> \param subsys_section ...
477!> \param cell_ref ...
478!> \param label ...
479!> \date    02.06.2000
480!> \par     History
481!>    - 11.2008 Teodoro Laino [tlaino] - rewrite and enabling user driven units
482!> \author  Matthias Krack
483!> \version 1.0
484! **************************************************************************************************
485   RECURSIVE SUBROUTINE write_cell(cell, subsys_section, cell_ref, label)
486
487      TYPE(cell_type), POINTER                           :: cell
488      TYPE(section_vals_type), POINTER                   :: subsys_section
489      TYPE(cell_type), OPTIONAL, POINTER                 :: cell_ref
490      CHARACTER(LEN=*), INTENT(IN), OPTIONAL             :: label
491
492      CHARACTER(len=*), PARAMETER :: routineN = 'write_cell', routineP = moduleN//':'//routineN
493
494      CHARACTER(LEN=default_string_length)               :: my_label, unit_str
495      INTEGER                                            :: output_unit
496      REAL(KIND=dp)                                      :: alpha, beta, gamma, val
497      REAL(KIND=dp), DIMENSION(3)                        :: abc
498      TYPE(cp_logger_type), POINTER                      :: logger
499      TYPE(enumeration_type), POINTER                    :: enum
500      TYPE(keyword_type), POINTER                        :: keyword
501      TYPE(section_type), POINTER                        :: section
502
503      NULLIFY (enum)
504      NULLIFY (keyword)
505      NULLIFY (logger)
506      NULLIFY (section)
507      logger => cp_get_default_logger()
508      my_label = "CELL|"
509      IF (PRESENT(label)) my_label = TRIM(label)
510      output_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%CELL", &
511                                         extension=".Log")
512      CALL section_vals_val_get(subsys_section, "PRINT%CELL%UNIT", c_val=unit_str)
513      IF (output_unit > 0) THEN
514         CALL get_cell(cell=cell, abc=abc, alpha=alpha, beta=beta, gamma=gamma)
515         WRITE (UNIT=output_unit, FMT='( )')
516         val = cp_unit_from_cp2k(cell%deth, TRIM(unit_str)//"^3")
517         WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.3)") &
518            TRIM(my_label)//" Volume ["//TRIM(unit_str)//"^3]:", val
519         val = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
520         WRITE (UNIT=output_unit, FMT="(T2,A,T30,3F10.3,4X,A6,F11.3)") &
521            TRIM(my_label)//" Vector a ["//TRIM(unit_str)//"]:", cell%hmat(:, 1)*val, &
522            "|a| = ", abc(1)*val, &
523            TRIM(my_label)//" Vector b ["//TRIM(unit_str)//"]:", cell%hmat(:, 2)*val, &
524            "|b| = ", abc(2)*val, &
525            TRIM(my_label)//" Vector c ["//TRIM(unit_str)//"]:", cell%hmat(:, 3)*val, &
526            "|c| = ", abc(3)*val
527         WRITE (UNIT=output_unit, FMT="(T2,A,T70,F11.3)") &
528            TRIM(my_label)//" Angle (b,c), alpha [degree]: ", alpha, &
529            TRIM(my_label)//" Angle (a,c), beta  [degree]: ", beta, &
530            TRIM(my_label)//" Angle (a,b), gamma [degree]: ", gamma
531         IF (cell%symmetry_id /= cell_sym_none) THEN
532            CALL create_cell_section(section)
533            keyword => section_get_keyword(section, "SYMMETRY")
534            CALL keyword_get(keyword, enum=enum)
535            WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
536               TRIM(my_label)//" Requested initial symmetry: ", &
537               ADJUSTR(TRIM(enum_i2c(enum, cell%symmetry_id)))
538            CALL section_release(section)
539         END IF
540         IF (cell%orthorhombic) THEN
541            WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
542               TRIM(my_label)//" Numerically orthorhombic: ", "YES"
543         ELSE
544            WRITE (UNIT=output_unit, FMT="(T2,A,T78,A3)") &
545               TRIM(my_label)//" Numerically orthorhombic: ", " NO"
546         END IF
547      END IF
548      CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
549                                        "PRINT%CELL")
550
551      IF (PRESENT(cell_ref)) THEN
552         CALL write_cell(cell_ref, subsys_section, label="CELL_REF|")
553      END IF
554
555   END SUBROUTINE write_cell
556
557END MODULE cell_methods
558