LCOV - code coverage report
Current view: top level - src/motion - cell_opt_utils.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:37c9bd6) Lines: 167 181 92.3 %
Date: 2023-03-30 11:55:16 Functions: 6 6 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2023 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief contains a functional that calculates the energy and its derivatives
      10             : !>      for the geometry optimizer
      11             : !> \par History
      12             : !>      03.2008 - Teodoro Laino [tlaino] - University of Zurich - Cell Optimization
      13             : ! **************************************************************************************************
      14             : MODULE cell_opt_utils
      15             :    USE cell_types,                      ONLY: &
      16             :         cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, cell_sym_monoclinic_gamma_ab, &
      17             :         cell_sym_orthorhombic, cell_sym_rhombohedral, cell_sym_tetragonal_ab, &
      18             :         cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type, &
      19             :         get_cell_param, set_cell_param
      20             :    USE cp_files,                        ONLY: close_file,&
      21             :                                               open_file
      22             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23             :                                               cp_logger_create,&
      24             :                                               cp_logger_get_default_unit_nr,&
      25             :                                               cp_logger_release,&
      26             :                                               cp_logger_set,&
      27             :                                               cp_logger_type,&
      28             :                                               cp_to_string
      29             :    USE cp_units,                        ONLY: cp_units_rad
      30             :    USE input_constants,                 ONLY: fix_none,&
      31             :                                               fix_x,&
      32             :                                               fix_xy,&
      33             :                                               fix_xz,&
      34             :                                               fix_y,&
      35             :                                               fix_yz,&
      36             :                                               fix_z
      37             :    USE input_cp2k_global,               ONLY: create_global_section
      38             :    USE input_enumeration_types,         ONLY: enum_i2c,&
      39             :                                               enumeration_type
      40             :    USE input_keyword_types,             ONLY: keyword_get,&
      41             :                                               keyword_type
      42             :    USE input_section_types,             ONLY: section_get_keyword,&
      43             :                                               section_release,&
      44             :                                               section_type,&
      45             :                                               section_vals_type,&
      46             :                                               section_vals_val_get,&
      47             :                                               section_vals_val_set
      48             :    USE kinds,                           ONLY: default_path_length,&
      49             :                                               default_string_length,&
      50             :                                               dp
      51             :    USE mathconstants,                   ONLY: sqrt3
      52             :    USE mathlib,                         ONLY: angle
      53             :    USE message_passing,                 ONLY: mp_para_env_type
      54             : #include "../base/base_uses.f90"
      55             : 
      56             :    IMPLICIT NONE
      57             :    PRIVATE
      58             : 
      59             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      60             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_opt_utils'
      61             : 
      62             :    PUBLIC :: get_ut_cell_matrix, get_dg_dh, gopt_new_logger_create, &
      63             :              gopt_new_logger_release, read_external_press_tensor, &
      64             :              apply_cell_constraints
      65             : 
      66             : CONTAINS
      67             : ! **************************************************************************************************
      68             : !> \brief Transform a general cell matrix into an upper triangular one
      69             : !> \param cell ...
      70             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
      71             : ! **************************************************************************************************
      72         192 :    SUBROUTINE get_ut_cell_matrix(cell)
      73             :       TYPE(cell_type), POINTER                           :: cell
      74             : 
      75             :       INTEGER, DIMENSION(3)                              :: periodic
      76             :       REAL(KIND=dp), DIMENSION(3)                        :: cell_angle, cell_length
      77             : 
      78             : ! orthorhombic cells are diagonal, and should stay exactly like this
      79             : 
      80         192 :       IF (.NOT. cell%orthorhombic) THEN
      81          56 :          CALL get_cell_param(cell, cell_length, cell_angle, cp_units_rad, periodic=periodic)
      82             :          CALL set_cell_param(cell, cell_length, cell_angle, periodic=periodic, &
      83          56 :                              do_init_cell=.TRUE.)
      84             :       END IF
      85             : 
      86         192 :    END SUBROUTINE get_ut_cell_matrix
      87             : 
      88             : ! **************************************************************************************************
      89             : !> \brief creates a new logger used for cell optimization algorithm
      90             : !> \param new_logger ...
      91             : !> \param root_section ...
      92             : !> \param para_env ...
      93             : !> \param project_name ...
      94             : !> \param id_run ...
      95             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
      96             : ! **************************************************************************************************
      97         850 :    SUBROUTINE gopt_new_logger_create(new_logger, root_section, para_env, project_name, &
      98             :                                      id_run)
      99             :       TYPE(cp_logger_type), POINTER                      :: new_logger
     100             :       TYPE(section_vals_type), POINTER                   :: root_section
     101             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     102             :       CHARACTER(len=default_string_length), INTENT(OUT)  :: project_name
     103             :       INTEGER, INTENT(IN)                                :: id_run
     104             : 
     105             :       CHARACTER(len=default_path_length)                 :: c_val, input_file_path, output_file_path
     106             :       CHARACTER(len=default_string_length)               :: label
     107             :       INTEGER                                            :: i, lp, unit_nr
     108             :       TYPE(cp_logger_type), POINTER                      :: logger
     109             :       TYPE(enumeration_type), POINTER                    :: enum
     110             :       TYPE(keyword_type), POINTER                        :: keyword
     111             :       TYPE(section_type), POINTER                        :: section
     112             : 
     113         850 :       NULLIFY (new_logger, logger, enum, keyword, section)
     114         850 :       logger => cp_get_default_logger()
     115             : 
     116         850 :       CALL create_global_section(section)
     117         850 :       keyword => section_get_keyword(section, "RUN_TYPE")
     118         850 :       CALL keyword_get(keyword, enum=enum)
     119         850 :       label = TRIM(enum_i2c(enum, id_run))
     120         850 :       CALL section_release(section)
     121             : 
     122             :       ! Redirecting output of the sub_calculation to a different file
     123         850 :       CALL section_vals_val_get(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
     124         850 :       input_file_path = TRIM(project_name)
     125         850 :       lp = LEN_TRIM(input_file_path)
     126         850 :       i = logger%iter_info%iteration(logger%iter_info%n_rlevel)
     127         850 :       input_file_path(lp + 1:LEN(input_file_path)) = "-"//TRIM(label)//"-"//ADJUSTL(cp_to_string(i))
     128         850 :       lp = LEN_TRIM(input_file_path)
     129         850 :       CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=input_file_path(1:lp))
     130         850 :       CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
     131             : 
     132             :       ! Redirecting output into a new file
     133         850 :       output_file_path = input_file_path(1:lp)//".out"
     134         850 :       IF (para_env%is_source()) THEN
     135             :          CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
     136         425 :                         file_action="WRITE", file_position="APPEND", unit_number=unit_nr)
     137             :       ELSE
     138         425 :          unit_nr = -1
     139             :       END IF
     140             :       CALL cp_logger_create(new_logger, para_env=para_env, default_global_unit_nr=unit_nr, &
     141         850 :                             close_global_unit_on_dealloc=.FALSE.)
     142         850 :       CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=c_val)
     143         850 :       IF (c_val /= "") THEN
     144         850 :          CALL cp_logger_set(new_logger, local_filename=TRIM(c_val)//"_localLog")
     145             :       END IF
     146         850 :       new_logger%iter_info%project_name = c_val
     147             :       CALL section_vals_val_get(root_section, "GLOBAL%PRINT_LEVEL", &
     148         850 :                                 i_val=new_logger%iter_info%print_level)
     149             : 
     150         850 :    END SUBROUTINE gopt_new_logger_create
     151             : 
     152             : ! **************************************************************************************************
     153             : !> \brief releases a new logger used for cell optimization algorithm
     154             : !> \param new_logger ...
     155             : !> \param root_section ...
     156             : !> \param para_env ...
     157             : !> \param project_name ...
     158             : !> \param id_run ...
     159             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     160             : ! **************************************************************************************************
     161         850 :    SUBROUTINE gopt_new_logger_release(new_logger, root_section, para_env, project_name, id_run)
     162             :       TYPE(cp_logger_type), POINTER                      :: new_logger
     163             :       TYPE(section_vals_type), POINTER                   :: root_section
     164             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     165             :       CHARACTER(len=default_string_length), INTENT(IN)   :: project_name
     166             :       INTEGER, INTENT(IN)                                :: id_run
     167             : 
     168             :       INTEGER                                            :: unit_nr
     169             : 
     170         850 :       IF (para_env%is_source()) THEN
     171         425 :          unit_nr = cp_logger_get_default_unit_nr(new_logger)
     172         425 :          CALL close_file(unit_number=unit_nr)
     173             :       END IF
     174         850 :       CALL cp_logger_release(new_logger)
     175         850 :       CALL section_vals_val_set(root_section, "GLOBAL%RUN_TYPE", i_val=id_run)
     176         850 :       CALL section_vals_val_set(root_section, "GLOBAL%PROJECT_NAME", c_val=project_name)
     177             : 
     178         850 :    END SUBROUTINE gopt_new_logger_release
     179             : 
     180             : ! **************************************************************************************************
     181             : !> \brief Reads the external pressure tensor
     182             : !> \param geo_section ...
     183             : !> \param cell ...
     184             : !> \param pres_ext ...
     185             : !> \param mtrx ...
     186             : !> \param rot ...
     187             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     188             : ! **************************************************************************************************
     189         192 :    SUBROUTINE read_external_press_tensor(geo_section, cell, pres_ext, mtrx, rot)
     190             :       TYPE(section_vals_type), POINTER                   :: geo_section
     191             :       TYPE(cell_type), POINTER                           :: cell
     192             :       REAL(KIND=dp), INTENT(OUT)                         :: pres_ext
     193             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: mtrx
     194             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: rot
     195             : 
     196             :       INTEGER                                            :: i, ind, j
     197             :       LOGICAL                                            :: check
     198             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: pres_ext_tens
     199         192 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pvals
     200             : 
     201         192 :       NULLIFY (pvals)
     202         192 :       mtrx = 0.0_dp
     203         192 :       pres_ext_tens = 0.0_dp
     204         192 :       pres_ext = 0.0_dp
     205         192 :       CALL section_vals_val_get(geo_section, "EXTERNAL_PRESSURE", r_vals=pvals)
     206         192 :       check = (SIZE(pvals) == 1) .OR. (SIZE(pvals) == 9)
     207         192 :       IF (.NOT. check) &
     208           0 :          CPABORT("EXTERNAL_PRESSURE can have 1 or 9 components only!")
     209             : 
     210         192 :       IF (SIZE(pvals) == 9) THEN
     211             :          ind = 0
     212         392 :          DO i = 1, 3
     213        1274 :             DO j = 1, 3
     214         882 :                ind = ind + 1
     215        1176 :                pres_ext_tens(j, i) = pvals(ind)
     216             :             END DO
     217             :          END DO
     218             :          ! Also the pressure tensor must be oriented in the same canonical directions
     219             :          ! of the simulation cell
     220        3920 :          pres_ext_tens = MATMUL(TRANSPOSE(rot), pres_ext_tens)
     221         392 :          DO i = 1, 3
     222         392 :             pres_ext = pres_ext + pres_ext_tens(i, i)
     223             :          END DO
     224          98 :          pres_ext = pres_ext/3.0_dp
     225         392 :          DO i = 1, 3
     226         392 :             pres_ext_tens(i, i) = pres_ext_tens(i, i) - pres_ext
     227             :          END DO
     228             :       ELSE
     229          94 :          pres_ext = pvals(1)
     230             :       END IF
     231             : 
     232        2496 :       IF (ANY(pres_ext_tens > 1.0E-5_dp)) THEN
     233           0 :          mtrx = cell%deth*MATMUL(cell%h_inv, MATMUL(pres_ext_tens, TRANSPOSE(cell%h_inv)))
     234             :       END IF
     235             : 
     236         192 :    END SUBROUTINE read_external_press_tensor
     237             : 
     238             : ! **************************************************************************************************
     239             : !> \brief Computes the derivatives for the cell
     240             : !> \param gradient ...
     241             : !> \param av_ptens ...
     242             : !> \param pres_ext ...
     243             : !> \param cell ...
     244             : !> \param mtrx ...
     245             : !> \param keep_angles ...
     246             : !> \param keep_symmetry ...
     247             : !> \param pres_int ...
     248             : !> \param pres_constr ...
     249             : !> \param constraint_id ...
     250             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     251             : ! **************************************************************************************************
     252        3377 :    SUBROUTINE get_dg_dh(gradient, av_ptens, pres_ext, cell, mtrx, keep_angles, &
     253             :                         keep_symmetry, pres_int, pres_constr, constraint_id)
     254             : 
     255             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gradient
     256             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: av_ptens
     257             :       REAL(KIND=dp), INTENT(IN)                          :: pres_ext
     258             :       TYPE(cell_type), POINTER                           :: cell
     259             :       REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN)         :: mtrx
     260             :       LOGICAL, INTENT(IN), OPTIONAL                      :: keep_angles, keep_symmetry
     261             :       REAL(KIND=dp), INTENT(OUT)                         :: pres_int, pres_constr
     262             :       INTEGER, INTENT(IN), OPTIONAL                      :: constraint_id
     263             : 
     264             :       INTEGER                                            :: i, my_constraint_id
     265             :       LOGICAL                                            :: my_keep_angles, my_keep_symmetry
     266             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: correction, pten_hinv_old, ptens
     267             : 
     268        3377 :       my_keep_angles = .FALSE.
     269        3377 :       IF (PRESENT(keep_angles)) my_keep_angles = keep_angles
     270        3377 :       my_keep_symmetry = .FALSE.
     271        3377 :       IF (PRESENT(keep_symmetry)) my_keep_symmetry = keep_symmetry
     272       23639 :       gradient = 0.0_dp
     273        3377 :       IF (PRESENT(constraint_id)) THEN
     274        3377 :          my_constraint_id = constraint_id
     275             :       ELSE
     276           0 :          my_constraint_id = fix_none
     277             :       END IF
     278             : 
     279       23639 :       gradient = 0.0_dp
     280             : 
     281        3377 :       ptens = av_ptens
     282             : 
     283             :       ! Evaluating the internal pressure
     284        3377 :       pres_int = 0.0_dp
     285       13508 :       DO i = 1, 3
     286       13508 :          pres_int = pres_int + ptens(i, i)
     287             :       END DO
     288        3377 :       pres_int = pres_int/3.0_dp
     289             : 
     290           0 :       SELECT CASE (my_constraint_id)
     291             :       CASE (fix_x)
     292           0 :          pres_constr = ptens(2, 2) + ptens(3, 3)
     293             :       CASE (fix_y)
     294           0 :          pres_constr = ptens(1, 1) + ptens(3, 3)
     295             :       CASE (fix_z)
     296           5 :          pres_constr = ptens(1, 1) + ptens(2, 2)
     297             :       CASE (fix_xy)
     298          10 :          pres_constr = ptens(3, 3)
     299             :       CASE (fix_xz)
     300           0 :          pres_constr = ptens(2, 2)
     301             :       CASE (fix_yz)
     302           0 :          pres_constr = ptens(1, 1)
     303             :       CASE (fix_none)
     304        3377 :          pres_constr = ptens(1, 1) + ptens(2, 2) + ptens(3, 3)
     305             :       END SELECT
     306        3377 :       pres_constr = pres_constr/3.0_dp
     307             : 
     308        3377 :       ptens(1, 1) = av_ptens(1, 1) - pres_ext
     309        3377 :       ptens(2, 2) = av_ptens(2, 2) - pres_ext
     310        3377 :       ptens(3, 3) = av_ptens(3, 3) - pres_ext
     311             : 
     312      175604 :       pten_hinv_old = cell%deth*MATMUL(cell%h_inv, ptens)
     313      135080 :       correction = MATMUL(mtrx, cell%hmat)
     314             : 
     315        3377 :       gradient(1) = pten_hinv_old(1, 1) - correction(1, 1)
     316        3377 :       gradient(2) = pten_hinv_old(2, 1) - correction(2, 1)
     317        3377 :       gradient(3) = pten_hinv_old(2, 2) - correction(2, 2)
     318        3377 :       gradient(4) = pten_hinv_old(3, 1) - correction(3, 1)
     319        3377 :       gradient(5) = pten_hinv_old(3, 2) - correction(3, 2)
     320        3377 :       gradient(6) = pten_hinv_old(3, 3) - correction(3, 3)
     321             : 
     322        3377 :       CALL apply_cell_constraints(gradient, cell, my_keep_angles, my_keep_symmetry, my_constraint_id)
     323             : 
     324       23639 :       gradient = -gradient
     325             : 
     326        3377 :    END SUBROUTINE get_dg_dh
     327             : 
     328             : ! **************************************************************************************************
     329             : !> \brief Apply cell constraints
     330             : !> \param gradient ...
     331             : !> \param cell ...
     332             : !> \param keep_angles ...
     333             : !> \param keep_symmetry ...
     334             : !> \param constraint_id ...
     335             : !> \author Matthias Krack (October 26, 2017, MK)
     336             : ! **************************************************************************************************
     337        3377 :    SUBROUTINE apply_cell_constraints(gradient, cell, keep_angles, keep_symmetry, constraint_id)
     338             : 
     339             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: gradient
     340             :       TYPE(cell_type), POINTER                           :: cell
     341             :       LOGICAL, INTENT(IN)                                :: keep_angles, keep_symmetry
     342             :       INTEGER, INTENT(IN)                                :: constraint_id
     343             : 
     344             :       REAL(KIND=dp)                                      :: a, a_length, ab_length, b_length, cosa, &
     345             :                                                             cosah, cosgamma, deriv_gamma, g, &
     346             :                                                             gamma, norm, norm_b, norm_c, sina, &
     347             :                                                             sinah, singamma
     348             : 
     349        3377 :       IF (keep_angles) THEN
     350             :          ! If we want to keep the angles constant we have to project out the
     351             :          ! components of the cell angles
     352         400 :          norm_b = DOT_PRODUCT(cell%hmat(:, 2), cell%hmat(:, 2))
     353         300 :          norm = DOT_PRODUCT(cell%hmat(1:2, 2), gradient(2:3))
     354         500 :          gradient(2:3) = cell%hmat(1:2, 2)/norm_b*norm
     355         400 :          norm_c = DOT_PRODUCT(cell%hmat(:, 3), cell%hmat(:, 3))
     356         400 :          norm = DOT_PRODUCT(cell%hmat(1:3, 3), gradient(4:6))
     357         700 :          gradient(4:6) = cell%hmat(1:3, 3)/norm_c*norm
     358             :          ! Retain an exact orthorhombic cell
     359             :          ! (off-diagonal elements must remain zero identically to keep QS fast)
     360         100 :          IF (cell%orthorhombic) THEN
     361          70 :             gradient(2) = 0.0_dp
     362          70 :             gradient(4) = 0.0_dp
     363          70 :             gradient(5) = 0.0_dp
     364             :          END IF
     365             :       END IF
     366             : 
     367        3377 :       IF (keep_symmetry) THEN
     368         568 :          SELECT CASE (cell%symmetry_id)
     369             :          CASE (cell_sym_cubic, &
     370             :                cell_sym_tetragonal_ab, &
     371             :                cell_sym_tetragonal_ac, &
     372             :                cell_sym_tetragonal_bc, &
     373             :                cell_sym_orthorhombic)
     374         100 :             SELECT CASE (cell%symmetry_id)
     375             :             CASE (cell_sym_cubic)
     376         100 :                g = (gradient(1) + gradient(3) + gradient(6))/3.0_dp
     377         100 :                gradient(1) = g
     378         100 :                gradient(3) = g
     379         465 :                gradient(6) = g
     380             :             CASE (cell_sym_tetragonal_ab, &
     381             :                   cell_sym_tetragonal_ac, &
     382             :                   cell_sym_tetragonal_bc)
     383         186 :                SELECT CASE (cell%symmetry_id)
     384             :                CASE (cell_sym_tetragonal_ab)
     385         186 :                   g = 0.5_dp*(gradient(1) + gradient(3))
     386         186 :                   gradient(1) = g
     387         186 :                   gradient(3) = g
     388             :                CASE (cell_sym_tetragonal_ac)
     389          93 :                   g = 0.5_dp*(gradient(1) + gradient(6))
     390          93 :                   gradient(1) = g
     391          93 :                   gradient(6) = g
     392             :                CASE (cell_sym_tetragonal_bc)
     393          86 :                   g = 0.5_dp*(gradient(3) + gradient(6))
     394          86 :                   gradient(3) = g
     395          86 :                   gradient(6) = g
     396             :                END SELECT
     397             :             CASE (cell_sym_orthorhombic)
     398             :                ! Nothing else to do
     399             :             END SELECT
     400         568 :             gradient(2) = 0.0_dp
     401         568 :             gradient(4) = 0.0_dp
     402         568 :             gradient(5) = 0.0_dp
     403             :          CASE (cell_sym_hexagonal)
     404         118 :             g = 0.5_dp*(gradient(1) + 0.5_dp*(gradient(2) + sqrt3*gradient(3)))
     405         118 :             gradient(1) = g
     406         118 :             gradient(2) = 0.5_dp*g
     407         118 :             gradient(3) = sqrt3*gradient(2)
     408         118 :             gradient(4) = 0.0_dp
     409         118 :             gradient(5) = 0.0_dp
     410             :          CASE (cell_sym_rhombohedral)
     411             :             a = (angle(cell%hmat(:, 3), cell%hmat(:, 2)) + &
     412             :                  angle(cell%hmat(:, 1), cell%hmat(:, 3)) + &
     413          83 :                  angle(cell%hmat(:, 1), cell%hmat(:, 2)))/3.0_dp
     414          83 :             cosa = COS(a)
     415          83 :             sina = SIN(a)
     416          83 :             cosah = COS(0.5_dp*a)
     417          83 :             sinah = SIN(0.5_dp*a)
     418          83 :             norm = cosa/cosah
     419          83 :             norm_c = SQRT(1.0_dp - norm*norm)
     420             :             g = (gradient(1) + gradient(2)*cosa + gradient(3)*sina + &
     421          83 :                  gradient(4)*cosah*norm + gradient(5)*sinah*norm + gradient(6)*norm_c)/3.0_dp
     422          83 :             gradient(1) = g
     423          83 :             gradient(2) = g*cosa
     424          83 :             gradient(3) = g*sina
     425          83 :             gradient(4) = g*cosah*norm
     426          83 :             gradient(5) = g*sinah*norm
     427          83 :             gradient(6) = g*norm_c
     428             :          CASE (cell_sym_monoclinic)
     429         126 :             gradient(2) = 0.0_dp
     430         126 :             gradient(5) = 0.0_dp
     431             :          CASE (cell_sym_monoclinic_gamma_ab)
     432             :             ! Cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg
     433             :             a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1) + &
     434             :                             cell%hmat(2, 1)*cell%hmat(2, 1) + &
     435         106 :                             cell%hmat(3, 1)*cell%hmat(3, 1))
     436             :             b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2) + &
     437             :                             cell%hmat(2, 2)*cell%hmat(2, 2) + &
     438         106 :                             cell%hmat(3, 2)*cell%hmat(3, 2))
     439         106 :             ab_length = 0.5_dp*(a_length + b_length)
     440         106 :             gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
     441         106 :             cosgamma = COS(gamma)
     442         106 :             singamma = SIN(gamma)
     443             :             ! Here, g is the average derivative of the cell vector length ab_length, and deriv_gamma is the derivative of the angle gamma
     444         106 :             g = 0.5_dp*(gradient(1) + cosgamma*gradient(2) + singamma*gradient(3))
     445         106 :             deriv_gamma = (gradient(3)*cosgamma - gradient(2)*singamma)/b_length
     446         106 :             gradient(1) = g
     447         106 :             gradient(2) = g*cosgamma - ab_length*singamma*deriv_gamma
     448         106 :             gradient(3) = g*singamma + ab_length*cosgamma*deriv_gamma
     449         106 :             gradient(4) = 0.0_dp
     450        1223 :             gradient(5) = 0.0_dp
     451             :          CASE (cell_sym_triclinic)
     452             :             ! Nothing to do
     453             :          END SELECT
     454             :       END IF
     455             : 
     456        3377 :       SELECT CASE (constraint_id)
     457             :       CASE (fix_x)
     458           0 :          gradient(1:2) = 0.0_dp
     459           0 :          gradient(4) = 0.0_dp
     460             :       CASE (fix_y)
     461           0 :          gradient(2:3) = 0.0_dp
     462           0 :          gradient(5) = 0.0_dp
     463             :       CASE (fix_z)
     464          20 :          gradient(4:6) = 0.0_dp
     465             :       CASE (fix_xy)
     466          60 :          gradient(1:5) = 0.0_dp
     467             :       CASE (fix_xz)
     468           0 :          gradient(1:2) = 0.0_dp
     469           0 :          gradient(4:6) = 0.0_dp
     470             :       CASE (fix_yz)
     471        3377 :          gradient(2:6) = 0.0_dp
     472             :       CASE (fix_none)
     473             :          ! Nothing to do
     474             :       END SELECT
     475             : 
     476        3377 :    END SUBROUTINE apply_cell_constraints
     477             : 
     478             : END MODULE cell_opt_utils

Generated by: LCOV version 1.15