LCOV - code coverage report
Current view: top level - src/subsys - cell_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:20fe009) Lines: 267 283 94.3 %
Date: 2022-07-05 19:56:53 Functions: 18 20 90.0 %

          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

Generated by: LCOV version 1.15