|
CP2K 2.4 (Revision 12889)
|
00001 !-----------------------------------------------------------------------------! 00002 ! CP2K: A general program to perform molecular dynamics simulations ! 00003 ! Copyright (C) 2000 - 2013 CP2K developers group ! 00004 !-----------------------------------------------------------------------------! 00005 00006 ! ***************************************************************************** 00012 MODULE cell_types 00013 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 00014 cp_print_key_unit_nr 00015 USE cp_para_types, ONLY: cp_para_env_type 00016 USE cp_parser_methods, ONLY: parser_get_next_line 00017 USE cp_parser_types, ONLY: cp_parser_type,& 00018 parser_create,& 00019 parser_release 00020 USE cp_units, ONLY: cp_unit_from_cp2k,& 00021 cp_unit_to_cp2k,& 00022 cp_units_rad 00023 USE f77_blas 00024 USE input_constants, ONLY: & 00025 cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, & 00026 cell_sym_none, cell_sym_orthorhombic, cell_sym_rhombohedral, & 00027 cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, & 00028 cell_sym_tetragonal_bc, cell_sym_triclinic, do_cell_cp2k, do_cell_xsc, & 00029 use_perd_none, use_perd_x, use_perd_xy, use_perd_xyz, use_perd_xz, & 00030 use_perd_y, use_perd_yz, use_perd_z 00031 USE input_cp2k, ONLY: parsed_cp2k_input 00032 USE input_cp2k_subsys, ONLY: create_cell_section 00033 USE input_enumeration_types, ONLY: enum_i2c,& 00034 enumeration_type 00035 USE input_keyword_types, ONLY: keyword_get,& 00036 keyword_type 00037 USE input_section_types, ONLY: section_get_keyword,& 00038 section_release,& 00039 section_type,& 00040 section_vals_get_subs_vals,& 00041 section_vals_type,& 00042 section_vals_val_get,& 00043 section_vals_val_set,& 00044 section_vals_val_unset 00045 USE kinds, ONLY: default_path_length,& 00046 default_string_length,& 00047 dp 00048 USE mathconstants, ONLY: degree,& 00049 sqrt3 00050 USE mathlib, ONLY: angle,& 00051 det_3x3,& 00052 inv_3x3 00053 USE termination, ONLY: stop_program 00054 00055 USE, INTRINSIC :: ISO_C_BINDING, ONLY: C_PTR, C_F_POINTER 00056 #include "cp_common_uses.h" 00057 00058 IMPLICIT NONE 00059 00060 PRIVATE 00061 00062 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cell_types' 00063 00064 INTEGER, SAVE, PRIVATE :: last_cell_id = 0 00065 00066 ! ***************************************************************************** 00070 TYPE cell_type 00071 INTEGER :: ref_count,id_nr 00072 INTEGER :: symmetry_id 00073 LOGICAL :: orthorhombic 00074 REAL(KIND = dp) :: deth 00075 INTEGER, DIMENSION(3) :: perd 00076 REAL(KIND = dp), DIMENSION(3,3) :: hmat,h_inv 00077 END TYPE cell_type 00078 00079 TYPE cell_p_type 00080 TYPE(cell_type),POINTER :: cell 00081 END TYPE cell_p_type 00082 00083 ! Public subroutines 00084 PUBLIC :: get_cell, get_cell_param, init_cell, read_cell,& 00085 write_cell, cell_create, cell_retain, cell_release,& 00086 cell_clone, cell_copy, compare_cells, parse_cell_line, set_cell_param, & 00087 pbc_cp2k_plumed, pbc_cp2k_plumed_getset_cell !for plumed 00088 00089 ! Public functions 00090 PUBLIC :: plane_distance, pbc, real_to_scaled, scaled_to_real 00091 00092 ! Public data types 00093 PUBLIC :: cell_type, cell_p_type 00094 00095 INTERFACE pbc 00096 MODULE PROCEDURE pbc1,pbc2,pbc3 00097 END INTERFACE 00098 00099 CONTAINS 00100 00101 ! ***************************************************************************** 00102 FUNCTION compare_cells (cell_1, cell_2, error) RESULT(compare) 00103 00104 TYPE(cell_type), POINTER :: cell_1, cell_2 00105 TYPE(cp_error_type), INTENT(INOUT) :: error 00106 LOGICAL :: compare 00107 00108 CHARACTER(len=*), PARAMETER :: routineN = 'compare_cells', 00109 routineP = moduleN//':'//routineN 00110 00111 LOGICAL :: failure 00112 00113 failure = .FALSE. 00114 CPPrecondition(ASSOCIATED(cell_1),cp_failure_level,routineP,error,failure) 00115 CPPrecondition(ASSOCIATED(cell_2),cp_failure_level,routineP,error,failure) 00116 00117 compare = ((cell_2%deth == cell_1%deth).AND.& 00118 (ALL(cell_2%perd == cell_1%perd)).AND.& 00119 (ALL(cell_2%hmat == cell_1%hmat)).AND.& 00120 (ALL(cell_2%h_inv == cell_1%h_inv)).AND.& 00121 (cell_2%orthorhombic.EQV.cell_1%orthorhombic).AND.& 00122 (cell_2%symmetry_id == cell_1%symmetry_id)) 00123 00124 END FUNCTION compare_cells 00125 00126 ! ***************************************************************************** 00127 SUBROUTINE cell_clone (cell_in, cell_out, error) 00128 00129 TYPE(cell_type), POINTER :: cell_in, cell_out 00130 TYPE(cp_error_type), INTENT(INOUT) :: error 00131 00132 CHARACTER(len=*), PARAMETER :: routineN = 'cell_clone', 00133 routineP = moduleN//':'//routineN 00134 00135 LOGICAL :: failure 00136 00137 failure = .FALSE. 00138 CPPrecondition(ASSOCIATED(cell_in),cp_failure_level,routineP,error,failure) 00139 CPPrecondition(ASSOCIATED(cell_out),cp_failure_level,routineP,error,failure) 00140 00141 cell_out%deth = cell_in%deth 00142 cell_out%perd = cell_in%perd 00143 cell_out%hmat = cell_in%hmat 00144 cell_out%h_inv = cell_in%h_inv 00145 cell_out%orthorhombic = cell_in%orthorhombic 00146 cell_out%symmetry_id = cell_in%symmetry_id 00147 cell_out%ref_count = 1 00148 last_cell_id = last_cell_id + 1 00149 cell_out%id_nr = last_cell_id 00150 00151 END SUBROUTINE cell_clone 00152 00153 ! ***************************************************************************** 00154 SUBROUTINE cell_copy (cell_in, cell_out, error) 00155 00156 TYPE(cell_type), POINTER :: cell_in, cell_out 00157 TYPE(cp_error_type), INTENT(INOUT) :: error 00158 00159 CHARACTER(len=*), PARAMETER :: routineN = 'cell_copy', 00160 routineP = moduleN//':'//routineN 00161 00162 LOGICAL :: failure 00163 00164 failure = .FALSE. 00165 CPPrecondition(ASSOCIATED(cell_in),cp_failure_level,routineP,error,failure) 00166 CPPrecondition(ASSOCIATED(cell_out),cp_failure_level,routineP,error,failure) 00167 00168 cell_out%deth = cell_in%deth 00169 cell_out%perd = cell_in%perd 00170 cell_out%hmat = cell_in%hmat 00171 cell_out%h_inv = cell_in%h_inv 00172 cell_out%orthorhombic = cell_in%orthorhombic 00173 cell_out%symmetry_id = cell_in%symmetry_id 00174 00175 END SUBROUTINE cell_copy 00176 00177 ! ***************************************************************************** 00183 SUBROUTINE parse_cell_line(input_line, cell_itimes, cell_time, h, vol, error) 00184 CHARACTER(LEN=*), INTENT(IN) :: input_line 00185 INTEGER, INTENT(OUT) :: cell_itimes 00186 REAL(KIND=dp), INTENT(OUT) :: cell_time 00187 REAL(KIND=dp), DIMENSION(3, 3), 00188 INTENT(OUT) :: h 00189 REAL(KIND=dp), INTENT(OUT) :: vol 00190 TYPE(cp_error_type), INTENT(INOUT) :: error 00191 00192 CHARACTER(len=*), PARAMETER :: routineN = 'parse_cell_line', 00193 routineP = moduleN//':'//routineN 00194 00195 INTEGER :: i, j 00196 LOGICAL :: failure 00197 00198 failure = .FALSE. 00199 READ (input_line,*) cell_itimes, cell_time,& 00200 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 00201 DO i = 1, 3 00202 DO j = 1, 3 00203 h(j,i) = cp_unit_to_cp2k(h(j,i), "angstrom", error=error) 00204 END DO 00205 END DO 00206 00207 END SUBROUTINE parse_cell_line 00208 00209 ! ***************************************************************************** 00215 SUBROUTINE get_cell(cell,alpha,beta,gamma,deth,orthorhombic,abc,periodic,& 00216 h,h_inv,id_nr,symmetry_id) 00217 00218 TYPE(cell_type), POINTER :: cell 00219 REAL(KIND=dp), INTENT(OUT), OPTIONAL :: alpha, beta, gamma, deth 00220 LOGICAL, INTENT(OUT), OPTIONAL :: orthorhombic 00221 REAL(KIND=dp), DIMENSION(3), 00222 INTENT(OUT), OPTIONAL :: abc 00223 INTEGER, DIMENSION(3), INTENT(OUT), 00224 OPTIONAL :: periodic 00225 REAL(KIND=dp), DIMENSION(3, 3), 00226 INTENT(OUT), OPTIONAL :: h, h_inv 00227 INTEGER, INTENT(out), OPTIONAL :: id_nr, symmetry_id 00228 00229 IF (PRESENT(deth)) deth = cell%deth 00230 IF (PRESENT(orthorhombic)) orthorhombic = cell%orthorhombic 00231 IF (PRESENT(periodic)) periodic(:) = cell%perd(:) 00232 IF (PRESENT(h)) h(:,:) = cell%hmat(:,:) 00233 IF (PRESENT(h_inv)) h_inv(:,:) = cell%h_inv(:,:) 00234 00235 ! Calculate the lengths of the cell vectors a, b, and c 00236 IF (PRESENT(abc)) THEN 00237 abc(1) = SQRT(cell%hmat(1,1)*cell%hmat(1,1) +& 00238 cell%hmat(2,1)*cell%hmat(2,1) +& 00239 cell%hmat(3,1)*cell%hmat(3,1)) 00240 abc(2) = SQRT(cell%hmat(1,2)*cell%hmat(1,2) +& 00241 cell%hmat(2,2)*cell%hmat(2,2) +& 00242 cell%hmat(3,2)*cell%hmat(3,2)) 00243 abc(3) = SQRT(cell%hmat(1,3)*cell%hmat(1,3) +& 00244 cell%hmat(2,3)*cell%hmat(2,3) +& 00245 cell%hmat(3,3)*cell%hmat(3,3)) 00246 END IF 00247 00248 ! Angles between the cell vectors a, b, and c 00249 ! alpha = <(b,c) 00250 IF (PRESENT(alpha)) alpha = angle(cell%hmat(:,2),cell%hmat(:,3))*degree 00251 ! beta = <(a,c) 00252 IF (PRESENT(beta)) beta = angle(cell%hmat(:,1),cell%hmat(:,3))*degree 00253 ! gamma = <(a,b) 00254 IF (PRESENT(gamma)) gamma = angle(cell%hmat(:,1),cell%hmat(:,2))*degree 00255 IF (PRESENT(id_nr)) id_nr = cell%id_nr 00256 IF (PRESENT(symmetry_id)) symmetry_id = cell%symmetry_id 00257 00258 END SUBROUTINE get_cell 00259 00260 ! ***************************************************************************** 00266 SUBROUTINE get_cell_param(cell,cell_length,cell_angle,units_angle,periodic,& 00267 error) 00268 00269 TYPE(cell_type), POINTER :: cell 00270 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: cell_length 00271 REAL(KIND=dp), DIMENSION(3), 00272 INTENT(OUT), OPTIONAL :: cell_angle 00273 INTEGER, INTENT(IN), OPTIONAL :: units_angle 00274 INTEGER, DIMENSION(3), INTENT(OUT), 00275 OPTIONAL :: periodic 00276 TYPE(cp_error_type), INTENT(INOUT) :: error 00277 00278 CHARACTER(len=*), PARAMETER :: routineN = 'get_cell_param', 00279 routineP = moduleN//':'//routineN 00280 00281 REAL(KIND=dp) :: alpha, beta, gamma 00282 00283 CALL get_cell(cell=cell,abc=cell_length) 00284 00285 IF (PRESENT(cell_angle)) THEN 00286 CALL get_cell(cell=cell,alpha=alpha,beta=beta,gamma=gamma) 00287 cell_angle(:) = (/alpha,beta,gamma/) 00288 IF (PRESENT(units_angle)) THEN 00289 IF (units_angle == cp_units_rad) cell_angle = cell_angle/degree 00290 END IF 00291 END IF 00292 00293 IF (PRESENT(periodic)) CALL get_cell(cell=cell,periodic=periodic) 00294 00295 END SUBROUTINE get_cell_param 00296 00297 ! ***************************************************************************** 00305 SUBROUTINE set_cell_param(cell,cell_length,cell_angle,periodic,do_init_cell,& 00306 error) 00307 00308 TYPE(cell_type), POINTER :: cell 00309 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: cell_length, cell_angle 00310 INTEGER, DIMENSION(3), INTENT(IN), 00311 OPTIONAL :: periodic 00312 LOGICAL, INTENT(IN) :: do_init_cell 00313 TYPE(cp_error_type), INTENT(INOUT) :: error 00314 00315 CHARACTER(len=*), PARAMETER :: routineN = 'set_cell_param', 00316 routineP = moduleN//':'//routineN 00317 00318 LOGICAL :: failure 00319 REAL(KIND=dp) :: cos_alpha, cos_beta, 00320 cos_gamma, sin_gamma 00321 00322 failure = .FALSE. 00323 CPPrecondition(ALL(cell_angle/=0.0_dp),cp_failure_level,routineP,error,failure) 00324 cos_gamma = COS(cell_angle(3)); IF (ABS(cos_gamma)<EPSILON(0.0_dp)) cos_gamma = 0.0_dp 00325 sin_gamma = SIN(cell_angle(3)); IF (ABS(sin_gamma)<EPSILON(0.0_dp)) sin_gamma = 0.0_dp 00326 cos_beta = COS(cell_angle(2)); IF (ABS(cos_beta )<EPSILON(0.0_dp)) cos_beta = 0.0_dp 00327 cos_alpha = COS(cell_angle(1)); IF (ABS(cos_alpha)<EPSILON(0.0_dp)) cos_alpha = 0.0_dp 00328 00329 cell%hmat(:,1) = (/ 1.0_dp, 0.0_dp, 0.0_dp/) 00330 cell%hmat(:,2) = (/cos_gamma, sin_gamma, 0.0_dp/) 00331 cell%hmat(:,3) = (/ cos_beta, (cos_alpha-cos_gamma*cos_beta)/sin_gamma, 0.0_dp/) 00332 cell%hmat(3,3) = SQRT(1.0_dp-cell%hmat(1,3)**2-cell%hmat(2,3)**2) 00333 00334 cell%hmat(:,1) = cell%hmat(:,1)*cell_length(1) 00335 cell%hmat(:,2) = cell%hmat(:,2)*cell_length(2) 00336 cell%hmat(:,3) = cell%hmat(:,3)*cell_length(3) 00337 00338 IF (do_init_cell) THEN 00339 IF (PRESENT(periodic)) THEN 00340 CALL init_cell(cell=cell,periodic=periodic) 00341 ELSE 00342 CALL init_cell(cell=cell) 00343 END IF 00344 END IF 00345 00346 END SUBROUTINE set_cell_param 00347 00348 ! ***************************************************************************** 00354 SUBROUTINE init_cell(cell,hmat,periodic) 00355 00356 TYPE(cell_type), POINTER :: cell 00357 REAL(KIND=dp), DIMENSION(3, 3), 00358 INTENT(IN), OPTIONAL :: hmat 00359 INTEGER, DIMENSION(3), INTENT(IN), 00360 OPTIONAL :: periodic 00361 00362 CHARACTER(len=*), PARAMETER :: routineN = 'init_cell', 00363 routineP = moduleN//':'//routineN 00364 00365 INTEGER :: dim 00366 REAL(KIND=dp) :: a, acosa, acosah, alpha, 00367 asina, asinah, beta, norm, 00368 norm_c 00369 REAL(KIND=dp), DIMENSION(3) :: abc 00370 00371 IF (PRESENT(hmat)) cell%hmat(:,:) = hmat(:,:) 00372 IF (PRESENT(periodic)) cell%perd(:) = periodic(:) 00373 00374 cell%deth = ABS(det_3x3(cell%hmat)) 00375 00376 IF (cell%deth < 1.0E-10_dp) THEN 00377 CALL stop_program(routineN,moduleN,__LINE__,& 00378 "An invalid set of cell vectors was specified. "//& 00379 "The determinant det(h) is too small") 00380 END IF 00381 cell%h_inv = inv_3x3(cell%hmat) 00382 00383 SELECT CASE (cell%symmetry_id) 00384 CASE (cell_sym_cubic,& 00385 cell_sym_tetragonal_ab,& 00386 cell_sym_tetragonal_ac,& 00387 cell_sym_tetragonal_bc,& 00388 cell_sym_orthorhombic) 00389 CALL get_cell(cell=cell,abc=abc) 00390 abc(2) = plane_distance(0,1,0,cell=cell) 00391 abc(3) = plane_distance(0,0,1,cell=cell) 00392 SELECT CASE (cell%symmetry_id) 00393 CASE (cell_sym_cubic) 00394 abc(1:3) = SUM(abc(1:3))/3.0_dp 00395 CASE (cell_sym_tetragonal_ab,& 00396 cell_sym_tetragonal_ac,& 00397 cell_sym_tetragonal_bc) 00398 SELECT CASE (cell%symmetry_id) 00399 CASE (cell_sym_tetragonal_ab) 00400 a = 0.5_dp*(abc(1) + abc(2)) 00401 abc(1) = a 00402 abc(2) = a 00403 CASE (cell_sym_tetragonal_ac) 00404 a = 0.5_dp*(abc(1) + abc(3)) 00405 abc(1) = a 00406 abc(3) = a 00407 CASE (cell_sym_tetragonal_bc) 00408 a = 0.5_dp*(abc(2) + abc(3)) 00409 abc(2) = a 00410 abc(3) = a 00411 END SELECT 00412 END SELECT 00413 cell%hmat(1,1) = abc(1); cell%hmat(1,2) = 0.0_dp; cell%hmat(1,3) = 0.0_dp 00414 cell%hmat(2,1) = 0.0_dp; cell%hmat(2,2) = abc(2); cell%hmat(2,3) = 0.0_dp 00415 cell%hmat(3,1) = 0.0_dp; cell%hmat(3,2) = 0.0_dp; cell%hmat(3,3) = abc(3) 00416 CASE (cell_sym_hexagonal) 00417 CALL get_cell(cell=cell,abc=abc) 00418 a = 0.5_dp*(abc(1) + abc(2)) 00419 acosa = 0.5_dp*a 00420 asina = sqrt3*acosa 00421 cell%hmat(1,1) = a ; cell%hmat(1,2) = acosa ; cell%hmat(1,3) = 0.0_dp 00422 cell%hmat(2,1) = 0.0_dp; cell%hmat(2,2) = asina ; cell%hmat(2,3) = 0.0_dp 00423 cell%hmat(3,1) = 0.0_dp; cell%hmat(3,2) = 0.0_dp; cell%hmat(3,3) = abc(3) 00424 CASE (cell_sym_rhombohedral) 00425 CALL get_cell(cell=cell,abc=abc) 00426 a = SUM(abc(1:3))/3.0_dp 00427 alpha = (angle(cell%hmat(:,3),cell%hmat(:,2)) +& 00428 angle(cell%hmat(:,1),cell%hmat(:,3)) +& 00429 angle(cell%hmat(:,1),cell%hmat(:,2)))/3.0_dp 00430 acosa = a*COS(alpha) 00431 asina = a*SIN(alpha) 00432 acosah = a*COS(0.5_dp*alpha) 00433 asinah = a*SIN(0.5_dp*alpha) 00434 norm = acosa/acosah 00435 norm_c = SQRT(1.0_dp - norm*norm) 00436 cell%hmat(1,1) = a ; cell%hmat(1,2) = acosa ; cell%hmat(1,3) = acosah*norm 00437 cell%hmat(2,1) = 0.0_dp; cell%hmat(2,2) = asina ; cell%hmat(2,3) = asinah*norm 00438 cell%hmat(3,1) = 0.0_dp; cell%hmat(3,2) = 0.0_dp; cell%hmat(3,3) = a*norm_c 00439 CASE (cell_sym_monoclinic) 00440 CALL get_cell(cell=cell,abc=abc) 00441 beta = angle(cell%hmat(:,1),cell%hmat(:,3)) 00442 cell%hmat(1,1) = abc(1); cell%hmat(1,2) = 0.0_dp; cell%hmat(1,3) = abc(3)*COS(beta) 00443 cell%hmat(2,1) = 0.0_dp; cell%hmat(2,2) = abc(2); cell%hmat(2,3) = 0.0_dp 00444 cell%hmat(3,1) = 0.0_dp; cell%hmat(3,2) = 0.0_dp; cell%hmat(3,3) = abc(3)*SIN(beta) 00445 CASE (cell_sym_triclinic) 00446 ! Nothing to do 00447 END SELECT 00448 00449 IF ((cell%hmat(1,2) == 0.0_dp).AND.(cell%hmat(1,3) == 0.0_dp).AND.& 00450 (cell%hmat(2,1) == 0.0_dp).AND.(cell%hmat(2,3) == 0.0_dp).AND.& 00451 (cell%hmat(3,1) == 0.0_dp).AND.(cell%hmat(3,2) == 0.0_dp)) THEN 00452 cell%orthorhombic = .TRUE. 00453 ELSE 00454 cell%orthorhombic = .FALSE. 00455 END IF 00456 00457 dim = COUNT(cell%perd == 1) 00458 IF ((dim /= 3).AND.(.NOT.cell%orthorhombic)) THEN 00459 CALL stop_program(routineN,moduleN,__LINE__,& 00460 "Non-orthorhombic and not periodic") 00461 END IF 00462 00463 END SUBROUTINE init_cell 00464 00465 ! ***************************************************************************** 00472 FUNCTION plane_distance(h,k,l,cell) RESULT(distance) 00473 00474 INTEGER, INTENT(IN) :: h, k, l 00475 TYPE(cell_type), POINTER :: cell 00476 REAL(KIND=dp) :: distance 00477 00478 REAL(KIND=dp) :: a, alpha, b, beta, c, cosa, 00479 cosb, cosg, d, gamma, x, y, z 00480 REAL(KIND=dp), DIMENSION(3) :: abc 00481 00482 x = REAL(h,KIND=dp) 00483 y = REAL(k,KIND=dp) 00484 z = REAL(l,KIND=dp) 00485 00486 CALL get_cell(cell=cell,abc=abc) 00487 00488 a = abc(1) 00489 b = abc(2) 00490 c = abc(3) 00491 00492 IF (cell%orthorhombic) THEN 00493 00494 d = (x/a)**2 + (y/b)**2 + (z/c)**2 00495 00496 ELSE 00497 00498 CALL get_cell(cell=cell,& 00499 alpha=alpha,& 00500 beta=beta,& 00501 gamma=gamma) 00502 00503 alpha = alpha/degree 00504 beta = beta/degree 00505 gamma = gamma/degree 00506 00507 cosa = COS(alpha) 00508 cosb = COS(beta) 00509 cosg = COS(gamma) 00510 00511 d = ((x*b*c*SIN(alpha))**2 +& 00512 (y*c*a*SIN(beta))**2 +& 00513 (z*a*b*SIN(gamma))**2 +& 00514 2.0_dp*a*b*c*(x*y*c*(cosa*cosb - cosg) +& 00515 z*x*b*(cosg*cosa - cosb) +& 00516 y*z*a*(cosb*cosg - cosa)))/& 00517 ((a*b*c)**2*(1.0_dp - cosa**2 - cosb**2 - cosg**2 +& 00518 2.0_dp*cosa*cosb*cosg)) 00519 00520 END IF 00521 00522 distance = 1.0_dp/SQRT(d) 00523 00524 END FUNCTION plane_distance 00525 00526 ! ***************************************************************************** 00533 FUNCTION pbc1(r,cell) RESULT(r_pbc) 00534 00535 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 00536 TYPE(cell_type), POINTER :: cell 00537 REAL(KIND=dp), DIMENSION(3) :: r_pbc 00538 00539 REAL(KIND=dp), DIMENSION(3) :: s 00540 00541 IF (cell%orthorhombic) THEN 00542 r_pbc(1) = r(1) - cell%hmat(1,1)*cell%perd(1)*ANINT(cell%h_inv(1,1)*r(1)) 00543 r_pbc(2) = r(2) - cell%hmat(2,2)*cell%perd(2)*ANINT(cell%h_inv(2,2)*r(2)) 00544 r_pbc(3) = r(3) - cell%hmat(3,3)*cell%perd(3)*ANINT(cell%h_inv(3,3)*r(3)) 00545 ELSE 00546 s(1) = cell%h_inv(1,1)*r(1) + cell%h_inv(1,2)*r(2) + cell%h_inv(1,3)*r(3) 00547 s(2) = cell%h_inv(2,1)*r(1) + cell%h_inv(2,2)*r(2) + cell%h_inv(2,3)*r(3) 00548 s(3) = cell%h_inv(3,1)*r(1) + cell%h_inv(3,2)*r(2) + cell%h_inv(3,3)*r(3) 00549 s(1) = s(1) - cell%perd(1)*ANINT(s(1)) 00550 s(2) = s(2) - cell%perd(2)*ANINT(s(2)) 00551 s(3) = s(3) - cell%perd(3)*ANINT(s(3)) 00552 r_pbc(1) = cell%hmat(1,1)*s(1) + cell%hmat(1,2)*s(2) + cell%hmat(1,3)*s(3) 00553 r_pbc(2) = cell%hmat(2,1)*s(1) + cell%hmat(2,2)*s(2) + cell%hmat(2,3)*s(3) 00554 r_pbc(3) = cell%hmat(3,1)*s(1) + cell%hmat(3,2)*s(2) + cell%hmat(3,3)*s(3) 00555 END IF 00556 00557 END FUNCTION pbc1 00558 00559 ! ***************************************************************************** 00566 FUNCTION pbc2(r,cell,nl) RESULT(r_pbc) 00567 00568 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 00569 TYPE(cell_type), POINTER :: cell 00570 INTEGER, DIMENSION(3), INTENT(IN) :: nl 00571 REAL(KIND=dp), DIMENSION(3) :: r_pbc 00572 00573 REAL(KIND=dp), DIMENSION(3) :: s 00574 00575 IF (cell%orthorhombic) THEN 00576 r_pbc(1) = r(1) - cell%hmat(1,1)*cell%perd(1)*& 00577 REAL(NINT(cell%h_inv(1,1)*r(1)) - nl(1),dp) 00578 r_pbc(2) = r(2) - cell%hmat(2,2)*cell%perd(2)*& 00579 REAL(NINT(cell%h_inv(2,2)*r(2)) - nl(2),dp) 00580 r_pbc(3) = r(3) - cell%hmat(3,3)*cell%perd(3)*& 00581 REAL(NINT(cell%h_inv(3,3)*r(3)) - nl(3),dp) 00582 ELSE 00583 s(1) = cell%h_inv(1,1)*r(1) + cell%h_inv(1,2)*r(2) + cell%h_inv(1,3)*r(3) 00584 s(2) = cell%h_inv(2,1)*r(1) + cell%h_inv(2,2)*r(2) + cell%h_inv(2,3)*r(3) 00585 s(3) = cell%h_inv(3,1)*r(1) + cell%h_inv(3,2)*r(2) + cell%h_inv(3,3)*r(3) 00586 s(1) = s(1) - cell%perd(1)*REAL(NINT(s(1)) - nl(1),dp) 00587 s(2) = s(2) - cell%perd(2)*REAL(NINT(s(2)) - nl(2),dp) 00588 s(3) = s(3) - cell%perd(3)*REAL(NINT(s(3)) - nl(3),dp) 00589 r_pbc(1) = cell%hmat(1,1)*s(1) + cell%hmat(1,2)*s(2) + cell%hmat(1,3)*s(3) 00590 r_pbc(2) = cell%hmat(2,1)*s(1) + cell%hmat(2,2)*s(2) + cell%hmat(2,3)*s(3) 00591 r_pbc(3) = cell%hmat(3,1)*s(1) + cell%hmat(3,2)*s(2) + cell%hmat(3,3)*s(3) 00592 END IF 00593 00594 END FUNCTION pbc2 00595 00596 ! ***************************************************************************** 00603 FUNCTION pbc3(ra,rb,cell) RESULT(rab_pbc) 00604 00605 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: ra, rb 00606 TYPE(cell_type), POINTER :: cell 00607 REAL(KIND=dp), DIMENSION(3) :: rab_pbc 00608 00609 INTEGER :: icell, jcell, kcell 00610 INTEGER, DIMENSION(3) :: periodic 00611 REAL(KIND=dp) :: rab2, rab2_pbc 00612 REAL(KIND=dp), DIMENSION(3) :: r, ra_pbc, rab, rb_image, 00613 rb_pbc, s2r 00614 00615 CALL get_cell(cell=cell,periodic=periodic) 00616 00617 ra_pbc(:) = pbc(ra(:),cell) 00618 rb_pbc(:) = pbc(rb(:),cell) 00619 00620 rab2_pbc = HUGE(1.0_dp) 00621 00622 DO icell=-periodic(1),periodic(1) 00623 DO jcell=-periodic(2),periodic(2) 00624 DO kcell=-periodic(3),periodic(3) 00625 r = REAL((/icell,jcell,kcell/),dp) 00626 CALL scaled_to_real(s2r,r,cell) 00627 rb_image(:) = rb_pbc(:) + s2r 00628 rab(:) = rb_image(:) - ra_pbc(:) 00629 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) 00630 IF (rab2 < rab2_pbc) THEN 00631 rab2_pbc = rab2 00632 rab_pbc(:) = rab(:) 00633 END IF 00634 END DO 00635 END DO 00636 END DO 00637 00638 END FUNCTION pbc3 00639 00640 ! ***************************************************************************** 00647 RECURSIVE SUBROUTINE read_cell( cell, cell_ref, use_ref_cell, cell_section,& 00648 00649 check_for_ref, para_env, error) 00650 TYPE(cell_type), POINTER :: cell, cell_ref 00651 LOGICAL, INTENT(OUT), OPTIONAL :: use_ref_cell 00652 TYPE(section_vals_type), OPTIONAL, 00653 POINTER :: cell_section 00654 LOGICAL, INTENT(IN), OPTIONAL :: check_for_ref 00655 TYPE(cp_para_env_type), POINTER :: para_env 00656 TYPE(cp_error_type), INTENT(INOUT) :: error 00657 00658 CHARACTER(len=*), PARAMETER :: routineN = 'read_cell', 00659 routineP = moduleN//':'//routineN 00660 00661 INTEGER :: my_per 00662 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell 00663 LOGICAL :: cell_read_a, cell_read_abc, cell_read_b, cell_read_c, 00664 cell_read_file, check, failure, my_check 00665 REAL(KIND=dp), DIMENSION(:), POINTER :: cell_angles, cell_par 00666 TYPE(section_vals_type), POINTER :: cell_ref_section 00667 00668 failure = .FALSE. 00669 my_check = .TRUE. 00670 NULLIFY(cell_ref_section, cell_par, multiple_unit_cell) 00671 IF (.NOT.failure) THEN 00672 IF (.NOT.ASSOCIATED(cell)) CALL cell_create(cell,error=error) 00673 IF (.NOT.ASSOCIATED(cell_ref)) CALL cell_create(cell_ref,error=error) 00674 IF (PRESENT(check_for_ref)) my_check = check_for_ref 00675 00676 cell%deth = 0.0_dp 00677 cell%orthorhombic = .FALSE. 00678 cell%perd(:) = 1 00679 cell%symmetry_id = cell_sym_none 00680 cell%hmat(:,:) = 0.0_dp 00681 cell%h_inv(:,:) = 0.0_dp 00682 cell_read_file = .FALSE. 00683 cell_read_a = .FALSE. 00684 cell_read_b = .FALSE. 00685 cell_read_c = .FALSE. 00686 ! Trying to read cell info from file 00687 CALL section_vals_val_get(cell_section,"CELL_FILE_NAME",explicit=cell_read_file,error=error) 00688 IF (cell_read_file) CALL read_cell_from_external_file(cell_section, para_env, error) 00689 00690 ! Trying to read cell info from the separate A, B, C vectors 00691 ! If cell information is provided through file A,B,C contain the file information.. 00692 ! a print warning is shown on screen.. 00693 CALL section_vals_val_get(cell_section,"A",explicit=cell_read_a,error=error) 00694 IF (cell_read_a) THEN 00695 CALL section_vals_val_get(cell_section,"A",r_vals=cell_par,error=error) 00696 cell%hmat(:,1) = cell_par(:) 00697 END IF 00698 CALL section_vals_val_get(cell_section,"B",explicit=cell_read_b,error=error) 00699 IF (cell_read_b) THEN 00700 CALL section_vals_val_get(cell_section,"B",r_vals=cell_par,error=error) 00701 cell%hmat(:,2) = cell_par(:) 00702 END IF 00703 CALL section_vals_val_get(cell_section,"C",explicit=cell_read_c,error=error) 00704 IF (cell_read_c) THEN 00705 CALL section_vals_val_get(cell_section,"C",r_vals=cell_par,error=error) 00706 cell%hmat(:,3) = cell_par(:) 00707 END IF 00708 check = ((cell_read_a.EQV.cell_read_b).AND.(cell_read_b.EQV.cell_read_c)) 00709 CALL cp_assert(check,cp_warning_level,cp_assertion_failed,routineP,& 00710 "Cell Information provided through vectors A, B and C. Not all three "//& 00711 "vectors were provided! Cell setup may be incomplete!"//& 00712 CPSourceFileRef,& 00713 only_ionode=.TRUE.) 00714 00715 ! Very last option.. Trying to read cell info from ABC keyword 00716 CALL section_vals_val_get(cell_section,"ABC",explicit=cell_read_abc,error=error) 00717 IF (cell_read_abc) THEN 00718 check = (cell_read_a.OR.cell_read_b.OR.cell_read_c) 00719 CALL cp_assert(.NOT.check,cp_warning_level,cp_assertion_failed,routineP,& 00720 "Cell Information provided through vectors A, B and C in conjunction with ABC."//& 00721 " The definition of the ABC keyword will override the one provided by A,B and C."//& 00722 CPSourceFileRef,& 00723 only_ionode=.TRUE.) 00724 cell%hmat = 0.0_dp 00725 CALL section_vals_val_get(cell_section,"ABC",r_vals=cell_par,error=error) 00726 CALL section_vals_val_get(cell_section,"ALPHA_BETA_GAMMA",r_vals=cell_angles,error=error) 00727 CALL set_cell_param(cell,cell_par,cell_angles,do_init_cell=.FALSE.,error=error) 00728 END IF 00729 00730 ! Multiple unit cell 00731 CALL section_vals_val_get(cell_section,"MULTIPLE_UNIT_CELL",i_vals=multiple_unit_cell,error=error) 00732 IF (ANY(multiple_unit_cell/=1)) CALL set_multiple_unit_cell(cell, multiple_unit_cell, error) 00733 00734 CALL section_vals_val_get(cell_section,"PERIODIC",i_val=my_per,error=error) 00735 SELECT CASE(my_per) 00736 CASE(use_perd_x) 00737 cell%perd = (/1,0,0/) 00738 CASE(use_perd_y) 00739 cell%perd = (/0,1,0/) 00740 CASE(use_perd_z) 00741 cell%perd = (/0,0,1/) 00742 CASE(use_perd_xy) 00743 cell%perd = (/1,1,0/) 00744 CASE(use_perd_xz) 00745 cell%perd = (/1,0,1/) 00746 CASE(use_perd_yz) 00747 cell%perd = (/0,1,1/) 00748 CASE(use_perd_xyz) 00749 cell%perd = (/1,1,1/) 00750 CASE(use_perd_none) 00751 cell%perd = (/0,0,0/) 00752 CASE DEFAULT 00753 CPPostcondition(.FALSE.,cp_failure_level,routineP,error,failure) 00754 END SELECT 00755 00756 ! Load requested cell symmetry 00757 CALL section_vals_val_get(cell_section,"SYMMETRY",i_val=cell%symmetry_id,error=error) 00758 00759 ! Initialize cell 00760 CALL init_cell(cell) 00761 00762 IF (.NOT.my_check) RETURN 00763 cell_ref_section => section_vals_get_subs_vals(cell_section,& 00764 "CELL_REF",error=error) 00765 IF (parsed_cp2k_input(cell_ref_section,check_this_section=.TRUE.,error=error)) THEN 00766 IF(PRESENT(use_ref_cell) ) use_ref_cell = .TRUE. 00767 CALL read_cell(cell_ref, cell_ref, use_ref_cell, cell_section=cell_ref_section,& 00768 check_for_ref=.FALSE., para_env=para_env, error=error) 00769 ELSE 00770 CALL cell_clone (cell, cell_ref, error) 00771 IF ( PRESENT ( use_ref_cell ) ) use_ref_cell = .FALSE. 00772 END IF 00773 END IF 00774 00775 END SUBROUTINE read_cell 00776 00777 ! ***************************************************************************** 00783 SUBROUTINE set_multiple_unit_cell(cell, multiple_unit_cell, error) 00784 00785 TYPE(cell_type), POINTER :: cell 00786 INTEGER, DIMENSION(:), POINTER :: multiple_unit_cell 00787 TYPE(cp_error_type), INTENT(INOUT) :: error 00788 00789 CHARACTER(len=*), PARAMETER :: routineN = 'set_multiple_unit_cell', 00790 routineP = moduleN//':'//routineN 00791 00792 LOGICAL :: failure 00793 00794 failure = .FALSE. 00795 00796 ! Fail is one of the value is set to zero.. 00797 CALL cp_assert(ALL(multiple_unit_cell>0),cp_fatal_level,cp_assertion_failed,routineP,& 00798 "CELL%MULTIPLE_UNIT_CELL accepts only integer values larger than 0! "//& 00799 "A value of 0 or negative is meaningless!"//& 00800 CPSourceFileRef) 00801 00802 ! scale abc accordingly user request 00803 cell%hmat(:,1) = cell%hmat(:,1)*multiple_unit_cell(1) 00804 cell%hmat(:,2) = cell%hmat(:,2)*multiple_unit_cell(2) 00805 cell%hmat(:,3) = cell%hmat(:,3)*multiple_unit_cell(3) 00806 00807 END SUBROUTINE set_multiple_unit_cell 00808 00809 ! ***************************************************************************** 00815 SUBROUTINE read_cell_from_external_file(cell_section, para_env, error) 00816 00817 TYPE(section_vals_type), POINTER :: cell_section 00818 TYPE(cp_para_env_type), POINTER :: para_env 00819 TYPE(cp_error_type), INTENT(INOUT) :: error 00820 00821 CHARACTER(len=*), PARAMETER :: routineN = 'read_cell_from_external_file', 00822 routineP = moduleN//':'//routineN 00823 00824 CHARACTER(LEN=default_path_length) :: cell_file_name 00825 INTEGER :: i, idum, j, my_format, n_rep, 00826 stat 00827 LOGICAL :: explicit, failure, my_end 00828 REAL(KIND=dp) :: xdum 00829 REAL(KIND=dp), DIMENSION(3, 3) :: hmat 00830 REAL(KIND=dp), DIMENSION(:), POINTER :: cell_par 00831 TYPE(cp_parser_type), POINTER :: parser 00832 00833 failure = .FALSE. 00834 NULLIFY(parser) 00835 CALL section_vals_val_get(cell_section,"CELL_FILE_NAME",c_val=cell_file_name,error=error) 00836 CALL section_vals_val_get(cell_section,"CELL_FILE_FORMAT",i_val=my_format, error=error) 00837 CALL parser_create(parser,cell_file_name, para_env=para_env,error=error) 00838 CALL parser_get_next_line(parser,1,error=error) 00839 SELECT CASE(my_format) 00840 CASE (do_cell_cp2k) 00841 my_end = .FALSE. 00842 DO WHILE (.NOT.my_end) 00843 READ(parser%input_line,*)idum,xdum,hmat(:,1),hmat(:,2),hmat(:,3) 00844 CALL parser_get_next_line(parser,1,at_end=my_end,error=error) 00845 END DO 00846 CASE (do_cell_xsc) 00847 READ(parser%input_line,*)idum,hmat(:,1),hmat(:,2),hmat(:,3) 00848 END SELECT 00849 CALL parser_release(parser,error=error) 00850 CALL section_vals_val_unset(cell_section,"CELL_FILE_NAME",error=error) 00851 CALL section_vals_val_unset(cell_section,"CELL_FILE_FORMAT",error=error) 00852 ! Conver to CP2K units 00853 DO i = 1, 3 00854 DO j = 1, 3 00855 hmat(j,i) = cp_unit_to_cp2k(hmat(j,i), "angstrom", error=error) 00856 END DO 00857 END DO 00858 ! Check if the cell was already defined 00859 explicit = .FALSE. 00860 CALL section_vals_val_get(cell_section,"A",n_rep_val=n_rep,error=error) 00861 explicit = explicit .OR. (n_rep==1) 00862 CALL section_vals_val_get(cell_section,"B",n_rep_val=n_rep,error=error) 00863 explicit = explicit .OR. (n_rep==1) 00864 CALL section_vals_val_get(cell_section,"C",n_rep_val=n_rep,error=error) 00865 explicit = explicit .OR. (n_rep==1) 00866 CALL section_vals_val_get(cell_section,"ABC",n_rep_val=n_rep,error=error) 00867 explicit = explicit .OR. (n_rep==1) 00868 ! Possibly print a warning 00869 CALL cp_assert(.NOT.explicit,cp_warning_level,cp_assertion_failed,routineP,& 00870 "Cell specification (A,B,C or ABC) provided together with the external "//& 00871 "cell setup! Ignoring (A,B,C or ABC) and proceeding with info read from the "//& 00872 "external file! "//& 00873 CPSourceFileRef,& 00874 only_ionode=.TRUE.) 00875 ! Copy cell information in the A, B, C fields..(we may need them later on..) 00876 ALLOCATE(cell_par(3), stat=stat) 00877 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00878 cell_par = hmat(:,1) 00879 CALL section_vals_val_set(cell_section,"A",r_vals_ptr=cell_par,error=error) 00880 ALLOCATE(cell_par(3), stat=stat) 00881 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00882 cell_par = hmat(:,2) 00883 CALL section_vals_val_set(cell_section,"B",r_vals_ptr=cell_par,error=error) 00884 ALLOCATE(cell_par(3), stat=stat) 00885 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 00886 cell_par = hmat(:,3) 00887 CALL section_vals_val_set(cell_section,"C",r_vals_ptr=cell_par,error=error) 00888 ! Unset possible keywords 00889 CALL section_vals_val_unset(cell_section,"ABC",error=error) 00890 CALL section_vals_val_unset(cell_section,"ALPHA_BETA_GAMMA",error=error) 00891 00892 END SUBROUTINE read_cell_from_external_file 00893 00894 ! ***************************************************************************** 00901 SUBROUTINE real_to_scaled(s,r,cell) 00902 00903 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: s 00904 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: r 00905 TYPE(cell_type), POINTER :: cell 00906 00907 IF (cell%orthorhombic) THEN 00908 s(1) = cell%h_inv(1,1)*r(1) 00909 s(2) = cell%h_inv(2,2)*r(2) 00910 s(3) = cell%h_inv(3,3)*r(3) 00911 ELSE 00912 s(1) = cell%h_inv(1,1)*r(1) + cell%h_inv(1,2)*r(2) + cell%h_inv(1,3)*r(3) 00913 s(2) = cell%h_inv(2,1)*r(1) + cell%h_inv(2,2)*r(2) + cell%h_inv(2,3)*r(3) 00914 s(3) = cell%h_inv(3,1)*r(1) + cell%h_inv(3,2)*r(2) + cell%h_inv(3,3)*r(3) 00915 END IF 00916 00917 END SUBROUTINE real_to_scaled 00918 00919 ! ***************************************************************************** 00926 SUBROUTINE scaled_to_real(r,s,cell) 00927 00928 REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: r 00929 REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: s 00930 TYPE(cell_type), POINTER :: cell 00931 00932 IF (cell%orthorhombic) THEN 00933 r(1) = cell%hmat(1,1)*s(1) 00934 r(2) = cell%hmat(2,2)*s(2) 00935 r(3) = cell%hmat(3,3)*s(3) 00936 ELSE 00937 r(1) = cell%hmat(1,1)*s(1) + cell%hmat(1,2)*s(2) + cell%hmat(1,3)*s(3) 00938 r(2) = cell%hmat(2,1)*s(1) + cell%hmat(2,2)*s(2) + cell%hmat(2,3)*s(3) 00939 r(3) = cell%hmat(3,1)*s(1) + cell%hmat(3,2)*s(2) + cell%hmat(3,3)*s(3) 00940 END IF 00941 00942 END SUBROUTINE scaled_to_real 00943 00944 ! ***************************************************************************** 00952 RECURSIVE SUBROUTINE write_cell(cell,subsys_section,cell_ref,label,error) 00953 00954 TYPE(cell_type), POINTER :: cell 00955 TYPE(section_vals_type), POINTER :: subsys_section 00956 TYPE(cell_type), OPTIONAL, POINTER :: cell_ref 00957 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: label 00958 TYPE(cp_error_type), INTENT(inout) :: error 00959 00960 CHARACTER(len=*), PARAMETER :: routineN = 'write_cell', 00961 routineP = moduleN//':'//routineN 00962 00963 CHARACTER(LEN=default_string_length) :: my_label, unit_str 00964 INTEGER :: output_unit 00965 REAL(KIND=dp) :: alpha, beta, gamma, val 00966 REAL(KIND=dp), DIMENSION(3) :: abc 00967 TYPE(cp_logger_type), POINTER :: logger 00968 TYPE(enumeration_type), POINTER :: enum 00969 TYPE(keyword_type), POINTER :: keyword 00970 TYPE(section_type), POINTER :: section 00971 00972 NULLIFY (enum) 00973 NULLIFY (keyword) 00974 NULLIFY (logger) 00975 NULLIFY (section) 00976 logger => cp_error_get_logger(error) 00977 my_label = "CELL|" 00978 IF (PRESENT(label)) my_label = TRIM(label) 00979 output_unit = cp_print_key_unit_nr(logger,subsys_section,"PRINT%CELL",& 00980 extension=".Log",error=error) 00981 CALL section_vals_val_get(subsys_section,"PRINT%CELL%UNIT",c_val=unit_str,error=error) 00982 IF (output_unit > 0) THEN 00983 CALL get_cell(cell=cell,abc=abc,alpha=alpha,beta=beta,gamma=gamma) 00984 WRITE (UNIT=output_unit, FMT='( )') 00985 val = cp_unit_from_cp2k(cell%deth,TRIM(unit_str)//"^3",error=error) 00986 WRITE (UNIT=output_unit,FMT="(T2,A,T61,F20.3)")& 00987 TRIM(my_label)//" Volume ["//TRIM(unit_str)//"^3]:",val 00988 val = cp_unit_from_cp2k(1.0_dp,TRIM(unit_str),error=error) 00989 WRITE (UNIT=output_unit,FMT="(T2,A,T30,3F10.3,4X,A6,F11.3)")& 00990 TRIM(my_label)//" Vector a ["//TRIM(unit_str)//"]:",cell%hmat(:,1)*val,& 00991 "|a| = ",abc(1)*val,& 00992 TRIM(my_label)//" Vector b ["//TRIM(unit_str)//"]:",cell%hmat(:,2)*val,& 00993 "|b| = ",abc(2)*val,& 00994 TRIM(my_label)//" Vector c ["//TRIM(unit_str)//"]:",cell%hmat(:,3)*val,& 00995 "|c| = ",abc(3)*val 00996 WRITE (UNIT=output_unit,FMT="(T2,A,T70,F11.3)")& 00997 TRIM(my_label)//" Angle (b,c), alpha [degree]: ",alpha,& 00998 TRIM(my_label)//" Angle (a,c), beta [degree]: ",beta,& 00999 TRIM(my_label)//" Angle (a,b), gamma [degree]: ",gamma 01000 IF (cell%symmetry_id /= cell_sym_none) THEN 01001 CALL create_cell_section(section,error=error) 01002 keyword => section_get_keyword(section,"SYMMETRY",error=error) 01003 CALL keyword_get(keyword,enum=enum,error=error) 01004 WRITE (UNIT=output_unit,FMT="(T2,A,T61,A20)")& 01005 TRIM(my_label)//" Requested initial symmetry: ",& 01006 ADJUSTR(TRIM(enum_i2c(enum,cell%symmetry_id,error=error))) 01007 CALL section_release(section,error=error) 01008 ELSE 01009 IF (cell%orthorhombic) THEN 01010 WRITE (UNIT=output_unit,FMT="(T2,A,T78,A3)")& 01011 TRIM(my_label)//" Orthorhombic: ","YES" 01012 ELSE 01013 WRITE (UNIT=output_unit,FMT="(T2,A,T78,A3)")& 01014 TRIM(my_label)//" Orthorhombic: "," NO" 01015 END IF 01016 END IF 01017 END IF 01018 CALL cp_print_key_finished_output(output_unit,logger,subsys_section,& 01019 "PRINT%CELL",error=error) 01020 01021 IF (PRESENT(cell_ref)) THEN 01022 CALL write_cell(cell_ref, subsys_section, label="CELL_REF|", error=error) 01023 END IF 01024 01025 END SUBROUTINE write_cell 01026 01027 ! ***************************************************************************** 01038 SUBROUTINE cell_create(cell,hmat,periodic,error) 01039 01040 TYPE(cell_type), POINTER :: cell 01041 REAL(KIND=dp), DIMENSION(3, 3), 01042 INTENT(IN), OPTIONAL :: hmat 01043 INTEGER, DIMENSION(3), INTENT(IN), 01044 OPTIONAL :: periodic 01045 TYPE(cp_error_type), INTENT(inout) :: error 01046 01047 CHARACTER(len=*), PARAMETER :: routineN = 'cell_create', 01048 routineP = moduleN//':'//routineN 01049 01050 INTEGER :: stat 01051 LOGICAL :: failure 01052 01053 failure = .FALSE. 01054 01055 CPPrecondition(.NOT.ASSOCIATED(cell),cp_failure_level,routineP,error,failure) 01056 ALLOCATE (cell,stat=stat) 01057 CPPostcondition(stat==0,cp_failure_level,routineP,error,failure) 01058 last_cell_id = last_cell_id + 1 01059 cell%id_nr = last_cell_id 01060 cell%ref_count = 1 01061 IF (PRESENT(periodic)) THEN 01062 cell%perd = periodic 01063 ELSE 01064 cell%perd = 1 01065 END IF 01066 cell%orthorhombic = .FALSE. 01067 cell%symmetry_id = cell_sym_none 01068 IF (.NOT.PRESENT(hmat)) RETURN 01069 IF (.NOT.failure) THEN 01070 CALL init_cell(cell,hmat,periodic) 01071 END IF 01072 01073 END SUBROUTINE cell_create 01074 01075 ! ***************************************************************************** 01084 SUBROUTINE cell_retain(cell,error) 01085 01086 TYPE(cell_type), POINTER :: cell 01087 TYPE(cp_error_type), INTENT(inout) :: error 01088 01089 CHARACTER(len=*), PARAMETER :: routineN = 'cell_retain', 01090 routineP = moduleN//':'//routineN 01091 01092 LOGICAL :: failure 01093 01094 failure = .FALSE. 01095 01096 CPPrecondition(ASSOCIATED(cell),cp_failure_level,routineP,error,failure) 01097 IF (.NOT. failure) THEN 01098 CPPreconditionNoFail(cell%ref_count>0,cp_failure_level,routineP,error) 01099 cell%ref_count=cell%ref_count+1 01100 END IF 01101 01102 END SUBROUTINE cell_retain 01103 01104 ! ***************************************************************************** 01113 SUBROUTINE cell_release(cell,error) 01114 01115 TYPE(cell_type), POINTER :: cell 01116 TYPE(cp_error_type), INTENT(inout) :: error 01117 01118 CHARACTER(len=*), PARAMETER :: routineN = 'cell_release', 01119 routineP = moduleN//':'//routineN 01120 01121 INTEGER :: stat 01122 LOGICAL :: failure 01123 01124 failure = .FALSE. 01125 01126 IF (ASSOCIATED(cell)) THEN 01127 CPPreconditionNoFail(cell%ref_count>0,cp_failure_level,routineP,error) 01128 cell%ref_count=cell%ref_count-1 01129 IF (cell%ref_count==0) THEN 01130 DEALLOCATE(cell,stat=stat) 01131 CPPostconditionNoFail(stat==0,cp_warning_level,routineP,error) 01132 END IF 01133 NULLIFY(cell) 01134 END IF 01135 01136 END SUBROUTINE cell_release 01137 01138 ! ***************************************************************************** 01145 SUBROUTINE pbc_cp2k_plumed(x,y,z) 01146 REAL(KIND=dp), INTENT(INOUT) :: x, y, z 01147 01148 REAL(KIND=dp), DIMENSION(3) :: r, r2 01149 TYPE(cell_type), POINTER :: cell 01150 01151 CALL pbc_cp2k_plumed_getset_cell(cell, set=.FALSE.) !retrieve the cell pointer 01152 01153 r(1)=x 01154 r(2)=y 01155 r(3)=z 01156 01157 r2 = pbc(r,cell) 01158 01159 x=r2(1) 01160 y=r2(2) 01161 z=r2(3) 01162 01163 END SUBROUTINE pbc_cp2k_plumed 01164 01165 ! ***************************************************************************** 01173 SUBROUTINE pbc_cp2k_plumed_getset_cell(cell, set) 01174 TYPE(cell_type), POINTER :: cell 01175 LOGICAL :: set 01176 01177 TYPE(cell_type), POINTER, SAVE :: stored_cell 01178 01179 IF (set .EQV. .TRUE.) THEN 01180 stored_cell => cell 01181 ELSE 01182 cell => stored_cell 01183 END IF 01184 01185 END SUBROUTINE pbc_cp2k_plumed_getset_cell 01186 END MODULE cell_types
1.7.3