(git:07a6c39)
Loading...
Searching...
No Matches
cell_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Handles all functions related to the CELL
10!> \par History
11!> 11.2008 Teodoro Laino [tlaino] - deeply cleaning cell_type from units
12!> 10.2014 Moved many routines from cell_types.F here.
13!> \author Matthias KracK (16.01.2002, based on a earlier version of CJM, JGH)
14! **************************************************************************************************
15MODULE cell_types
16 USE cp_units, ONLY: cp_unit_to_cp2k
17 USE kinds, ONLY: dp
18 USE mathconstants, ONLY: degree
19 USE mathlib, ONLY: angle
20#include "../base/base_uses.f90"
21
22 IMPLICIT NONE
23
24 PRIVATE
25
26 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_types'
27
28 ! Impose cell symmetry
29 INTEGER, PARAMETER, PUBLIC :: cell_sym_none = 0, &
30 cell_sym_triclinic = 1, &
31 cell_sym_monoclinic = 2, &
32 cell_sym_monoclinic_gamma_ab = 3, &
33 cell_sym_orthorhombic = 4, &
34 cell_sym_tetragonal_ab = 5, &
35 cell_sym_tetragonal_ac = 6, &
36 cell_sym_tetragonal_bc = 7, &
37 cell_sym_rhombohedral = 8, &
38 cell_sym_hexagonal_gamma_60 = 9, &
39 cell_sym_hexagonal_gamma_120 = 10, &
40 cell_sym_cubic = 11
41
42 INTEGER, PARAMETER, PUBLIC :: use_perd_none = 0, &
43 use_perd_x = 1, &
44 use_perd_y = 2, &
45 use_perd_z = 3, &
46 use_perd_xy = 4, &
47 use_perd_xz = 5, &
48 use_perd_yz = 6, &
49 use_perd_xyz = 7
50
51 CHARACTER(LEN=3), DIMENSION(7), &
52 PARAMETER, PUBLIC :: periodicity_string = [" X", " Y", " Z", &
53 " XY", " XZ", " YZ", &
54 "XYZ"]
55
56! **************************************************************************************************
57!> \brief Type defining parameters related to the simulation cell
58!> \version 1.0
59! **************************************************************************************************
60 TYPE cell_type
61 CHARACTER(LEN=12) :: tag = "CELL"
62 INTEGER :: ref_count = -1, &
63 symmetry_id = use_perd_none
64 LOGICAL :: orthorhombic = .false. ! actually means a diagonal hmat
65 LOGICAL :: input_cell_canonicalized = .false.
66 REAL(kind=dp) :: deth = 0.0_dp
67 INTEGER, DIMENSION(3) :: perd = -1
68 REAL(kind=dp), DIMENSION(3, 3) :: hmat = 0.0_dp, &
69 h_inv = 0.0_dp, &
70 input_hmat = 0.0_dp, &
71 input_to_canonical = 0.0_dp, &
72 input_recip_to_canonical = 0.0_dp
73 END TYPE cell_type
74
75 TYPE cell_p_type
76 TYPE(cell_type), POINTER :: cell => null()
77 END TYPE cell_p_type
78
79 ! Public data types
80 PUBLIC :: cell_type, &
81 cell_p_type
82
83 ! Public subroutines
84 PUBLIC :: cell_clone, &
85 cell_copy, &
86 cell_transform_input_cartesian, &
87 cell_transform_input_reciprocal, &
88 cell_release, &
89 cell_retain, &
90 get_cell, &
91 parse_cell_line
92
93#if defined (__PLUMED2)
94 PUBLIC :: pbc_cp2k_plumed_getset_cell
95#endif
96
97 ! Public functions
98 PUBLIC :: plane_distance, &
99 pbc, &
100 pbc_stable, &
101 real_to_scaled, &
102 scaled_to_real
103
104 INTERFACE pbc
105 MODULE PROCEDURE pbc1, pbc2, pbc3, pbc4
106 END INTERFACE
107
108CONTAINS
109
110! **************************************************************************************************
111!> \brief Select a stable periodic image index close to half-cell boundaries.
112!> \param s Scaled coordinate
113!> \return image_shift ...
114! **************************************************************************************************
115 PURE ELEMENTAL FUNCTION pbc_image_shift(s) RESULT(image_shift)
116
117 REAL(KIND=dp), INTENT(IN) :: s
118 REAL(KIND=dp) :: image_shift
119
120 REAL(KIND=dp) :: half_boundary, tolerance
121
122 image_shift = anint(s)
123 half_boundary = anint(s - 0.5_dp) + 0.5_dp
124 tolerance = min(1.0e-8_dp, 64.0_dp*epsilon(1.0_dp)*max(1.0_dp, abs(s)))
125 IF (abs(s - half_boundary) <= tolerance) THEN
126 ! Select the lower side for every lattice-equivalent half-cell boundary.
127 image_shift = half_boundary + 0.5_dp
128 END IF
129
130 END FUNCTION pbc_image_shift
131
132! **************************************************************************************************
133!> \brief Clone cell variable
134!> \param cell_in Cell variable to be clone
135!> \param cell_out Cloned cell variable
136!> \param tag Optional new tag for cloned cell variable
137!> \par History
138!> - Optional tag added (17.05.2023, MK)
139! **************************************************************************************************
140 SUBROUTINE cell_clone(cell_in, cell_out, tag)
141
142 TYPE(cell_type), POINTER :: cell_in, cell_out
143 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
144
145 cell_out = cell_in
146 cell_out%ref_count = 1
147 IF (PRESENT(tag)) cell_out%tag = tag
148
149 END SUBROUTINE cell_clone
150
151! **************************************************************************************************
152!> \brief Copy cell variable
153!> \param cell_in Cell variable to be copied
154!> \param cell_out Copy of cell variable
155!> \param tag Optional new tag
156!> \par History
157!> - Optional tag added (17.05.2023, MK)
158! **************************************************************************************************
159 SUBROUTINE cell_copy(cell_in, cell_out, tag)
160
161 TYPE(cell_type), POINTER :: cell_in, cell_out
162 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
163
164 cell_out%deth = cell_in%deth
165 cell_out%perd = cell_in%perd
166 cell_out%hmat = cell_in%hmat
167 cell_out%h_inv = cell_in%h_inv
168 cell_out%input_cell_canonicalized = cell_in%input_cell_canonicalized
169 cell_out%input_hmat = cell_in%input_hmat
170 cell_out%input_to_canonical = cell_in%input_to_canonical
171 cell_out%input_recip_to_canonical = cell_in%input_recip_to_canonical
172 cell_out%orthorhombic = cell_in%orthorhombic
173 cell_out%symmetry_id = cell_in%symmetry_id
174 IF (PRESENT(tag)) THEN
175 cell_out%tag = tag
176 ELSE
177 cell_out%tag = cell_in%tag
178 END IF
179
180 END SUBROUTINE cell_copy
181
182! **************************************************************************************************
183!> \brief Read cell info from a line (parsed from a file)
184!> \param input_line ...
185!> \param cell_itimes ...
186!> \param cell_time ...
187!> \param h ...
188!> \param vol ...
189!> \date 19.02.2008
190!> \author Teodoro Laino [tlaino] - University of Zurich
191!> \version 1.0
192! **************************************************************************************************
193 SUBROUTINE parse_cell_line(input_line, cell_itimes, cell_time, h, vol)
194
195 CHARACTER(LEN=*), INTENT(IN) :: input_line
196 INTEGER, INTENT(OUT) :: cell_itimes
197 REAL(kind=dp), INTENT(OUT) :: cell_time
198 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT) :: h
199 REAL(kind=dp), INTENT(OUT) :: vol
200
201 INTEGER :: i, j
202
203 READ (input_line, *) cell_itimes, cell_time, &
204 h(1, 1), h(2, 1), h(3, 1), h(1, 2), h(2, 2), h(3, 2), h(1, 3), h(2, 3), h(3, 3), vol
205 DO i = 1, 3
206 DO j = 1, 3
207 h(j, i) = cp_unit_to_cp2k(h(j, i), "angstrom")
208 END DO
209 END DO
210
211 END SUBROUTINE parse_cell_line
212
213! **************************************************************************************************
214!> \brief Get informations about a simulation cell.
215!> \param cell ...
216!> \param alpha ...
217!> \param beta ...
218!> \param gamma ...
219!> \param deth ...
220!> \param orthorhombic ...
221!> \param abc ...
222!> \param periodic ...
223!> \param h ...
224!> \param h_inv ...
225!> \param symmetry_id ...
226!> \param tag ...
227!> \date 16.01.2002
228!> \author Matthias Krack
229!> \version 1.0
230! **************************************************************************************************
231 SUBROUTINE get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, &
232 h, h_inv, symmetry_id, tag)
233
234 TYPE(cell_type), POINTER :: cell
235 REAL(kind=dp), INTENT(OUT), OPTIONAL :: alpha, beta, gamma, deth
236 LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic
237 REAL(kind=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: abc
238 INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: periodic
239 REAL(kind=dp), DIMENSION(3, 3), INTENT(OUT), &
240 OPTIONAL :: h, h_inv
241 INTEGER, INTENT(OUT), OPTIONAL :: symmetry_id
242 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: tag
243
244 cpassert(ASSOCIATED(cell))
245
246 IF (PRESENT(deth)) deth = cell%deth ! the volume
247 IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic
248 IF (PRESENT(periodic)) periodic(:) = cell%perd(:)
249 IF (PRESENT(h)) h(:, :) = cell%hmat(:, :)
250 IF (PRESENT(h_inv)) h_inv(:, :) = cell%h_inv(:, :)
251
252 ! Calculate the lengths of the cell vectors a, b, and c
253 IF (PRESENT(abc)) THEN
254 abc(1) = sqrt(cell%hmat(1, 1)*cell%hmat(1, 1) + &
255 cell%hmat(2, 1)*cell%hmat(2, 1) + &
256 cell%hmat(3, 1)*cell%hmat(3, 1))
257 abc(2) = sqrt(cell%hmat(1, 2)*cell%hmat(1, 2) + &
258 cell%hmat(2, 2)*cell%hmat(2, 2) + &
259 cell%hmat(3, 2)*cell%hmat(3, 2))
260 abc(3) = sqrt(cell%hmat(1, 3)*cell%hmat(1, 3) + &
261 cell%hmat(2, 3)*cell%hmat(2, 3) + &
262 cell%hmat(3, 3)*cell%hmat(3, 3))
263 END IF
264
265 ! Angles between the cell vectors a, b, and c
266 ! alpha = <(b,c)
267 IF (PRESENT(alpha)) alpha = angle(cell%hmat(:, 2), cell%hmat(:, 3))*degree
268 ! beta = <(a,c)
269 IF (PRESENT(beta)) beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))*degree
270 ! gamma = <(a,b)
271 IF (PRESENT(gamma)) gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))*degree
272 IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id
273 IF (PRESENT(tag)) tag = cell%tag
274
275 END SUBROUTINE get_cell
276
277! **************************************************************************************************
278!> \brief Transform a Cartesian real-space vector from the user input cell frame
279!> into CP2K's canonical internal cell frame.
280!> \param cell ...
281!> \param vector ...
282! **************************************************************************************************
283 SUBROUTINE cell_transform_input_cartesian(cell, vector)
284
285 TYPE(cell_type), POINTER :: cell
286 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: vector
287
288 cpassert(ASSOCIATED(cell))
289
290 IF (cell%input_cell_canonicalized) vector = matmul(cell%input_to_canonical, vector)
291
292 END SUBROUTINE cell_transform_input_cartesian
293
294! **************************************************************************************************
295!> \brief Transform a Cartesian reciprocal-space vector from the user input cell
296!> frame into CP2K's canonical internal cell frame.
297!> \param cell ...
298!> \param vector ...
299! **************************************************************************************************
300 SUBROUTINE cell_transform_input_reciprocal(cell, vector)
301
302 TYPE(cell_type), POINTER :: cell
303 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: vector
304
305 cpassert(ASSOCIATED(cell))
306
307 IF (cell%input_cell_canonicalized) vector = matmul(cell%input_recip_to_canonical, vector)
308
309 END SUBROUTINE cell_transform_input_reciprocal
310
311! **************************************************************************************************
312!> \brief Calculate the distance between two lattice planes as defined by
313!> a triple of Miller indices (hkl).
314!> \param h ...
315!> \param k ...
316!> \param l ...
317!> \param cell ...
318!> \return ...
319!> \date 18.11.2004
320!> \author Matthias Krack
321!> \version 1.0
322! **************************************************************************************************
323 FUNCTION plane_distance(h, k, l, cell) RESULT(distance)
324
325 INTEGER, INTENT(IN) :: h, k, l
326 TYPE(cell_type), POINTER :: cell
327 REAL(kind=dp) :: distance
328
329 REAL(kind=dp) :: a, alpha, b, beta, c, cosa, cosb, cosg, &
330 d, gamma, x, y, z
331 REAL(kind=dp), DIMENSION(3) :: abc
332
333 x = real(h, kind=dp)
334 y = real(k, kind=dp)
335 z = real(l, kind=dp)
336
337 CALL get_cell(cell=cell, abc=abc)
338
339 a = abc(1)
340 b = abc(2)
341 c = abc(3)
342
343 IF (cell%orthorhombic) THEN
344
345 d = (x/a)**2 + (y/b)**2 + (z/c)**2
346
347 ELSE
348
349 CALL get_cell(cell=cell, &
350 alpha=alpha, &
351 beta=beta, &
352 gamma=gamma)
353
354 alpha = alpha/degree
355 beta = beta/degree
356 gamma = gamma/degree
357
358 cosa = cos(alpha)
359 cosb = cos(beta)
360 cosg = cos(gamma)
361
362 d = ((x*b*c*sin(alpha))**2 + &
363 (y*c*a*sin(beta))**2 + &
364 (z*a*b*sin(gamma))**2 + &
365 2.0_dp*a*b*c*(x*y*c*(cosa*cosb - cosg) + &
366 z*x*b*(cosg*cosa - cosb) + &
367 y*z*a*(cosb*cosg - cosa)))/ &
368 ((a*b*c)**2*(1.0_dp - cosa**2 - cosb**2 - cosg**2 + &
369 2.0_dp*cosa*cosb*cosg))
370
371 END IF
372
373 distance = 1.0_dp/sqrt(d)
374
375 END FUNCTION plane_distance
376
377! **************************************************************************************************
378!> \brief Apply the periodic boundary conditions defined by a simulation
379!> cell to a position vector r.
380!> \param r ...
381!> \param cell ...
382!> \return ...
383!> \date 16.01.2002
384!> \author Matthias Krack
385!> \version 1.0
386! **************************************************************************************************
387 FUNCTION pbc1(r, cell) RESULT(r_pbc)
388
389 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: r
390 TYPE(cell_type), POINTER :: cell
391 REAL(kind=dp), DIMENSION(3) :: r_pbc
392
393 REAL(kind=dp), DIMENSION(3) :: s
394
395 cpassert(ASSOCIATED(cell))
396
397 IF (cell%orthorhombic) THEN
398 r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*anint(cell%h_inv(1, 1)*r(1))
399 r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*anint(cell%h_inv(2, 2)*r(2))
400 r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*anint(cell%h_inv(3, 3)*r(3))
401 ELSE
402 s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
403 s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
404 s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
405 s(1) = s(1) - cell%perd(1)*anint(s(1))
406 s(2) = s(2) - cell%perd(2)*anint(s(2))
407 s(3) = s(3) - cell%perd(3)*anint(s(3))
408 r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
409 r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
410 r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
411 END IF
412
413 END FUNCTION pbc1
414
415! **************************************************************************************************
416!> \brief Apply a stable periodic-image convention for k-point Bloch gauges.
417!> \param r ...
418!> \param cell ...
419!> \return r_pbc ...
420! **************************************************************************************************
421 FUNCTION pbc_stable(r, cell) RESULT(r_pbc)
422
423 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: r
424 TYPE(cell_type), POINTER :: cell
425 REAL(kind=dp), DIMENSION(3) :: r_pbc
426
427 REAL(kind=dp), DIMENSION(3) :: s
428
429 cpassert(ASSOCIATED(cell))
430
431 IF (cell%orthorhombic) THEN
432 r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
433 pbc_image_shift(cell%h_inv(1, 1)*r(1))
434 r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
435 pbc_image_shift(cell%h_inv(2, 2)*r(2))
436 r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
437 pbc_image_shift(cell%h_inv(3, 3)*r(3))
438 ELSE
439 s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
440 s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
441 s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
442 s(1) = s(1) - cell%perd(1)*pbc_image_shift(s(1))
443 s(2) = s(2) - cell%perd(2)*pbc_image_shift(s(2))
444 s(3) = s(3) - cell%perd(3)*pbc_image_shift(s(3))
445 r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
446 r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
447 r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
448 END IF
449
450 END FUNCTION pbc_stable
451
452! **************************************************************************************************
453!> \brief Apply the periodic boundary conditions defined by a simulation
454!> cell to a position vector r subtracting nl from the periodic images
455!> \param r ...
456!> \param cell ...
457!> \param nl ...
458!> \return ...
459!> \date 16.01.2002
460!> \author Matthias Krack
461!> \version 1.0
462! **************************************************************************************************
463 FUNCTION pbc2(r, cell, nl) RESULT(r_pbc)
464
465 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: r
466 TYPE(cell_type), POINTER :: cell
467 INTEGER, DIMENSION(3), INTENT(IN) :: nl
468 REAL(kind=dp), DIMENSION(3) :: r_pbc
469
470 REAL(kind=dp), DIMENSION(3) :: s
471
472 cpassert(ASSOCIATED(cell))
473
474 IF (cell%orthorhombic) THEN
475 r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
476 REAL(nint(cell%h_inv(1, 1)*r(1)) - nl(1), dp)
477 r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
478 REAL(nint(cell%h_inv(2, 2)*r(2)) - nl(2), dp)
479 r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
480 REAL(nint(cell%h_inv(3, 3)*r(3)) - nl(3), dp)
481 ELSE
482 s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
483 s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
484 s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
485 s(1) = s(1) - cell%perd(1)*real(nint(s(1)) - nl(1), dp)
486 s(2) = s(2) - cell%perd(2)*real(nint(s(2)) - nl(2), dp)
487 s(3) = s(3) - cell%perd(3)*real(nint(s(3)) - nl(3), dp)
488 r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
489 r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
490 r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
491 END IF
492
493 END FUNCTION pbc2
494
495! **************************************************************************************************
496!> \brief Apply the periodic boundary conditions defined by the simulation
497!> cell cell to the vector pointing from atom a to atom b.
498!> \param ra ...
499!> \param rb ...
500!> \param cell ...
501!> \return ...
502!> \date 11.03.2004
503!> \author Matthias Krack
504!> \version 1.0
505! **************************************************************************************************
506 FUNCTION pbc3(ra, rb, cell) RESULT(rab_pbc)
507
508 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: ra, rb
509 TYPE(cell_type), POINTER :: cell
510 REAL(kind=dp), DIMENSION(3) :: rab_pbc
511
512 INTEGER :: icell, jcell, kcell
513 INTEGER, DIMENSION(3) :: periodic
514 REAL(kind=dp) :: rab2, rab2_pbc
515 REAL(kind=dp), DIMENSION(3) :: r, ra_pbc, rab, rb_image, rb_pbc, s2r
516
517 CALL get_cell(cell=cell, periodic=periodic)
518
519 ra_pbc(:) = pbc(ra(:), cell)
520 rb_pbc(:) = pbc(rb(:), cell)
521
522 rab2_pbc = huge(1.0_dp)
523
524 DO icell = -periodic(1), periodic(1)
525 DO jcell = -periodic(2), periodic(2)
526 DO kcell = -periodic(3), periodic(3)
527 r = real([icell, jcell, kcell], dp)
528 CALL scaled_to_real(s2r, r, cell)
529 rb_image(:) = rb_pbc(:) + s2r
530 rab(:) = rb_image(:) - ra_pbc(:)
531 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
532 IF (rab2 < rab2_pbc) THEN
533 rab2_pbc = rab2
534 rab_pbc(:) = rab(:)
535 END IF
536 END DO
537 END DO
538 END DO
539
540 END FUNCTION pbc3
541
542 !if positive_range == true, r(i) (or s(i)) in range [0, hmat(i,i)],
543 !else, r(i) (s(i)) in range [-hmat(i,i)/2, hmat(i,i)/2]
544! **************************************************************************************************
545!> \brief ...
546!> \param r ...
547!> \param cell ...
548!> \param positive_range ...
549!> \return ...
550! **************************************************************************************************
551 FUNCTION pbc4(r, cell, positive_range) RESULT(r_pbc)
552
553 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: r
554 TYPE(cell_type), POINTER :: cell
555 LOGICAL :: positive_range
556 REAL(kind=dp), DIMENSION(3) :: r_pbc
557
558 REAL(kind=dp), DIMENSION(3) :: s
559
560 cpassert(ASSOCIATED(cell))
561
562 IF (positive_range) THEN
563 IF (cell%orthorhombic) THEN
564 r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*floor(cell%h_inv(1, 1)*r(1))
565 r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*floor(cell%h_inv(2, 2)*r(2))
566 r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*floor(cell%h_inv(3, 3)*r(3))
567 ELSE
568 s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
569 s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
570 s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
571 s(1) = s(1) - cell%perd(1)*floor(s(1))
572 s(2) = s(2) - cell%perd(2)*floor(s(2))
573 s(3) = s(3) - cell%perd(3)*floor(s(3))
574 r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
575 r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
576 r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
577 END IF
578 ELSE
579 r_pbc = pbc1(r, cell)
580 END IF
581
582 END FUNCTION pbc4
583
584! **************************************************************************************************
585!> \brief Transform real to scaled cell coordinates.
586!> s=h_inv*r
587!> \param s ...
588!> \param r ...
589!> \param cell ...
590!> \date 16.01.2002
591!> \author Matthias Krack
592!> \version 1.0
593! **************************************************************************************************
594 SUBROUTINE real_to_scaled(s, r, cell)
595
596 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: s
597 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: r
598 TYPE(cell_type), POINTER :: cell
599
600 cpassert(ASSOCIATED(cell))
601
602 IF (cell%orthorhombic) THEN
603 s(1) = cell%h_inv(1, 1)*r(1)
604 s(2) = cell%h_inv(2, 2)*r(2)
605 s(3) = cell%h_inv(3, 3)*r(3)
606 ELSE
607 s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
608 s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
609 s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
610 END IF
611
612 END SUBROUTINE real_to_scaled
613
614! **************************************************************************************************
615!> \brief Transform scaled cell coordinates real coordinates.
616!> r=h*s
617!> \param r ...
618!> \param s ...
619!> \param cell ...
620!> \date 16.01.2002
621!> \author Matthias Krack
622!> \version 1.0
623! **************************************************************************************************
624 SUBROUTINE scaled_to_real(r, s, cell)
625
626 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: r
627 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: s
628 TYPE(cell_type), POINTER :: cell
629
630 cpassert(ASSOCIATED(cell))
631
632 IF (cell%orthorhombic) THEN
633 r(1) = cell%hmat(1, 1)*s(1)
634 r(2) = cell%hmat(2, 2)*s(2)
635 r(3) = cell%hmat(3, 3)*s(3)
636 ELSE
637 r(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
638 r(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
639 r(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
640 END IF
641
642 END SUBROUTINE scaled_to_real
643! **************************************************************************************************
644!> \brief retains the given cell (see doc/ReferenceCounting.html)
645!> \param cell the cell to retain
646!> \par History
647!> 09.2003 created [fawzi]
648!> \author Fawzi Mohamed
649! **************************************************************************************************
650 SUBROUTINE cell_retain(cell)
651
652 TYPE(cell_type), POINTER :: cell
653
654 cpassert(ASSOCIATED(cell))
655 cpassert(cell%ref_count > 0)
656 cell%ref_count = cell%ref_count + 1
657
658 END SUBROUTINE cell_retain
659
660! **************************************************************************************************
661!> \brief releases the given cell (see doc/ReferenceCounting.html)
662!> \param cell the cell to release
663!> \par History
664!> 09.2003 created [fawzi]
665!> \author Fawzi Mohamed
666! **************************************************************************************************
667 SUBROUTINE cell_release(cell)
668
669 TYPE(cell_type), POINTER :: cell
670
671 IF (ASSOCIATED(cell)) THEN
672 cpassert(cell%ref_count > 0)
673 cell%ref_count = cell%ref_count - 1
674 IF (cell%ref_count == 0) THEN
675 DEALLOCATE (cell)
676 END IF
677 NULLIFY (cell)
678 END IF
679
680 END SUBROUTINE cell_release
681
682#if defined (__PLUMED2)
683! **************************************************************************************************
684!> \brief For the interface with plumed, pass a cell pointer and retrieve it
685!> later. It's a hack, but avoids passing the cell back and forth
686!> across the Fortran/C++ interface
687!> \param cell ...
688!> \param set ...
689!> \date 28.02.2013
690!> \author RK
691!> \version 1.0
692! **************************************************************************************************
693 SUBROUTINE pbc_cp2k_plumed_getset_cell(cell, set)
694
695 TYPE(cell_type), POINTER :: cell
696 LOGICAL :: set
697
698 TYPE(cell_type), POINTER, SAVE :: stored_cell
699
700 IF (set) THEN
701 stored_cell => cell
702 ELSE
703 cell => stored_cell
704 END IF
705
706 END SUBROUTINE pbc_cp2k_plumed_getset_cell
707#endif
708
709END MODULE cell_types
cell_types::pbc
Definition cell_types.F:104
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::cell_transform_input_reciprocal
subroutine, public cell_transform_input_reciprocal(cell, vector)
Transform a Cartesian reciprocal-space vector from the user input cell frame into CP2K's canonical in...
Definition cell_types.F:301
cell_types::scaled_to_real
subroutine, public scaled_to_real(r, s, cell)
Transform scaled cell coordinates real coordinates. r=h*s.
Definition cell_types.F:625
cell_types::use_perd_xyz
integer, parameter, public use_perd_xyz
Definition cell_types.F:42
cell_types::parse_cell_line
subroutine, public parse_cell_line(input_line, cell_itimes, cell_time, h, vol)
Read cell info from a line (parsed from a file)
Definition cell_types.F:194
cell_types::cell_sym_monoclinic
integer, parameter, public cell_sym_monoclinic
Definition cell_types.F:29
cell_types::use_perd_y
integer, parameter, public use_perd_y
Definition cell_types.F:42
cell_types::cell_sym_triclinic
integer, parameter, public cell_sym_triclinic
Definition cell_types.F:29
cell_types::cell_sym_tetragonal_ab
integer, parameter, public cell_sym_tetragonal_ab
Definition cell_types.F:29
cell_types::use_perd_xz
integer, parameter, public use_perd_xz
Definition cell_types.F:42
cell_types::cell_sym_rhombohedral
integer, parameter, public cell_sym_rhombohedral
Definition cell_types.F:29
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:595
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:668
cell_types::use_perd_x
integer, parameter, public use_perd_x
Definition cell_types.F:42
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:141
cell_types::cell_sym_tetragonal_ac
integer, parameter, public cell_sym_tetragonal_ac
Definition cell_types.F:29
cell_types::use_perd_z
integer, parameter, public use_perd_z
Definition cell_types.F:42
cell_types::use_perd_yz
integer, parameter, public use_perd_yz
Definition cell_types.F:42
cell_types::get_cell
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:233
cell_types::use_perd_none
integer, parameter, public use_perd_none
Definition cell_types.F:42
cell_types::cell_retain
subroutine, public cell_retain(cell)
retains the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:651
cell_types::cell_sym_hexagonal_gamma_60
integer, parameter, public cell_sym_hexagonal_gamma_60
Definition cell_types.F:29
cell_types::periodicity_string
character(len=3), dimension(7), parameter, public periodicity_string
Definition cell_types.F:51
cell_types::cell_sym_orthorhombic
integer, parameter, public cell_sym_orthorhombic
Definition cell_types.F:29
cell_types::cell_sym_none
integer, parameter, public cell_sym_none
Definition cell_types.F:29
cell_types::cell_sym_hexagonal_gamma_120
integer, parameter, public cell_sym_hexagonal_gamma_120
Definition cell_types.F:29
cell_types::cell_copy
subroutine, public cell_copy(cell_in, cell_out, tag)
Copy cell variable.
Definition cell_types.F:160
cell_types::cell_sym_monoclinic_gamma_ab
integer, parameter, public cell_sym_monoclinic_gamma_ab
Definition cell_types.F:29
cell_types::cell_sym_cubic
integer, parameter, public cell_sym_cubic
Definition cell_types.F:29
cell_types::use_perd_xy
integer, parameter, public use_perd_xy
Definition cell_types.F:42
cell_types::cell_transform_input_cartesian
subroutine, public cell_transform_input_cartesian(cell, vector)
Transform a Cartesian real-space vector from the user input cell frame into CP2K's canonical internal...
Definition cell_types.F:284
cell_types::pbc_stable
real(kind=dp) function, dimension(3), public pbc_stable(r, cell)
Apply a stable periodic-image convention for k-point Bloch gauges.
Definition cell_types.F:422
cell_types::cell_sym_tetragonal_bc
integer, parameter, public cell_sym_tetragonal_bc
Definition cell_types.F:29
cell_types::plane_distance
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:324
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_to_cp2k
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:1222
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::degree
real(kind=dp), parameter, public degree
Definition mathconstants.F:139
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::angle
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:184
cell_types::cell_p_type
Definition cell_types.F:75
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60