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