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 : minimum_image_weights, &
92 : parse_cell_line
93 :
94 : #if defined (__PLUMED2)
95 : PUBLIC :: pbc_cp2k_plumed_getset_cell
96 : #endif
97 :
98 : ! Public functions
99 : PUBLIC :: plane_distance, &
100 : pbc, &
101 : pbc_stable, &
102 : real_to_scaled, &
103 : scaled_to_real
104 :
105 : INTERFACE pbc
106 : MODULE PROCEDURE pbc1, pbc2, pbc3, pbc4
107 : END INTERFACE
108 :
109 : CONTAINS
110 :
111 : ! **************************************************************************************************
112 : !> \brief Select a stable periodic image index close to half-cell boundaries.
113 : !> \param s Scaled coordinate
114 : !> \return image_shift ...
115 : ! **************************************************************************************************
116 428049 : PURE ELEMENTAL FUNCTION pbc_image_shift(s) RESULT(image_shift)
117 :
118 : REAL(KIND=dp), INTENT(IN) :: s
119 : REAL(KIND=dp) :: image_shift
120 :
121 : REAL(KIND=dp) :: half_boundary, tolerance
122 :
123 428049 : image_shift = ANINT(s)
124 428049 : half_boundary = ANINT(s - 0.5_dp) + 0.5_dp
125 428049 : tolerance = MIN(1.0e-8_dp, 64.0_dp*EPSILON(1.0_dp)*MAX(1.0_dp, ABS(s)))
126 428049 : IF (ABS(s - half_boundary) <= tolerance) THEN
127 : ! Select the lower side for every lattice-equivalent half-cell boundary.
128 91973 : image_shift = half_boundary + 0.5_dp
129 : END IF
130 :
131 428049 : END FUNCTION pbc_image_shift
132 :
133 : ! **************************************************************************************************
134 : !> \brief Clone cell variable
135 : !> \param cell_in Cell variable to be clone
136 : !> \param cell_out Cloned cell variable
137 : !> \param tag Optional new tag for cloned cell variable
138 : !> \par History
139 : !> - Optional tag added (17.05.2023, MK)
140 : ! **************************************************************************************************
141 70253 : SUBROUTINE cell_clone(cell_in, cell_out, tag)
142 :
143 : TYPE(cell_type), POINTER :: cell_in, cell_out
144 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
145 :
146 70253 : cell_out = cell_in
147 70253 : cell_out%ref_count = 1
148 11941 : IF (PRESENT(tag)) cell_out%tag = tag
149 :
150 70253 : END SUBROUTINE cell_clone
151 :
152 : ! **************************************************************************************************
153 : !> \brief Copy cell variable
154 : !> \param cell_in Cell variable to be copied
155 : !> \param cell_out Copy of cell variable
156 : !> \param tag Optional new tag
157 : !> \par History
158 : !> - Optional tag added (17.05.2023, MK)
159 : ! **************************************************************************************************
160 233337 : SUBROUTINE cell_copy(cell_in, cell_out, tag)
161 :
162 : TYPE(cell_type), POINTER :: cell_in, cell_out
163 : CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: tag
164 :
165 233337 : cell_out%deth = cell_in%deth
166 1866696 : cell_out%perd = cell_in%perd
167 6066762 : cell_out%hmat = cell_in%hmat
168 6066762 : cell_out%h_inv = cell_in%h_inv
169 233337 : cell_out%input_cell_canonicalized = cell_in%input_cell_canonicalized
170 6066762 : cell_out%input_hmat = cell_in%input_hmat
171 6066762 : cell_out%input_to_canonical = cell_in%input_to_canonical
172 6066762 : cell_out%input_recip_to_canonical = cell_in%input_recip_to_canonical
173 233337 : cell_out%orthorhombic = cell_in%orthorhombic
174 233337 : cell_out%symmetry_id = cell_in%symmetry_id
175 233337 : IF (PRESENT(tag)) THEN
176 13677 : cell_out%tag = tag
177 : ELSE
178 219660 : cell_out%tag = cell_in%tag
179 : END IF
180 :
181 233337 : END SUBROUTINE cell_copy
182 :
183 : ! **************************************************************************************************
184 : !> \brief Read cell info from a line (parsed from a file)
185 : !> \param input_line ...
186 : !> \param cell_itimes ...
187 : !> \param cell_time ...
188 : !> \param h ...
189 : !> \param vol ...
190 : !> \date 19.02.2008
191 : !> \author Teodoro Laino [tlaino] - University of Zurich
192 : !> \version 1.0
193 : ! **************************************************************************************************
194 344 : SUBROUTINE parse_cell_line(input_line, cell_itimes, cell_time, h, vol)
195 :
196 : CHARACTER(LEN=*), INTENT(IN) :: input_line
197 : INTEGER, INTENT(OUT) :: cell_itimes
198 : REAL(KIND=dp), INTENT(OUT) :: cell_time
199 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: h
200 : REAL(KIND=dp), INTENT(OUT) :: vol
201 :
202 : INTEGER :: i, j
203 :
204 344 : READ (input_line, *) cell_itimes, cell_time, &
205 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
206 1376 : DO i = 1, 3
207 4472 : DO j = 1, 3
208 4128 : h(j, i) = cp_unit_to_cp2k(h(j, i), "angstrom")
209 : END DO
210 : END DO
211 :
212 344 : END SUBROUTINE parse_cell_line
213 :
214 : ! **************************************************************************************************
215 : !> \brief Get informations about a simulation cell.
216 : !> \param cell ...
217 : !> \param alpha ...
218 : !> \param beta ...
219 : !> \param gamma ...
220 : !> \param deth ...
221 : !> \param orthorhombic ...
222 : !> \param abc ...
223 : !> \param periodic ...
224 : !> \param h ...
225 : !> \param h_inv ...
226 : !> \param symmetry_id ...
227 : !> \param tag ...
228 : !> \date 16.01.2002
229 : !> \author Matthias Krack
230 : !> \version 1.0
231 : ! **************************************************************************************************
232 142934429 : SUBROUTINE get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, &
233 : h, h_inv, symmetry_id, tag)
234 :
235 : TYPE(cell_type), POINTER :: cell
236 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha, beta, gamma, deth
237 : LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic
238 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: abc
239 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: periodic
240 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
241 : OPTIONAL :: h, h_inv
242 : INTEGER, INTENT(OUT), OPTIONAL :: symmetry_id
243 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: tag
244 :
245 0 : CPASSERT(ASSOCIATED(cell))
246 :
247 142934429 : IF (PRESENT(deth)) deth = cell%deth ! the volume
248 142934429 : IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic
249 548618954 : IF (PRESENT(periodic)) periodic(:) = cell%perd(:)
250 143236985 : IF (PRESENT(h)) h(:, :) = cell%hmat(:, :)
251 142935317 : IF (PRESENT(h_inv)) h_inv(:, :) = cell%h_inv(:, :)
252 :
253 : ! Calculate the lengths of the cell vectors a, b, and c
254 142934429 : IF (PRESENT(abc)) THEN
255 : abc(1) = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
256 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
257 7488575 : cell%hmat(3, 1)*cell%hmat(3, 1))
258 : abc(2) = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
259 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
260 7488575 : cell%hmat(3, 2)*cell%hmat(3, 2))
261 : abc(3) = SQRT(cell%hmat(1, 3)*cell%hmat(1, 3) + &
262 : cell%hmat(2, 3)*cell%hmat(2, 3) + &
263 7488575 : cell%hmat(3, 3)*cell%hmat(3, 3))
264 : END IF
265 :
266 : ! Angles between the cell vectors a, b, and c
267 : ! alpha = <(b,c)
268 142934429 : IF (PRESENT(alpha)) alpha = angle(cell%hmat(:, 2), cell%hmat(:, 3))*degree
269 : ! beta = <(a,c)
270 142934429 : IF (PRESENT(beta)) beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))*degree
271 : ! gamma = <(a,b)
272 142934429 : IF (PRESENT(gamma)) gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))*degree
273 142934429 : IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id
274 142934429 : IF (PRESENT(tag)) tag = cell%tag
275 :
276 142934429 : END SUBROUTINE get_cell
277 :
278 : ! **************************************************************************************************
279 : !> \brief Transform a Cartesian real-space vector from the user input cell frame
280 : !> into CP2K's canonical internal cell frame.
281 : !> \param cell ...
282 : !> \param vector ...
283 : ! **************************************************************************************************
284 787869 : SUBROUTINE cell_transform_input_cartesian(cell, vector)
285 :
286 : TYPE(cell_type), POINTER :: cell
287 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: vector
288 :
289 787869 : CPASSERT(ASSOCIATED(cell))
290 :
291 954261 : IF (cell%input_cell_canonicalized) vector = MATMUL(cell%input_to_canonical, vector)
292 :
293 787869 : END SUBROUTINE cell_transform_input_cartesian
294 :
295 : ! **************************************************************************************************
296 : !> \brief Transform a Cartesian reciprocal-space vector from the user input cell
297 : !> frame into CP2K's canonical internal cell frame.
298 : !> \param cell ...
299 : !> \param vector ...
300 : ! **************************************************************************************************
301 10 : SUBROUTINE cell_transform_input_reciprocal(cell, vector)
302 :
303 : TYPE(cell_type), POINTER :: cell
304 : REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: vector
305 :
306 10 : CPASSERT(ASSOCIATED(cell))
307 :
308 10 : IF (cell%input_cell_canonicalized) vector = MATMUL(cell%input_recip_to_canonical, vector)
309 :
310 10 : END SUBROUTINE cell_transform_input_reciprocal
311 :
312 : ! **************************************************************************************************
313 : !> \brief Calculate the distance between two lattice planes as defined by
314 : !> a triple of Miller indices (hkl).
315 : !> \param h ...
316 : !> \param k ...
317 : !> \param l ...
318 : !> \param cell ...
319 : !> \return ...
320 : !> \date 18.11.2004
321 : !> \author Matthias Krack
322 : !> \version 1.0
323 : ! **************************************************************************************************
324 7380269 : FUNCTION plane_distance(h, k, l, cell) RESULT(distance)
325 :
326 : INTEGER, INTENT(IN) :: h, k, l
327 : TYPE(cell_type), POINTER :: cell
328 : REAL(KIND=dp) :: distance
329 :
330 : REAL(KIND=dp) :: a, alpha, b, beta, c, cosa, cosb, cosg, &
331 : d, gamma, x, y, z
332 : REAL(KIND=dp), DIMENSION(3) :: abc
333 :
334 7380269 : x = REAL(h, KIND=dp)
335 7380269 : y = REAL(k, KIND=dp)
336 7380269 : z = REAL(l, KIND=dp)
337 :
338 7380269 : CALL get_cell(cell=cell, abc=abc)
339 :
340 7380269 : a = abc(1)
341 7380269 : b = abc(2)
342 7380269 : c = abc(3)
343 :
344 7380269 : IF (cell%orthorhombic) THEN
345 :
346 7165642 : d = (x/a)**2 + (y/b)**2 + (z/c)**2
347 :
348 : ELSE
349 :
350 : CALL get_cell(cell=cell, &
351 : alpha=alpha, &
352 : beta=beta, &
353 214627 : gamma=gamma)
354 :
355 214627 : alpha = alpha/degree
356 214627 : beta = beta/degree
357 214627 : gamma = gamma/degree
358 :
359 214627 : cosa = COS(alpha)
360 214627 : cosb = COS(beta)
361 214627 : cosg = COS(gamma)
362 :
363 : d = ((x*b*c*SIN(alpha))**2 + &
364 : (y*c*a*SIN(beta))**2 + &
365 : (z*a*b*SIN(gamma))**2 + &
366 : 2.0_dp*a*b*c*(x*y*c*(cosa*cosb - cosg) + &
367 : z*x*b*(cosg*cosa - cosb) + &
368 : y*z*a*(cosb*cosg - cosa)))/ &
369 : ((a*b*c)**2*(1.0_dp - cosa**2 - cosb**2 - cosg**2 + &
370 214627 : 2.0_dp*cosa*cosb*cosg))
371 :
372 : END IF
373 :
374 7380269 : distance = 1.0_dp/SQRT(d)
375 :
376 7380269 : END FUNCTION plane_distance
377 :
378 : ! **************************************************************************************************
379 : !> \brief Apply the periodic boundary conditions defined by a simulation
380 : !> cell to a position vector r.
381 : !> \param r ...
382 : !> \param cell ...
383 : !> \return ...
384 : !> \date 16.01.2002
385 : !> \author Matthias Krack
386 : !> \version 1.0
387 : ! **************************************************************************************************
388 411309350 : FUNCTION pbc1(r, cell) RESULT(r_pbc)
389 :
390 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
391 : TYPE(cell_type), POINTER :: cell
392 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
393 :
394 : REAL(KIND=dp), DIMENSION(3) :: s
395 :
396 411309350 : CPASSERT(ASSOCIATED(cell))
397 :
398 411309350 : IF (cell%orthorhombic) THEN
399 385325278 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*ANINT(cell%h_inv(1, 1)*r(1))
400 385325278 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*ANINT(cell%h_inv(2, 2)*r(2))
401 385325278 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*ANINT(cell%h_inv(3, 3)*r(3))
402 : ELSE
403 25984072 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
404 25984072 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
405 25984072 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
406 25984072 : s(1) = s(1) - cell%perd(1)*ANINT(s(1))
407 25984072 : s(2) = s(2) - cell%perd(2)*ANINT(s(2))
408 25984072 : s(3) = s(3) - cell%perd(3)*ANINT(s(3))
409 25984072 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
410 25984072 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
411 25984072 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
412 : END IF
413 :
414 411309350 : END FUNCTION pbc1
415 :
416 : ! **************************************************************************************************
417 : !> \brief Distribute an atom-pair contribution equally among its nearest candidate images.
418 : !> Images within 1.0E-6 of the minimum distance receive normalized equal weights.
419 : !> \param ra Position of the first atom.
420 : !> \param rb Position of the second atom.
421 : !> \param cell Simulation cell used to wrap the positions and translate the images.
422 : !> \param index_to_cell Candidate lattice translations, indexed (direction, image).
423 : !> \param weights Normalized image weights; non-nearest images receive zero.
424 : ! **************************************************************************************************
425 42641 : SUBROUTINE minimum_image_weights(ra, rb, cell, index_to_cell, weights)
426 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
427 : TYPE(cell_type), INTENT(IN), POINTER :: cell
428 : INTEGER, DIMENSION(:, :), INTENT(IN) :: index_to_cell
429 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: weights
430 :
431 : INTEGER :: i_image
432 : REAL(KIND=dp), DIMENSION(3) :: ra_pbc, rab, rb_pbc
433 85282 : REAL(KIND=dp), DIMENSION(SIZE(index_to_cell, 2)) :: distances
434 :
435 42641 : CPASSERT(SIZE(index_to_cell, 1) == 3)
436 42641 : CPASSERT(SIZE(index_to_cell, 2) > 0)
437 42641 : CPASSERT(SIZE(weights) == SIZE(index_to_cell, 2))
438 :
439 42641 : ra_pbc = pbc(ra, cell)
440 42641 : rb_pbc = pbc(rb, cell)
441 355326 : DO i_image = 1, SIZE(index_to_cell, 2)
442 6253700 : rab = ra_pbc - (rb_pbc + MATMUL(cell%hmat, REAL(index_to_cell(:, i_image), KIND=dp)))
443 1293381 : distances(i_image) = NORM2(rab)
444 : END DO
445 :
446 668011 : weights = MERGE(1.0_dp, 0.0_dp, distances <= MINVAL(distances) + 1.0E-6_dp)
447 668011 : weights = weights/SUM(weights)
448 :
449 42641 : END SUBROUTINE minimum_image_weights
450 :
451 : ! **************************************************************************************************
452 : !> \brief Apply a stable periodic-image convention for k-point Bloch gauges.
453 : !> \param r ...
454 : !> \param cell ...
455 : !> \return r_pbc ...
456 : ! **************************************************************************************************
457 142683 : FUNCTION pbc_stable(r, cell) RESULT(r_pbc)
458 :
459 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
460 : TYPE(cell_type), POINTER :: cell
461 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
462 :
463 : REAL(KIND=dp), DIMENSION(3) :: s
464 :
465 142683 : CPASSERT(ASSOCIATED(cell))
466 :
467 142683 : IF (cell%orthorhombic) THEN
468 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
469 118006 : pbc_image_shift(cell%h_inv(1, 1)*r(1))
470 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
471 118006 : pbc_image_shift(cell%h_inv(2, 2)*r(2))
472 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
473 118006 : pbc_image_shift(cell%h_inv(3, 3)*r(3))
474 : ELSE
475 24677 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
476 24677 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
477 24677 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
478 24677 : s(1) = s(1) - cell%perd(1)*pbc_image_shift(s(1))
479 24677 : s(2) = s(2) - cell%perd(2)*pbc_image_shift(s(2))
480 24677 : s(3) = s(3) - cell%perd(3)*pbc_image_shift(s(3))
481 24677 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
482 24677 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
483 24677 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
484 : END IF
485 :
486 142683 : END FUNCTION pbc_stable
487 :
488 : ! **************************************************************************************************
489 : !> \brief Apply the periodic boundary conditions defined by a simulation
490 : !> cell to a position vector r subtracting nl from the periodic images
491 : !> \param r ...
492 : !> \param cell ...
493 : !> \param nl ...
494 : !> \return ...
495 : !> \date 16.01.2002
496 : !> \author Matthias Krack
497 : !> \version 1.0
498 : ! **************************************************************************************************
499 0 : FUNCTION pbc2(r, cell, nl) RESULT(r_pbc)
500 :
501 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
502 : TYPE(cell_type), POINTER :: cell
503 : INTEGER, DIMENSION(3), INTENT(IN) :: nl
504 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
505 :
506 : REAL(KIND=dp), DIMENSION(3) :: s
507 :
508 0 : CPASSERT(ASSOCIATED(cell))
509 :
510 0 : IF (cell%orthorhombic) THEN
511 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
512 0 : REAL(NINT(cell%h_inv(1, 1)*r(1)) - nl(1), dp)
513 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
514 0 : REAL(NINT(cell%h_inv(2, 2)*r(2)) - nl(2), dp)
515 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
516 0 : REAL(NINT(cell%h_inv(3, 3)*r(3)) - nl(3), dp)
517 : ELSE
518 0 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
519 0 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
520 0 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
521 0 : s(1) = s(1) - cell%perd(1)*REAL(NINT(s(1)) - nl(1), dp)
522 0 : s(2) = s(2) - cell%perd(2)*REAL(NINT(s(2)) - nl(2), dp)
523 0 : s(3) = s(3) - cell%perd(3)*REAL(NINT(s(3)) - nl(3), dp)
524 0 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
525 0 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
526 0 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
527 : END IF
528 :
529 0 : END FUNCTION pbc2
530 :
531 : ! **************************************************************************************************
532 : !> \brief Apply the periodic boundary conditions defined by the simulation
533 : !> cell cell to the vector pointing from atom a to atom b.
534 : !> \param ra ...
535 : !> \param rb ...
536 : !> \param cell ...
537 : !> \return ...
538 : !> \date 11.03.2004
539 : !> \author Matthias Krack
540 : !> \version 1.0
541 : ! **************************************************************************************************
542 134996946 : FUNCTION pbc3(ra, rb, cell) RESULT(rab_pbc)
543 :
544 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
545 : TYPE(cell_type), POINTER :: cell
546 : REAL(KIND=dp), DIMENSION(3) :: rab_pbc
547 :
548 : INTEGER :: icell, jcell, kcell
549 : INTEGER, DIMENSION(3) :: periodic
550 : REAL(KIND=dp) :: rab2, rab2_pbc
551 : REAL(KIND=dp), DIMENSION(3) :: r, ra_pbc, rab, rb_image, rb_pbc, s2r
552 :
553 134996946 : CALL get_cell(cell=cell, periodic=periodic)
554 :
555 134996946 : ra_pbc(:) = pbc(ra(:), cell)
556 134996946 : rb_pbc(:) = pbc(rb(:), cell)
557 :
558 134996946 : rab2_pbc = HUGE(1.0_dp)
559 :
560 533373680 : DO icell = -periodic(1), periodic(1)
561 1721891996 : DO jcell = -periodic(2), periodic(2)
562 5145821168 : DO kcell = -periodic(3), periodic(3)
563 14235704472 : r = REAL([icell, jcell, kcell], dp)
564 3558926118 : CALL scaled_to_real(s2r, r, cell)
565 14235704472 : rb_image(:) = rb_pbc(:) + s2r
566 14235704472 : rab(:) = rb_image(:) - ra_pbc(:)
567 3558926118 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
568 4747444434 : IF (rab2 < rab2_pbc) THEN
569 2797605292 : rab2_pbc = rab2
570 2797605292 : rab_pbc(:) = rab(:)
571 : END IF
572 : END DO
573 : END DO
574 : END DO
575 :
576 134996946 : END FUNCTION pbc3
577 :
578 : !if positive_range == true, r(i) (or s(i)) in range [0, hmat(i,i)],
579 : !else, r(i) (s(i)) in range [-hmat(i,i)/2, hmat(i,i)/2]
580 : ! **************************************************************************************************
581 : !> \brief ...
582 : !> \param r ...
583 : !> \param cell ...
584 : !> \param positive_range ...
585 : !> \return ...
586 : ! **************************************************************************************************
587 307505 : FUNCTION pbc4(r, cell, positive_range) RESULT(r_pbc)
588 :
589 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
590 : TYPE(cell_type), POINTER :: cell
591 : LOGICAL :: positive_range
592 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
593 :
594 : REAL(KIND=dp), DIMENSION(3) :: s
595 :
596 307505 : CPASSERT(ASSOCIATED(cell))
597 :
598 307505 : IF (positive_range) THEN
599 307505 : IF (cell%orthorhombic) THEN
600 75229 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*FLOOR(cell%h_inv(1, 1)*r(1))
601 75229 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*FLOOR(cell%h_inv(2, 2)*r(2))
602 75229 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*FLOOR(cell%h_inv(3, 3)*r(3))
603 : ELSE
604 232276 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
605 232276 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
606 232276 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
607 232276 : s(1) = s(1) - cell%perd(1)*FLOOR(s(1))
608 232276 : s(2) = s(2) - cell%perd(2)*FLOOR(s(2))
609 232276 : s(3) = s(3) - cell%perd(3)*FLOOR(s(3))
610 232276 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
611 232276 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
612 232276 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
613 : END IF
614 : ELSE
615 0 : r_pbc = pbc1(r, cell)
616 : END IF
617 :
618 307505 : END FUNCTION pbc4
619 :
620 : ! **************************************************************************************************
621 : !> \brief Transform real to scaled cell coordinates.
622 : !> s=h_inv*r
623 : !> \param s ...
624 : !> \param r ...
625 : !> \param cell ...
626 : !> \date 16.01.2002
627 : !> \author Matthias Krack
628 : !> \version 1.0
629 : ! **************************************************************************************************
630 111320013 : SUBROUTINE real_to_scaled(s, r, cell)
631 :
632 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: s
633 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
634 : TYPE(cell_type), POINTER :: cell
635 :
636 111320013 : CPASSERT(ASSOCIATED(cell))
637 :
638 111320013 : IF (cell%orthorhombic) THEN
639 100324462 : s(1) = cell%h_inv(1, 1)*r(1)
640 100324462 : s(2) = cell%h_inv(2, 2)*r(2)
641 100324462 : s(3) = cell%h_inv(3, 3)*r(3)
642 : ELSE
643 10995551 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
644 10995551 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
645 10995551 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
646 : END IF
647 :
648 111320013 : END SUBROUTINE real_to_scaled
649 :
650 : ! **************************************************************************************************
651 : !> \brief Transform scaled cell coordinates real coordinates.
652 : !> r=h*s
653 : !> \param r ...
654 : !> \param s ...
655 : !> \param cell ...
656 : !> \date 16.01.2002
657 : !> \author Matthias Krack
658 : !> \version 1.0
659 : ! **************************************************************************************************
660 3730327966 : SUBROUTINE scaled_to_real(r, s, cell)
661 :
662 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: r
663 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s
664 : TYPE(cell_type), POINTER :: cell
665 :
666 3730327966 : CPASSERT(ASSOCIATED(cell))
667 :
668 3730327966 : IF (cell%orthorhombic) THEN
669 3454027005 : r(1) = cell%hmat(1, 1)*s(1)
670 3454027005 : r(2) = cell%hmat(2, 2)*s(2)
671 3454027005 : r(3) = cell%hmat(3, 3)*s(3)
672 : ELSE
673 276300961 : r(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
674 276300961 : r(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
675 276300961 : r(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
676 : END IF
677 :
678 3730327966 : END SUBROUTINE scaled_to_real
679 : ! **************************************************************************************************
680 : !> \brief retains the given cell (see doc/ReferenceCounting.html)
681 : !> \param cell the cell to retain
682 : !> \par History
683 : !> 09.2003 created [fawzi]
684 : !> \author Fawzi Mohamed
685 : ! **************************************************************************************************
686 80768 : SUBROUTINE cell_retain(cell)
687 :
688 : TYPE(cell_type), POINTER :: cell
689 :
690 80768 : CPASSERT(ASSOCIATED(cell))
691 80768 : CPASSERT(cell%ref_count > 0)
692 80768 : cell%ref_count = cell%ref_count + 1
693 :
694 80768 : END SUBROUTINE cell_retain
695 :
696 : ! **************************************************************************************************
697 : !> \brief releases the given cell (see doc/ReferenceCounting.html)
698 : !> \param cell the cell to release
699 : !> \par History
700 : !> 09.2003 created [fawzi]
701 : !> \author Fawzi Mohamed
702 : ! **************************************************************************************************
703 216135 : SUBROUTINE cell_release(cell)
704 :
705 : TYPE(cell_type), POINTER :: cell
706 :
707 216135 : IF (ASSOCIATED(cell)) THEN
708 165160 : CPASSERT(cell%ref_count > 0)
709 165160 : cell%ref_count = cell%ref_count - 1
710 165160 : IF (cell%ref_count == 0) THEN
711 84393 : DEALLOCATE (cell)
712 : END IF
713 165160 : NULLIFY (cell)
714 : END IF
715 :
716 216135 : END SUBROUTINE cell_release
717 :
718 : #if defined (__PLUMED2)
719 : ! **************************************************************************************************
720 : !> \brief For the interface with plumed, pass a cell pointer and retrieve it
721 : !> later. It's a hack, but avoids passing the cell back and forth
722 : !> across the Fortran/C++ interface
723 : !> \param cell ...
724 : !> \param set ...
725 : !> \date 28.02.2013
726 : !> \author RK
727 : !> \version 1.0
728 : ! **************************************************************************************************
729 2 : SUBROUTINE pbc_cp2k_plumed_getset_cell(cell, set)
730 :
731 : TYPE(cell_type), POINTER :: cell
732 : LOGICAL :: set
733 :
734 : TYPE(cell_type), POINTER, SAVE :: stored_cell
735 :
736 2 : IF (set) THEN
737 2 : stored_cell => cell
738 : ELSE
739 0 : cell => stored_cell
740 : END IF
741 :
742 2 : END SUBROUTINE pbc_cp2k_plumed_getset_cell
743 : #endif
744 :
745 312685 : END MODULE cell_types
|
