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