Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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_x = 0, &
43 : use_perd_y = 1, &
44 : use_perd_z = 2, &
45 : use_perd_xy = 3, &
46 : use_perd_xz = 4, &
47 : use_perd_yz = 5, &
48 : use_perd_xyz = 6, &
49 : use_perd_none = 7
50 :
51 : ! **************************************************************************************************
52 : !> \brief Type defining parameters related to the simulation cell
53 : !> \version 1.0
54 : ! **************************************************************************************************
55 : TYPE cell_type
56 : CHARACTER(LEN=12) :: tag = "CELL"
57 : INTEGER :: ref_count = -1, &
58 : symmetry_id = use_perd_none
59 : LOGICAL :: orthorhombic = .FALSE. ! actually means a diagonal hmat
60 : REAL(KIND=dp) :: deth = 0.0_dp
61 : INTEGER, DIMENSION(3) :: perd = -1
62 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat = 0.0_dp, &
63 : h_inv = 0.0_dp
64 : END TYPE cell_type
65 :
66 : TYPE cell_p_type
67 : TYPE(cell_type), POINTER :: cell => NULL()
68 : END TYPE cell_p_type
69 :
70 : ! Public data types
71 : PUBLIC :: cell_type, &
72 : cell_p_type
73 :
74 : ! Public subroutines
75 : PUBLIC :: cell_clone, &
76 : cell_copy, &
77 : cell_release, &
78 : cell_retain, &
79 : get_cell, &
80 : parse_cell_line
81 :
82 : #if defined (__PLUMED2)
83 : PUBLIC :: pbc_cp2k_plumed_getset_cell
84 : #endif
85 :
86 : ! Public functions
87 : PUBLIC :: plane_distance, &
88 : pbc, &
89 : real_to_scaled, &
90 : scaled_to_real
91 :
92 : INTERFACE pbc
93 : MODULE PROCEDURE pbc1, pbc2, pbc3, pbc4
94 : END INTERFACE
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief Clone cell variable
100 : !> \param cell_in Cell variable to be clone
101 : !> \param cell_out Cloned cell variable
102 : !> \param tag Optional new tag for cloned cell variable
103 : !> \par History
104 : !> - Optional tag added (17.05.2023, MK)
105 : ! **************************************************************************************************
106 31336 : SUBROUTINE cell_clone(cell_in, cell_out, tag)
107 :
108 : TYPE(cell_type), POINTER :: cell_in, cell_out
109 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
110 :
111 31336 : cell_out = cell_in
112 31336 : cell_out%ref_count = 1
113 20097 : IF (PRESENT(tag)) cell_out%tag = tag
114 :
115 31336 : END SUBROUTINE cell_clone
116 :
117 : ! **************************************************************************************************
118 : !> \brief Copy cell variable
119 : !> \param cell_in Cell variable to be copied
120 : !> \param cell_out Copy of cell variable
121 : !> \param tag Optional new tag
122 : !> \par History
123 : !> - Optional tag added (17.05.2023, MK)
124 : ! **************************************************************************************************
125 228562 : SUBROUTINE cell_copy(cell_in, cell_out, tag)
126 :
127 : TYPE(cell_type), POINTER :: cell_in, cell_out
128 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
129 :
130 228562 : cell_out%deth = cell_in%deth
131 1828496 : cell_out%perd = cell_in%perd
132 5942612 : cell_out%hmat = cell_in%hmat
133 5942612 : cell_out%h_inv = cell_in%h_inv
134 228562 : cell_out%orthorhombic = cell_in%orthorhombic
135 228562 : cell_out%symmetry_id = cell_in%symmetry_id
136 228562 : IF (PRESENT(tag)) THEN
137 8902 : cell_out%tag = tag
138 : ELSE
139 219660 : cell_out%tag = cell_in%tag
140 : END IF
141 :
142 228562 : END SUBROUTINE cell_copy
143 :
144 : ! **************************************************************************************************
145 : !> \brief Read cell info from a line (parsed from a file)
146 : !> \param input_line ...
147 : !> \param cell_itimes ...
148 : !> \param cell_time ...
149 : !> \param h ...
150 : !> \param vol ...
151 : !> \date 19.02.2008
152 : !> \author Teodoro Laino [tlaino] - University of Zurich
153 : !> \version 1.0
154 : ! **************************************************************************************************
155 368 : SUBROUTINE parse_cell_line(input_line, cell_itimes, cell_time, h, vol)
156 :
157 : CHARACTER(LEN=*), INTENT(IN) :: input_line
158 : INTEGER, INTENT(OUT) :: cell_itimes
159 : REAL(KIND=dp), INTENT(OUT) :: cell_time
160 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: h
161 : REAL(KIND=dp), INTENT(OUT) :: vol
162 :
163 : INTEGER :: i, j
164 :
165 368 : READ (input_line, *) cell_itimes, cell_time, &
166 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
167 1472 : DO i = 1, 3
168 4784 : DO j = 1, 3
169 4416 : h(j, i) = cp_unit_to_cp2k(h(j, i), "angstrom")
170 : END DO
171 : END DO
172 :
173 368 : END SUBROUTINE parse_cell_line
174 :
175 : ! **************************************************************************************************
176 : !> \brief Get informations about a simulation cell.
177 : !> \param cell ...
178 : !> \param alpha ...
179 : !> \param beta ...
180 : !> \param gamma ...
181 : !> \param deth ...
182 : !> \param orthorhombic ...
183 : !> \param abc ...
184 : !> \param periodic ...
185 : !> \param h ...
186 : !> \param h_inv ...
187 : !> \param symmetry_id ...
188 : !> \param tag ...
189 : !> \date 16.01.2002
190 : !> \author Matthias Krack
191 : !> \version 1.0
192 : ! **************************************************************************************************
193 133026866 : SUBROUTINE get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, &
194 : h, h_inv, symmetry_id, tag)
195 :
196 : TYPE(cell_type), POINTER :: cell
197 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha, beta, gamma, deth
198 : LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic
199 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: abc
200 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: periodic
201 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
202 : OPTIONAL :: h, h_inv
203 : INTEGER, INTENT(OUT), OPTIONAL :: symmetry_id
204 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: tag
205 :
206 0 : CPASSERT(ASSOCIATED(cell))
207 :
208 133026866 : IF (PRESENT(deth)) deth = cell%deth ! the volume
209 133026866 : IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic
210 510577631 : IF (PRESENT(periodic)) periodic(:) = cell%perd(:)
211 133456070 : IF (PRESENT(h)) h(:, :) = cell%hmat(:, :)
212 133026962 : IF (PRESENT(h_inv)) h_inv(:, :) = cell%h_inv(:, :)
213 :
214 : ! Calculate the lengths of the cell vectors a, b, and c
215 133026866 : IF (PRESENT(abc)) THEN
216 : abc(1) = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
217 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
218 7024900 : cell%hmat(3, 1)*cell%hmat(3, 1))
219 : abc(2) = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
220 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
221 7024900 : cell%hmat(3, 2)*cell%hmat(3, 2))
222 : abc(3) = SQRT(cell%hmat(1, 3)*cell%hmat(1, 3) + &
223 : cell%hmat(2, 3)*cell%hmat(2, 3) + &
224 7024900 : cell%hmat(3, 3)*cell%hmat(3, 3))
225 : END IF
226 :
227 : ! Angles between the cell vectors a, b, and c
228 : ! alpha = <(b,c)
229 133026866 : IF (PRESENT(alpha)) alpha = angle(cell%hmat(:, 2), cell%hmat(:, 3))*degree
230 : ! beta = <(a,c)
231 133026866 : IF (PRESENT(beta)) beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))*degree
232 : ! gamma = <(a,b)
233 133026866 : IF (PRESENT(gamma)) gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))*degree
234 133026866 : IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id
235 133026866 : IF (PRESENT(tag)) tag = cell%tag
236 :
237 133026866 : END SUBROUTINE get_cell
238 :
239 : ! **************************************************************************************************
240 : !> \brief Calculate the distance between two lattice planes as defined by
241 : !> a triple of Miller indices (hkl).
242 : !> \param h ...
243 : !> \param k ...
244 : !> \param l ...
245 : !> \param cell ...
246 : !> \return ...
247 : !> \date 18.11.2004
248 : !> \author Matthias Krack
249 : !> \version 1.0
250 : ! **************************************************************************************************
251 6924651 : FUNCTION plane_distance(h, k, l, cell) RESULT(distance)
252 :
253 : INTEGER, INTENT(IN) :: h, k, l
254 : TYPE(cell_type), POINTER :: cell
255 : REAL(KIND=dp) :: distance
256 :
257 : REAL(KIND=dp) :: a, alpha, b, beta, c, cosa, cosb, cosg, &
258 : d, gamma, x, y, z
259 : REAL(KIND=dp), DIMENSION(3) :: abc
260 :
261 6924651 : x = REAL(h, KIND=dp)
262 6924651 : y = REAL(k, KIND=dp)
263 6924651 : z = REAL(l, KIND=dp)
264 :
265 6924651 : CALL get_cell(cell=cell, abc=abc)
266 :
267 6924651 : a = abc(1)
268 6924651 : b = abc(2)
269 6924651 : c = abc(3)
270 :
271 6924651 : IF (cell%orthorhombic) THEN
272 :
273 6790946 : d = (x/a)**2 + (y/b)**2 + (z/c)**2
274 :
275 : ELSE
276 :
277 : CALL get_cell(cell=cell, &
278 : alpha=alpha, &
279 : beta=beta, &
280 133705 : gamma=gamma)
281 :
282 133705 : alpha = alpha/degree
283 133705 : beta = beta/degree
284 133705 : gamma = gamma/degree
285 :
286 133705 : cosa = COS(alpha)
287 133705 : cosb = COS(beta)
288 133705 : cosg = COS(gamma)
289 :
290 : d = ((x*b*c*SIN(alpha))**2 + &
291 : (y*c*a*SIN(beta))**2 + &
292 : (z*a*b*SIN(gamma))**2 + &
293 : 2.0_dp*a*b*c*(x*y*c*(cosa*cosb - cosg) + &
294 : z*x*b*(cosg*cosa - cosb) + &
295 : y*z*a*(cosb*cosg - cosa)))/ &
296 : ((a*b*c)**2*(1.0_dp - cosa**2 - cosb**2 - cosg**2 + &
297 133705 : 2.0_dp*cosa*cosb*cosg))
298 :
299 : END IF
300 :
301 6924651 : distance = 1.0_dp/SQRT(d)
302 :
303 6924651 : END FUNCTION plane_distance
304 :
305 : ! **************************************************************************************************
306 : !> \brief Apply the periodic boundary conditions defined by a simulation
307 : !> cell to a position vector r.
308 : !> \param r ...
309 : !> \param cell ...
310 : !> \return ...
311 : !> \date 16.01.2002
312 : !> \author Matthias Krack
313 : !> \version 1.0
314 : ! **************************************************************************************************
315 373796861 : FUNCTION pbc1(r, cell) RESULT(r_pbc)
316 :
317 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
318 : TYPE(cell_type), POINTER :: cell
319 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
320 :
321 : REAL(KIND=dp), DIMENSION(3) :: s
322 :
323 373796861 : CPASSERT(ASSOCIATED(cell))
324 :
325 373796861 : IF (cell%orthorhombic) THEN
326 348268825 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*ANINT(cell%h_inv(1, 1)*r(1))
327 348268825 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*ANINT(cell%h_inv(2, 2)*r(2))
328 348268825 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*ANINT(cell%h_inv(3, 3)*r(3))
329 : ELSE
330 25528036 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
331 25528036 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
332 25528036 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
333 25528036 : s(1) = s(1) - cell%perd(1)*ANINT(s(1))
334 25528036 : s(2) = s(2) - cell%perd(2)*ANINT(s(2))
335 25528036 : s(3) = s(3) - cell%perd(3)*ANINT(s(3))
336 25528036 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
337 25528036 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
338 25528036 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
339 : END IF
340 :
341 373796861 : END FUNCTION pbc1
342 :
343 : ! **************************************************************************************************
344 : !> \brief Apply the periodic boundary conditions defined by a simulation
345 : !> cell to a position vector r subtracting nl from the periodic images
346 : !> \param r ...
347 : !> \param cell ...
348 : !> \param nl ...
349 : !> \return ...
350 : !> \date 16.01.2002
351 : !> \author Matthias Krack
352 : !> \version 1.0
353 : ! **************************************************************************************************
354 142966578 : FUNCTION pbc2(r, cell, nl) RESULT(r_pbc)
355 :
356 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
357 : TYPE(cell_type), POINTER :: cell
358 : INTEGER, DIMENSION(3), INTENT(IN) :: nl
359 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
360 :
361 : REAL(KIND=dp), DIMENSION(3) :: s
362 :
363 142966578 : CPASSERT(ASSOCIATED(cell))
364 :
365 142966578 : IF (cell%orthorhombic) THEN
366 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
367 132353154 : REAL(NINT(cell%h_inv(1, 1)*r(1)) - nl(1), dp)
368 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
369 132353154 : REAL(NINT(cell%h_inv(2, 2)*r(2)) - nl(2), dp)
370 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
371 132353154 : REAL(NINT(cell%h_inv(3, 3)*r(3)) - nl(3), dp)
372 : ELSE
373 10613424 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
374 10613424 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
375 10613424 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
376 10613424 : s(1) = s(1) - cell%perd(1)*REAL(NINT(s(1)) - nl(1), dp)
377 10613424 : s(2) = s(2) - cell%perd(2)*REAL(NINT(s(2)) - nl(2), dp)
378 10613424 : s(3) = s(3) - cell%perd(3)*REAL(NINT(s(3)) - nl(3), dp)
379 10613424 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
380 10613424 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
381 10613424 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
382 : END IF
383 :
384 142966578 : END FUNCTION pbc2
385 :
386 : ! **************************************************************************************************
387 : !> \brief Apply the periodic boundary conditions defined by the simulation
388 : !> cell cell to the vector pointing from atom a to atom b.
389 : !> \param ra ...
390 : !> \param rb ...
391 : !> \param cell ...
392 : !> \return ...
393 : !> \date 11.03.2004
394 : !> \author Matthias Krack
395 : !> \version 1.0
396 : ! **************************************************************************************************
397 125680225 : FUNCTION pbc3(ra, rb, cell) RESULT(rab_pbc)
398 :
399 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
400 : TYPE(cell_type), POINTER :: cell
401 : REAL(KIND=dp), DIMENSION(3) :: rab_pbc
402 :
403 : INTEGER :: icell, jcell, kcell
404 : INTEGER, DIMENSION(3) :: periodic
405 : REAL(KIND=dp) :: rab2, rab2_pbc
406 : REAL(KIND=dp), DIMENSION(3) :: r, ra_pbc, rab, rb_image, rb_pbc, s2r
407 :
408 125680225 : CALL get_cell(cell=cell, periodic=periodic)
409 :
410 125680225 : ra_pbc(:) = pbc(ra(:), cell)
411 125680225 : rb_pbc(:) = pbc(rb(:), cell)
412 :
413 125680225 : rab2_pbc = HUGE(1.0_dp)
414 :
415 498347250 : DO icell = -periodic(1), periodic(1)
416 1611977171 : DO jcell = -periodic(2), periodic(2)
417 4822809099 : DO kcell = -periodic(3), periodic(3)
418 13346048612 : r = REAL((/icell, jcell, kcell/), dp)
419 3336512153 : CALL scaled_to_real(s2r, r, cell)
420 13346048612 : rb_image(:) = rb_pbc(:) + s2r
421 13346048612 : rab(:) = rb_image(:) - ra_pbc(:)
422 3336512153 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
423 4450142074 : IF (rab2 < rab2_pbc) THEN
424 2610675568 : rab2_pbc = rab2
425 2610675568 : rab_pbc(:) = rab(:)
426 : END IF
427 : END DO
428 : END DO
429 : END DO
430 :
431 125680225 : END FUNCTION pbc3
432 :
433 : !if positive_range == true, r(i) (or s(i)) in range [0, hmat(i,i)],
434 : !else, r(i) (s(i)) in range [-hmat(i,i)/2, hmat(i,i)/2]
435 : ! **************************************************************************************************
436 : !> \brief ...
437 : !> \param r ...
438 : !> \param cell ...
439 : !> \param positive_range ...
440 : !> \return ...
441 : ! **************************************************************************************************
442 238 : FUNCTION pbc4(r, cell, positive_range) RESULT(r_pbc)
443 :
444 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
445 : TYPE(cell_type), POINTER :: cell
446 : LOGICAL :: positive_range
447 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
448 :
449 : REAL(KIND=dp), DIMENSION(3) :: s
450 :
451 238 : CPASSERT(ASSOCIATED(cell))
452 :
453 238 : IF (positive_range) THEN
454 238 : IF (cell%orthorhombic) THEN
455 238 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*FLOOR(cell%h_inv(1, 1)*r(1))
456 238 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*FLOOR(cell%h_inv(2, 2)*r(2))
457 238 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*FLOOR(cell%h_inv(3, 3)*r(3))
458 : ELSE
459 0 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
460 0 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
461 0 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
462 0 : s(1) = s(1) - cell%perd(1)*FLOOR(s(1))
463 0 : s(2) = s(2) - cell%perd(2)*FLOOR(s(2))
464 0 : s(3) = s(3) - cell%perd(3)*FLOOR(s(3))
465 0 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
466 0 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
467 0 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
468 : END IF
469 : ELSE
470 0 : r_pbc = pbc1(r, cell)
471 : END IF
472 :
473 238 : END FUNCTION pbc4
474 :
475 : ! **************************************************************************************************
476 : !> \brief Transform real to scaled cell coordinates.
477 : !> s=h_inv*r
478 : !> \param s ...
479 : !> \param r ...
480 : !> \param cell ...
481 : !> \date 16.01.2002
482 : !> \author Matthias Krack
483 : !> \version 1.0
484 : ! **************************************************************************************************
485 92508800 : SUBROUTINE real_to_scaled(s, r, cell)
486 :
487 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: s
488 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
489 : TYPE(cell_type), POINTER :: cell
490 :
491 92508800 : CPASSERT(ASSOCIATED(cell))
492 :
493 92508800 : IF (cell%orthorhombic) THEN
494 86249442 : s(1) = cell%h_inv(1, 1)*r(1)
495 86249442 : s(2) = cell%h_inv(2, 2)*r(2)
496 86249442 : s(3) = cell%h_inv(3, 3)*r(3)
497 : ELSE
498 6259358 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
499 6259358 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
500 6259358 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
501 : END IF
502 :
503 92508800 : END SUBROUTINE real_to_scaled
504 :
505 : ! **************************************************************************************************
506 : !> \brief Transform scaled cell coordinates real coordinates.
507 : !> r=h*s
508 : !> \param r ...
509 : !> \param s ...
510 : !> \param cell ...
511 : !> \date 16.01.2002
512 : !> \author Matthias Krack
513 : !> \version 1.0
514 : ! **************************************************************************************************
515 3468576630 : SUBROUTINE scaled_to_real(r, s, cell)
516 :
517 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: r
518 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s
519 : TYPE(cell_type), POINTER :: cell
520 :
521 3468576630 : CPASSERT(ASSOCIATED(cell))
522 :
523 3468576630 : IF (cell%orthorhombic) THEN
524 3196184843 : r(1) = cell%hmat(1, 1)*s(1)
525 3196184843 : r(2) = cell%hmat(2, 2)*s(2)
526 3196184843 : r(3) = cell%hmat(3, 3)*s(3)
527 : ELSE
528 272391787 : r(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
529 272391787 : r(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
530 272391787 : r(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
531 : END IF
532 :
533 3468576630 : END SUBROUTINE scaled_to_real
534 : ! **************************************************************************************************
535 : !> \brief retains the given cell (see doc/ReferenceCounting.html)
536 : !> \param cell the cell to retain
537 : !> \par History
538 : !> 09.2003 created [fawzi]
539 : !> \author Fawzi Mohamed
540 : ! **************************************************************************************************
541 56341 : SUBROUTINE cell_retain(cell)
542 :
543 : TYPE(cell_type), POINTER :: cell
544 :
545 56341 : CPASSERT(ASSOCIATED(cell))
546 56341 : CPASSERT(cell%ref_count > 0)
547 56341 : cell%ref_count = cell%ref_count + 1
548 :
549 56341 : END SUBROUTINE cell_retain
550 :
551 : ! **************************************************************************************************
552 : !> \brief releases the given cell (see doc/ReferenceCounting.html)
553 : !> \param cell the cell to release
554 : !> \par History
555 : !> 09.2003 created [fawzi]
556 : !> \author Fawzi Mohamed
557 : ! **************************************************************************************************
558 150439 : SUBROUTINE cell_release(cell)
559 :
560 : TYPE(cell_type), POINTER :: cell
561 :
562 150439 : IF (ASSOCIATED(cell)) THEN
563 116805 : CPASSERT(cell%ref_count > 0)
564 116805 : cell%ref_count = cell%ref_count - 1
565 116805 : IF (cell%ref_count == 0) THEN
566 60464 : DEALLOCATE (cell)
567 : END IF
568 116805 : NULLIFY (cell)
569 : END IF
570 :
571 150439 : END SUBROUTINE cell_release
572 :
573 : #if defined (__PLUMED2)
574 : ! **************************************************************************************************
575 : !> \brief For the interface with plumed, pass a cell pointer and retrieve it
576 : !> later. It's a hack, but avoids passing the cell back and forth
577 : !> across the Fortran/C++ interface
578 : !> \param cell ...
579 : !> \param set ...
580 : !> \date 28.02.2013
581 : !> \author RK
582 : !> \version 1.0
583 : ! **************************************************************************************************
584 2 : SUBROUTINE pbc_cp2k_plumed_getset_cell(cell, set)
585 :
586 : TYPE(cell_type), POINTER :: cell
587 : LOGICAL :: set
588 :
589 : TYPE(cell_type), POINTER, SAVE :: stored_cell
590 :
591 2 : IF (set) THEN
592 2 : stored_cell => cell
593 : ELSE
594 0 : cell => stored_cell
595 : END IF
596 :
597 2 : END SUBROUTINE pbc_cp2k_plumed_getset_cell
598 : #endif
599 :
600 0 : END MODULE cell_types
|
