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 : real_to_scaled, &
101 : scaled_to_real
102 :
103 : INTERFACE pbc
104 : MODULE PROCEDURE pbc1, pbc2, pbc3, pbc4
105 : END INTERFACE
106 :
107 : CONTAINS
108 :
109 : ! **************************************************************************************************
110 : !> \brief Clone cell variable
111 : !> \param cell_in Cell variable to be clone
112 : !> \param cell_out Cloned cell variable
113 : !> \param tag Optional new tag for cloned cell variable
114 : !> \par History
115 : !> - Optional tag added (17.05.2023, MK)
116 : ! **************************************************************************************************
117 66245 : SUBROUTINE cell_clone(cell_in, cell_out, tag)
118 :
119 : TYPE(cell_type), POINTER :: cell_in, cell_out
120 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
121 :
122 66245 : cell_out = cell_in
123 66245 : cell_out%ref_count = 1
124 11110 : IF (PRESENT(tag)) cell_out%tag = tag
125 :
126 66245 : END SUBROUTINE cell_clone
127 :
128 : ! **************************************************************************************************
129 : !> \brief Copy cell variable
130 : !> \param cell_in Cell variable to be copied
131 : !> \param cell_out Copy of cell variable
132 : !> \param tag Optional new tag
133 : !> \par History
134 : !> - Optional tag added (17.05.2023, MK)
135 : ! **************************************************************************************************
136 232654 : SUBROUTINE cell_copy(cell_in, cell_out, tag)
137 :
138 : TYPE(cell_type), POINTER :: cell_in, cell_out
139 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
140 :
141 232654 : cell_out%deth = cell_in%deth
142 1861232 : cell_out%perd = cell_in%perd
143 6049004 : cell_out%hmat = cell_in%hmat
144 6049004 : cell_out%h_inv = cell_in%h_inv
145 232654 : cell_out%input_cell_canonicalized = cell_in%input_cell_canonicalized
146 6049004 : cell_out%input_hmat = cell_in%input_hmat
147 6049004 : cell_out%input_to_canonical = cell_in%input_to_canonical
148 6049004 : cell_out%input_recip_to_canonical = cell_in%input_recip_to_canonical
149 232654 : cell_out%orthorhombic = cell_in%orthorhombic
150 232654 : cell_out%symmetry_id = cell_in%symmetry_id
151 232654 : IF (PRESENT(tag)) THEN
152 12994 : cell_out%tag = tag
153 : ELSE
154 219660 : cell_out%tag = cell_in%tag
155 : END IF
156 :
157 232654 : END SUBROUTINE cell_copy
158 :
159 : ! **************************************************************************************************
160 : !> \brief Read cell info from a line (parsed from a file)
161 : !> \param input_line ...
162 : !> \param cell_itimes ...
163 : !> \param cell_time ...
164 : !> \param h ...
165 : !> \param vol ...
166 : !> \date 19.02.2008
167 : !> \author Teodoro Laino [tlaino] - University of Zurich
168 : !> \version 1.0
169 : ! **************************************************************************************************
170 368 : SUBROUTINE parse_cell_line(input_line, cell_itimes, cell_time, h, vol)
171 :
172 : CHARACTER(LEN=*), INTENT(IN) :: input_line
173 : INTEGER, INTENT(OUT) :: cell_itimes
174 : REAL(KIND=dp), INTENT(OUT) :: cell_time
175 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: h
176 : REAL(KIND=dp), INTENT(OUT) :: vol
177 :
178 : INTEGER :: i, j
179 :
180 368 : READ (input_line, *) cell_itimes, cell_time, &
181 736 : 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
182 1472 : DO i = 1, 3
183 4784 : DO j = 1, 3
184 4416 : h(j, i) = cp_unit_to_cp2k(h(j, i), "angstrom")
185 : END DO
186 : END DO
187 :
188 368 : END SUBROUTINE parse_cell_line
189 :
190 : ! **************************************************************************************************
191 : !> \brief Get informations about a simulation cell.
192 : !> \param cell ...
193 : !> \param alpha ...
194 : !> \param beta ...
195 : !> \param gamma ...
196 : !> \param deth ...
197 : !> \param orthorhombic ...
198 : !> \param abc ...
199 : !> \param periodic ...
200 : !> \param h ...
201 : !> \param h_inv ...
202 : !> \param symmetry_id ...
203 : !> \param tag ...
204 : !> \date 16.01.2002
205 : !> \author Matthias Krack
206 : !> \version 1.0
207 : ! **************************************************************************************************
208 136235612 : SUBROUTINE get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, &
209 : h, h_inv, symmetry_id, tag)
210 :
211 : TYPE(cell_type), POINTER :: cell
212 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha, beta, gamma, deth
213 : LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic
214 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: abc
215 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: periodic
216 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
217 : OPTIONAL :: h, h_inv
218 : INTEGER, INTENT(OUT), OPTIONAL :: symmetry_id
219 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: tag
220 :
221 0 : CPASSERT(ASSOCIATED(cell))
222 :
223 136235612 : IF (PRESENT(deth)) deth = cell%deth ! the volume
224 136235612 : IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic
225 522349709 : IF (PRESENT(periodic)) periodic(:) = cell%perd(:)
226 136454108 : IF (PRESENT(h)) h(:, :) = cell%hmat(:, :)
227 136236116 : IF (PRESENT(h_inv)) h_inv(:, :) = cell%h_inv(:, :)
228 :
229 : ! Calculate the lengths of the cell vectors a, b, and c
230 136235612 : IF (PRESENT(abc)) THEN
231 : abc(1) = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
232 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
233 7339427 : cell%hmat(3, 1)*cell%hmat(3, 1))
234 : abc(2) = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
235 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
236 7339427 : cell%hmat(3, 2)*cell%hmat(3, 2))
237 : abc(3) = SQRT(cell%hmat(1, 3)*cell%hmat(1, 3) + &
238 : cell%hmat(2, 3)*cell%hmat(2, 3) + &
239 7339427 : cell%hmat(3, 3)*cell%hmat(3, 3))
240 : END IF
241 :
242 : ! Angles between the cell vectors a, b, and c
243 : ! alpha = <(b,c)
244 136235612 : IF (PRESENT(alpha)) alpha = angle(cell%hmat(:, 2), cell%hmat(:, 3))*degree
245 : ! beta = <(a,c)
246 136235612 : IF (PRESENT(beta)) beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))*degree
247 : ! gamma = <(a,b)
248 136235612 : IF (PRESENT(gamma)) gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))*degree
249 136235612 : IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id
250 136235612 : IF (PRESENT(tag)) tag = cell%tag
251 :
252 136235612 : END SUBROUTINE get_cell
253 :
254 : ! **************************************************************************************************
255 : !> \brief Transform a Cartesian real-space vector from the user input cell frame
256 : !> into CP2K's canonical internal cell frame.
257 : !> \param cell ...
258 : !> \param vector ...
259 : ! **************************************************************************************************
260 784157 : SUBROUTINE cell_transform_input_cartesian(cell, vector)
261 :
262 : TYPE(cell_type), POINTER :: cell
263 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: vector
264 :
265 784157 : CPASSERT(ASSOCIATED(cell))
266 :
267 950549 : IF (cell%input_cell_canonicalized) vector = MATMUL(cell%input_to_canonical, vector)
268 :
269 784157 : END SUBROUTINE cell_transform_input_cartesian
270 :
271 : ! **************************************************************************************************
272 : !> \brief Transform a Cartesian reciprocal-space vector from the user input cell
273 : !> frame into CP2K's canonical internal cell frame.
274 : !> \param cell ...
275 : !> \param vector ...
276 : ! **************************************************************************************************
277 0 : SUBROUTINE cell_transform_input_reciprocal(cell, vector)
278 :
279 : TYPE(cell_type), POINTER :: cell
280 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: vector
281 :
282 0 : CPASSERT(ASSOCIATED(cell))
283 :
284 0 : IF (cell%input_cell_canonicalized) vector = MATMUL(cell%input_recip_to_canonical, vector)
285 :
286 0 : END SUBROUTINE cell_transform_input_reciprocal
287 :
288 : ! **************************************************************************************************
289 : !> \brief Calculate the distance between two lattice planes as defined by
290 : !> a triple of Miller indices (hkl).
291 : !> \param h ...
292 : !> \param k ...
293 : !> \param l ...
294 : !> \param cell ...
295 : !> \return ...
296 : !> \date 18.11.2004
297 : !> \author Matthias Krack
298 : !> \version 1.0
299 : ! **************************************************************************************************
300 7236380 : FUNCTION plane_distance(h, k, l, cell) RESULT(distance)
301 :
302 : INTEGER, INTENT(IN) :: h, k, l
303 : TYPE(cell_type), POINTER :: cell
304 : REAL(KIND=dp) :: distance
305 :
306 : REAL(KIND=dp) :: a, alpha, b, beta, c, cosa, cosb, cosg, &
307 : d, gamma, x, y, z
308 : REAL(KIND=dp), DIMENSION(3) :: abc
309 :
310 7236380 : x = REAL(h, KIND=dp)
311 7236380 : y = REAL(k, KIND=dp)
312 7236380 : z = REAL(l, KIND=dp)
313 :
314 7236380 : CALL get_cell(cell=cell, abc=abc)
315 :
316 7236380 : a = abc(1)
317 7236380 : b = abc(2)
318 7236380 : c = abc(3)
319 :
320 7236380 : IF (cell%orthorhombic) THEN
321 :
322 7042433 : d = (x/a)**2 + (y/b)**2 + (z/c)**2
323 :
324 : ELSE
325 :
326 : CALL get_cell(cell=cell, &
327 : alpha=alpha, &
328 : beta=beta, &
329 193947 : gamma=gamma)
330 :
331 193947 : alpha = alpha/degree
332 193947 : beta = beta/degree
333 193947 : gamma = gamma/degree
334 :
335 193947 : cosa = COS(alpha)
336 193947 : cosb = COS(beta)
337 193947 : cosg = COS(gamma)
338 :
339 : d = ((x*b*c*SIN(alpha))**2 + &
340 : (y*c*a*SIN(beta))**2 + &
341 : (z*a*b*SIN(gamma))**2 + &
342 : 2.0_dp*a*b*c*(x*y*c*(cosa*cosb - cosg) + &
343 : z*x*b*(cosg*cosa - cosb) + &
344 : y*z*a*(cosb*cosg - cosa)))/ &
345 : ((a*b*c)**2*(1.0_dp - cosa**2 - cosb**2 - cosg**2 + &
346 193947 : 2.0_dp*cosa*cosb*cosg))
347 :
348 : END IF
349 :
350 7236380 : distance = 1.0_dp/SQRT(d)
351 :
352 7236380 : END FUNCTION plane_distance
353 :
354 : ! **************************************************************************************************
355 : !> \brief Apply the periodic boundary conditions defined by a simulation
356 : !> cell to a position vector r.
357 : !> \param r ...
358 : !> \param cell ...
359 : !> \return ...
360 : !> \date 16.01.2002
361 : !> \author Matthias Krack
362 : !> \version 1.0
363 : ! **************************************************************************************************
364 395239065 : FUNCTION pbc1(r, cell) RESULT(r_pbc)
365 :
366 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
367 : TYPE(cell_type), POINTER :: cell
368 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
369 :
370 : REAL(KIND=dp), DIMENSION(3) :: s
371 :
372 395239065 : CPASSERT(ASSOCIATED(cell))
373 :
374 395239065 : IF (cell%orthorhombic) THEN
375 369675300 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*ANINT(cell%h_inv(1, 1)*r(1))
376 369675300 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*ANINT(cell%h_inv(2, 2)*r(2))
377 369675300 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*ANINT(cell%h_inv(3, 3)*r(3))
378 : ELSE
379 25563765 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
380 25563765 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
381 25563765 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
382 25563765 : s(1) = s(1) - cell%perd(1)*ANINT(s(1))
383 25563765 : s(2) = s(2) - cell%perd(2)*ANINT(s(2))
384 25563765 : s(3) = s(3) - cell%perd(3)*ANINT(s(3))
385 25563765 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
386 25563765 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
387 25563765 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
388 : END IF
389 :
390 395239065 : END FUNCTION pbc1
391 :
392 : ! **************************************************************************************************
393 : !> \brief Apply the periodic boundary conditions defined by a simulation
394 : !> cell to a position vector r subtracting nl from the periodic images
395 : !> \param r ...
396 : !> \param cell ...
397 : !> \param nl ...
398 : !> \return ...
399 : !> \date 16.01.2002
400 : !> \author Matthias Krack
401 : !> \version 1.0
402 : ! **************************************************************************************************
403 142966578 : FUNCTION pbc2(r, cell, nl) RESULT(r_pbc)
404 :
405 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
406 : TYPE(cell_type), POINTER :: cell
407 : INTEGER, DIMENSION(3), INTENT(IN) :: nl
408 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
409 :
410 : REAL(KIND=dp), DIMENSION(3) :: s
411 :
412 142966578 : CPASSERT(ASSOCIATED(cell))
413 :
414 142966578 : IF (cell%orthorhombic) THEN
415 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
416 132353154 : REAL(NINT(cell%h_inv(1, 1)*r(1)) - nl(1), dp)
417 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
418 132353154 : REAL(NINT(cell%h_inv(2, 2)*r(2)) - nl(2), dp)
419 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
420 132353154 : REAL(NINT(cell%h_inv(3, 3)*r(3)) - nl(3), dp)
421 : ELSE
422 10613424 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
423 10613424 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
424 10613424 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
425 10613424 : s(1) = s(1) - cell%perd(1)*REAL(NINT(s(1)) - nl(1), dp)
426 10613424 : s(2) = s(2) - cell%perd(2)*REAL(NINT(s(2)) - nl(2), dp)
427 10613424 : s(3) = s(3) - cell%perd(3)*REAL(NINT(s(3)) - nl(3), dp)
428 10613424 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
429 10613424 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
430 10613424 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
431 : END IF
432 :
433 142966578 : END FUNCTION pbc2
434 :
435 : ! **************************************************************************************************
436 : !> \brief Apply the periodic boundary conditions defined by the simulation
437 : !> cell cell to the vector pointing from atom a to atom b.
438 : !> \param ra ...
439 : !> \param rb ...
440 : !> \param cell ...
441 : !> \return ...
442 : !> \date 11.03.2004
443 : !> \author Matthias Krack
444 : !> \version 1.0
445 : ! **************************************************************************************************
446 128490057 : FUNCTION pbc3(ra, rb, cell) RESULT(rab_pbc)
447 :
448 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
449 : TYPE(cell_type), POINTER :: cell
450 : REAL(KIND=dp), DIMENSION(3) :: rab_pbc
451 :
452 : INTEGER :: icell, jcell, kcell
453 : INTEGER, DIMENSION(3) :: periodic
454 : REAL(KIND=dp) :: rab2, rab2_pbc
455 : REAL(KIND=dp), DIMENSION(3) :: r, ra_pbc, rab, rb_image, rb_pbc, s2r
456 :
457 128490057 : CALL get_cell(cell=cell, periodic=periodic)
458 :
459 128490057 : ra_pbc(:) = pbc(ra(:), cell)
460 128490057 : rb_pbc(:) = pbc(rb(:), cell)
461 :
462 128490057 : rab2_pbc = HUGE(1.0_dp)
463 :
464 509570884 : DO icell = -periodic(1), periodic(1)
465 1648425681 : DO jcell = -periodic(2), periodic(2)
466 4932107303 : DO kcell = -periodic(3), periodic(3)
467 13648686716 : r = REAL([icell, jcell, kcell], dp)
468 3412171679 : CALL scaled_to_real(s2r, r, cell)
469 13648686716 : rb_image(:) = rb_pbc(:) + s2r
470 13648686716 : rab(:) = rb_image(:) - ra_pbc(:)
471 3412171679 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
472 4551026476 : IF (rab2 < rab2_pbc) THEN
473 2678646096 : rab2_pbc = rab2
474 2678646096 : rab_pbc(:) = rab(:)
475 : END IF
476 : END DO
477 : END DO
478 : END DO
479 :
480 128490057 : END FUNCTION pbc3
481 :
482 : !if positive_range == true, r(i) (or s(i)) in range [0, hmat(i,i)],
483 : !else, r(i) (s(i)) in range [-hmat(i,i)/2, hmat(i,i)/2]
484 : ! **************************************************************************************************
485 : !> \brief ...
486 : !> \param r ...
487 : !> \param cell ...
488 : !> \param positive_range ...
489 : !> \return ...
490 : ! **************************************************************************************************
491 19091 : FUNCTION pbc4(r, cell, positive_range) RESULT(r_pbc)
492 :
493 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
494 : TYPE(cell_type), POINTER :: cell
495 : LOGICAL :: positive_range
496 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
497 :
498 : REAL(KIND=dp), DIMENSION(3) :: s
499 :
500 19091 : CPASSERT(ASSOCIATED(cell))
501 :
502 19091 : IF (positive_range) THEN
503 19091 : IF (cell%orthorhombic) THEN
504 11556 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*FLOOR(cell%h_inv(1, 1)*r(1))
505 11556 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*FLOOR(cell%h_inv(2, 2)*r(2))
506 11556 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*FLOOR(cell%h_inv(3, 3)*r(3))
507 : ELSE
508 7535 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
509 7535 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
510 7535 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
511 7535 : s(1) = s(1) - cell%perd(1)*FLOOR(s(1))
512 7535 : s(2) = s(2) - cell%perd(2)*FLOOR(s(2))
513 7535 : s(3) = s(3) - cell%perd(3)*FLOOR(s(3))
514 7535 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
515 7535 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
516 7535 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
517 : END IF
518 : ELSE
519 0 : r_pbc = pbc1(r, cell)
520 : END IF
521 :
522 19091 : END FUNCTION pbc4
523 :
524 : ! **************************************************************************************************
525 : !> \brief Transform real to scaled cell coordinates.
526 : !> s=h_inv*r
527 : !> \param s ...
528 : !> \param r ...
529 : !> \param cell ...
530 : !> \date 16.01.2002
531 : !> \author Matthias Krack
532 : !> \version 1.0
533 : ! **************************************************************************************************
534 110624578 : SUBROUTINE real_to_scaled(s, r, cell)
535 :
536 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: s
537 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
538 : TYPE(cell_type), POINTER :: cell
539 :
540 110624578 : CPASSERT(ASSOCIATED(cell))
541 :
542 110624578 : IF (cell%orthorhombic) THEN
543 99964564 : s(1) = cell%h_inv(1, 1)*r(1)
544 99964564 : s(2) = cell%h_inv(2, 2)*r(2)
545 99964564 : s(3) = cell%h_inv(3, 3)*r(3)
546 : ELSE
547 10660014 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
548 10660014 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
549 10660014 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
550 : END IF
551 :
552 110624578 : END SUBROUTINE real_to_scaled
553 :
554 : ! **************************************************************************************************
555 : !> \brief Transform scaled cell coordinates real coordinates.
556 : !> r=h*s
557 : !> \param r ...
558 : !> \param s ...
559 : !> \param cell ...
560 : !> \date 16.01.2002
561 : !> \author Matthias Krack
562 : !> \version 1.0
563 : ! **************************************************************************************************
564 3576745176 : SUBROUTINE scaled_to_real(r, s, cell)
565 :
566 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: r
567 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s
568 : TYPE(cell_type), POINTER :: cell
569 :
570 3576745176 : CPASSERT(ASSOCIATED(cell))
571 :
572 3576745176 : IF (cell%orthorhombic) THEN
573 3302903255 : r(1) = cell%hmat(1, 1)*s(1)
574 3302903255 : r(2) = cell%hmat(2, 2)*s(2)
575 3302903255 : r(3) = cell%hmat(3, 3)*s(3)
576 : ELSE
577 273841921 : r(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
578 273841921 : r(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
579 273841921 : r(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
580 : END IF
581 :
582 3576745176 : END SUBROUTINE scaled_to_real
583 : ! **************************************************************************************************
584 : !> \brief retains the given cell (see doc/ReferenceCounting.html)
585 : !> \param cell the cell to retain
586 : !> \par History
587 : !> 09.2003 created [fawzi]
588 : !> \author Fawzi Mohamed
589 : ! **************************************************************************************************
590 75688 : SUBROUTINE cell_retain(cell)
591 :
592 : TYPE(cell_type), POINTER :: cell
593 :
594 75688 : CPASSERT(ASSOCIATED(cell))
595 75688 : CPASSERT(cell%ref_count > 0)
596 75688 : cell%ref_count = cell%ref_count + 1
597 :
598 75688 : END SUBROUTINE cell_retain
599 :
600 : ! **************************************************************************************************
601 : !> \brief releases the given cell (see doc/ReferenceCounting.html)
602 : !> \param cell the cell to release
603 : !> \par History
604 : !> 09.2003 created [fawzi]
605 : !> \author Fawzi Mohamed
606 : ! **************************************************************************************************
607 203034 : SUBROUTINE cell_release(cell)
608 :
609 : TYPE(cell_type), POINTER :: cell
610 :
611 203034 : IF (ASSOCIATED(cell)) THEN
612 155306 : CPASSERT(cell%ref_count > 0)
613 155306 : cell%ref_count = cell%ref_count - 1
614 155306 : IF (cell%ref_count == 0) THEN
615 79618 : DEALLOCATE (cell)
616 : END IF
617 155306 : NULLIFY (cell)
618 : END IF
619 :
620 203034 : END SUBROUTINE cell_release
621 :
622 : #if defined (__PLUMED2)
623 : ! **************************************************************************************************
624 : !> \brief For the interface with plumed, pass a cell pointer and retrieve it
625 : !> later. It's a hack, but avoids passing the cell back and forth
626 : !> across the Fortran/C++ interface
627 : !> \param cell ...
628 : !> \param set ...
629 : !> \date 28.02.2013
630 : !> \author RK
631 : !> \version 1.0
632 : ! **************************************************************************************************
633 2 : SUBROUTINE pbc_cp2k_plumed_getset_cell(cell, set)
634 :
635 : TYPE(cell_type), POINTER :: cell
636 : LOGICAL :: set
637 :
638 : TYPE(cell_type), POINTER, SAVE :: stored_cell
639 :
640 2 : IF (set) THEN
641 2 : stored_cell => cell
642 : ELSE
643 0 : cell => stored_cell
644 : END IF
645 :
646 2 : END SUBROUTINE pbc_cp2k_plumed_getset_cell
647 : #endif
648 :
649 0 : END MODULE cell_types
|
