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