Line data Source code
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 : ! **************************************************************************************************
15 : MODULE 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 :
108 : CONTAINS
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 138846 : 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 138846 : image_shift = ANINT(s)
123 138846 : half_boundary = ANINT(s - 0.5_dp) + 0.5_dp
124 138846 : tolerance = MIN(1.0e-8_dp, 64.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, ABS(s)))
125 138846 : IF (ABS(s - half_boundary) <= tolerance) THEN
126 : ! Select the lower side for every lattice-equivalent half-cell boundary.
127 15540 : image_shift = half_boundary + 0.5_dp
128 : END IF
129 :
130 138846 : 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 67816 : 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 67816 : cell_out = cell_in
146 67816 : cell_out%ref_count = 1
147 11422 : IF (PRESENT(tag)) cell_out%tag = tag
148 :
149 67816 : 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 232714 : 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 232714 : cell_out%deth = cell_in%deth
165 1861712 : cell_out%perd = cell_in%perd
166 6050564 : cell_out%hmat = cell_in%hmat
167 6050564 : cell_out%h_inv = cell_in%h_inv
168 232714 : cell_out%input_cell_canonicalized = cell_in%input_cell_canonicalized
169 6050564 : cell_out%input_hmat = cell_in%input_hmat
170 6050564 : cell_out%input_to_canonical = cell_in%input_to_canonical
171 6050564 : cell_out%input_recip_to_canonical = cell_in%input_recip_to_canonical
172 232714 : cell_out%orthorhombic = cell_in%orthorhombic
173 232714 : cell_out%symmetry_id = cell_in%symmetry_id
174 232714 : IF (PRESENT(tag)) THEN
175 13054 : cell_out%tag = tag
176 : ELSE
177 219660 : cell_out%tag = cell_in%tag
178 : END IF
179 :
180 232714 : 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 344 : 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 344 : READ (input_line, *) cell_itimes, cell_time, &
204 688 : 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 1376 : DO i = 1, 3
206 4472 : DO j = 1, 3
207 4128 : h(j, i) = cp_unit_to_cp2k(h(j, i), "angstrom")
208 : END DO
209 : END DO
210 :
211 344 : 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 140250216 : 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 0 : CPASSERT(ASSOCIATED(cell))
245 :
246 140250216 : IF (PRESENT(deth)) deth = cell%deth ! the volume
247 140250216 : IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic
248 538133355 : IF (PRESENT(periodic)) periodic(:) = cell%perd(:)
249 140609412 : IF (PRESENT(h)) h(:, :) = cell%hmat(:, :)
250 140250720 : IF (PRESENT(h_inv)) h_inv(:, :) = cell%h_inv(:, :)
251 :
252 : ! Calculate the lengths of the cell vectors a, b, and c
253 140250216 : 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 7404072 : 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 7404072 : 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 7404072 : 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 140250216 : IF (PRESENT(alpha)) alpha = angle(cell%hmat(:, 2), cell%hmat(:, 3))*degree
268 : ! beta = <(a,c)
269 140250216 : IF (PRESENT(beta)) beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))*degree
270 : ! gamma = <(a,b)
271 140250216 : IF (PRESENT(gamma)) gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))*degree
272 140250216 : IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id
273 140250216 : IF (PRESENT(tag)) tag = cell%tag
274 :
275 140250216 : 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 786695 : 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 786695 : CPASSERT(ASSOCIATED(cell))
289 :
290 953087 : IF (cell%input_cell_canonicalized) vector = MATMUL(cell%input_to_canonical, vector)
291 :
292 786695 : 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 0 : 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 0 : CPASSERT(ASSOCIATED(cell))
306 :
307 0 : IF (cell%input_cell_canonicalized) vector = MATMUL(cell%input_recip_to_canonical, vector)
308 :
309 0 : 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 7295178 : 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 7295178 : x = REAL(h, KIND=dp)
334 7295178 : y = REAL(k, KIND=dp)
335 7295178 : z = REAL(l, KIND=dp)
336 :
337 7295178 : CALL get_cell(cell=cell, abc=abc)
338 :
339 7295178 : a = abc(1)
340 7295178 : b = abc(2)
341 7295178 : c = abc(3)
342 :
343 7295178 : IF (cell%orthorhombic) THEN
344 :
345 7086623 : 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 208555 : gamma=gamma)
353 :
354 208555 : alpha = alpha/degree
355 208555 : beta = beta/degree
356 208555 : gamma = gamma/degree
357 :
358 208555 : cosa = COS(alpha)
359 208555 : cosb = COS(beta)
360 208555 : 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 208555 : 2.0_dp*cosa*cosb*cosg))
370 :
371 : END IF
372 :
373 7295178 : distance = 1.0_dp/SQRT(d)
374 :
375 7295178 : 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 412187116 : 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 412187116 : CPASSERT(ASSOCIATED(cell))
396 :
397 412187116 : IF (cell%orthorhombic) THEN
398 386381562 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*ANINT(cell%h_inv(1, 1)*r(1))
399 386381562 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*ANINT(cell%h_inv(2, 2)*r(2))
400 386381562 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*ANINT(cell%h_inv(3, 3)*r(3))
401 : ELSE
402 25805554 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
403 25805554 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
404 25805554 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
405 25805554 : s(1) = s(1) - cell%perd(1)*ANINT(s(1))
406 25805554 : s(2) = s(2) - cell%perd(2)*ANINT(s(2))
407 25805554 : s(3) = s(3) - cell%perd(3)*ANINT(s(3))
408 25805554 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
409 25805554 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
410 25805554 : 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 412187116 : 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 46282 : 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 46282 : CPASSERT(ASSOCIATED(cell))
430 :
431 46282 : IF (cell%orthorhombic) THEN
432 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
433 36998 : pbc_image_shift(cell%h_inv(1, 1)*r(1))
434 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
435 36998 : pbc_image_shift(cell%h_inv(2, 2)*r(2))
436 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
437 36998 : pbc_image_shift(cell%h_inv(3, 3)*r(3))
438 : ELSE
439 9284 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
440 9284 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
441 9284 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
442 9284 : s(1) = s(1) - cell%perd(1)*pbc_image_shift(s(1))
443 9284 : s(2) = s(2) - cell%perd(2)*pbc_image_shift(s(2))
444 9284 : s(3) = s(3) - cell%perd(3)*pbc_image_shift(s(3))
445 9284 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
446 9284 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
447 9284 : 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 46282 : 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 0 : 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 0 : CPASSERT(ASSOCIATED(cell))
473 :
474 0 : IF (cell%orthorhombic) THEN
475 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
476 0 : 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 0 : 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 0 : REAL(NINT(cell%h_inv(3, 3)*r(3)) - nl(3), dp)
481 : ELSE
482 0 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
483 0 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
484 0 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
485 0 : s(1) = s(1) - cell%perd(1)*REAL(NINT(s(1)) - nl(1), dp)
486 0 : s(2) = s(2) - cell%perd(2)*REAL(NINT(s(2)) - nl(2), dp)
487 0 : s(3) = s(3) - cell%perd(3)*REAL(NINT(s(3)) - nl(3), dp)
488 0 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
489 0 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
490 0 : 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 0 : 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 132406589 : 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 132406589 : CALL get_cell(cell=cell, periodic=periodic)
518 :
519 132406589 : ra_pbc(:) = pbc(ra(:), cell)
520 132406589 : rb_pbc(:) = pbc(rb(:), cell)
521 :
522 132406589 : rab2_pbc = HUGE(1.0_dp)
523 :
524 523099872 : DO icell = -periodic(1), periodic(1)
525 1688654957 : DO jcell = -periodic(2), periodic(2)
526 5046383795 : DO kcell = -periodic(3), periodic(3)
527 13960541708 : r = REAL([icell, jcell, kcell], dp)
528 3490135427 : CALL scaled_to_real(s2r, r, cell)
529 13960541708 : rb_image(:) = rb_pbc(:) + s2r
530 13960541708 : rab(:) = rb_image(:) - ra_pbc(:)
531 3490135427 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
532 4655690512 : IF (rab2 < rab2_pbc) THEN
533 2736867664 : rab2_pbc = rab2
534 2736867664 : rab_pbc(:) = rab(:)
535 : END IF
536 : END DO
537 : END DO
538 : END DO
539 :
540 132406589 : 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 304169 : 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 304169 : CPASSERT(ASSOCIATED(cell))
561 :
562 304169 : IF (positive_range) THEN
563 304169 : IF (cell%orthorhombic) THEN
564 72102 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*FLOOR(cell%h_inv(1, 1)*r(1))
565 72102 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*FLOOR(cell%h_inv(2, 2)*r(2))
566 72102 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*FLOOR(cell%h_inv(3, 3)*r(3))
567 : ELSE
568 232067 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
569 232067 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
570 232067 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
571 232067 : s(1) = s(1) - cell%perd(1)*FLOOR(s(1))
572 232067 : s(2) = s(2) - cell%perd(2)*FLOOR(s(2))
573 232067 : s(3) = s(3) - cell%perd(3)*FLOOR(s(3))
574 232067 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
575 232067 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
576 232067 : 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 0 : r_pbc = pbc1(r, cell)
580 : END IF
581 :
582 304169 : 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 111002306 : 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 111002306 : CPASSERT(ASSOCIATED(cell))
601 :
602 111002306 : IF (cell%orthorhombic) THEN
603 100093976 : s(1) = cell%h_inv(1, 1)*r(1)
604 100093976 : s(2) = cell%h_inv(2, 2)*r(2)
605 100093976 : s(3) = cell%h_inv(3, 3)*r(3)
606 : ELSE
607 10908330 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
608 10908330 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
609 10908330 : 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 111002306 : 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 3659191210 : 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 3659191210 : CPASSERT(ASSOCIATED(cell))
631 :
632 3659191210 : IF (cell%orthorhombic) THEN
633 3383846503 : r(1) = cell%hmat(1, 1)*s(1)
634 3383846503 : r(2) = cell%hmat(2, 2)*s(2)
635 3383846503 : r(3) = cell%hmat(3, 3)*s(3)
636 : ELSE
637 275344707 : r(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
638 275344707 : r(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
639 275344707 : 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 3659191210 : 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 77552 : SUBROUTINE cell_retain(cell)
651 :
652 : TYPE(cell_type), POINTER :: cell
653 :
654 77552 : CPASSERT(ASSOCIATED(cell))
655 77552 : CPASSERT(cell%ref_count > 0)
656 77552 : cell%ref_count = cell%ref_count + 1
657 :
658 77552 : 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 207785 : SUBROUTINE cell_release(cell)
668 :
669 : TYPE(cell_type), POINTER :: cell
670 :
671 207785 : IF (ASSOCIATED(cell)) THEN
672 158811 : CPASSERT(cell%ref_count > 0)
673 158811 : cell%ref_count = cell%ref_count - 1
674 158811 : IF (cell%ref_count == 0) THEN
675 81259 : DEALLOCATE (cell)
676 : END IF
677 158811 : NULLIFY (cell)
678 : END IF
679 :
680 207785 : 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 2 : 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 2 : IF (set) THEN
701 2 : stored_cell => cell
702 : ELSE
703 0 : cell => stored_cell
704 : END IF
705 :
706 2 : END SUBROUTINE pbc_cp2k_plumed_getset_cell
707 : #endif
708 :
709 0 : END MODULE cell_types
|
