Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2022 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 : cp_units_rad
18 : USE kinds, ONLY: dp
19 : USE mathconstants, ONLY: degree,&
20 : sqrt3
21 : USE mathlib, ONLY: angle,&
22 : det_3x3,&
23 : inv_3x3
24 : #include "../base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_types'
31 :
32 : INTEGER, SAVE, PRIVATE :: last_cell_id = 0
33 :
34 : ! Impose cell symmetry
35 : INTEGER, PARAMETER, PUBLIC :: cell_sym_none = 0, &
36 : cell_sym_triclinic = 1, &
37 : cell_sym_monoclinic = 2, &
38 : cell_sym_monoclinic_gamma_ab = 3, &
39 : cell_sym_orthorhombic = 4, &
40 : cell_sym_tetragonal_ab = 5, &
41 : cell_sym_tetragonal_ac = 6, &
42 : cell_sym_tetragonal_bc = 7, &
43 : cell_sym_rhombohedral = 8, &
44 : cell_sym_hexagonal = 9, &
45 : cell_sym_cubic = 10
46 :
47 : INTEGER, PARAMETER, PUBLIC :: use_perd_x = 0, &
48 : use_perd_y = 1, &
49 : use_perd_z = 2, &
50 : use_perd_xy = 3, &
51 : use_perd_xz = 4, &
52 : use_perd_yz = 5, &
53 : use_perd_xyz = 6, &
54 : use_perd_none = 7
55 :
56 : ! **************************************************************************************************
57 : !> \brief Type defining parameters related to the simulation cell
58 : !> \version 1.0
59 : ! **************************************************************************************************
60 : TYPE cell_type
61 : INTEGER :: id_nr, ref_count, symmetry_id
62 : LOGICAL :: orthorhombic ! actually means a diagonal hmat
63 : REAL(KIND=dp) :: deth
64 : INTEGER, DIMENSION(3) :: perd
65 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat, h_inv
66 : END TYPE cell_type
67 :
68 : TYPE cell_p_type
69 : TYPE(cell_type), POINTER :: cell
70 : END TYPE cell_p_type
71 :
72 : ! Public data types
73 : PUBLIC :: cell_type, cell_p_type
74 :
75 : ! Public subroutines
76 : PUBLIC :: get_cell, get_cell_param, init_cell, &
77 : cell_create, cell_retain, cell_release, &
78 : cell_clone, cell_copy, parse_cell_line, set_cell_param
79 :
80 : #if defined (__PLUMED2)
81 : PUBLIC :: pbc_cp2k_plumed_getset_cell
82 : #endif
83 :
84 : ! Public functions
85 : PUBLIC :: plane_distance, pbc, real_to_scaled, scaled_to_real
86 :
87 : INTERFACE pbc
88 : MODULE PROCEDURE pbc1, pbc2, pbc3, pbc4
89 : END INTERFACE
90 :
91 : CONTAINS
92 :
93 : ! **************************************************************************************************
94 : !> \brief ...
95 : !> \param cell_in ...
96 : !> \param cell_out ...
97 : ! **************************************************************************************************
98 27401 : SUBROUTINE cell_clone(cell_in, cell_out)
99 :
100 : TYPE(cell_type), INTENT(IN) :: cell_in
101 : TYPE(cell_type), INTENT(OUT) :: cell_out
102 :
103 27401 : cell_out%deth = cell_in%deth
104 109604 : cell_out%perd = cell_in%perd
105 356213 : cell_out%hmat = cell_in%hmat
106 356213 : cell_out%h_inv = cell_in%h_inv
107 27401 : cell_out%orthorhombic = cell_in%orthorhombic
108 27401 : cell_out%symmetry_id = cell_in%symmetry_id
109 27401 : cell_out%ref_count = 1
110 27401 : last_cell_id = last_cell_id + 1
111 27401 : cell_out%id_nr = last_cell_id
112 :
113 27401 : END SUBROUTINE cell_clone
114 :
115 : ! **************************************************************************************************
116 : !> \brief ...
117 : !> \param cell_in ...
118 : !> \param cell_out ...
119 : ! **************************************************************************************************
120 227034 : SUBROUTINE cell_copy(cell_in, cell_out)
121 :
122 : TYPE(cell_type), INTENT(IN) :: cell_in
123 : TYPE(cell_type), INTENT(INOUT) :: cell_out
124 :
125 227034 : cell_out%deth = cell_in%deth
126 908136 : cell_out%perd = cell_in%perd
127 2951442 : cell_out%hmat = cell_in%hmat
128 2951442 : cell_out%h_inv = cell_in%h_inv
129 227034 : cell_out%orthorhombic = cell_in%orthorhombic
130 227034 : cell_out%symmetry_id = cell_in%symmetry_id
131 :
132 227034 : END SUBROUTINE cell_copy
133 :
134 : ! **************************************************************************************************
135 : !> \brief Read cell info from a line (parsed from a file)
136 : !> \param input_line ...
137 : !> \param cell_itimes ...
138 : !> \param cell_time ...
139 : !> \param h ...
140 : !> \param vol ...
141 : !> \date 19.02.2008
142 : !> \author Teodoro Laino [tlaino] - University of Zurich
143 : !> \version 1.0
144 : ! **************************************************************************************************
145 344 : SUBROUTINE parse_cell_line(input_line, cell_itimes, cell_time, h, vol)
146 : CHARACTER(LEN=*), INTENT(IN) :: input_line
147 : INTEGER, INTENT(OUT) :: cell_itimes
148 : REAL(KIND=dp), INTENT(OUT) :: cell_time
149 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: h
150 : REAL(KIND=dp), INTENT(OUT) :: vol
151 :
152 : INTEGER :: i, j
153 :
154 344 : READ (input_line, *) cell_itimes, cell_time, &
155 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
156 1376 : DO i = 1, 3
157 4472 : DO j = 1, 3
158 4128 : h(j, i) = cp_unit_to_cp2k(h(j, i), "angstrom")
159 : END DO
160 : END DO
161 :
162 344 : END SUBROUTINE parse_cell_line
163 :
164 : ! **************************************************************************************************
165 : !> \brief Get informations about a simulation cell.
166 : !> \param cell ...
167 : !> \param alpha ...
168 : !> \param beta ...
169 : !> \param gamma ...
170 : !> \param deth ...
171 : !> \param orthorhombic ...
172 : !> \param abc ...
173 : !> \param periodic ...
174 : !> \param h ...
175 : !> \param h_inv ...
176 : !> \param id_nr ...
177 : !> \param symmetry_id ...
178 : !> \date 16.01.2002
179 : !> \author Matthias Krack
180 : !> \version 1.0
181 : ! **************************************************************************************************
182 124334530 : SUBROUTINE get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, &
183 : h, h_inv, id_nr, symmetry_id)
184 :
185 : TYPE(cell_type), INTENT(IN) :: cell
186 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha, beta, gamma, deth
187 : LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic
188 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: abc
189 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: periodic
190 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT), &
191 : OPTIONAL :: h, h_inv
192 : INTEGER, INTENT(OUT), OPTIONAL :: id_nr, symmetry_id
193 :
194 124334530 : IF (PRESENT(deth)) deth = cell%deth ! the volume
195 124334530 : IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic
196 477578146 : IF (PRESENT(periodic)) periodic(:) = cell%perd(:)
197 124349914 : IF (PRESENT(h)) h(:, :) = cell%hmat(:, :)
198 124334626 : IF (PRESENT(h_inv)) h_inv(:, :) = cell%h_inv(:, :)
199 :
200 : ! Calculate the lengths of the cell vectors a, b, and c
201 124334530 : IF (PRESENT(abc)) THEN
202 : abc(1) = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
203 : cell%hmat(2, 1)*cell%hmat(2, 1) + &
204 6547009 : cell%hmat(3, 1)*cell%hmat(3, 1))
205 : abc(2) = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
206 : cell%hmat(2, 2)*cell%hmat(2, 2) + &
207 6547009 : cell%hmat(3, 2)*cell%hmat(3, 2))
208 : abc(3) = SQRT(cell%hmat(1, 3)*cell%hmat(1, 3) + &
209 : cell%hmat(2, 3)*cell%hmat(2, 3) + &
210 6547009 : cell%hmat(3, 3)*cell%hmat(3, 3))
211 : END IF
212 :
213 : ! Angles between the cell vectors a, b, and c
214 : ! alpha = <(b,c)
215 124334530 : IF (PRESENT(alpha)) alpha = angle(cell%hmat(:, 2), cell%hmat(:, 3))*degree
216 : ! beta = <(a,c)
217 124334530 : IF (PRESENT(beta)) beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))*degree
218 : ! gamma = <(a,b)
219 124334530 : IF (PRESENT(gamma)) gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))*degree
220 124334530 : IF (PRESENT(id_nr)) id_nr = cell%id_nr
221 124334530 : IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id
222 :
223 124334530 : END SUBROUTINE get_cell
224 :
225 : ! **************************************************************************************************
226 : !> \brief Access internal type variables
227 : !> \param cell ...
228 : !> \param cell_length ...
229 : !> \param cell_angle ...
230 : !> \param units_angle ...
231 : !> \param periodic ...
232 : !> \date 04.04.2002
233 : !> \author Matthias Krack
234 : !> \version 1.0
235 : ! **************************************************************************************************
236 56 : SUBROUTINE get_cell_param(cell, cell_length, cell_angle, units_angle, periodic)
237 :
238 : TYPE(cell_type), INTENT(IN) :: cell
239 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: cell_length
240 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT), OPTIONAL :: cell_angle
241 : INTEGER, INTENT(IN), OPTIONAL :: units_angle
242 : INTEGER, DIMENSION(3), INTENT(OUT), OPTIONAL :: periodic
243 :
244 : REAL(KIND=dp) :: alpha, beta, gamma
245 :
246 56 : CALL get_cell(cell=cell, abc=cell_length)
247 :
248 56 : IF (PRESENT(cell_angle)) THEN
249 56 : CALL get_cell(cell=cell, alpha=alpha, beta=beta, gamma=gamma)
250 224 : cell_angle(:) = (/alpha, beta, gamma/)
251 56 : IF (PRESENT(units_angle)) THEN
252 280 : IF (units_angle == cp_units_rad) cell_angle = cell_angle/degree
253 : END IF
254 : END IF
255 :
256 56 : IF (PRESENT(periodic)) CALL get_cell(cell=cell, periodic=periodic)
257 :
258 56 : END SUBROUTINE get_cell_param
259 :
260 : ! **************************************************************************************************
261 : !> \brief Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma)
262 : !> using the convention: a parallel to the x axis, b in the x-y plane and
263 : !> and c univoquely determined; gamma is the angle between a and b; beta
264 : !> is the angle between c and a and alpha is the angle between c and b
265 : !> \param cell ...
266 : !> \param cell_length ...
267 : !> \param cell_angle ...
268 : !> \param periodic ...
269 : !> \param do_init_cell ...
270 : !> \date 03.2008
271 : !> \author Teodoro Laino
272 : ! **************************************************************************************************
273 14477 : SUBROUTINE set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
274 :
275 : TYPE(cell_type), INTENT(INOUT) :: cell
276 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: cell_length, cell_angle
277 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
278 : LOGICAL, INTENT(IN) :: do_init_cell
279 :
280 : REAL(KIND=dp) :: cos_alpha, cos_beta, cos_gamma, eps, &
281 : sin_gamma
282 :
283 57908 : CPASSERT(ALL(cell_angle /= 0.0_dp))
284 14477 : eps = EPSILON(0.0_dp)
285 14477 : cos_gamma = COS(cell_angle(3)); IF (ABS(cos_gamma) < eps) cos_gamma = 0.0_dp
286 14477 : IF (ABS(ABS(cos_gamma) - 1.0_dp) < eps) cos_gamma = SIGN(1.0_dp, cos_gamma)
287 14477 : sin_gamma = SIN(cell_angle(3)); IF (ABS(sin_gamma) < eps) sin_gamma = 0.0_dp
288 14477 : IF (ABS(ABS(sin_gamma) - 1.0_dp) < eps) sin_gamma = SIGN(1.0_dp, sin_gamma)
289 14477 : cos_beta = COS(cell_angle(2)); IF (ABS(cos_beta) < eps) cos_beta = 0.0_dp
290 14477 : IF (ABS(ABS(cos_beta) - 1.0_dp) < eps) cos_beta = SIGN(1.0_dp, cos_beta)
291 14477 : cos_alpha = COS(cell_angle(1)); IF (ABS(cos_alpha) < eps) cos_alpha = 0.0_dp
292 14477 : IF (ABS(ABS(cos_alpha) - 1.0_dp) < eps) cos_alpha = SIGN(1.0_dp, cos_alpha)
293 :
294 57908 : cell%hmat(:, 1) = (/1.0_dp, 0.0_dp, 0.0_dp/)
295 57908 : cell%hmat(:, 2) = (/cos_gamma, sin_gamma, 0.0_dp/)
296 57908 : cell%hmat(:, 3) = (/cos_beta, (cos_alpha - cos_gamma*cos_beta)/sin_gamma, 0.0_dp/)
297 14477 : cell%hmat(3, 3) = SQRT(1.0_dp - cell%hmat(1, 3)**2 - cell%hmat(2, 3)**2)
298 :
299 57908 : cell%hmat(:, 1) = cell%hmat(:, 1)*cell_length(1)
300 57908 : cell%hmat(:, 2) = cell%hmat(:, 2)*cell_length(2)
301 57908 : cell%hmat(:, 3) = cell%hmat(:, 3)*cell_length(3)
302 :
303 14477 : IF (do_init_cell) THEN
304 66 : IF (PRESENT(periodic)) THEN
305 66 : CALL init_cell(cell=cell, periodic=periodic)
306 : ELSE
307 0 : CALL init_cell(cell=cell)
308 : END IF
309 : END IF
310 :
311 14477 : END SUBROUTINE set_cell_param
312 :
313 : ! **************************************************************************************************
314 : !> \brief Initialise/readjust a simulation cell after hmat has been changed
315 : !> \param cell ...
316 : !> \param hmat ...
317 : !> \param periodic ...
318 : !> \date 16.01.2002
319 : !> \author Matthias Krack
320 : !> \version 1.0
321 : ! **************************************************************************************************
322 337359 : SUBROUTINE init_cell(cell, hmat, periodic)
323 :
324 : TYPE(cell_type), INTENT(INOUT) :: cell
325 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
326 : OPTIONAL :: hmat
327 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
328 :
329 : REAL(KIND=dp), PARAMETER :: eps_hmat = 1.0E-14_dp
330 :
331 : INTEGER :: dim
332 : REAL(KIND=dp) :: a, acosa, acosah, acosgamma, alpha, &
333 : asina, asinah, asingamma, beta, gamma, &
334 : norm, norm_c
335 : REAL(KIND=dp), DIMENSION(3) :: abc
336 :
337 419799 : IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
338 337557 : IF (PRESENT(periodic)) cell%perd(:) = periodic(:)
339 :
340 337359 : cell%deth = ABS(det_3x3(cell%hmat))
341 :
342 337359 : IF (cell%deth < 1.0E-10_dp) THEN
343 : CALL cp_abort(__LOCATION__, &
344 : "An invalid set of cell vectors was specified. "// &
345 0 : "The determinant det(h) is too small")
346 : END IF
347 :
348 341397 : SELECT CASE (cell%symmetry_id)
349 : CASE (cell_sym_cubic, &
350 : cell_sym_tetragonal_ab, &
351 : cell_sym_tetragonal_ac, &
352 : cell_sym_tetragonal_bc, &
353 : cell_sym_orthorhombic)
354 4038 : CALL get_cell(cell=cell, abc=abc)
355 4038 : abc(2) = plane_distance(0, 1, 0, cell=cell)
356 4038 : abc(3) = plane_distance(0, 0, 1, cell=cell)
357 4038 : SELECT CASE (cell%symmetry_id)
358 : CASE (cell_sym_cubic)
359 7579 : abc(1:3) = SUM(abc(1:3))/3.0_dp
360 : CASE (cell_sym_tetragonal_ab, &
361 : cell_sym_tetragonal_ac, &
362 : cell_sym_tetragonal_bc)
363 1322 : SELECT CASE (cell%symmetry_id)
364 : CASE (cell_sym_tetragonal_ab)
365 1322 : a = 0.5_dp*(abc(1) + abc(2))
366 1322 : abc(1) = a
367 1322 : abc(2) = a
368 : CASE (cell_sym_tetragonal_ac)
369 661 : a = 0.5_dp*(abc(1) + abc(3))
370 661 : abc(1) = a
371 661 : abc(3) = a
372 : CASE (cell_sym_tetragonal_bc)
373 612 : a = 0.5_dp*(abc(2) + abc(3))
374 612 : abc(2) = a
375 612 : abc(3) = a
376 : END SELECT
377 : END SELECT
378 4038 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = 0.0_dp
379 4038 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
380 4038 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
381 : CASE (cell_sym_hexagonal)
382 842 : CALL get_cell(cell=cell, abc=abc)
383 842 : a = 0.5_dp*(abc(1) + abc(2))
384 842 : acosa = 0.5_dp*a
385 842 : asina = sqrt3*acosa
386 842 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = 0.0_dp
387 842 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = 0.0_dp
388 842 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
389 : CASE (cell_sym_rhombohedral)
390 593 : CALL get_cell(cell=cell, abc=abc)
391 2372 : a = SUM(abc(1:3))/3.0_dp
392 : alpha = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
393 : angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
394 593 : angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
395 593 : acosa = a*COS(alpha)
396 593 : asina = a*SIN(alpha)
397 593 : acosah = a*COS(0.5_dp*alpha)
398 593 : asinah = a*SIN(0.5_dp*alpha)
399 593 : norm = acosa/acosah
400 593 : norm_c = SQRT(1.0_dp - norm*norm)
401 593 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosa; cell%hmat(1, 3) = acosah*norm
402 593 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asina; cell%hmat(2, 3) = asinah*norm
403 593 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = a*norm_c
404 : CASE (cell_sym_monoclinic)
405 886 : CALL get_cell(cell=cell, abc=abc)
406 886 : beta = angle(cell%hmat(:, 1), cell%hmat(:, 3))
407 886 : cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)
408 886 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
409 886 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)
410 : CASE (cell_sym_monoclinic_gamma_ab)
411 : ! Cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg
412 746 : CALL get_cell(cell=cell, abc=abc)
413 746 : a = 0.5_dp*(abc(1) + abc(2))
414 746 : gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
415 746 : acosgamma = a*COS(gamma)
416 746 : asingamma = a*SIN(gamma)
417 746 : cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosgamma; cell%hmat(1, 3) = 0.0_dp
418 746 : cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asingamma; cell%hmat(2, 3) = 0.0_dp
419 338105 : cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
420 : CASE (cell_sym_triclinic)
421 : ! Nothing to do
422 : END SELECT
423 :
424 : ! Do we have an (almost) orthorhombic cell?
425 : IF ((ABS(cell%hmat(1, 2)) < eps_hmat) .AND. (ABS(cell%hmat(1, 3)) < eps_hmat) .AND. &
426 : (ABS(cell%hmat(2, 1)) < eps_hmat) .AND. (ABS(cell%hmat(2, 3)) < eps_hmat) .AND. &
427 337359 : (ABS(cell%hmat(3, 1)) < eps_hmat) .AND. (ABS(cell%hmat(3, 2)) < eps_hmat)) THEN
428 315790 : cell%orthorhombic = .TRUE.
429 : ELSE
430 21569 : cell%orthorhombic = .FALSE.
431 : END IF
432 :
433 : ! Retain an exact orthorhombic cell
434 : ! (off-diagonal elements must remain zero identically to keep QS fast)
435 337359 : IF (cell%orthorhombic) THEN
436 315790 : cell%hmat(1, 2) = 0.0_dp
437 315790 : cell%hmat(1, 3) = 0.0_dp
438 315790 : cell%hmat(2, 1) = 0.0_dp
439 315790 : cell%hmat(2, 3) = 0.0_dp
440 315790 : cell%hmat(3, 1) = 0.0_dp
441 315790 : cell%hmat(3, 2) = 0.0_dp
442 : END IF
443 :
444 1349436 : dim = COUNT(cell%perd == 1)
445 337359 : IF ((dim == 1) .AND. (.NOT. cell%orthorhombic)) THEN
446 0 : CPABORT("Non-orthorhombic and not periodic")
447 : END IF
448 :
449 : ! Update deth and hmat_inv with enforced symmetry
450 337359 : cell%deth = ABS(det_3x3(cell%hmat))
451 337359 : IF (cell%deth < 1.0E-10_dp) THEN
452 : CALL cp_abort(__LOCATION__, &
453 : "An invalid set of cell vectors was obtained after applying "// &
454 0 : "the requested cell symmetry. The determinant det(h) is too small")
455 : END IF
456 4385667 : cell%h_inv = inv_3x3(cell%hmat)
457 :
458 337359 : END SUBROUTINE init_cell
459 :
460 : ! **************************************************************************************************
461 : !> \brief Calculate the distance between two lattice planes as defined by
462 : !> a triple of Miller indices (hkl).
463 : !> \param h ...
464 : !> \param k ...
465 : !> \param l ...
466 : !> \param cell ...
467 : !> \return ...
468 : !> \date 18.11.2004
469 : !> \author Matthias Krack
470 : !> \version 1.0
471 : ! **************************************************************************************************
472 6450540 : FUNCTION plane_distance(h, k, l, cell) RESULT(distance)
473 :
474 : INTEGER, INTENT(IN) :: h, k, l
475 : TYPE(cell_type), INTENT(IN) :: cell
476 : REAL(KIND=dp) :: distance
477 :
478 : REAL(KIND=dp) :: a, alpha, b, beta, c, cosa, cosb, cosg, &
479 : d, gamma, x, y, z
480 : REAL(KIND=dp), DIMENSION(3) :: abc
481 :
482 6450540 : x = REAL(h, KIND=dp)
483 6450540 : y = REAL(k, KIND=dp)
484 6450540 : z = REAL(l, KIND=dp)
485 :
486 6450540 : CALL get_cell(cell=cell, abc=abc)
487 :
488 6450540 : a = abc(1)
489 6450540 : b = abc(2)
490 6450540 : c = abc(3)
491 :
492 6450540 : IF (cell%orthorhombic) THEN
493 :
494 6397042 : d = (x/a)**2 + (y/b)**2 + (z/c)**2
495 :
496 : ELSE
497 :
498 : CALL get_cell(cell=cell, &
499 : alpha=alpha, &
500 : beta=beta, &
501 53498 : gamma=gamma)
502 :
503 53498 : alpha = alpha/degree
504 53498 : beta = beta/degree
505 53498 : gamma = gamma/degree
506 :
507 53498 : cosa = COS(alpha)
508 53498 : cosb = COS(beta)
509 53498 : cosg = COS(gamma)
510 :
511 : d = ((x*b*c*SIN(alpha))**2 + &
512 : (y*c*a*SIN(beta))**2 + &
513 : (z*a*b*SIN(gamma))**2 + &
514 : 2.0_dp*a*b*c*(x*y*c*(cosa*cosb - cosg) + &
515 : z*x*b*(cosg*cosa - cosb) + &
516 : y*z*a*(cosb*cosg - cosa)))/ &
517 : ((a*b*c)**2*(1.0_dp - cosa**2 - cosb**2 - cosg**2 + &
518 53498 : 2.0_dp*cosa*cosb*cosg))
519 :
520 : END IF
521 :
522 6450540 : distance = 1.0_dp/SQRT(d)
523 :
524 6450540 : END FUNCTION plane_distance
525 :
526 : ! **************************************************************************************************
527 : !> \brief Apply the periodic boundary conditions defined by a simulation
528 : !> cell to a position vector r.
529 : !> \param r ...
530 : !> \param cell ...
531 : !> \return ...
532 : !> \date 16.01.2002
533 : !> \author Matthias Krack
534 : !> \version 1.0
535 : ! **************************************************************************************************
536 332486867 : FUNCTION pbc1(r, cell) RESULT(r_pbc)
537 :
538 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
539 : TYPE(cell_type), INTENT(IN) :: cell
540 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
541 :
542 : REAL(KIND=dp), DIMENSION(3) :: s
543 :
544 332486867 : IF (cell%orthorhombic) THEN
545 306514128 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*ANINT(cell%h_inv(1, 1)*r(1))
546 306514128 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*ANINT(cell%h_inv(2, 2)*r(2))
547 306514128 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*ANINT(cell%h_inv(3, 3)*r(3))
548 : ELSE
549 25972739 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
550 25972739 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
551 25972739 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
552 25972739 : s(1) = s(1) - cell%perd(1)*ANINT(s(1))
553 25972739 : s(2) = s(2) - cell%perd(2)*ANINT(s(2))
554 25972739 : s(3) = s(3) - cell%perd(3)*ANINT(s(3))
555 25972739 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
556 25972739 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
557 25972739 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
558 : END IF
559 :
560 332486867 : END FUNCTION pbc1
561 :
562 : ! **************************************************************************************************
563 : !> \brief Apply the periodic boundary conditions defined by a simulation
564 : !> cell to a position vector r subtracting nl from the periodic images
565 : !> \param r ...
566 : !> \param cell ...
567 : !> \param nl ...
568 : !> \return ...
569 : !> \date 16.01.2002
570 : !> \author Matthias Krack
571 : !> \version 1.0
572 : ! **************************************************************************************************
573 142304514 : FUNCTION pbc2(r, cell, nl) RESULT(r_pbc)
574 :
575 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
576 : TYPE(cell_type), INTENT(IN) :: cell
577 : INTEGER, DIMENSION(3), INTENT(IN) :: nl
578 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
579 :
580 : REAL(KIND=dp), DIMENSION(3) :: s
581 :
582 142304514 : IF (cell%orthorhombic) THEN
583 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)* &
584 132353154 : REAL(NINT(cell%h_inv(1, 1)*r(1)) - nl(1), dp)
585 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)* &
586 132353154 : REAL(NINT(cell%h_inv(2, 2)*r(2)) - nl(2), dp)
587 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)* &
588 132353154 : REAL(NINT(cell%h_inv(3, 3)*r(3)) - nl(3), dp)
589 : ELSE
590 9951360 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
591 9951360 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
592 9951360 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
593 9951360 : s(1) = s(1) - cell%perd(1)*REAL(NINT(s(1)) - nl(1), dp)
594 9951360 : s(2) = s(2) - cell%perd(2)*REAL(NINT(s(2)) - nl(2), dp)
595 9951360 : s(3) = s(3) - cell%perd(3)*REAL(NINT(s(3)) - nl(3), dp)
596 9951360 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
597 9951360 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
598 9951360 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
599 : END IF
600 :
601 142304514 : END FUNCTION pbc2
602 :
603 : ! **************************************************************************************************
604 : !> \brief Apply the periodic boundary conditions defined by the simulation
605 : !> cell cell to the vector pointing from atom a to atom b.
606 : !> \param ra ...
607 : !> \param rb ...
608 : !> \param cell ...
609 : !> \return ...
610 : !> \date 11.03.2004
611 : !> \author Matthias Krack
612 : !> \version 1.0
613 : ! **************************************************************************************************
614 117615083 : FUNCTION pbc3(ra, rb, cell) RESULT(rab_pbc)
615 :
616 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb
617 : TYPE(cell_type), INTENT(IN) :: cell
618 : REAL(KIND=dp), DIMENSION(3) :: rab_pbc
619 :
620 : INTEGER :: icell, jcell, kcell
621 : INTEGER, DIMENSION(3) :: periodic
622 : REAL(KIND=dp) :: rab2, rab2_pbc
623 : REAL(KIND=dp), DIMENSION(3) :: r, ra_pbc, rab, rb_image, rb_pbc, s2r
624 :
625 117615083 : CALL get_cell(cell=cell, periodic=periodic)
626 :
627 117615083 : ra_pbc(:) = pbc(ra(:), cell)
628 117615083 : rb_pbc(:) = pbc(rb(:), cell)
629 :
630 117615083 : rab2_pbc = HUGE(1.0_dp)
631 :
632 466114640 : DO icell = -periodic(1), periodic(1)
633 1507268375 : DO jcell = -periodic(2), periodic(2)
634 4508769657 : DO kcell = -periodic(3), periodic(3)
635 12476465460 : r = REAL((/icell, jcell, kcell/), dp)
636 3119116365 : CALL scaled_to_real(s2r, r, cell)
637 12476465460 : rb_image(:) = rb_pbc(:) + s2r
638 12476465460 : rab(:) = rb_image(:) - ra_pbc(:)
639 3119116365 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
640 4160270100 : IF (rab2 < rab2_pbc) THEN
641 2427952528 : rab2_pbc = rab2
642 2427952528 : rab_pbc(:) = rab(:)
643 : END IF
644 : END DO
645 : END DO
646 : END DO
647 :
648 117615083 : END FUNCTION pbc3
649 :
650 : !if positive_range == true, r(i) (or s(i)) in range [0, hmat(i,i)],
651 : !else, r(i) (s(i)) in range [-hmat(i,i)/2, hmat(i,i)/2]
652 : ! **************************************************************************************************
653 : !> \brief ...
654 : !> \param r ...
655 : !> \param cell ...
656 : !> \param positive_range ...
657 : !> \return ...
658 : ! **************************************************************************************************
659 238 : FUNCTION pbc4(r, cell, positive_range) RESULT(r_pbc)
660 :
661 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
662 : TYPE(cell_type), INTENT(IN) :: cell
663 : LOGICAL :: positive_range
664 : REAL(KIND=dp), DIMENSION(3) :: r_pbc
665 :
666 : REAL(KIND=dp), DIMENSION(3) :: s
667 :
668 238 : IF (positive_range) THEN
669 238 : IF (cell%orthorhombic) THEN
670 238 : r_pbc(1) = r(1) - cell%hmat(1, 1)*cell%perd(1)*FLOOR(cell%h_inv(1, 1)*r(1))
671 238 : r_pbc(2) = r(2) - cell%hmat(2, 2)*cell%perd(2)*FLOOR(cell%h_inv(2, 2)*r(2))
672 238 : r_pbc(3) = r(3) - cell%hmat(3, 3)*cell%perd(3)*FLOOR(cell%h_inv(3, 3)*r(3))
673 : ELSE
674 0 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
675 0 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
676 0 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
677 0 : s(1) = s(1) - cell%perd(1)*FLOOR(s(1))
678 0 : s(2) = s(2) - cell%perd(2)*FLOOR(s(2))
679 0 : s(3) = s(3) - cell%perd(3)*FLOOR(s(3))
680 0 : r_pbc(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
681 0 : r_pbc(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
682 0 : r_pbc(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
683 : END IF
684 : ELSE
685 0 : r_pbc = pbc1(r, cell)
686 : END IF
687 :
688 238 : END FUNCTION pbc4
689 :
690 : ! **************************************************************************************************
691 : !> \brief Transform real to scaled cell coordinates.
692 : !> s=h_inv*r
693 : !> \param s ...
694 : !> \param r ...
695 : !> \param cell ...
696 : !> \date 16.01.2002
697 : !> \author Matthias Krack
698 : !> \version 1.0
699 : ! **************************************************************************************************
700 59249454 : SUBROUTINE real_to_scaled(s, r, cell)
701 :
702 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: s
703 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r
704 : TYPE(cell_type), INTENT(IN) :: cell
705 :
706 59249454 : IF (cell%orthorhombic) THEN
707 55923662 : s(1) = cell%h_inv(1, 1)*r(1)
708 55923662 : s(2) = cell%h_inv(2, 2)*r(2)
709 55923662 : s(3) = cell%h_inv(3, 3)*r(3)
710 : ELSE
711 3325792 : s(1) = cell%h_inv(1, 1)*r(1) + cell%h_inv(1, 2)*r(2) + cell%h_inv(1, 3)*r(3)
712 3325792 : s(2) = cell%h_inv(2, 1)*r(1) + cell%h_inv(2, 2)*r(2) + cell%h_inv(2, 3)*r(3)
713 3325792 : s(3) = cell%h_inv(3, 1)*r(1) + cell%h_inv(3, 2)*r(2) + cell%h_inv(3, 3)*r(3)
714 : END IF
715 :
716 59249454 : END SUBROUTINE real_to_scaled
717 :
718 : ! **************************************************************************************************
719 : !> \brief Transform scaled cell coordinates real coordinates.
720 : !> r=h*s
721 : !> \param r ...
722 : !> \param s ...
723 : !> \param cell ...
724 : !> \date 16.01.2002
725 : !> \author Matthias Krack
726 : !> \version 1.0
727 : ! **************************************************************************************************
728 3230782142 : SUBROUTINE scaled_to_real(r, s, cell)
729 :
730 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: r
731 : REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s
732 : TYPE(cell_type), INTENT(IN) :: cell
733 :
734 3230782142 : IF (cell%orthorhombic) THEN
735 2965431883 : r(1) = cell%hmat(1, 1)*s(1)
736 2965431883 : r(2) = cell%hmat(2, 2)*s(2)
737 2965431883 : r(3) = cell%hmat(3, 3)*s(3)
738 : ELSE
739 265350259 : r(1) = cell%hmat(1, 1)*s(1) + cell%hmat(1, 2)*s(2) + cell%hmat(1, 3)*s(3)
740 265350259 : r(2) = cell%hmat(2, 1)*s(1) + cell%hmat(2, 2)*s(2) + cell%hmat(2, 3)*s(3)
741 265350259 : r(3) = cell%hmat(3, 1)*s(1) + cell%hmat(3, 2)*s(2) + cell%hmat(3, 3)*s(3)
742 : END IF
743 :
744 3230782142 : END SUBROUTINE scaled_to_real
745 :
746 : ! **************************************************************************************************
747 : !> \brief allocates and initializes a cell
748 : !> \param cell the cell to initialize
749 : !> \param hmat the h matrix that defines the cell
750 : !> \param periodic periodicity of the cell
751 : !> \par History
752 : !> 09.2003 created [fawzi]
753 : !> \author Fawzi Mohamed
754 : ! **************************************************************************************************
755 51828 : SUBROUTINE cell_create(cell, hmat, periodic)
756 :
757 : TYPE(cell_type), POINTER :: cell
758 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
759 : OPTIONAL :: hmat
760 : INTEGER, DIMENSION(3), INTENT(IN), OPTIONAL :: periodic
761 :
762 51828 : CPASSERT(.NOT. ASSOCIATED(cell))
763 51828 : ALLOCATE (cell)
764 51828 : last_cell_id = last_cell_id + 1
765 51828 : cell%id_nr = last_cell_id
766 51828 : cell%ref_count = 1
767 51828 : IF (PRESENT(periodic)) THEN
768 1336 : cell%perd = periodic
769 : ELSE
770 205976 : cell%perd = 1
771 : END IF
772 51828 : cell%orthorhombic = .FALSE.
773 51828 : cell%symmetry_id = cell_sym_none
774 51828 : IF (PRESENT(hmat)) CALL init_cell(cell, hmat)
775 :
776 51828 : END SUBROUTINE cell_create
777 :
778 : ! **************************************************************************************************
779 : !> \brief retains the given cell (see doc/ReferenceCounting.html)
780 : !> \param cell the cell to retain
781 : !> \par History
782 : !> 09.2003 created [fawzi]
783 : !> \author Fawzi Mohamed
784 : ! **************************************************************************************************
785 47674 : SUBROUTINE cell_retain(cell)
786 :
787 : TYPE(cell_type), POINTER :: cell
788 :
789 47674 : CPASSERT(ASSOCIATED(cell))
790 47674 : CPASSERT(cell%ref_count > 0)
791 47674 : cell%ref_count = cell%ref_count + 1
792 :
793 47674 : END SUBROUTINE cell_retain
794 :
795 : ! **************************************************************************************************
796 : !> \brief releases the given cell (see doc/ReferenceCounting.html)
797 : !> \param cell the cell to release
798 : !> \par History
799 : !> 09.2003 created [fawzi]
800 : !> \author Fawzi Mohamed
801 : ! **************************************************************************************************
802 128175 : SUBROUTINE cell_release(cell)
803 :
804 : TYPE(cell_type), POINTER :: cell
805 :
806 128175 : IF (ASSOCIATED(cell)) THEN
807 99514 : CPASSERT(cell%ref_count > 0)
808 99514 : cell%ref_count = cell%ref_count - 1
809 99514 : IF (cell%ref_count == 0) THEN
810 51840 : DEALLOCATE (cell)
811 : END IF
812 99514 : NULLIFY (cell)
813 : END IF
814 :
815 128175 : END SUBROUTINE cell_release
816 :
817 : #if defined (__PLUMED2)
818 : ! **************************************************************************************************
819 : !> \brief For the interface with plumed, pass a cell pointer and retrieve it
820 : !> later. It's a hack, but avoids passing the cell back and forth
821 : !> across the Fortran/C++ interface
822 : !> \param cell ...
823 : !> \param set ...
824 : !> \date 28.02.2013
825 : !> \author RK
826 : !> \version 1.0
827 : ! **************************************************************************************************
828 2 : SUBROUTINE pbc_cp2k_plumed_getset_cell(cell, set)
829 : TYPE(cell_type), POINTER :: cell
830 : LOGICAL :: set
831 :
832 : TYPE(cell_type), POINTER, SAVE :: stored_cell
833 :
834 2 : IF (set .EQV. .TRUE.) THEN
835 2 : stored_cell => cell
836 : ELSE
837 0 : cell => stored_cell
838 : END IF
839 :
840 2 : END SUBROUTINE pbc_cp2k_plumed_getset_cell
841 : #endif
842 :
843 0 : END MODULE cell_types
|
